Loading [MathJax]/jax/input/TeX/config.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_j.cc
Go to the documentation of this file.
1// Copyright (C) 2010--2018, 2021, 2022, 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
19static char help[] =
20 "\nSSA_TESTJ\n"
21 " Testing program for the finite element implementation of the SSA.\n"
22 " Does a time-independent calculation. Does not run IceModel or a derived\n"
23 " class thereof. Uses verification test J. Also may be used in a PISM\n"
24 " software (regression) test.\n\n";
25
26#include <petscsys.h>
27
28#include "pism/stressbalance/ssa/SSATestCase.hh"
29#include "pism/util/Mask.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
39std::shared_ptr<Grid> ssa_test_j_grid(std::shared_ptr<Context> ctx, int Mx, int My) {
40 return SSATestCase::grid(ctx, Mx, My, 300e3, 300e3, grid::CELL_CENTER, grid::XY_PERIODIC);
41}
42
44{
45public:
46 SSATestCaseJ(std::shared_ptr<SSA> ssa) : SSATestCase(ssa) {
47
49 // 0.01 water fraction
50 m_ice_enthalpy.set(EC.enthalpy(273.15, 0.01, 0.0));
51 }
52
53protected:
54 virtual void initializeSSACoefficients();
55
56 virtual void exactSolution(int i, int j,
57 double x, double y, double *u, double *v);
58};
59
61 m_tauc.set(0.0); // irrelevant for test J
62 m_geometry.bed_elevation.set(-1000.0); // assures shelf is floating (maximum ice thickness is 770 m)
64
65 /* use Ritz et al (2001) value of 30 MPa year for typical vertically-averaged viscosity */
66 double ocean_rho = m_config->get_number("constants.sea_water.density"),
67 ice_rho = m_config->get_number("constants.ice.density");
68 const double nu0 = units::convert(m_sys, 30.0, "MPa year", "Pa s"); /* = 9.45e14 Pa s */
69 const double H0 = 500.; /* 500 m typical thickness */
70
71 // Test J has a viscosity that is independent of velocity. So we force a
72 // constant viscosity by settting the strength_extension
73 // thickness larger than the given ice thickness. (max = 770m).
74 m_ssa->strength_extension->set_notional_strength(nu0 * H0);
75 m_ssa->strength_extension->set_min_thickness(800);
76
78
79 for (auto p = m_grid->points(); p; p.next()) {
80 const int i = p.i(), j = p.j();
81
82 const double myx = m_grid->x(i), myy = m_grid->y(j);
83
84 // set H,h on regular grid
85 struct TestJParameters J_parameters = exactJ(myx, myy);
86
87 m_geometry.ice_thickness(i,j) = J_parameters.H;
88 m_geometry.ice_surface_elevation(i,j) = (1.0 - ice_rho / ocean_rho) * J_parameters.H; // FIXME issue #15
89
90 // special case at center point: here we set bc_values at (i,j) by
91 // setting bc_mask and bc_values appropriately
92 if ((i == ((int)m_grid->Mx()) / 2) and
93 (j == ((int)m_grid->My()) / 2)) {
94 m_bc_mask(i,j) = 1;
95 m_bc_values(i,j).u = J_parameters.u;
96 m_bc_values(i,j).v = J_parameters.v;
97 }
98 }
99
100 // communicate what we have set
105}
106
107void SSATestCaseJ::exactSolution(int /*i*/, int /*j*/,
108 double x, double y,
109 double *u, double *v) {
110 struct TestJParameters J_parameters = exactJ(x, y);
111 *u = J_parameters.u;
112 *v = J_parameters.v;
113}
114
115} // end of namespace stressbalance
116} // end of namespace pism
117
118int main(int argc, char *argv[]) {
119
120 using namespace pism;
121 using namespace pism::stressbalance;
122
123 MPI_Comm com = MPI_COMM_WORLD;
124 petsc::Initializer petsc(argc, argv, help);
125
126 com = PETSC_COMM_WORLD;
127
128 /* This explicit scoping forces destructors to be called before PetscFinalize() */
129 try {
130 std::shared_ptr<Context> ctx = context_from_options(com, "ssa_testj");
131 Config::Ptr config = ctx->config();
132
133 std::string usage = "\n"
134 "usage of SSA_TESTJ:\n"
135 " run ssafe_test -Mx <number> -My <number> -ssa_method <fd|fem>\n"
136 "\n";
137
138 bool stop = show_usage_check_req_opts(*ctx->log(), "ssa_testj", {}, usage);
139
140 if (stop) {
141 return 0;
142 }
143
144 // Parameters that can be overridden by command line options
145
146 unsigned int Mx = config->get_number("grid.Mx");
147 unsigned int My = config->get_number("grid.My");
148
149 auto method = config->get_string("stress_balance.ssa.method");
150 auto output_file = config->get_string("output.file");
151
152 bool write_output = config->get_string("output.size") != "none";
153
154 config->set_flag("basal_resistance.pseudo_plastic.enabled", false);
155 config->set_string("stress_balance.ssa.flow_law", "isothermal_glen");
156
157 auto grid = ssa_test_j_grid(ctx, Mx, My);
158 SSATestCaseJ testcase(SSATestCase::solver(grid, method));
159 testcase.init();
160 testcase.run();
161 testcase.report("J");
162 if (write_output) {
163 testcase.write(output_file);
164 }
165 }
166 catch (...) {
167 handle_fatal_errors(com);
168 return 1;
169 }
170
171 return 0;
172}
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.
array::Scalar2 ice_surface_elevation
Definition Geometry.hh:57
array::CellType2 cell_type
Definition Geometry.hh:55
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 initializeSSACoefficients()
Set up the coefficient variables as appropriate for the test case.
Definition ssa_test_j.cc:60
SSATestCaseJ(std::shared_ptr< SSA > ssa)
Definition ssa_test_j.cc:46
virtual void exactSolution(int i, int j, double x, double y, double *u, double *v)
std::shared_ptr< const pism::Grid > m_grid
virtual void init()
Initialize the test case at the start of a run.
const units::System::Ptr m_sys
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
std::shared_ptr< SSA > m_ssa
virtual void write(const std::string &filename)
Save the computation and data to a file.
#define H0
Definition exactTestM.c:39
struct TestJParameters exactJ(const double x, const double y)
@ XY_PERIODIC
Definition Grid.hh:54
@ CELL_CENTER
Definition Grid.hh:56
std::shared_ptr< Grid > ssa_test_j_grid(std::shared_ptr< Context > ctx, int Mx, int My)
Definition ssa_test_j.cc:39
Stress balance models and related diagnostics.
Definition Isochrones.hh:31
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition Units.cc:70
@ MASK_FLOATING
Definition Mask.hh:34
int main(int argc, char *argv[])
static char help[]
Definition ssa_test_j.cc:19