tssa_testj.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tssa_testj.cc (6281B)
---
1 // Copyright (C) 2010--2018 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/iceModelVec.hh"
35 #include "pism/util/io/File.hh"
36 #include "pism/util/petscwrappers/PetscInitializer.hh"
37 #include "pism/util/pism_utilities.hh"
38 #include "pism/util/pism_options.hh"
39 #include "pism/verification/tests/exactTestsIJ.h"
40
41 namespace pism {
42 namespace stressbalance {
43
44 class SSATestCaseJ: public SSATestCase
45 {
46 public:
47 SSATestCaseJ(Context::Ptr ctx, int Mx, int My, SSAFactory ssafactory)
48 : SSATestCase(ctx, Mx, My, 300e3, 300e3, CELL_CENTER, XY_PERIODIC) {
49 m_config->set_flag("basal_resistance.pseudo_plastic.enabled", false);
50
51 m_enthalpyconverter = EnthalpyConverter::Ptr(new EnthalpyConverter(*m_config));
52 m_config->set_string("stress_balance.ssa.flow_law", "isothermal_glen");
53
54 m_ssa = ssafactory(m_grid);
55 }
56
57 protected:
58 virtual void initializeSSACoefficients();
59
60 virtual void exactSolution(int i, int j,
61 double x, double y, double *u, double *v);
62 };
63
64 void SSATestCaseJ::initializeSSACoefficients() {
65 m_tauc.set(0.0); // irrelevant for test J
66 m_geometry.bed_elevation.set(-1000.0); // assures shelf is floating (maximum ice thickness is 770 m)
67 m_geometry.cell_type.set(MASK_FLOATING);
68
69 double enth0 = m_enthalpyconverter->enthalpy(273.15, 0.01, 0.0); // 0.01 water fraction
70 m_ice_enthalpy.set(enth0);
71
72 /* use Ritz et al (2001) value of 30 MPa year for typical vertically-averaged viscosity */
73 double ocean_rho = m_config->get_number("constants.sea_water.density"),
74 ice_rho = m_config->get_number("constants.ice.density");
75 const double nu0 = units::convert(m_sys, 30.0, "MPa year", "Pa s"); /* = 9.45e14 Pa s */
76 const double H0 = 500.; /* 500 m typical thickness */
77
78 // Test J has a viscosity that is independent of velocity. So we force a
79 // constant viscosity by settting the strength_extension
80 // thickness larger than the given ice thickness. (max = 770m).
81 m_ssa->strength_extension->set_notional_strength(nu0 * H0);
82 m_ssa->strength_extension->set_min_thickness(800);
83
84 IceModelVec::AccessList list{&m_geometry.ice_thickness, &m_geometry.ice_surface_elevation, &m_bc_mask, &m_bc_values};
85
86 for (Points p(*m_grid); p; p.next()) {
87 const int i = p.i(), j = p.j();
88
89 const double myx = m_grid->x(i), myy = m_grid->y(j);
90
91 // set H,h on regular grid
92 struct TestJParameters J_parameters = exactJ(myx, myy);
93
94 m_geometry.ice_thickness(i,j) = J_parameters.H;
95 m_geometry.ice_surface_elevation(i,j) = (1.0 - ice_rho / ocean_rho) * J_parameters.H; // FIXME issue #15
96
97 // special case at center point: here we set bc_values at (i,j) by
98 // setting bc_mask and bc_values appropriately
99 if ((i == ((int)m_grid->Mx()) / 2) and
100 (j == ((int)m_grid->My()) / 2)) {
101 m_bc_mask(i,j) = 1;
102 m_bc_values(i,j).u = J_parameters.u;
103 m_bc_values(i,j).v = J_parameters.v;
104 }
105 }
106
107 // communicate what we have set
108 m_geometry.ice_surface_elevation.update_ghosts();
109 m_geometry.ice_thickness.update_ghosts();
110 m_bc_mask.update_ghosts();
111 m_bc_values.update_ghosts();
112 }
113
114 void SSATestCaseJ::exactSolution(int /*i*/, int /*j*/,
115 double x, double y,
116 double *u, double *v) {
117 struct TestJParameters J_parameters = exactJ(x, y);
118 *u = J_parameters.u;
119 *v = J_parameters.v;
120 }
121
122 } // end of namespace stressbalance
123 } // end of namespace pism
124
125 int main(int argc, char *argv[]) {
126
127 using namespace pism;
128 using namespace pism::stressbalance;
129
130 MPI_Comm com = MPI_COMM_WORLD;
131 petsc::Initializer petsc(argc, argv, help);
132
133 com = PETSC_COMM_WORLD;
134
135 /* This explicit scoping forces destructors to be called before PetscFinalize() */
136 try {
137 Context::Ptr ctx = context_from_options(com, "ssa_testj");
138 Config::Ptr config = ctx->config();
139
140 std::string usage = "\n"
141 "usage of SSA_TESTJ:\n"
142 " run ssafe_test -Mx <number> -My <number> -ssa_method <fd|fem>\n"
143 "\n";
144
145 bool stop = show_usage_check_req_opts(*ctx->log(), "ssa_testj", {}, usage);
146
147 if (stop) {
148 return 0;
149 }
150
151 // Parameters that can be overridden by command line options
152
153 unsigned int Mx = config->get_number("grid.Mx");
154 unsigned int My = config->get_number("grid.My");
155
156 auto method = config->get_string("stress_balance.ssa.method");
157 auto output = config->get_string("output.file_name");
158
159 // Determine the kind of solver to use.
160 SSAFactory ssafactory = NULL;
161 if (method == "fem") {
162 ssafactory = SSAFEMFactory;
163 } else if (method == "fd") {
164 ssafactory = SSAFDFactory;
165 } else {
166 /* can't happen */
167 }
168
169 SSATestCaseJ testcase(ctx, Mx, My, ssafactory);
170 testcase.init();
171 testcase.run();
172 testcase.report("J");
173 testcase.write(output);
174 }
175 catch (...) {
176 handle_fatal_errors(com);
177 return 1;
178 }
179
180 return 0;
181 }