25#include "pism/stressbalance/blatter/Blatter.hh"
26#include "pism/util/error_handling.hh"
27#include "pism/util/Vector2d.hh"
29#include "pism/stressbalance/blatter/util/DataAccess.hh"
30#include "pism/stressbalance/blatter/util/grid_hierarchy.hh"
31#include "pism/util/node_types.hh"
33#include "pism/rheology/FlowLawFactory.hh"
35#include "pism/geometry/Geometry.hh"
36#include "pism/stressbalance/StressBalance.hh"
37#include "pism/util/array/Array3D.hh"
38#include "pism/util/pism_options.hh"
39#include "pism/util/pism_utilities.hh"
40#include "pism/util/fem/Quadrature.hh"
43namespace stressbalance {
63 int ierr = DMDAGetLocalInfo(
m_da, &info);
PISM_CHK(ierr,
"DMDAGetLocalInfo");
75 for (
int j = info.gys; j < info.gys + info.gym - 1; j++) {
76 for (
int i = info.gxs; i < info.gxs + info.gxm - 1; i++) {
85 if (p[
k].thickness < min_thickness) {
94 node_type(ii, jj) +=
static_cast<double>(interior);
103 for (
int j = info.ys; j < info.ys + info.ym; j++) {
104 for (
int i = info.xs; i < info.xs + info.xm; i++) {
106 switch ((
int)node_type(i, j)) {
150 for (
int n = 0;
n < N; ++
n) {
176 for (
int n = 0;
n < N; ++
n) {
184 return grounded and floating;
203 for (
int n = 0;
n < N; ++
n) {
205 if (z[
k] > sea_level[
k]) {
212 return above and below;
227 const int *node_type,
228 const double *ice_bottom,
229 const double *sea_level) {
236 for (
int n = 0;
n < N; ++
n) {
244 for (
int n = 0;
n < N; ++
n) {
245 if (ice_bottom[nodes[
n]] < sea_level[nodes[
n]]) {
261 m_parameters(grid,
"bp_input_parameters", array::WITH_GHOSTS),
262 m_face4(grid->dx(), grid->dy(), fem::Q1Quadrature4()),
263 m_face100(grid->dx(), grid->dy(), fem::Q1QuadratureN(10))
269 auto pism_da =
grid->get_dm(1, 0);
271 int ierr =
setup(*pism_da,
grid->periodicity(), Mz, coarsening_factor,
"bp_");
274 "Failed to allocate a Blatter solver instance");
278 std::vector<double> sigma(Mz);
279 double dz = 1.0 / (Mz - 1.0);
280 for (
int i = 0; i < Mz; ++i) {
287 .long_name(
"u velocity component on the sigma grid")
292 .long_name(
"v velocity component on the sigma grid")
295 std::map<std::string,std::string> z_attrs =
297 {
"long_name",
"scaled Z-coordinate in the ice (z_base=0, z_surface=1)"},
301 m_u_sigma->metadata(0).z().set_name(
"z_sigma");
302 m_v_sigma->metadata(0).z().set_name(
"z_sigma");
304 for (
const auto &z_attr : z_attrs) {
305 m_u_sigma->metadata(0).z().set_string(z_attr.first, z_attr.second);
306 m_v_sigma->metadata(0).z().set_string(z_attr.first, z_attr.second);
317 double g =
m_config->get_number(
"constants.standard_gravity");
326 double E =
m_config->get_number(
"stress_balance.blatter.enhancement_factor");
329 double softness =
m_config->get_number(
"flow_law.isothermal_Glen.ice_softness"),
346 Mat mrestrict, Vec rscale, Mat inject,
347 DM coarse,
void *ctx) {
354 PetscErrorCode ierr =
restrict_data(fine, coarse,
"3D_DM"); CHKERRQ(ierr);
365 ierr = DMGetOptionsPrefix(dm_fine, &prefix); CHKERRQ(ierr);
370 ierr =
setup_level(dm_coarse, mg_levels); CHKERRQ(ierr);
384 int coarsening_factor,
385 const std::string &prefix) {
387 PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pism_da, &comm); CHKERRQ(ierr);
391 auto option =
pism::printf(
"-%spc_mg_levels", prefix.c_str());
399 int c = coarsening_factor;
404 if (((mz - 1) / c) * c != mz - 1) {
405 int N = std::pow(c, (
int)mg_levels - 1);
406 auto message =
pism::printf(
"Blatter stress balance solver: settings\n"
407 "stress_balance.blatter.Mz = %d,\n"
408 "stress_balance.blatter.coarsening_factor = %d,\n"
409 "and '%s %d' are not compatible.\n"
410 "To use N = %d multigrid levels with the coarsening factor C = %d\n"
411 "stress_balance.blatter.Mz has to be equal to A * C^(N - 1) + 1\n"
412 "for some positive integer A, e.g. %d, %d, %d, ...",
413 Mz, c, option.c_str(), mg_levels, mg_levels, c,
414 N + 1, 2*N + 1, 3*N + 1);
417 mz = (mz - 1) / c + 1;
430 ierr = DMDAGetInfo(pism_da,
445 ierr = DMDAGetOwnershipRanges(pism_da, &lx, &ly, NULL); CHKERRQ(ierr);
453 bx = (periodicity &
grid::X_PERIODIC) != 0 ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE,
454 by = (periodicity &
grid::Y_PERIODIC) != 0 ? DM_BOUNDARY_PERIODIC : DM_BOUNDARY_NONE,
455 bz = DM_BOUNDARY_NONE;
457 ierr = DMDACreate3d(comm,
467 ierr = DMSetOptionsPrefix(
m_da, prefix.c_str()); CHKERRQ(ierr);
470 ierr = DMDASetRefinementFactor(
m_da, coarsening_factor, 1, 1); CHKERRQ(ierr);
472 ierr = DMSetFromOptions(
m_da); CHKERRQ(ierr);
474 ierr = DMSetUp(
m_da); CHKERRQ(ierr);
486 ierr = DMCreateGlobalVector(
m_da,
m_x.
rawptr()); CHKERRQ(ierr);
488 ierr = VecSetOptionsPrefix(
m_x, prefix.c_str()); CHKERRQ(ierr);
490 ierr = VecSetFromOptions(
m_x); CHKERRQ(ierr);
497 ierr = SNESCreate(comm,
m_snes.
rawptr()); CHKERRQ(ierr);
499 ierr = SNESSetOptionsPrefix(
m_snes, prefix.c_str()); CHKERRQ(ierr);
501 ierr = SNESSetDM(
m_snes,
m_da); CHKERRQ(ierr);
503 ierr = DMDASNESSetFunctionLocal(
m_da, INSERT_VALUES,
504#
if PETSC_VERSION_LT(3,21,0)
509 this); CHKERRQ(ierr);
511 ierr = DMDASNESSetJacobianLocal(
m_da,
512#
if PETSC_VERSION_LT(3,21,0)
517 this); CHKERRQ(ierr);
519 ierr = SNESSetFromOptions(
m_snes); CHKERRQ(ierr);
522 PetscBool ksp_use_ew = PETSC_FALSE;
523 ierr = SNESKSPGetUseEW(
m_snes, &ksp_use_ew); CHKERRQ(ierr);
536 ice_density =
m_config->get_number(
"constants.ice.density"),
537 water_density =
m_config->get_number(
"constants.sea_water.density"),
538 alpha = ice_density / water_density;
549 for (
auto p =
m_grid->points(); p; p.next()) {
550 const int i = p.i(), j = p.j();
553 b_grounded = b(i, j),
554 b_floating = sea_level(i, j) - alpha * H(i, j),
555 s_grounded = b(i, j) + H(i, j),
556 s_floating = sea_level(i, j) + (1.0 - alpha) * H(i, j);
561 m_parameters(i, j).bed = std::max(b_grounded, b_floating);
563 m_parameters(i, j).floatation = s_floating - s_grounded;
570 double H_min =
m_config->get_number(
"stress_balance.ice_free_thickness_standard");
582 const auto *enthalpy = inputs.
enthalpy;
584 const auto &zlevels = enthalpy->
levels();
585 auto Mz = zlevels.size();
591 int ierr = DMDAGetLocalInfo(da, &info);
PISM_CHK(ierr,
"DMDAGetLocalInfo");
601 for (
auto p =
m_grid->points(); p; p.next()) {
602 const int i = p.i(), j = p.j();
604 double H = ice_thickness(i, j);
606 const double *E = enthalpy->get_column(i, j);
608 for (
int k = 0;
k < Mz_sigma; ++
k) {
610 z =
grid_z(0.0, H, Mz_sigma,
k),
612 pressure =
m_EC->pressure(depth),
615 auto k0 =
m_grid->kBelowHeight(z);
618 double lambda = (z - zlevels[k0]) / (zlevels[k0 + 1] - zlevels[k0]);
620 E_local = (1.0 - lambda) * E[k0] + lambda * E[k0 + 1];
625 hardness[j][i][
k] =
m_flow_law->hardness(E_local, pressure);
641 double *bottom_elevation,
642 double *ice_thickness,
643 double *surface_elevation,
644 double *sea_level)
const {
649 auto p = P[I.j][I.i];
651 node_type[
n] =
static_cast<int>(p.node_type);
652 bottom_elevation[
n] = p.bed;
653 ice_thickness[
n] = p.thickness;
655 if (surface_elevation !=
nullptr) {
656 surface_elevation[
n] = p.bed + p.thickness;
659 if (sea_level !=
nullptr) {
660 sea_level[
n] = p.sea_level;
666 m_log->message(2,
"* Initializing the Blatter stress balance...\n");
674 unsigned int start = input_file.
nrecords() - 1;
676 if (u_sigma_found and v_sigma_found) {
677 m_log->message(3,
"Reading uvel_sigma and vvel_sigma...\n");
685 "uvel_sigma and vvel_sigma not found");
688 int ierr = VecSet(
m_x, 0.0);
PISM_CHK(ierr,
"VecSet");
705 int ierr = DMDAGetLocalInfo(
m_da, &info);
PISM_CHK(ierr,
"DMDAGetLocalInfo");
712 double R_min = 1e16, R_max = 0.0, R_avg = 0.0;
714 double n_cells = 0.0;
717 for (
int j = info.ys; j < info.ys + info.ym - 1; j++) {
718 for (
int i = info.xs; i < info.xs + info.xm - 1; i++) {
725 for (
int k = 0;
k < 4; ++
k) {
736 for (
int k = 0;
k < 4; ++
k) {
737 dz_max = std::max(P[
k].thickness / (info.mz - 1), dz_max);
739 double R = dz_max / dxy;
741 R_min = std::min(R, R_min);
742 R_max = std::max(R, R_max);
749 R_avg /= std::max(n_cells, 1.0);
755 "Blatter solver: %d * (%d - 1) = %d active elements\n",
756 (
int)n_cells, (
int)info.mz, (
int)(n_cells * (info.mz - 1)));
760 " Vertical spacing (m): min = %3.3f, max = %3.3f, avg = %3.3f\n",
761 R_min * dxy, R_max * dxy, R_avg * dxy);
763 " Aspect ratios: min = %3.3f, max = %3.3f, avg = %3.3f, max/min = %3.3f\n",
764 R_min, R_max, R_avg, R_max / R_min);
779 PISM_CHK(ierr,
"SNESGetConvergedReason");
782 PISM_CHK(ierr,
"SNESGetIterationNumber");
784 ierr = SNESGetLinearSolveIterations(
m_snes, &result.
ksp_it);
785 PISM_CHK(ierr,
"SNESGetLinearSolveIterations");
788 ierr = SNESGetKSP(
m_snes, &ksp);
791 ierr = KSPGetConvergedReason(ksp, &result.
ksp_reason);
792 PISM_CHK(ierr,
"KSPGetConvergedReason");
795 ierr = KSPGetPC(ksp, &pc);
799 ierr = PCGetType(pc, &pc_type);
802 if (std::string(pc_type) == PCMG) {
804 ierr = PCMGGetCoarseSolve(pc, &coarse_ksp);
805 PISM_CHK(ierr,
"PCMGGetCoarseSolve");
808 PISM_CHK(ierr,
"KSPGetIterationNumber");
821 int snes_total_it = 0;
822 int ksp_total_it = 0;
828 schoofLen =
m_config->get_number(
"flow_law.Schoof_regularizing_length",
"m"),
829 schoofVel =
m_config->get_number(
"flow_law.Schoof_regularizing_velocity",
"m second-1"),
831 eps = pow(schoofVel / schoofLen, 2.0),
833 gamma = std::floor(std::log10(eps)),
848 PetscInt snes_max_it = 0;
849 ierr = SNESGetTolerances(
m_snes, NULL, NULL, NULL, &snes_max_it, NULL);
850 PISM_CHK(ierr,
"SNESGetTolerances");
852 double lambda = lambda_min;
853 double delta = delta0;
860 "Blatter solver: Starting parameter continuation with lambda = %f\n",
863 for (
int N = 0; N < Nc; ++N) {
867 m_log->message(2,
"Blatter solver: step %d with lambda = %f, eps = %e\n",
874 ksp_total_it += info.
ksp_it;
877 m_log->message(2,
"Blatter solver continuation step #%d: %s\n"
878 " lambda = %f, eps = %e\n"
879 " SNES: %d, KSP: %d\n",
884 " Coarse MG level KSP (last iteration): %d\n",
894 info.
ksp_it = ksp_total_it;
904 double F = (snes_max_it - info.
snes_it) / (
double)snes_max_it;
905 delta *= 1.0 + A *
F *
F;
908 delta = std::min(delta, delta_max);
911 if (lambda + delta > lambda_max) {
912 delta = lambda_max - lambda;
915 m_log->message(2,
" Advancing lambda from %f to %f (delta = %f)\n",
916 lambda, lambda + delta, delta);
919 }
else if (info.
snes_reason == SNES_DIVERGED_LINE_SEARCH or
929 if (lambda < lambda_min) {
930 m_log->message(2,
"Blatter solver: Parameter continuation failed at step 1\n");
934 if (std::fabs(delta - delta_min) < 1e-6) {
935 m_log->message(2,
"Blatter solver: cannot reduce the continuation step\n");
942 delta =
pism::clip(delta, delta_min, delta_max);
948 m_log->message(2,
" Back-tracking to lambda = %f using delta = %f\n",
955 m_log->message(2,
"Blatter solver failed after %d parameter continuation steps\n", Nc);
966 ierr = SNESKSPSetUseEW(snes, PETSC_FALSE);
970 ierr = SNESGetKSP(snes, &ksp);
973 ierr = KSPSetTolerances(ksp, rtol, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
980 ierr = SNESKSPSetUseEW(snes, PETSC_TRUE);
PISM_CHK(ierr,
"SNESKSPSetUseEW");
989 schoofLen =
m_config->get_number(
"flow_law.Schoof_regularizing_length",
"m"),
990 schoofVel =
m_config->get_number(
"flow_law.Schoof_regularizing_velocity",
"m second-1");
1005 int snes_total_it = 0;
1006 int ksp_total_it = 0;
1008 double ksp_rtol = 1e-5;
1012 ierr = VecNorm(
m_x, NORM_INFINITY, &norm);
PISM_CHK(ierr,
"VecNorm");
1019 "Blatter solver: zero initial guess\n"
1020 " Disabling the Eisenstat-Walker method of adjusting solver tolerances\n");
1028 snes_total_it += info.
snes_it;
1029 ksp_total_it += info.
ksp_it;
1034 m_log->message(2,
"Blatter solver: %s\n", SNESConvergedReasons[info.
snes_reason]);
1038 m_log->message(2,
" Trying without the Eisenstat-Walker method of adjusting solver tolerances\n");
1041 if (not (info.
snes_reason == SNES_DIVERGED_LINE_SEARCH or
1056 snes_total_it += info.
snes_it;
1057 ksp_total_it += info.
ksp_it;
1062 m_log->message(2,
"Blatter solver: %s\n", SNESConvergedReasons[info.
snes_reason]);
1068 snes_total_it += info.
snes_it;
1069 ksp_total_it += info.
ksp_it;
1073 m_log->message(2,
"Blatter solver: %s\n", SNESConvergedReasons[info.
snes_reason]);
1081 "Blatter solver: %s. Done.\n"
1082 " SNES: %d, KSP: %d\n",
1084 snes_total_it, ksp_total_it);
1087 " Level 0 KSP (last iteration): %d\n",
1106 int ierr = DMDAVecGetArray(
m_da,
m_x, &x);
PISM_CHK(ierr,
"DMDAVecGetArray");
1112 for (
auto p =
m_grid->points(); p; p.next()) {
1113 const int i = p.i(), j = p.j();
1118 for (
int k = 0;
k < Mz; ++
k) {
1119 u[
k] = x[j][i][
k].
u;
1120 v[
k] = x[j][i][
k].
v;
1124 ierr = DMDAVecRestoreArray(
m_da,
m_x, &x);
PISM_CHK(ierr,
"DMDAVecRestoreArray");
1129 int ierr = DMDAVecGetArray(
m_da,
m_x, &x);
PISM_CHK(ierr,
"DMDAVecGetArray");
1133 for (
auto p =
m_grid->points(); p; p.next()) {
1134 const int i = p.i(), j = p.j();
1136 result(i, j) = x[j][i][0];
1139 ierr = DMDAVecRestoreArray(
m_da,
m_x, &x);
PISM_CHK(ierr,
"DMDAVecRestoreArray");
1146 int ierr = DMDAVecGetArray(
m_da,
m_x, &x);
PISM_CHK(ierr,
"DMDAVecGetArray");
1152 for (
auto p =
m_grid->points(); p; p.next()) {
1153 const int i = p.i(), j = p.j();
1158 for (
int k = 0;
k < Mz; ++
k) {
1159 x[j][i][
k].
u = u[
k];
1160 x[j][i][
k].
v = v[
k];
1164 ierr = DMDAVecRestoreArray(
m_da,
m_x, &x);
PISM_CHK(ierr,
"DMDAVecRestoreArray");
1168 PetscErrorCode ierr;
1171 ierr = DMDAVecGetArray(
m_da,
m_x, &x);
PISM_CHK(ierr,
"DMDAVecGetArray");
1177 for (
auto p =
m_grid->points(); p; p.next()) {
1178 const int i = p.i(), j = p.j();
1186 double dz = H / (Mz - 1);
1187 for (
int k = 0;
k < Mz - 1; ++
k) {
1188 V += x[j][i][
k] + x[j][i][
k + 1];
1190 V *= (0.5 * dz) / H;
1196 ierr = DMDAVecRestoreArray(
m_da,
m_x, &x);
PISM_CHK(ierr,
"DMDAVecRestoreArray");
#define ICE_GOLDSBY_KOHLSTEDT
std::shared_ptr< const Grid > grid() const
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
unsigned int nrecords() const
Get the number of records. Uses the length of an unlimited dimension.
bool variable_exists(const std::string &short_name) const
Checks if a variable exists.
High-level PISM I/O class.
array::Scalar1 sea_level_elevation
array::CellType2 cell_type
array::Scalar2 ice_thickness
array::Scalar2 bed_elevation
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.
double * get_column(int i, int j)
A virtual class collecting methods common to ice and bedrock 3D fields.
const std::vector< double > & levels() const
void set(double c)
Result: v[j] <- c for all j.
void update_ghosts()
Updates ghost points.
void local_to_global(unsigned int k, int &i, int &j) const
Convert a local degree of freedom index k to a global degree of freedom index (i,j).
void reset(int i, int j)
Initialize the Element to element (i, j) for the purposes of inserting into global residual and Jacob...
void nodal_values(const array::Scalar &x_global, int *result) const
Get nodal values of an integer mask.
GlobalIndex local_to_global(int i, int j, int k, unsigned int n) const
Q1 element with sides parallel to X and Y axes.
unsigned int n_pts() const
Number of quadrature points.
The 1-point Gaussian quadrature on the square [-1,1]*[-1,1].
std::shared_ptr< FlowLaw > create()
void remove(const std::string &name)
void define_model_state_impl(const File &output) const
The default (empty implementation).
void get_basal_velocity(array::Vector &result)
virtual Vector2d u_bc(double x, double y, double z) const
void compute_node_type(double min_thickness)
void write_model_state_impl(const File &output) const
The default (empty implementation).
virtual bool marine_boundary(int face, const int *node_type, const double *ice_bottom, const double *sea_level)
fem::Q1Element3Face m_face4
static bool partially_submerged_face(int face, const double *z, const double *sea_level)
virtual void init_2d_parameters(const Inputs &inputs)
std::shared_ptr< array::Array3D > m_u_sigma
std::shared_ptr< array::Array3D > velocity_v_sigma() const
SolutionInfo parameter_continuation()
void update(const Inputs &inputs, bool)
PetscErrorCode setup(DM pism_da, grid::Periodicity p, int Mz, int coarsening_factor, const std::string &prefix)
void init_ice_hardness(const Inputs &inputs, const petsc::DM &da)
static PetscErrorCode jacobian_callback(DMDALocalInfo *info, const Vector2d ***x, Mat A, Mat J, Blatter *solver)
static PetscErrorCode function_callback(DMDALocalInfo *info, const Vector2d ***x, Vector2d ***f, Blatter *solver)
fem::Q1Element3Face m_face100
void compute_averaged_velocity(array::Vector &result)
std::shared_ptr< array::Array3D > velocity_u_sigma() const
Blatter(std::shared_ptr< const Grid > grid, int Mz, int coarsening_factor)
static bool exterior_element(const int *node_type)
static bool grounding_line(const double *F)
array::Array2D< Parameters > m_parameters
virtual bool dirichlet_node(const DMDALocalInfo &info, const fem::Element3::GlobalIndex &I)
std::shared_ptr< array::Array3D > m_v_sigma
virtual void nodal_parameter_values(const fem::Q1Element3 &element, Parameters **P, int i, int j, int *node_type, double *bottom, double *thickness, double *surface, double *sea_level) const
void set_initial_guess(const array::Array3D &u_sigma, const array::Array3D &v_sigma)
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.
array::Scalar m_basal_frictional_heating
std::shared_ptr< rheology::FlowLaw > m_flow_law
EnthalpyConverter::Ptr m_EC
array::Vector2 m_velocity
Shallow stress balance (such as the SSA).
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
const int n_chi
Number of shape functions on a Q1 element.
const unsigned int incident_nodes[n_faces][4]
@ PISM_READONLY
open an existing file for reading only
static void disable_ew(::SNES snes, double rtol)
static void enable_ew(::SNES snes)
static PetscErrorCode blatter_restriction_hook(DM fine, Mat mrestrict, Vec rscale, Mat inject, DM coarse, void *ctx)
static PetscErrorCode blatter_coarsening_hook(DM dm_fine, DM dm_coarse, void *ctx)
double grid_z(double b, double H, int Mz, int k)
PetscErrorCode setup_level(DM dm, int mg_levels)
static double F(double SL, double B, double H, double alpha)
InputOptions process_input_options(MPI_Comm com, Config::ConstPtr config)
DMDALocalInfo grid_transpose(const DMDALocalInfo &input)
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
std::string printf(const char *format,...)
PetscErrorCode restrict_data(DM fine, DM coarse, const char *dm_name)
Restrict model parameters from the fine grid onto the coarse grid.
void GlobalMin(MPI_Comm comm, double *local, double *result, int count)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
PetscErrorCode create_restriction(DM fine, DM coarse, const char *dm_name)
Create the restriction matrix.
PetscInt mg_coarse_ksp_it
SNESConvergedReason snes_reason
KSPConvergedReason ksp_reason