19#include "pism/util/Grid.hh"
20#include "pism/stressbalance/ssa/SSAFEM.hh"
21#include "pism/util/fem/Quadrature.hh"
22#include "pism/util/fem/DirichletData.hh"
23#include "pism/util/Mask.hh"
24#include "pism/basalstrength/basal_resistance.hh"
25#include "pism/rheology/FlowLaw.hh"
26#include "pism/util/pism_options.hh"
27#include "pism/util/error_handling.hh"
28#include "pism/util/Vars.hh"
29#include "pism/stressbalance/StressBalance.hh"
30#include "pism/geometry/Geometry.hh"
32#include "pism/util/node_types.hh"
33#include "pism/util/pism_utilities.hh"
34#include "pism/util/Interpolation1D.hh"
35#include "pism/util/petscwrappers/DM.hh"
36#include "pism/util/petscwrappers/Vec.hh"
37#include "pism/util/petscwrappers/Viewer.hh"
40namespace stressbalance {
49 m_bc_mask(grid,
"bc_mask"),
50 m_bc_values(grid,
"_bc"),
52 m_coefficients(grid,
"ssa_coefficients", array::WITH_GHOSTS, 1),
53 m_node_type(m_grid,
"node_type"),
54 m_boundary_integral(m_grid,
"boundary_integral"),
55 m_element_index(*grid),
56 m_q1_element(*grid, fem::Q1Quadrature4()) {
60 const double ice_density =
m_config->get_number(
"constants.ice.density");
61 m_alpha = 1 - ice_density /
m_config->get_number(
"constants.sea_water.density");
62 m_rho_g = ice_density *
m_config->get_number(
"constants.standard_gravity");
81 ierr = DMDASNESSetFunctionLocal(*dm, INSERT_VALUES,
82#
if PETSC_VERSION_LT(3,21,0)
88 PISM_CHK(ierr,
"DMDASNESSetFunctionLocal");
90 ierr = DMDASNESSetJacobianLocal(*dm,
91#
if PETSC_VERSION_LT(3,21,0)
97 PISM_CHK(ierr,
"DMDASNESSetJacobianLocal");
99 ierr = DMSetMatType(*dm,
"baij");
103 PISM_CHK(ierr,
"DMSetApplicationContext");
105 ierr = SNESSetDM(
m_snes, *dm);
109 int snes_max_it = 200;
110 ierr = SNESSetTolerances(
m_snes, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, snes_max_it,
112 PISM_CHK(ierr,
"SNESSetTolerances");
114 ierr = SNESSetFromOptions(
m_snes);
115 PISM_CHK(ierr,
"SNESSetFromOptions");
118 "node types: interior, boundary, exterior");
123 "residual contribution from lateral boundaries");
132 if (
m_grid->variables().is_available(
"ssa_driving_stress_x") and
133 m_grid->variables().is_available(
"ssa_driving_stress_y")) {
138 m_log->message(2,
" [using the SNES-based finite element method implementation]\n");
144 "Enforce Dirichlet conditions with this additional scaling",
173 if (reason->failed()) {
175 "SSAFEM solve failed to converge (SNES reason %s)",
176 reason->description().c_str());
179 if (
m_log->get_threshold() > 2) {
180 m_stdout_ssa +=
"SSAFEM converged (SNES reason " + reason->description() +
")";
202 ierr = PetscViewerASCIIOpen(
m_grid->com, filename->c_str(), viewer.
rawptr());
203 PISM_CHK(ierr,
"PetscViewerASCIIOpen");
205 ierr = PetscViewerASCIIPrintf(viewer,
"SNES before SSASolve_FE\n");
206 PISM_CHK(ierr,
"PetscViewerASCIIPrintf");
208 ierr = SNESView(
m_snes, viewer);
211 ierr = PetscViewerASCIIPrintf(viewer,
"solution vector before SSASolve_FE\n");
212 PISM_CHK(ierr,
"PetscViewerASCIIPrintf");
219 if (
m_log->get_threshold() >= 2) {
228 SNESConvergedReason snes_reason;
229 ierr = SNESGetConvergedReason(
m_snes, &snes_reason);
PISM_CHK(ierr,
"SNESGetConvergedReason");
232 if (not reason->failed()) {
237 bool view_solution =
options::Bool(
"-ssa_view_solution",
"view solution of the SSA system");
240 ierr = PetscViewerASCIIOpen(
m_grid->com, filename->c_str(), viewer.
rawptr());
241 PISM_CHK(ierr,
"PetscViewerASCIIOpen");
243 ierr = PetscViewerASCIIPrintf(viewer,
"solution vector after SSASolve\n");
244 PISM_CHK(ierr,
"PetscViewerASCIIPrintf");
278 const std::vector<double> &z =
m_grid->z();
288 if (use_explicit_driving_stress) {
294 for (
auto p =
m_grid->points(); p; p.next()) {
295 const int i = p.i(), j = p.j();
300 if (use_explicit_driving_stress) {
301 tau_d.
u = (*m_driving_stress_x)(i, j);
302 tau_d.
v = (*m_driving_stress_y)(i, j);
310 m_grid->kBelowHeight(thickness),
330 const bool use_cfbc =
m_config->get_flag(
"stress_balance.calving_front_stress_bc");
334 m_config->get_number(
"stress_balance.ice_free_thickness_standard"),
350 double *hardness)
const {
352 unsigned int n = E.
n_pts();
354 for (
unsigned int q = 0; q <
n; q++) {
363 for (
int k = 0;
k < E.
n_chi();
k++) {
373 mask[q] =
m_gc.
mask(sea_level, bed, thickness[q]);
382 const unsigned int n = E.
n_pts();
384 for (
unsigned int q = 0; q <
n; q++) {
387 for (
int k = 0;
k < E.
n_chi();
k++) {
442 const unsigned int n = E.
n_pts();
444 for (
unsigned int q = 0; q <
n; q++) {
456 for (
int k = 0;
k < E.
n_chi();
k++) {
470 const int M =
m_gc.
mask(sea_level, b, H);
475 h_x = grounded ? b_x + H_x :
m_alpha * H_x,
476 h_y = grounded ? b_y + H_y :
m_alpha * H_y;
478 result[q].
u = - pressure * h_x;
479 result[q].
v = - pressure * h_y;
510 double *nuH,
double *dnuH,
511 double *beta,
double *dbeta) {
513 if (thickness < strength_extension->get_min_thickness()) {
515 if (dnuH !=
nullptr) {
524 if (dnuH !=
nullptr) {
538 if (dbeta !=
nullptr) {
558 const unsigned int Nq = 2;
559 double chi_b[Nq][Nq];
564 for (
int i = 0; i < 2; ++i) {
565 for (
int j = 0; j < 2; ++j) {
575 P1Element2 p1_element[Nk] = {P1Element2(*
m_grid, Q_p1, 0),
576 P1Element2(*
m_grid, Q_p1, 1),
577 P1Element2(*
m_grid, Q_p1, 2),
578 P1Element2(*
m_grid, Q_p1, 3)};
583 use_cfbc =
m_config->get_flag(
"stress_balance.calving_front_stress_bc");
586 ice_density =
m_config->get_number(
"constants.ice.density"),
587 ocean_density =
m_config->get_number(
"constants.sea_water.density"),
588 standard_gravity =
m_config->get_number(
"constants.standard_gravity");
613 for (
int j = ys; j < ys + ym; j++) {
614 for (
int i = xs; i < xs + xm; i++) {
632 E = &p1_element[type];
642 std::vector<Vector2d> I(Nk);
654 double psi[2] = {0.0, 0.0};
656 const unsigned int n_sides = E->
n_sides();
658 for (
unsigned int s = 0; s < n_sides; ++s) {
666 n1 = (n0 + 1) % n_sides;
674 for (
unsigned int q = 0; q < Nq; ++q) {
680 psi[0] = chi_b[0][q];
681 psi[1] = chi_b[1][q];
686 H = H_nodal[n0] * psi[0] + H_nodal[n1] * psi[1],
687 bed = b_nodal[n0] * psi[0] + b_nodal[n1] * psi[1],
688 sea_level = sl_nodal[n0] * psi[0] + sl_nodal[n1] * psi[1];
692 P_ice = 0.5 * ice_density * standard_gravity * H,
694 ice_density, ocean_density,
698 double dP = H * (P_ice - P_water);
710 I[n0] += W * (- psi[0] * dP) * E->
normal(s);
711 I[n1] += W * (- psi[1] * dP) * E->
normal(s);
740 const bool use_cfbc =
m_config->get_flag(
"stress_balance.calving_front_stress_bc");
747 P1Element2 p1_element[Nk] = {P1Element2(*
m_grid, Q_p1, 0),
748 P1Element2(*
m_grid, Q_p1, 1),
749 P1Element2(*
m_grid, Q_p1, 2),
750 P1Element2(*
m_grid, Q_p1, 3)};
757 for (
auto p =
m_grid->points(); p; p.next()) {
758 const int i = p.i(), j = p.j();
767 Vector2d U[Nq_max], U_x[Nq_max], U_y[Nq_max];
778 for (
int j = ys; j < ys + ym; j++) {
779 for (
int i = xs; i < xs + xm; i++) {
799 E = &p1_element[type];
810 const unsigned int Nq = E->
n_pts();
816 double thickness[Nq_max];
818 double hardness[Nq_max];
827 if (use_explicit_driving_stress) {
840 if (dirichlet_data) {
843 dirichlet_data.
enforce(*E, velocity_nodal);
855 for (
unsigned int k = 0;
k < Nk;
k++) {
860 for (
unsigned int q = 0; q < Nq; q++) {
864 double eta = 0.0, beta = 0.0;
866 U[q], U_x[q], U_y[q],
867 &eta, NULL, &beta, NULL);
870 const Vector2d tau_b = U[q] * (- beta);
875 u_y_plus_v_x = U_y[q].
u + U_x[q].
v;
878 for (
int k = 0;
k < E->
n_chi();
k++) {
881 residual[
k].
u += W * (eta * (psi.
dx * (4.0 * u_x + 2.0 * v_y) + psi.
dy * u_y_plus_v_x)
882 - psi.
val * (tau_b.
u + tau_d[q].
u));
883 residual[
k].
v += W * (eta * (psi.
dx * u_y_plus_v_x + psi.
dy * (2.0 * u_x + 4.0 * v_y))
884 - psi.
val * (tau_b.
v + tau_d[q].
v));
898 if (dirichlet_data) {
899 dirichlet_data.
fix_residual(velocity_global, residual_global);
914 Vector2d const *
const *
const residual_global) {
916 bool monitorFunction =
options::Bool(
"-ssa_monitor_function",
"monitor the SSA residual");
917 if (not monitorFunction) {
921 ierr = PetscPrintf(
m_grid->com,
922 "SSA Solution and Function values (pointwise residuals)\n");
927 for (
auto p =
m_grid->points(); p; p.next()) {
928 const int i = p.i(), j = p.j();
930 ierr = PetscSynchronizedPrintf(
m_grid->com,
931 "[%2d, %2d] u=(%12.10e, %12.10e) f=(%12.4e, %12.4e)\n",
933 velocity_global[j][i].
u, velocity_global[j][i].
v,
934 residual_global[j][i].
u, residual_global[j][i].
v);
935 PISM_CHK(ierr,
"PetscSynchronizedPrintf");
942 ierr = PetscSynchronizedFlush(
m_grid->com, NULL);
943 PISM_CHK(ierr,
"PetscSynchronizedFlush");
959 const bool use_cfbc =
m_config->get_flag(
"stress_balance.calving_front_stress_bc");
966 P1Element2 p1_element[Nk] = {P1Element2(*
m_grid, Q_p1, 0),
967 P1Element2(*
m_grid, Q_p1, 1),
968 P1Element2(*
m_grid, Q_p1, 2),
969 P1Element2(*
m_grid, Q_p1, 3)};
972 PetscErrorCode ierr = MatZeroEntries(Jac);
981 Vector2d U[Nq_max], U_x[Nq_max], U_y[Nq_max];
992 for (
int j = ys; j < ys + ym; j++) {
993 for (
int i = xs; i < xs + xm; i++) {
1013 E = &p1_element[type];
1029 double thickness[Nq_max];
1030 double tauc[Nq_max];
1031 double hardness[Nq_max];
1038 mask, thickness, tauc, hardness);
1048 if (dirichlet_data) {
1049 dirichlet_data.
enforce(*E, velocity_nodal);
1053 E->
evaluate(velocity_nodal, U, U_x, U_y);
1059 double K[2*Nk][2*Nk];
1061 ierr = PetscMemzero(K,
sizeof(K));
1064 for (
unsigned int q = 0; q < Nq; q++) {
1072 u_y_plus_v_x = U_y[q].
u + U_x[q].
v;
1074 double eta = 0.0, deta = 0.0, beta = 0.0, dbeta = 0.0;
1076 U[q], U_x[q], U_y[q],
1077 &eta, &deta, &beta, &dbeta);
1079 for (
unsigned int l = 0; l < n_chi; l++) {
1086 gamma_u = (2.0 * u_x + v_y) *
phi.dx + 0.5 * u_y_plus_v_x *
phi.dy,
1087 gamma_v = 0.5 * u_y_plus_v_x *
phi.dx + (u_x + 2.0 * v_y) *
phi.dy;
1091 eta_u = deta * gamma_u,
1092 eta_v = deta * gamma_v;
1096 taub_xu = -dbeta * u * u *
phi.val - beta *
phi.val,
1097 taub_xv = -dbeta * u * v *
phi.val,
1098 taub_yu = -dbeta * v * u *
phi.val,
1099 taub_yv = -dbeta * v * v *
phi.val - beta *
phi.val;
1101 for (
unsigned int k = 0;
k < n_chi;
k++) {
1108 ierr = PetscPrintf(PETSC_COMM_SELF,
"eta=0 i %d j %d q %d k %d\n", i, j, q,
k);
1113 K[
k*2 + 0][l*2 + 0] += W * (eta_u * (psi.
dx * (4 * u_x + 2 * v_y) + psi.
dy * u_y_plus_v_x)
1114 + eta * (4 * psi.
dx *
phi.dx + psi.
dy *
phi.dy) - psi.
val * taub_xu);
1116 K[
k*2 + 0][l*2 + 1] += W * (eta_v * (psi.
dx * (4 * u_x + 2 * v_y) + psi.
dy * u_y_plus_v_x)
1117 + eta * (2 * psi.
dx *
phi.dy + psi.
dy *
phi.dx) - psi.
val * taub_xv);
1119 K[
k*2 + 1][l*2 + 0] += W * (eta_u * (psi.
dx * u_y_plus_v_x + psi.
dy * (2 * u_x + 4 * v_y))
1120 + eta * (psi.
dx *
phi.dy + 2 * psi.
dy *
phi.dx) - psi.
val * taub_yu);
1122 K[
k*2 + 1][l*2 + 1] += W * (eta_v * (psi.
dx * u_y_plus_v_x + psi.
dy * (2 * u_x + 4 * v_y))
1123 + eta * (psi.
dx *
phi.dx + 4 * psi.
dy *
phi.dy) - psi.
val * taub_yv);
1142 if (dirichlet_data) {
1154 ierr = MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
1155 PISM_CHK(ierr,
"MatAssemblyBegin");
1157 ierr = MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
1160 ierr = MatSetOption(Jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
1163 ierr = MatSetOption(Jac, MAT_SYMMETRIC, PETSC_TRUE);
1170 PetscErrorCode ierr;
1171 bool mon_jac =
options::Bool(
"-ssa_monitor_jacobian",
"monitor the SSA Jacobian");
1180 ierr = SNESGetIterationNumber(
m_snes, &iter);
1181 PISM_CHK(ierr,
"SNESGetIterationNumber");
1183 auto file_name =
pism::printf(
"PISM_SSAFEM_J%d.m", (
int)iter);
1186 "writing Matlab-readable file for SSAFEM system A xsoln = rhs to file `%s' ...\n",
1191 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);
1192 PISM_CHK(ierr,
"PetscViewerSetType");
1194 ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
1195 PISM_CHK(ierr,
"PetscViewerPushFormat");
1197 ierr = PetscViewerFileSetName(viewer, file_name.c_str());
1198 PISM_CHK(ierr,
"PetscViewerFileSetName");
1200 ierr = PetscObjectSetName((PetscObject) Jac,
"A");
1201 PISM_CHK(ierr,
"PetscObjectSetName");
1203 ierr = MatView(Jac, viewer);
1206 ierr = PetscViewerPopFormat(viewer);
1207 PISM_CHK(ierr,
"PetscViewerPopFormat");
1212 Vector2d const *
const *
const velocity,
1219 MPI_Comm com = MPI_COMM_SELF;
1220 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)fe->
da, &com); CHKERRQ(ierr);
1222 SETERRQ(com, 1,
"A PISM callback failed");
1228 Vector2d const *
const *
const velocity,
1235 MPI_Comm com = MPI_COMM_SELF;
1236 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)fe->
da, &com); CHKERRQ(ierr);
1238 SETERRQ(com, 1,
"A PISM callback failed");
const units::System::Ptr m_sys
unit system used by this component
const Config::ConstPtr m_config
configuration database used by this component
const Logger::ConstPtr m_log
logger (for easy access)
const std::shared_ptr< const Grid > m_grid
grid used by this component
int mask(double sea_level, double bed, double thickness) const
array::Scalar1 sea_level_elevation
array::Scalar2 ice_thickness
array::Scalar2 bed_elevation
virtual void drag_with_derivative(double tauc, double vx, double vy, double *drag, double *ddrag) const
Compute the drag coefficient and its derivative with respect to .
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
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
void copy_from(const Array2D< T > &source)
void add(double alpha, const Array2D< T > &x)
double * get_column(int i, int j)
void set_interpolation_type(InterpolationType type)
void set(double c)
Result: v[j] <- c for all j.
std::shared_ptr< petsc::DM > dm() const
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
void fix_residual(Vector2d const *const *x_global, Vector2d **r_global)
void enforce(const Element2 &element, Vector2d *x_e)
void fix_residual_homogeneous(Vector2d **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...
double side_length(unsigned int side) const
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 array::Scalar &x_global, int *result) const
Get nodal values of an integer mask.
unsigned int n_sides() const
Vector2d normal(unsigned int side) const
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.
double weight(unsigned int q) const
Weight of the quadrature point q
const Germ & chi(unsigned int q, unsigned int k) const
unsigned int n_pts() const
Number of quadrature points.
The mapping from global to local degrees of freedom.
P1 element embedded in Q1Element2.
double weight(int k) const
QuadPoint point(int k) const
std::shared_ptr< TerminationReason > solve_with_reason(const Inputs &inputs)
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 ...
void driving_stress(const fem::Element &E, const Coefficients *x, Vector2d *result) const
const array::Scalar * m_driving_stress_x
void cache_inputs(const Inputs &inputs)
Initialize stored data from the coefficients in the SSA. Called by SSAFEM::solve.
static PetscErrorCode function_callback(DMDALocalInfo *info, Vector2d const *const *velocity, Vector2d **residual, CallbackData *fe)
SNES callbacks.
double m_beta_ice_free_bedrock
SSAFEM(std::shared_ptr< const Grid > g)
static void explicit_driving_stress(const fem::Element &E, const Coefficients *x, Vector2d *result)
CallbackData m_callback_data
void compute_local_jacobian(Vector2d const *const *velocity, Mat J)
Implements the callback for computing the Jacobian.
array::Vector1 m_boundary_integral
Boundary integral (CFBC contribution to the residual).
static PetscErrorCode jacobian_callback(DMDALocalInfo *info, Vector2d const *const *velocity, Mat A, Mat J, CallbackData *fe)
void monitor_jacobian(Mat Jac)
void monitor_function(Vector2d const *const *velocity_global, Vector2d const *const *residual_global)
virtual void solve(const Inputs &inputs)
fem::Q1Element2 m_q1_element
virtual void init_impl()
Initialize a generic regular-grid SSA solver.
std::shared_ptr< TerminationReason > solve_nocache()
void cache_residual_cfbc(const Inputs &inputs)
Compute and cache residual contributions from the integral over the lateral boundary.
void compute_local_function(Vector2d const *const *velocity, Vector2d **residual)
Implements the callback for computing the residual.
fem::ElementIterator m_element_index
void quad_point_values(const fem::Element &E, 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.
const array::Scalar * m_driving_stress_y
array::Scalar1 m_node_type
Storage for node types (interior, boundary, exterior).
array::Vector1 m_bc_values
array::Array2D< Coefficients > m_coefficients
double get_notional_strength() const
Returns strength = (viscosity times thickness).
array::Vector m_velocity_global
SSAStrengthExtension * strength_extension
virtual void init_impl()
Initialize a generic regular-grid SSA solver.
const array::Vector1 & velocity() const
Get the thickness-advective 2D velocity.
std::shared_ptr< rheology::FlowLaw > m_flow_law
IceBasalResistancePlasticLaw * m_basal_sliding_law
array::Vector2 m_velocity
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
Germ chi(unsigned int k, const QuadPoint &pt)
Linear basis functions on the interval [-1, -1].
const unsigned int MAX_QUADRATURE_SIZE
ElementType element_type(const int node_type[q1::n_chi])
bool ice_free_land(int M)
bool grounded(int M)
Grounded cell (grounded ice or ice-free).
bool ocean(int M)
An ocean cell (floating ice or ice-free).
bool Bool(const std::string &option, const std::string &description)
double averaged_hardness(const FlowLaw &ice, double thickness, unsigned int kbelowH, const double *zlevels, const double *enthalpy)
Computes vertical average of B(E, p) ice hardness, namely .
static double secondInvariant_2D(const Vector2d &U_x, const Vector2d &U_y)
void compute_node_types(const array::Scalar1 &ice_thickness, double thickness_threshold, array::Scalar &result)
double average_water_column_pressure(double ice_thickness, double bed, double floatation_level, double rho_ice, double rho_water, double g)
std::string printf(const char *format,...)
void handle_fatal_errors(MPI_Comm com)
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.
Adaptor for gluing SNESDAFormFunction callbacks to an SSAFEM.
double thickness
ice thickness
double tauc
basal yield stress
double sea_level
sea level
Vector2d driving_stress
prescribed gravitational driving stress
double hardness
ice hardness