PISM, A Parallel Ice Sheet Model
stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
|
Implements the forward problem of the map taking \(\tau_c\) to the corresponding solution of the SSA. More...
#include <IP_SSATaucForwardProblem.hh>
Public Types | |
typedef array::Scalar | DesignVec |
The function space for the design variable, i.e. \(\tau_c\). More... | |
typedef array::Scalar1 | DesignVecGhosted |
typedef array::Vector | StateVec |
The function space for the state variable, \(u_{\rm SSA}\). More... | |
typedef array::Vector1 | StateVec1 |
Public Member Functions | |
IP_SSATaucForwardProblem (std::shared_ptr< const Grid > g, IPDesignVariableParameterization &tp) | |
virtual | ~IP_SSATaucForwardProblem ()=default |
void | init () |
virtual void | set_tauc_fixed_locations (array::Scalar &locations) |
Selects nodes where \(\tau_c\) (more specifically \(\zeta\)) should not be adjusted. More... | |
virtual std::shared_ptr< array::Vector > | solution () |
Returns the last solution of the SSA as computed by linearize_at. More... | |
virtual IPDesignVariableParameterization & | tauc_param () |
Exposes the \(\tau_c\) parameterization being used. More... | |
virtual void | set_design (array::Scalar &zeta) |
Sets the current value of of the design paramter \(\zeta\). More... | |
virtual std::shared_ptr< TerminationReason > | linearize_at (array::Scalar &zeta) |
Sets the current value of the design variable \(\zeta\) and solves the SSA to find the associated \(u_{\rm SSA}\). More... | |
virtual void | assemble_residual (array::Vector &u, array::Vector &R) |
Computes the residual function \(\mathcal{R}(u, \zeta)\) as defined in the class-level documentation. More... | |
virtual void | assemble_residual (array::Vector &u, Vec R) |
Computes the residual function \(\mathcal{R}(u, \zeta)\) defined in the class-level documentation. More... | |
virtual void | assemble_jacobian_state (array::Vector &u, Mat J) |
Assembles the state Jacobian matrix. More... | |
virtual void | apply_jacobian_design (array::Vector &u, array::Scalar &dzeta, array::Vector &du) |
Applies the design Jacobian matrix to a perturbation of the design variable. More... | |
virtual void | apply_jacobian_design (array::Vector &u, array::Scalar &dzeta, Vec du) |
Applies the design Jacobian matrix to a perturbation of the design variable. More... | |
virtual void | apply_jacobian_design (array::Vector &u, array::Scalar &dzeta, Vector2d **du_a) |
Applies the design Jacobian matrix to a perturbation of the design variable. More... | |
virtual void | apply_jacobian_design_transpose (array::Vector &u, array::Vector &du, array::Scalar &dzeta) |
Applies the transpose of the design Jacobian matrix to a perturbation of the state variable. More... | |
virtual void | apply_jacobian_design_transpose (array::Vector &u, array::Vector &du, Vec dzeta) |
Applies the transpose of the design Jacobian matrix to a perturbation of the state variable. More... | |
virtual void | apply_jacobian_design_transpose (array::Vector &u, array::Vector &du, double **dzeta) |
Applies the transpose of the design Jacobian matrix to a perturbation of the state variable. More... | |
virtual void | apply_linearization (array::Scalar &dzeta, array::Vector &du) |
Applies the linearization of the forward map (i.e. the reduced gradient \(DF\) described in the class-level documentation.) More... | |
virtual void | apply_linearization_transpose (array::Vector &du, array::Scalar &dzeta) |
Applies the transpose of the linearization of the forward map (i.e. the transpose of the reduced gradient \(DF\) described in the class-level documentation.) More... | |
petsc::DM & | get_da () const |
Exposes the DMDA of the underlying grid for the benefit of TAO. More... | |
Public Member Functions inherited from pism::stressbalance::SSAFEM | |
SSAFEM (std::shared_ptr< const Grid > g) | |
virtual | ~SSAFEM ()=default |
Public Member Functions inherited from pism::stressbalance::SSA | |
SSA (std::shared_ptr< const Grid > g) | |
virtual | ~SSA () |
virtual void | update (const Inputs &inputs, bool full_update) |
Update the SSA solution. More... | |
void | set_initial_guess (const array::Vector &guess) |
Set the initial guess of the SSA velocity. More... | |
virtual std::string | stdout_report () const |
Produce a report string for the standard output. More... | |
const array::Vector & | driving_stress () const |
Public Member Functions inherited from pism::stressbalance::ShallowStressBalance | |
ShallowStressBalance (std::shared_ptr< const Grid > g) | |
virtual | ~ShallowStressBalance () |
void | init () |
const array::Vector1 & | velocity () const |
Get the thickness-advective 2D velocity. More... | |
const array::Scalar & | basal_frictional_heating () |
Get the basal frictional heating (for the adaptive energy time-stepping). More... | |
void | compute_basal_frictional_heating (const array::Vector &velocity, const array::Scalar &tauc, const array::CellType &mask, array::Scalar &result) const |
Compute the basal frictional heating. More... | |
std::shared_ptr< const rheology::FlowLaw > | flow_law () const |
EnthalpyConverter::Ptr | enthalpy_converter () const |
const IceBasalResistancePlasticLaw * | sliding_law () const |
double | flow_enhancement_factor () const |
Public Member Functions inherited from pism::Component | |
Component (std::shared_ptr< const Grid > grid) | |
virtual | ~Component ()=default |
DiagnosticList | diagnostics () const |
TSDiagnosticList | ts_diagnostics () const |
std::shared_ptr< const Grid > | grid () const |
const Time & | time () const |
const Profiling & | profiling () const |
void | define_model_state (const File &output) const |
Define model state variables in an output file. More... | |
void | write_model_state (const File &output) const |
Write model state variables to an output file. More... | |
MaxTimestep | max_timestep (double t) const |
Reports the maximum time-step the model can take at time t. More... | |
Additional Inherited Members | |
Public Attributes inherited from pism::stressbalance::SSA | |
SSAStrengthExtension * | strength_extension |
Protected Types inherited from pism::Component | |
enum | RegriddingFlag { REGRID_WITHOUT_REGRID_VARS , NO_REGRID_WITHOUT_REGRID_VARS } |
This flag determines whether a variable is read from the -regrid_file file even if it is not listed among variables in -regrid_vars . More... | |
Protected Member Functions inherited from pism::stressbalance::SSAFEM | |
virtual void | init_impl () |
Initialize a generic regular-grid SSA solver. More... | |
void | cache_inputs (const Inputs &inputs) |
Initialize stored data from the coefficients in the SSA. Called by SSAFEM::solve. More... | |
void | quad_point_values (const fem::Element &Q, const Coefficients *x, int *mask, double *thickness, double *tauc, double *hardness) const |
Compute quadrature point values of various coefficients given a quadrature Q and nodal values. More... | |
void | explicit_driving_stress (const fem::Element &E, const Coefficients *x, Vector2d *driving_stress) const |
void | driving_stress (const fem::Element &E, const Coefficients *x, Vector2d *driving_stress) const |
void | PointwiseNuHAndBeta (double thickness, double hardness, int mask, double tauc, const Vector2d &U, const Vector2d &U_x, const Vector2d &U_y, double *nuH, double *dnuH, double *beta, double *dbeta) |
Compute the "(regularized effective viscosity) x (ice thickness)" and effective viscous bed strength from the current solution, at a single quadrature point. More... | |
void | compute_local_function (Vector2d const *const *const velocity, Vector2d **residual) |
Implements the callback for computing the residual. More... | |
void | compute_local_jacobian (Vector2d const *const *const velocity, Mat J) |
Implements the callback for computing the Jacobian. More... | |
virtual void | solve (const Inputs &inputs) |
std::shared_ptr< TerminationReason > | solve_with_reason (const Inputs &inputs) |
std::shared_ptr< TerminationReason > | solve_nocache () |
Protected Member Functions inherited from pism::stressbalance::SSA | |
virtual void | define_model_state_impl (const File &output) const |
The default (empty implementation). More... | |
virtual void | write_model_state_impl (const File &output) const |
The default (empty implementation). More... | |
virtual DiagnosticList | diagnostics_impl () const |
virtual void | compute_driving_stress (const array::Scalar &ice_thickness, const array::Scalar1 &surface_elevation, const array::CellType1 &cell_type, const array::Scalar1 *no_model_mask, array::Vector &result) const |
Compute the gravitational driving stress. More... | |
void | extrapolate_velocity (const array::CellType1 &cell_type, array::Vector1 &velocity) const |
Protected Member Functions inherited from pism::Component | |
virtual MaxTimestep | max_timestep_impl (double t) const |
virtual TSDiagnosticList | ts_diagnostics_impl () const |
void | regrid (const std::string &module_name, array::Array &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS) |
Implements the forward problem of the map taking \(\tau_c\) to the corresponding solution of the SSA.
The class SSAFEM solves the SSA, and the solution depends on a large number of parameters. Considering all of these to be fixed except for \(\tau_c\), we obtain a map \(F_{\rm SSA}\) taking \(\tau_c\) to the corresponding solution \(u_{\rm SSA}\) of the SSA. This is a forward problem which we would like to invert; given \(u_{\rm SSA}\) determine \(\tau_c\) such that \(F_{\rm SSA}(\tau_c) = u_{\rm SSA}\). The forward problem actually implemented by IP_SSATaucForwardProblem is a mild variation \(F_{\rm SSA}\).
First, given the constraint \(\tau_c\ge 0\) it can be helpful to parameterize \(\tau_c\) by some other parameter \(\zeta\),
\[ \tau_c = g(\zeta), \]
where the function \(g\) is non-negative. The function \(g\) is specified by an instance of IPDesignVariableParameterization.
Second, there may be locations where the value of \(\tau_c\) (and hence \(\zeta\)) is known a-priori, and should not be adjusted. Let \(\pi\) be the map that replaces the values of \(\zeta\) with known values at these locations.
IP_SSATaucForwardProblem implements the forward problem
\[ F(\zeta) = F_{\rm SSA}(g(\pi(\zeta))). \]
Algorithms to solve the inverse problem make use of variations on the linearization of \(F\), which are explained below.
The solution of the SSA in SSAFEM is implemented using SNES to solve a nonlinear residual function of the form
\[ \mathcal{R}_{SSA}(u;\tau_c, \text{other parameters})= \vec 0, \]
usually using Newton's method. Let
\[ \mathcal{R}(u, \zeta) = \mathcal{R}_{SSA}(u; g(\pi(\zeta)), \text{other parameters}). \]
The derivative of \(\mathcal{R}\) with respect to \( u\) is called the state Jacobian, \(J_{\rm State}\). Specifically, if \(u=[U_j]\) then
\[ (J_{\rm State})_{ij} = \frac{d \mathcal{R}_i}{dU_j}. \]
This is exactly the same Jacobian that is computed by SSAFEM for solving the SSA via SNES. The matrix for the design Jacobian can be obtained with assemble_jacobian_state.
The derivative of \(\mathcal{R}\) with respect to \( \zeta\) is called the design Jacobian, \(J_{\rm Design}\). Specifically, if \(\zeta=[Z_k]\) then
\[ (J_{\rm Design})_{ik} = \frac{d \mathcal{R}_i}{dZ_k}. \]
The map \(J_{\rm Design}\) can be applied to a vector \(d\zeta\) with apply_jacobian_design. For inverse methods using adjoints, one also needs to be able to apply the transpose of \(J_{\rm Design}\), which is done using apply_jacobian_design_transpose.
The forward problem map \(F\) solves the implicit equation
\[ \mathcal{R}(F(\zeta), \zeta) = 0. \]
Its linearisation \(DF\) then satisfies
\[ J_{\rm State}\; DF\; d\zeta + J_{\rm Design}\; d\zeta = 0 \]
for any perturbation \(d\zeta\). Hence
\[ DF = -J_{\rm State}^{-1}\circ J_{\rm Design}. \]
This derivative is sometimes called the reduced gradient in the literature. To apply \(DF\) to a perturbation \(d\zeta\), use apply_linearization. Adjoint methods require the transpose of this map; to apply \(DF^t\) to \(du\) use apply_linearization_transpose.
Definition at line 103 of file IP_SSATaucForwardProblem.hh.