PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
Public Types | Public Member Functions | Protected Attributes | List of all members
pism::inverse::IP_SSATaucForwardProblem Class Reference

Implements the forward problem of the map taking \(\tau_c\) to the corresponding solution of the SSA. More...

#include <IP_SSATaucForwardProblem.hh>

+ Inheritance diagram for pism::inverse::IP_SSATaucForwardProblem:

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::Vectorsolution ()
 Returns the last solution of the SSA as computed by linearize_at. More...
 
virtual IPDesignVariableParameterizationtauc_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< TerminationReasonlinearize_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::DMget_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::Vectordriving_stress () const
 
- Public Member Functions inherited from pism::stressbalance::ShallowStressBalance
 ShallowStressBalance (std::shared_ptr< const Grid > g)
 
virtual ~ShallowStressBalance ()
 
void init ()
 
const array::Vector1velocity () const
 Get the thickness-advective 2D velocity. More...
 
const array::Scalarbasal_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::FlowLawflow_law () const
 
EnthalpyConverter::Ptr enthalpy_converter () const
 
const IceBasalResistancePlasticLawsliding_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 Gridgrid () const
 
const Timetime () const
 
const Profilingprofiling () 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...
 

Protected Attributes

array::Scalarm_zeta
 Current value of zeta, provided from caller. More...
 
array::Scalar1 m_dzeta_local
 Storage for d_zeta with ghosts, if needed when an argument d_zeta is ghost-less. More...
 
array::Scalar2 m_tauc_copy
 Storage for tauc (avoids modifying fields obtained via pism::Vars) More...
 
array::Scalarm_fixed_tauc_locations
 Locations where \(\tau_c\) should not be adjusted. More...
 
IPDesignVariableParameterizationm_tauc_param
 The function taking \(\zeta\) to \(\tau_c\). More...
 
std::shared_ptr< array::Vectorm_velocity_shared
 Copy of the velocity field managed using a shared pointer. More...
 
array::Vector m_du_global
 Temporary storage when state vectors need to be used without ghosts. More...
 
array::Vector1 m_du_local
 Temporary storage when state vectors need to be used with ghosts. More...
 
fem::ElementIterator m_element_index
 
fem::Q1Element2 m_element
 
petsc::KSP m_ksp
 KSP used in apply_linearization and apply_linearization_transpose. More...
 
petsc::Mat m_J_state
 Mat used in apply_linearization and apply_linearization_transpose. More...
 
SNESConvergedReason m_reason
 
bool m_rebuild_J_state
 Flag indicating that the state jacobian matrix needs rebuilding. More...
 
- Protected Attributes inherited from pism::stressbalance::SSAFEM
array::Scalar1 m_bc_mask
 
array::Vector1 m_bc_values
 
GeometryCalculator m_gc
 
double m_alpha
 
double m_rho_g
 
array::Array2D< Coefficientsm_coefficients
 
CallbackData m_callback_data
 
petsc::SNES m_snes
 
array::Scalar1 m_node_type
 Storage for node types (interior, boundary, exterior). More...
 
array::Vector1 m_boundary_integral
 Boundary integral (CFBC contribution to the residual). More...
 
double m_dirichletScale
 
double m_beta_ice_free_bedrock
 
double m_epsilon_ssa
 
fem::ElementIterator m_element_index
 
fem::Q1Element2 m_q1_element
 
const array::Scalarm_driving_stress_x
 
const array::Scalarm_driving_stress_y
 
- Protected Attributes inherited from pism::stressbalance::SSA
array::CellType2 m_mask
 
array::Vector m_taud
 
std::string m_stdout_ssa
 
std::shared_ptr< petsc::DMm_da
 
array::Vector m_velocity_global
 
int m_event_ssa
 
- Protected Attributes inherited from pism::stressbalance::ShallowStressBalance
IceBasalResistancePlasticLawm_basal_sliding_law
 
std::shared_ptr< rheology::FlowLawm_flow_law
 
EnthalpyConverter::Ptr m_EC
 
array::Vector2 m_velocity
 
array::Scalar m_basal_frictional_heating
 
double m_e_factor
 flow enhancement factor More...
 
- Protected Attributes inherited from pism::Component
const std::shared_ptr< const Gridm_grid
 grid used by this component More...
 
const Config::ConstPtr m_config
 configuration database used by this component More...
 
const units::System::Ptr m_sys
 unit system used by this component More...
 
const Logger::ConstPtr m_log
 logger (for easy access) More...
 

Additional Inherited Members

- Public Attributes inherited from pism::stressbalance::SSA
SSAStrengthExtensionstrength_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< TerminationReasonsolve_with_reason (const Inputs &inputs)
 
std::shared_ptr< TerminationReasonsolve_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)
 

Detailed Description

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.


The documentation for this class was generated from the following files: