Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ssa_test_i.cc
Go to the documentation of this file.
1// Copyright (C) 2010--2018, 2021, 2022, 2023, 2024 Ed Bueler, Constantine Khroulev, and David Maxwell
2//
3// This file is part of PISM.
4//
5// PISM is free software; you can redistribute it and/or modify it under the
6// terms of the GNU General Public License as published by the Free Software
7// Foundation; either version 3 of the License, or (at your option) any later
8// version.
9//
10// PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12// FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13// details.
14//
15// You should have received a copy of the GNU General Public License
16// along with PISM; if not, write to the Free Software
17// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19#include <memory>
20static char help[] =
21 "\nSSA_TESTI\n"
22 " Testing program for the finite element implementation of the SSA.\n"
23 " Does a time-independent calculation. Does not run IceModel or a derived\n"
24 " class thereof. Uses verification test I. Also may be used in a PISM\n"
25 " software (regression) test.\n\n";
26
27#include "pism/stressbalance/ssa/SSAFD.hh"
28#include "pism/stressbalance/ssa/SSAFEM.hh"
29#include "pism/stressbalance/ssa/SSATestCase.hh"
30#include "pism/util/Context.hh"
31#include "pism/util/error_handling.hh"
32#include "pism/util/petscwrappers/PetscInitializer.hh"
33#include "pism/util/pism_options.hh"
34#include "pism/verification/tests/exactTestsIJ.h"
35
36namespace pism {
37namespace stressbalance {
38
39const double m_schoof = 10; // (pure number)
40const double L_schoof = 40e3; // meters
41const double aspect_schoof = 0.05; // (pure)
43 // = 2000 m THICKNESS
44const double B_schoof = 3.7e8; // Pa s^{1/3}; hardness
45 // given on p. 239 of Schoof; why so big?
46
47std::shared_ptr<Grid> ssa_test_i_grid(std::shared_ptr<Context> ctx, int Mx, int My) {
48 return SSATestCase::grid(ctx, Mx, My,
49 std::max(60.0e3, ((Mx - 1) / 2) * (2.0 * (3.0 * L_schoof) / (My - 1))),
51}
52
54public:
55 SSATestCaseI(std::shared_ptr<SSA> ssa) : SSATestCase(ssa) {
57 // 0.01 water fraction
58 m_ice_enthalpy.set(EC.enthalpy(273.15, 0.01, 0.0));
59 }
60
61protected:
62 virtual void initializeSSACoefficients();
63
64 virtual void exactSolution(int i, int j,
65 double x, double y, double *u, double *v);
66
67};
68
70
71 m_bc_mask.set(0);
73
75
76 for (auto p = m_grid->points(); p; p.next()) {
77 const int i = p.i(), j = p.j();
78
79 // Evaluate the exact solution and yield stress. Exact u, v will only be used at the
80 // grid edge.
81 auto I = exactI(m_schoof, m_grid->x(i), m_grid->y(j));
82
83 m_geometry.bed_elevation(i, j) = I.bed;
84
85 m_tauc(i,j) = I.tauc;
86
87 if (grid::domain_edge(*m_grid, i, j)) {
88 m_bc_mask(i, j) = 1;
89 m_bc_values(i, j) = { I.u, I.v };
90 }
91 }
92
94
98}
99
100
101void SSATestCaseI::exactSolution(int /*i*/, int /*j*/,
102 double x, double y,
103 double *u, double *v) {
104 auto I = exactI(m_schoof, x, y);
105 *u = I.u;
106 *v = I.v;
107}
108
109} // end of namespace stressbalance
110} // end of namespace pism
111
112int main(int argc, char *argv[]) {
113
114 using namespace pism;
115 using namespace pism::stressbalance;
116
117 MPI_Comm com = MPI_COMM_WORLD;
118 petsc::Initializer petsc(argc, argv, help);
119
120 com = PETSC_COMM_WORLD;
121
122 /* This explicit scoping forces destructors to be called before PetscFinalize() */
123 try {
124 std::shared_ptr<Context> ctx = context_from_options(com, "ssa_testi");
125 Config::Ptr config = ctx->config();
126
127 std::string usage = "\n"
128 "usage of SSA_TESTi:\n"
129 " run ssa_testi -Mx <number> -My <number> -ssa_method <fd|fem>\n"
130 "\n";
131
132 bool stop = show_usage_check_req_opts(*ctx->log(), "ssa_testi", {}, usage);
133
134 if (stop) {
135 return 0;
136 }
137
138 // Parameters that can be overridden by command line options
139 unsigned int Mx = config->get_number("grid.Mx");
140 unsigned int My = config->get_number("grid.My");
141
142 auto method = config->get_string("stress_balance.ssa.method");
143 auto output_file = config->get_string("output.file");
144
145 bool write_output = config->get_string("output.size") != "none";
146
147 // we have to set parameters *before* `ssa` is allocated
148 config->set_flag("basal_resistance.pseudo_plastic.enabled", false);
149
150 config->set_string("stress_balance.ssa.flow_law", "isothermal_glen");
151 config->set_number("flow_law.isothermal_Glen.ice_softness",
152 pow(B_schoof, -config->get_number("stress_balance.ssa.Glen_exponent")));
153
154 // The finite difference code uses the following flag to treat the non-periodic grid correctly.
155 config->set_flag("stress_balance.ssa.compute_surface_gradient_inward", true);
156 config->set_number("stress_balance.ssa.epsilon", 0.0); // don't use this lower bound
157
158 auto grid = ssa_test_i_grid(ctx, Mx, My);
159 SSATestCaseI testcase(SSATestCase::solver(grid, method));
160 testcase.init();
161 testcase.run();
162 testcase.report("I");
163 if (write_output) {
164 testcase.write(output_file);
165 }
166 }
167 catch (...) {
168 handle_fatal_errors(com);
169 return 1;
170 }
171
172 return 0;
173}
std::shared_ptr< Config > Ptr
double enthalpy(double T, double omega, double P) const
Compute enthalpy from absolute temperature, liquid water fraction, and pressure.
Converts between specific enthalpy and temperature or liquid content.
void ensure_consistency(double ice_free_thickness_threshold)
Definition Geometry.cc:114
array::Scalar2 ice_surface_elevation
Definition Geometry.hh:57
array::Scalar2 ice_thickness
Definition Geometry.hh:51
array::Scalar2 bed_elevation
Definition Geometry.hh:47
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition Array.hh:64
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
void update_ghosts()
Updates ghost points.
Definition Array.cc:615
virtual void exactSolution(int i, int j, double x, double y, double *u, double *v)
virtual void initializeSSACoefficients()
Set up the coefficient variables as appropriate for the test case.
Definition ssa_test_i.cc:69
SSATestCaseI(std::shared_ptr< SSA > ssa)
Definition ssa_test_i.cc:55
std::shared_ptr< const pism::Grid > m_grid
virtual void init()
Initialize the test case at the start of a run.
static std::shared_ptr< Grid > grid(std::shared_ptr< Context > ctx, int Mx, int My, double Lx, double Ly, grid::Registration registration, grid::Periodicity periodicity)
virtual void report(const std::string &testname)
Report on the generated solution.
virtual void run()
Solve the SSA.
const Config::ConstPtr m_config
virtual void write(const std::string &filename)
Save the computation and data to a file.
struct TestIParameters exactI(const double m, const double x, const double y)
@ NOT_PERIODIC
Definition Grid.hh:54
@ CELL_CORNER
Definition Grid.hh:56
bool domain_edge(const Grid &grid, int i, int j)
Definition Grid.hh:409
const double B_schoof
Definition ssa_test_i.cc:44
const double m_schoof
Definition ssa_test_i.cc:39
const double aspect_schoof
Definition ssa_test_i.cc:41
const double L_schoof
Definition ssa_test_i.cc:40
std::shared_ptr< Grid > ssa_test_i_grid(std::shared_ptr< Context > ctx, int Mx, int My)
Definition ssa_test_i.cc:47
const double H0_schoof
Definition ssa_test_i.cc:42
Stress balance models and related diagnostics.
Definition Isochrones.hh:31
int main(int argc, char *argv[])
static char help[]
Definition ssa_test_i.cc:20