20 #include "pism/basalstrength/basal_resistance.hh"
21 #include "pism/util/IceGrid.hh"
22 #include "pism/util/Mask.hh"
23 #include "pism/util/ConfigInterface.hh"
24 #include "pism/util/Vars.hh"
25 #include "pism/util/error_handling.hh"
26 #include "pism/util/pism_utilities.hh"
27 #include "pism/rheology/FlowLaw.hh"
28 #include "pism/geometry/Geometry.hh"
29 #include "pism/stressbalance/StressBalance.hh"
30 #include "pism/util/petscwrappers/DM.hh"
31 #include "pism/util/petscwrappers/Vec.hh"
41 m_dzeta_local(m_grid,
"d_zeta_local",
WITH_GHOSTS, m_stencil_width),
42 m_fixed_design_locations(NULL),
44 m_du_global(m_grid,
"linearization work vector (sans ghosts)",
46 m_du_local(m_grid,
"linearization work vector (with ghosts)",
48 m_hardav(m_grid,
"hardav",
WITH_GHOSTS, m_stencil_width),
49 m_element_index(*m_grid),
50 m_element(*m_grid, fem::Q1Quadrature4()),
51 m_rebuild_J_state(true)
60 ierr = DMSetMatType(*
m_da, MATBAIJ);
69 double ksp_rtol = 1e-12;
70 ierr = KSPSetTolerances(
m_ksp, ksp_rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
74 ierr = KSPGetPC(
m_ksp, &pc);
77 ierr = PCSetType(pc, PCBJACOBI);
80 ierr = KSPSetFromOptions(
m_ksp);
107 inputs.
bc_mask =
m_grid->variables().get_2d_mask(
"vel_bc_mask");
136 const int i = p.i(), j = p.j();
244 dzeta_local = &dzeta;
249 list.
add(*dzeta_local);
253 const int i = p.i(), j = p.j();
265 Vector2 U[Nq_max], U_x[Nq_max], U_y[Nq_max];
289 for (
int j =ys; j<ys+ym; j++) {
290 for (
int i =xs; i<xs+xm; i++) {
293 for (
unsigned int k=0;
k<Nk;
k++) {
318 for (
unsigned int k=0;
k<Nk;
k++) {
324 double thickness[Nq_max];
329 double hardness[Nq_max];
334 mask, thickness, tauc, hardness);
337 for (
unsigned int q = 0; q < Nq; q++) {
339 double Duqq[3] = {U_x[q].
u, U_y[q].
v, 0.5 * (U_y[q].
u + U_x[q].
v)};
346 d_nuH *= (2.0 * thickness[q]);
351 for (
unsigned int k = 0;
k < Nk;
k++) {
353 du_e[
k].
u += W*d_nuH*(testqk.
dx*(2*Duqq[0] + Duqq[1]) + testqk.
dy*Duqq[2]);
354 du_e[
k].
v += W*d_nuH*(testqk.
dy*(2*Duqq[1] + Duqq[0]) + testqk.
dx*Duqq[2]);
433 Vector2 U[Nq_max], U_x[Nq_max], U_y[Nq_max];
452 const int i = p.i(), j = p.j();
465 for (
int j = ys; j < ys + ym; j++) {
466 for (
int i = xs; i < xs + xm; i++) {
485 for (
unsigned int k = 0;
k < Nk;
k++) {
489 double thickness[Nq_max];
494 double hardness[Nq_max];
499 mask, thickness, tauc, hardness);
502 for (
unsigned int q = 0; q < Nq; q++) {
504 double Duqq[3] = {U_x[q].
u, U_y[q].
v, 0.5 * (U_y[q].
u + U_x[q].
v)};
512 d_nuH_dB *= (2.0 * thickness[q]);
517 for (
unsigned int k = 0;
k < Nk;
k++) {
518 dzeta_e[
k] += W*d_nuH_dB*
m_element.
chi(q,
k).
val*((du_dx_q[q].
u*(2*Duqq[0] + Duqq[1]) +
519 du_dy_q[q].u*Duqq[2]) +
520 (du_dy_q[q].
v*(2*Duqq[1] + Duqq[0]) +
521 du_dx_q[q].v*Duqq[2]));
534 const int i = p.i(), j = p.j();
538 dzeta_a[j][i] *= dB_dzeta;
579 KSPConvergedReason reason;
580 ierr = KSPGetConvergedReason(
m_ksp, &reason);
581 PISM_CHK(ierr,
"KSPGetConvergedReason");
585 " failed to converge (KSP reason %s)",
586 KSPConvergedReasons[reason]);
589 "IP_SSAHardavForwardProblem::apply_linearization converged"
590 " (KSP reason %s)\n",
591 KSPConvergedReasons[reason]);
645 KSPConvergedReason reason;
646 ierr = KSPGetConvergedReason(
m_ksp, &reason);
647 PISM_CHK(ierr,
"KSPGetConvergedReason");
651 KSPConvergedReasons[reason]);
654 "IP_SSAHardavForwardProblem::apply_linearization converged (KSP reason %s)\n",
655 KSPConvergedReasons[reason]);
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
const Config::ConstPtr m_config
configuration database used by this component
const Logger::ConstPtr m_log
logger (for easy access)
const IceGrid::ConstPtr m_grid
grid used by this component
void ensure_consistency(double ice_free_thickness_threshold)
IceModelVec2S bed_elevation
IceModelVec2S ice_area_specific_volume
IceModelVec2S sea_level_elevation
IceModelVec2S ice_thickness
std::shared_ptr< const IceGrid > ConstPtr
A simple class "hiding" the fact that the mask is stored as floating-point scalars (instead of intege...
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.
unsigned int stencil_width() const
Get the stencil width of the current IceModelVec. Returns 0 if ghosts are not available.
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
void set(double c)
Result: v[j] <- c for all j.
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
void failed()
Indicates a failure of a parallel section.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
std::shared_ptr< TerminationReason > Ptr
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
void enforce_homogeneous(const Element2 &element, double *x_e)
void fix_residual_homogeneous(double **r_global)
void enforce(const Element2 &element, Vector2 *x_e)
void enforce_homogeneous(const Element2 &element, Vector2 *x_e)
void fix_residual_homogeneous(Vector2 **r)
void constrain(Element2 &element)
Constrain element, i.e. ensure that quadratures do not contribute to Dirichlet nodes by marking corre...
void evaluate(const T *x, T *vals, T *dx, T *dy)
Given nodal values, compute the values and partial derivatives at the quadrature points.
void reset(int i, int j)
Initialize the Element to element (i, j) for the purposes of inserting into global residual and Jacob...
void add_contribution(const T *local, T **y_global) const
Add the values of element-local contributions y to the global vector y_global.
void nodal_values(const IceModelVec2Int &x_global, int *result) const
Get nodal values of an integer mask.
int xm
total number of elements to loop over in the x-direction.
int ym
total number of elements to loop over in the y-direction.
int ys
y-coordinate of the first element to loop over.
int xs
x-coordinate of the first element to loop over.
const Germ & chi(unsigned int q, unsigned int k) const
int n_pts() const
Number of quadrature points.
double weight(unsigned int q) const
Weight of the quadrature point q
virtual void convertToDesignVariable(IceModelVec2S &zeta, IceModelVec2S &d, bool communicate=true)
Transforms a vector of values to a vector of values.
virtual void toDesignVariable(double zeta, double *value, double *derivative)=0
Converts from parameterization value to .
petsc::Mat m_J_state
Mat used in apply_linearization and apply_linearization_transpose.
IceModelVec2V::Ptr m_velocity_shared
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...
IPDesignVariableParameterization & m_design_param
The function taking to .
virtual TerminationReason::Ptr linearize_at(IceModelVec2S &zeta)
Sets the current value of the design variable and solves the SSA to find the associated .
IceModelVec2S * m_zeta
Current value of zeta, provided from caller.
virtual void set_design(IceModelVec2S &zeta)
Sets the current value of of the design paramter .
virtual void apply_jacobian_design(IceModelVec2V &u, IceModelVec2S &dzeta, IceModelVec2V &du)
Applies the design Jacobian matrix to a perturbation of the design variable.
virtual void assemble_residual(IceModelVec2V &u, IceModelVec2V &R)
Computes the residual function as defined in the class-level documentation.
IceModelVec2Int * m_fixed_design_locations
Locations where should not be adjusted.
IceModelVec2V m_du_global
Temporary storage when state vectors need to be used without ghosts.
fem::ElementIterator m_element_index
virtual void assemble_jacobian_state(IceModelVec2V &u, Mat J)
Assembles the state Jacobian matrix.
bool m_rebuild_J_state
Flag indicating that the state jacobian matrix needs rebuilding.
IP_SSAHardavForwardProblem(IceGrid::ConstPtr g, IPDesignVariableParameterization &tp)
Constructs from the same objects as SSAFEM, plus a specification of how is parameterized.
IceModelVec2S m_dzeta_local
Storage for d_zeta with ghosts, if needed when an argument d_zeta is ghost-less.
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...
IceModelVec2V m_du_local
Temporary storage when state vectors need to be used with ghosts.
IceModelVec2S m_hardav
Vertically-averaged ice hardness.
fem::Q1Element2 m_element
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.
petsc::KSP m_ksp
KSP used in apply_linearization and apply_linearization_transpose.
std::shared_ptr< Wrapper > Ptr
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.
void cache_inputs(const Inputs &inputs)
Initialize stored data from the coefficients in the SSA. Called by SSAFEM::solve.
IceModelVec2< Coefficients > m_coefficients
void compute_local_function(Vector2 const *const *const velocity, Vector2 **residual)
Implements the callback for computing the residual.
IceModelVec2Int m_bc_mask
TerminationReason::Ptr solve_nocache()
void compute_local_jacobian(Vector2 const *const *const velocity, Mat J)
Implements the callback for computing the Jacobian.
IceModelVec2V m_bc_values
double get_min_thickness() const
Returns minimum thickness to trigger use of extension.
SSAStrengthExtension * strength_extension
std::shared_ptr< petsc::DM > m_da
std::shared_ptr< rheology::FlowLaw > m_flow_law
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
const unsigned int MAX_QUADRATURE_SIZE
static double secondInvariant_2D(const Vector2 &U_x, const Vector2 &U_y)
double dy
Function derivative with respect to y.
double val
Function value.
double dx
Function deriviative with respect to x.
Struct for gathering the value and derivative of a function at a point.