tssa_testi.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tssa_testi.cc (6905B)
---
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_TESTI\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 I. 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/Context.hh"
31 #include "pism/util/VariableMetadata.hh"
32 #include "pism/util/error_handling.hh"
33 #include "pism/util/iceModelVec.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
43 const double m_schoof = 10; // (pure number)
44 const double L_schoof = 40e3; // meters
45 const double aspect_schoof = 0.05; // (pure)
46 const double H0_schoof = aspect_schoof * L_schoof;
47 // = 2000 m THICKNESS
48 const double B_schoof = 3.7e8; // Pa s^{1/3}; hardness
49 // given on p. 239 of Schoof; why so big?
50
51 class SSATestCaseI: public SSATestCase {
52 public:
53 SSATestCaseI(Context::Ptr ctx, int Mx, int My, SSAFactory ssafactory)
54 : SSATestCase(ctx,
55 Mx, My,
56 std::max(60.0e3, ((Mx - 1) / 2) * (2.0 * (3.0 * L_schoof) / (My - 1))),
57 3.0 * L_schoof,
58 CELL_CORNER,
59 NOT_PERIODIC) {
60 m_enthalpyconverter = EnthalpyConverter::Ptr(new EnthalpyConverter(*m_config));
61
62 m_config->set_flag("basal_resistance.pseudo_plastic.enabled", false);
63
64 m_config->set_string("stress_balance.ssa.flow_law", "isothermal_glen");
65 m_config->set_number("flow_law.isothermal_Glen.ice_softness", pow(B_schoof, -m_config->get_number("stress_balance.ssa.Glen_exponent")));
66
67 m_ssa = ssafactory(m_grid);
68 }
69
70 protected:
71 virtual void initializeSSACoefficients();
72
73 virtual void exactSolution(int i, int j,
74 double x, double y, double *u, double *v);
75
76 };
77
78 void SSATestCaseI::initializeSSACoefficients() {
79
80 m_bc_mask.set(0);
81 m_geometry.ice_thickness.set(H0_schoof);
82
83 double enth0 = m_enthalpyconverter->enthalpy(273.15, 0.01, 0.0); // 0.01 water fraction
84 m_ice_enthalpy.set(enth0);
85
86 // ssa->strength_extension->set_min_thickness(2*H0_schoof);
87
88 // The finite difference code uses the following flag to treat the non-periodic grid correctly.
89 m_config->set_flag("stress_balance.ssa.compute_surface_gradient_inward", true);
90 m_config->set_number("stress_balance.ssa.epsilon", 0.0); // don't use this lower bound
91
92 IceModelVec::AccessList list{&m_tauc, &m_bc_values, &m_bc_mask, &m_geometry.ice_surface_elevation, &m_geometry.bed_elevation};
93
94 double standard_gravity = m_config->get_number("constants.standard_gravity"),
95 ice_rho = m_config->get_number("constants.ice.density");
96
97 for (Points p(*m_grid); p; p.next()) {
98 const int i = p.i(), j = p.j();
99
100 const double y = m_grid->y(j);
101 const double theta = atan(0.001); /* a slope of 1/1000, a la Siple streams */
102 const double f = ice_rho * standard_gravity * H0_schoof * tan(theta);
103 m_tauc(i,j) = f * pow(fabs(y / L_schoof), m_schoof);
104 }
105 m_tauc.update_ghosts();
106
107 for (Points p(*m_grid); p; p.next()) {
108 const int i = p.i(), j = p.j();
109
110 const double myx = m_grid->x(i), myy=m_grid->y(j);
111 // eval exact solution; will only use exact vels if at edge
112 struct TestIParameters I_parameters = exactI(m_schoof, myx, myy);
113 m_geometry.bed_elevation(i, j) = I_parameters.bed;
114 m_geometry.ice_surface_elevation(i,j) = m_geometry.bed_elevation(i,j) + H0_schoof;
115
116 bool edge = ((j == 0) || (j == (int)m_grid->My() - 1) ||
117 (i == 0) || (i == (int)m_grid->Mx() - 1));
118 if (edge) {
119 m_bc_mask(i,j) = 1;
120 m_bc_values(i,j).u = I_parameters.u;
121 m_bc_values(i,j).v = I_parameters.v;
122 }
123 }
124
125 // communicate what we have set
126 m_geometry.ice_surface_elevation.update_ghosts();
127 m_geometry.bed_elevation.update_ghosts();
128 m_bc_mask.update_ghosts();
129 m_bc_values.update_ghosts();
130 }
131
132
133 void SSATestCaseI::exactSolution(int /*i*/, int /*j*/,
134 double x, double y,
135 double *u, double *v) {
136 struct TestIParameters I_parameters = exactI(m_schoof, x, y);
137 *u = I_parameters.u;
138 *v = I_parameters.v;
139 }
140
141 } // end of namespace stressbalance
142 } // end of namespace pism
143
144 int main(int argc, char *argv[]) {
145
146 using namespace pism;
147 using namespace pism::stressbalance;
148
149 MPI_Comm com = MPI_COMM_WORLD;
150 petsc::Initializer petsc(argc, argv, help);
151
152 com = PETSC_COMM_WORLD;
153
154 /* This explicit scoping forces destructors to be called before PetscFinalize() */
155 try {
156 Context::Ptr ctx = context_from_options(com, "ssa_testi");
157 Config::Ptr config = ctx->config();
158
159 std::string usage = "\n"
160 "usage of SSA_TESTi:\n"
161 " run ssa_testi -Mx <number> -My <number> -ssa_method <fd|fem>\n"
162 "\n";
163
164 bool stop = show_usage_check_req_opts(*ctx->log(), "ssa_testi", {}, usage);
165
166 if (stop) {
167 return 0;
168 }
169
170 // Parameters that can be overridden by command line options
171 unsigned int Mx = config->get_number("grid.Mx");
172 unsigned int My = config->get_number("grid.My");
173
174 auto method = config->get_string("stress_balance.ssa.method");
175 auto output_file = config->get_string("output.file_name");
176
177 // Determine the kind of solver to use.
178 SSAFactory ssafactory = NULL;
179 if (method == "fem") {
180 ssafactory = SSAFEMFactory;
181 } else if (method == "fd") {
182 ssafactory = SSAFDFactory;
183 } else {
184 /* can't happen */
185 }
186
187 SSATestCaseI testcase(ctx, Mx, My, ssafactory);
188 testcase.init();
189 testcase.run();
190 testcase.report("I");
191 testcase.write(output_file);
192 }
193 catch (...) {
194 handle_fatal_errors(com);
195 return 1;
196 }
197
198 return 0;
199 }