20 #include "pism/util/IceGrid.hh"
21 #include "pism/util/ConfigInterface.hh"
22 #include "pism/util/Context.hh"
39 : m_ssaforward(ssaforward),
40 m_dGlobal(d0.grid(),
"design variable (global)",
WITHOUT_GHOSTS, d0.stencil_width()),
42 m_dzeta(d0.grid(),
"dzeta",
WITH_GHOSTS, d0.stencil_width()),
43 m_u(d0.grid(),
"state variable",
WITH_GHOSTS, u_obs.stencil_width()),
44 m_du(d0.grid(),
"du",
WITH_GHOSTS, u_obs.stencil_width()),
47 m_d_Jdesign(d0.grid(),
"Jdesign design variable",
WITH_GHOSTS, d0.stencil_width()),
48 m_u_Jdesign(d0.grid(),
"Jdesign state variable",
WITH_GHOSTS, u_obs.stencil_width()),
49 m_designFunctional(designFunctional),
50 m_stateFunctional(stateFunctional)
56 double stressScale = grid->ctx()->config()->get_number(
"inverse.design.param_tauc_scale");
59 m_velocityScale = grid->ctx()->config()->get_number(
"inverse.ssa.velocity_scale",
"m second-1");
83 ierr = DMSetMatType(da, MATBAIJ);
89 int nLocalNodes = grid->xm()*grid->ym();
90 int nGlobalNodes = grid->Mx()*grid->My();
91 ierr = MatCreateShell(grid->com, 2*nLocalNodes, nLocalNodes, 2*nGlobalNodes, nGlobalNodes,
95 ierr = MatShellSetOperation(
m_Jdesign, MATOP_MULT,
97 PISM_CHK(ierr,
"MatShellSetOperation");
99 ierr = MatShellSetOperation(
m_Jdesign, MATOP_MULT_TRANSPOSE,
101 PISM_CHK(ierr,
"MatShellSetOperation");
119 m_x->scatterToA(
m_d->vec());
125 ierr = TaoSetStateDesignIS(tao,
126 m_x->blockBIndexSet() ,
127 m_x->blockAIndexSet() );
128 PISM_CHK(ierr,
"TaoSetStateDesignIS");
145 ierr = TaoGetSolutionStatus(tao, &its, NULL, NULL, NULL, NULL, NULL);
146 PISM_CHK(ierr,
"TaoGetSolutionStatus");
149 for (
int k = 0;
k < nListeners;
k++) {
161 double *value, Vec gradient) {
190 if (reason->failed()) {
240 *s = SAME_NONZERO_PATTERN;
261 PISM_CHK(ierr,
"DMGlobalToLocalBegin");
264 PISM_CHK(ierr,
"DMGlobalToLocalEnd");
278 ierr = DMGlobalToLocalBegin(*
m_du.
dm(), x, INSERT_VALUES,
m_du.
vec());
279 PISM_CHK(ierr,
"DMGlobalToLocalBegin");
281 ierr = DMGlobalToLocalEnd(*
m_du.
dm(), x, INSERT_VALUES,
m_du.
vec());
282 PISM_CHK(ierr,
"DMGlobalToLocalEnd");
295 PetscErrorCode ierr = MatShellGetContext(A,&ctx);
296 PISM_CHK(ierr,
"MatShellGetContext");
300 MPI_Comm com = MPI_COMM_SELF;
301 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)A, &com); CHKERRQ(ierr);
303 SETERRQ(com, 1,
"A PISM callback failed");
312 PetscErrorCode ierr = MatShellGetContext(A,&ctx);
313 PISM_CHK(ierr,
"MatShellGetContext");
317 MPI_Comm com = MPI_COMM_SELF;
318 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)A, &com); CHKERRQ(ierr);
320 SETERRQ(com, 1,
"A PISM callback failed");
static TerminationReason::Ptr success()
std::shared_ptr< const IceGrid > ConstPtr
void copy_from(const IceModelVec2S &source)
std::shared_ptr< IceModelVec2S > Ptr
std::shared_ptr< IceModelVec2V > Ptr
void copy_from(const IceModelVec2< T > &source)
unsigned int stencil_width() const
Get the stencil width of the current IceModelVec. Returns 0 if ghosts are not available.
std::shared_ptr< petsc::DM > dm() const
IceGrid::ConstPtr grid() const
std::shared_ptr< TerminationReason > Ptr
virtual void valueAt(IMVecType &x, double *OUTPUT)=0
Computes the value of the functional at the vector x.
virtual void gradientAt(IMVecType &x, IMVecType &gradient)=0
Computes the gradient of the functional at the vector x.
Abstract base class for functions from ice model vectors to .
virtual void set_design(IceModelVec2S &zeta)
Sets the current value of of the design paramter .
virtual void assemble_residual(IceModelVec2V &u, IceModelVec2V &R)
Computes the residual function as defined in the class-level documentation.
virtual void apply_jacobian_design(IceModelVec2V &u, IceModelVec2S &dzeta, IceModelVec2V &du)
Applies the design Jacobian matrix to a perturbation of the design variable.
petsc::DM & get_da() const
Exposes the DMDA of the underlying grid for the benefit of TAO.
virtual IceModelVec2V::Ptr solution()
Returns the last solution of the SSA as computed by linearize_at.
virtual void assemble_jacobian_state(IceModelVec2V &u, Mat J)
Assembles the state Jacobian matrix.
virtual void apply_jacobian_design_transpose(IceModelVec2V &u, IceModelVec2V &du, IceModelVec2S &dzeta)
Applies the transpose of the design Jacobian matrix to a perturbation of the state variable.
virtual TerminationReason::Ptr linearize_at(IceModelVec2S &zeta)
Sets the current value of the design variable and solves the SSA to find the associated .
Implements the forward problem of the map taking to the corresponding solution of the SSA.
virtual DesignVec::Ptr designSolution()
virtual void applyConstraintsJacobianDesignTranspose(Vec x, Vec y)
IceModelVec2V m_u_Jdesign
IPFunctional< IceModelVec2V > & m_stateFunctional
virtual void evaluateObjectiveAndGradient(Tao tao, Vec x, double *value, Vec gradient)
virtual StateVec::Ptr stateSolution()
static PetscErrorCode jacobian_design_transpose_callback(Mat A, Vec x, Vec y)
virtual void evaluateConstraints(Tao tao, Vec x, Vec r)
IPFunctional< IceModelVec2S > & m_designFunctional
IceModelVec2S m_d_Jdesign
virtual void setInitialGuess(DesignVec &d0)
std::vector< Listener::Ptr > m_listeners
virtual void evaluateConstraintsJacobianState(Tao tao, Vec x, Mat Jstate, Mat Jpc, Mat Jinv, MatStructure *s)
DesignVec::Ptr m_grad_design
IP_SSATaucTaoTikhonovProblemLCL(IP_SSATaucForwardProblem &ssaforward, DesignVec &d0, StateVec &u_obs, double eta, IPFunctional< DesignVec > &designFunctional, IPFunctional< StateVec > &stateFunctional)
std::unique_ptr< IPTwoBlockVec > m_x
IP_SSATaucForwardProblem & m_ssaforward
StateVec::Ptr m_constraints
virtual TerminationReason::Ptr formInitialGuess(Vec *x)
virtual void applyConstraintsJacobianDesign(Vec x, Vec y)
virtual void evaluateConstraintsJacobianDesign(Tao tao, Vec x, Mat Jdesign)
StateVec::Ptr m_grad_state
static PetscErrorCode jacobian_design_callback(Mat A, Vec x, Vec y)
double m_constraintsScale
Defines a Tikhonov minimization problem of determining from SSA velocities to be solved with a TaoBa...
static void connect(Tao tao, Problem &p, Vec c, Mat Jc, Mat Jd, Mat Jcpc=NULL, Mat Jcinv=NULL)
static void connect(Tao tao, Problem &p)
Adaptor to connect a TAO objective and gradient function callback to a C++ object method.
#define PISM_CHK(errcode, name)
void handle_fatal_errors(MPI_Comm com)