20 #include "pism/util/TerminationReason.hh"
21 #include "pism/util/pism_options.hh"
22 #include "pism/util/ConfigInterface.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/Context.hh"
25 #include "pism/util/petscwrappers/Vec.hh"
34 : m_design_stencil_width(d0.stencil_width()),
35 m_state_stencil_width(u_obs.stencil_width()),
36 m_ssaforward(ssaforward),
37 m_x(d0.grid(),
"x",
WITH_GHOSTS, m_design_stencil_width),
40 m_tmp_D1Local(d0.grid(),
"work vector",
WITH_GHOSTS, m_design_stencil_width),
41 m_tmp_D2Local(d0.grid(),
"work vector",
WITH_GHOSTS, m_design_stencil_width),
44 m_tmp_S1Local(d0.grid(),
"work vector",
WITH_GHOSTS, m_state_stencil_width),
45 m_tmp_S2Local(d0.grid(),
"work vector",
WITH_GHOSTS, m_state_stencil_width),
49 m_d_diff(d0.grid(),
"d_diff",
WITH_GHOSTS, m_design_stencil_width),
50 m_d_diff_lin(d0.grid(),
"d_diff linearized",
WITH_GHOSTS, m_design_stencil_width),
51 m_h(d0.grid(),
"h",
WITH_GHOSTS, m_design_stencil_width),
54 m_dh_dalpha(d0.grid(),
"dh_dalpha",
WITH_GHOSTS, m_design_stencil_width),
60 m_u_diff(d0.grid(),
"du",
WITH_GHOSTS, m_state_stencil_width),
62 m_designFunctional(designFunctional),
63 m_stateFunctional(stateFunctional),
75 ierr = KSPSetOptionsPrefix(
m_ksp,
"inv_gn_");
76 PISM_CHK(ierr,
"KSPSetOptionsPrefix");
78 double ksp_rtol = 1e-5;
79 ierr = KSPSetTolerances(
m_ksp, ksp_rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
82 ierr = KSPSetType(
m_ksp, KSPCG);
86 ierr = KSPGetPC(
m_ksp, &pc);
89 ierr = PCSetType(pc, PCNONE);
92 ierr = KSPSetFromOptions(
m_ksp);
95 int nLocalNodes = grid->xm()*grid->ym();
96 int nGlobalNodes = grid->Mx()*grid->My();
97 ierr = MatCreateShell(grid->com, nLocalNodes, nLocalNodes,
113 m_tikhonov_atol = grid->ctx()->config()->get_number(
"inverse.tikhonov.atol");
114 m_tikhonov_rtol = grid->ctx()->config()->get_number(
"inverse.tikhonov.rtol");
115 m_tikhonov_ptol = grid->ctx()->config()->get_number(
"inverse.tikhonov.ptol");
138 ierr = DMGlobalToLocalBegin(*
m_x.
dm(), x, INSERT_VALUES,
m_x.
vec());
139 PISM_CHK(ierr,
"DMGlobalToLocalBegin");
141 ierr = DMGlobalToLocalEnd(*
m_x.
dm(), x, INSERT_VALUES,
m_x.
vec());
142 PISM_CHK(ierr,
"DMGlobalToLocalEnd");
155 ierr = VecCopy(GNx.
vec(), y);
PISM_CHK(ierr,
"VecCopy");
182 KSPConvergedReason ksp_reason;
183 ierr = KSPGetConvergedReason(
m_ksp ,&ksp_reason);
184 PISM_CHK(ierr,
"KSPGetConvergedReason");
206 *value =
m_alpha*dValue + sValue;
212 double designNorm, stateNorm, sumNorm;
213 double dWeight, sWeight;
220 designNorm *= dWeight;
221 stateNorm *= sWeight;
226 "----------------------------------------------------------\n");
228 "IP_SSATaucTikhonovGNSolver Iteration %d: misfit %g; functional %g \n",
233 double relsum = (sumNorm/
std::max(designNorm,stateNorm));
235 "design norm %g stateNorm %g sum %g; relative difference %g\n",
236 designNorm, stateNorm, sumNorm, relsum);
265 if (reason->failed()) {
286 double valDesign, valState;
305 double descent_derivative;
312 if (descent_derivative >=0) {
313 printf(
"descent derivative: %g\n",descent_derivative);
322 if (step_reason->succeeded()) {
323 if (
m_value <= old_value + 1e-3*alpha*descent_derivative) {
328 printf(
"forward solve failed in linsearch. Shrinking.\n");
332 printf(
"alpha= %g; derivative = %g\n",alpha,descent_derivative);
345 " IP_SSATaucTikhonovGNSolver::solve.");
351 double dlogalpha = 0;
356 if (step_reason->failed()) {
358 reason->set_root_cause(step_reason);
365 if (reason->done()) {
375 if (step_reason->failed()) {
377 reason->set_root_cause(step_reason);
382 if (step_reason->failed()) {
385 reason->set_root_cause(step_reason);
391 if (step_reason->failed()) {
394 reason->set_root_cause(step_reason);
424 KSPConvergedReason ksp_reason;
425 ierr = KSPGetConvergedReason(
m_ksp,&ksp_reason);
426 PISM_CHK(ierr,
"KSPGetConvergedReason");
450 double ddisc_sq_dalpha;
454 if (ddisc_sq_dalpha <= 0) {
458 "Adaptive Tikhonov sanity check failed (dh/dalpha= %g <= 0)."
459 " Tighten inv_gn_ksp_rtol?\n",
466 double ddisc_sq_dalpha_a;
468 double ddisc_sq_dalpha_b;
470 ddisc_sq_dalpha = 2*
m_alpha*(ddisc_sq_dalpha_a+
m_alpha*ddisc_sq_dalpha_b);
473 "Adaptive Tikhonov sanity check recovery attempt: dh/dalpha= %g. \n",
487 if (fabs(*dlogalpha)> stepmax) {
488 double sgn = *dlogalpha > 0 ? 1 : -1;
489 *dlogalpha = stepmax*sgn;
static TerminationReason::Ptr keep_iterating()
static TerminationReason::Ptr success()
static TerminationReason::Ptr max_iter()
std::shared_ptr< const IceGrid > ConstPtr
void add(double alpha, const IceModelVec2S &x)
void copy_from(const IceModelVec2S &source)
void copy_from(const IceModelVec2< T > &source)
void add(double alpha, const IceModelVec2< T > &x)
void update_ghosts()
Updates ghost points.
std::shared_ptr< petsc::DM > dm() const
void set(double c)
Result: v[j] <- c for all j.
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
IceGrid::ConstPtr grid() const
std::vector< double > norm(int n) const
Computes the norm of all the components of an IceModelVec.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
std::shared_ptr< TerminationReason > Ptr
Abstract base class for IPFunctionals arising from an inner product.
virtual IceModelVec2V::Ptr solution()
Returns the last solution of the SSA as computed by linearize_at.
virtual void apply_linearization_transpose(IceModelVec2V &du, IceModelVec2S &dzeta)
Applies the transpose of the linearization of the forward map (i.e. the transpose of the reduced grad...
virtual void apply_linearization(IceModelVec2S &dzeta, IceModelVec2V &du)
Applies the linearization of the forward map (i.e. the reduced gradient described in the class-level...
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 void evaluateGNFunctional(DesignVec &h, double *value)
virtual TerminationReason::Ptr init()
virtual TerminationReason::Ptr check_convergence()
virtual void assemble_GN_rhs(DesignVec &out)
virtual TerminationReason::Ptr linesearch()
IPInnerProductFunctional< StateVec > & m_stateFunctional
virtual TerminationReason::Ptr evaluate_objective_and_gradient()
DesignVec m_dh_dalphaGlobal
IPInnerProductFunctional< DesignVec > & m_designFunctional
virtual TerminationReason::Ptr compute_dlogalpha(double *dalpha)
const unsigned int m_design_stencil_width
IP_SSATaucForwardProblem & m_ssaforward
virtual void apply_GN(IceModelVec2S &h, IceModelVec2S &out)
virtual TerminationReason::Ptr solve()
IP_SSATaucTikhonovGNSolver(IP_SSATaucForwardProblem &ssaforward, DesignVec &d0, StateVec &u_obs, double eta, IPInnerProductFunctional< DesignVec > &designFunctional, IPInnerProductFunctional< StateVec > &stateFunctional)
virtual TerminationReason::Ptr solve_linearized()
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
bool Bool(const std::string &option, const std::string &description)
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
std::string printf(const char *format,...)