PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
ssa_testj.cc
Go to the documentation of this file.
1 // Copyright (C) 2010--2018, 2021, 2022 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 static 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 "pism/basalstrength/basal_resistance.hh" // IceBasalResistancePlasticLaw
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/Mask.hh"
31 #include "pism/util/Context.hh"
32 #include "pism/util/VariableMetadata.hh"
33 #include "pism/util/error_handling.hh"
34 #include "pism/util/io/File.hh"
35 #include "pism/util/petscwrappers/PetscInitializer.hh"
36 #include "pism/util/pism_utilities.hh"
37 #include "pism/util/pism_options.hh"
38 #include "pism/verification/tests/exactTestsIJ.h"
39 
40 namespace pism {
41 namespace stressbalance {
42 
44 {
45 public:
46  SSATestCaseJ(std::shared_ptr<Context> ctx, int Mx, int My, SSAFactory ssafactory)
47  : SSATestCase(ctx, Mx, My, 300e3, 300e3, grid::CELL_CENTER, grid::XY_PERIODIC) {
48  m_config->set_flag("basal_resistance.pseudo_plastic.enabled", false);
49 
51  m_config->set_string("stress_balance.ssa.flow_law", "isothermal_glen");
52 
53  m_ssa = ssafactory(m_grid);
54  }
55 
56 protected:
57  virtual void initializeSSACoefficients();
58 
59  virtual void exactSolution(int i, int j,
60  double x, double y, double *u, double *v);
61 };
62 
64  m_tauc.set(0.0); // irrelevant for test J
65  m_geometry.bed_elevation.set(-1000.0); // assures shelf is floating (maximum ice thickness is 770 m)
67 
68  double enth0 = m_enthalpyconverter->enthalpy(273.15, 0.01, 0.0); // 0.01 water fraction
69  m_ice_enthalpy.set(enth0);
70 
71  /* use Ritz et al (2001) value of 30 MPa year for typical vertically-averaged viscosity */
72  double ocean_rho = m_config->get_number("constants.sea_water.density"),
73  ice_rho = m_config->get_number("constants.ice.density");
74  const double nu0 = units::convert(m_sys, 30.0, "MPa year", "Pa s"); /* = 9.45e14 Pa s */
75  const double H0 = 500.; /* 500 m typical thickness */
76 
77  // Test J has a viscosity that is independent of velocity. So we force a
78  // constant viscosity by settting the strength_extension
79  // thickness larger than the given ice thickness. (max = 770m).
82 
84 
85  for (auto p = m_grid->points(); p; p.next()) {
86  const int i = p.i(), j = p.j();
87 
88  const double myx = m_grid->x(i), myy = m_grid->y(j);
89 
90  // set H,h on regular grid
91  struct TestJParameters J_parameters = exactJ(myx, myy);
92 
93  m_geometry.ice_thickness(i,j) = J_parameters.H;
94  m_geometry.ice_surface_elevation(i,j) = (1.0 - ice_rho / ocean_rho) * J_parameters.H; // FIXME issue #15
95 
96  // special case at center point: here we set bc_values at (i,j) by
97  // setting bc_mask and bc_values appropriately
98  if ((i == ((int)m_grid->Mx()) / 2) and
99  (j == ((int)m_grid->My()) / 2)) {
100  m_bc_mask(i,j) = 1;
101  m_bc_values(i,j).u = J_parameters.u;
102  m_bc_values(i,j).v = J_parameters.v;
103  }
104  }
105 
106  // communicate what we have set
111 }
112 
113 void SSATestCaseJ::exactSolution(int /*i*/, int /*j*/,
114  double x, double y,
115  double *u, double *v) {
116  struct TestJParameters J_parameters = exactJ(x, y);
117  *u = J_parameters.u;
118  *v = J_parameters.v;
119 }
120 
121 } // end of namespace stressbalance
122 } // end of namespace pism
123 
124 int main(int argc, char *argv[]) {
125 
126  using namespace pism;
127  using namespace pism::stressbalance;
128 
129  MPI_Comm com = MPI_COMM_WORLD;
130  petsc::Initializer petsc(argc, argv, help);
131 
132  com = PETSC_COMM_WORLD;
133 
134  /* This explicit scoping forces destructors to be called before PetscFinalize() */
135  try {
136  std::shared_ptr<Context> ctx = context_from_options(com, "ssa_testj");
137  Config::Ptr config = ctx->config();
138 
139  std::string usage = "\n"
140  "usage of SSA_TESTJ:\n"
141  " run ssafe_test -Mx <number> -My <number> -ssa_method <fd|fem>\n"
142  "\n";
143 
144  bool stop = show_usage_check_req_opts(*ctx->log(), "ssa_testj", {}, usage);
145 
146  if (stop) {
147  return 0;
148  }
149 
150  // Parameters that can be overridden by command line options
151 
152  unsigned int Mx = config->get_number("grid.Mx");
153  unsigned int My = config->get_number("grid.My");
154 
155  auto method = config->get_string("stress_balance.ssa.method");
156  auto output = config->get_string("output.file");
157 
158  bool write_output = config->get_string("output.size") != "none";
159 
160  // Determine the kind of solver to use.
161  SSAFactory ssafactory = NULL;
162  if (method == "fem") {
163  ssafactory = SSAFEMFactory;
164  } else if (method == "fd") {
165  ssafactory = SSAFDFactory;
166  } else {
167  /* can't happen */
168  }
169 
170  SSATestCaseJ testcase(ctx, Mx, My, ssafactory);
171  testcase.init();
172  testcase.run();
173  testcase.report("J");
174  if (write_output) {
175  testcase.write(output);
176  }
177  }
178  catch (...) {
179  handle_fatal_errors(com);
180  return 1;
181  }
182 
183  return 0;
184 }
std::shared_ptr< Config > Ptr
std::shared_ptr< EnthalpyConverter > Ptr
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:65
void set(double c)
Result: v[j] <- c for all j.
Definition: Array.cc:707
void update_ghosts()
Updates ghost points.
Definition: Array.cc:693
void set_min_thickness(double my_min_thickness)
Set minimum thickness to trigger use of extension.
Definition: SSA.cc:50
void set_notional_strength(double my_nuH)
Set strength = (viscosity times thickness).
Definition: SSA.cc:41
virtual void initializeSSACoefficients()
Set up the coefficient variables as appropriate for the test case.
Definition: ssa_testj.cc:63
SSATestCaseJ(std::shared_ptr< Context > ctx, int Mx, int My, SSAFactory ssafactory)
Definition: ssa_testj.cc:46
virtual void exactSolution(int i, int j, double x, double y, double *u, double *v)
Definition: ssa_testj.cc:113
virtual void init()
Initialize the test case at the start of a run.
Definition: SSATestCase.cc:94
EnthalpyConverter::Ptr m_enthalpyconverter
Definition: SSATestCase.hh:102
const units::System::Ptr m_sys
Definition: SSATestCase.hh:96
std::shared_ptr< Grid > m_grid
Definition: SSATestCase.hh:94
virtual void report(const std::string &testname)
Report on the generated solution.
Definition: SSATestCase.cc:121
virtual void run()
Solve the SSA.
Definition: SSATestCase.cc:103
virtual void write(const std::string &filename)
Save the computation and data to a file.
Definition: SSATestCase.cc:296
SSAStrengthExtension * strength_extension
Definition: SSA.hh:115
const double H0
Definition: exactTestK.c:37
struct TestJParameters exactJ(const double x, const double y)
Definition: exactTestsIJ.c:80
@ XY_PERIODIC
Definition: Grid.hh:51
@ CELL_CENTER
Definition: Grid.hh:53
SSA * SSAFDFactory(std::shared_ptr< const Grid > g)
Constructs a new SSAFD.
Definition: SSAFD.cc:49
SSA *(* SSAFactory)(std::shared_ptr< const Grid >)
Definition: SSA.hh:77
SSA * SSAFEMFactory(std::shared_ptr< const Grid > g)
Factory function for constructing a new SSAFEM.
Definition: SSAFEM.cc:112
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
std::shared_ptr< Context > context_from_options(MPI_Comm com, const std::string &prefix, bool print)
Create a default context using options.
Definition: Context.cc:170
@ MASK_FLOATING
Definition: Mask.hh:34
void handle_fatal_errors(MPI_Comm com)
bool show_usage_check_req_opts(const Logger &log, const std::string &execname, const std::vector< std::string > &required_options, const std::string &usage)
In a single call a driver program can provide a usage string to the user and check if required option...
Definition: pism_options.cc:51
int main(int argc, char *argv[])
Definition: ssa_testj.cc:124
static char help[]
Definition: ssa_testj.cc:19