tssa_test_plug.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_test_plug.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 /* This file implements a test case for the ssa: plug flow. The geometry
20 consists of a constant surface slope in the positive x-direction, and the
21 ice is pinned on the y-boundaries. There is no basal shear stress, and hence
22 the the only nonzero terms in the SSA are the "p-laplacian" and the driving
23 stress.
24 */
25
26 static char help[] =
27 "\nSSA_TEST_PLUG\n"
28 " Testing program for the finite element implementation of the SSA.\n"
29 " Does a time-independent calculation. Does not run IceModel or a derived\n"
30 " class thereof.\n\n";
31
32 #include <cmath>
33
34 #include "pism/basalstrength/basal_resistance.hh" // IceBasalResistancePlasticLaw
35 #include "pism/stressbalance/ssa/SSAFD.hh"
36 #include "pism/stressbalance/ssa/SSAFEM.hh"
37 #include "pism/stressbalance/ssa/SSATestCase.hh"
38 #include "pism/util/Context.hh"
39 #include "pism/util/VariableMetadata.hh"
40 #include "pism/util/error_handling.hh"
41 #include "pism/util/iceModelVec.hh"
42 #include "pism/util/io/File.hh"
43 #include "pism/util/petscwrappers/PetscInitializer.hh"
44 #include "pism/util/pism_utilities.hh"
45 #include "pism/util/pism_options.hh"
46 #include "pism/verification/tests/exactTestsIJ.h"
47
48 namespace pism {
49 namespace stressbalance {
50
51 class SSATestCasePlug: public SSATestCase {
52 public:
53 SSATestCasePlug(Context::Ptr ctx, int Mx, int My,
54 double n, SSAFactory ssafactory)
55 : SSATestCase(ctx, Mx, My, 50e3, 50e3, CELL_CORNER, NOT_PERIODIC) {
56 H0 = 2000.; // m
57 L = 50.e3; // 50km half-width
58 dhdx = 0.001; // pure number, slope of surface & bed
59 tauc0 = 0.; // No basal shear stress
60
61 // Pa s^{1/3}; hardness given on p. 239 of Schoof; why so big?
62 B0 = 3.7e8;
63
64 this->glen_n = n;
65
66 // Basal sliding law parameters are irrelevant because tauc=0
67
68 // Enthalpy converter is irrelevant (but still required) for this test.
69 m_enthalpyconverter = EnthalpyConverter::Ptr(new EnthalpyConverter(*m_config));
70
71 // Use constant hardness
72 m_config->set_string("stress_balance.ssa.flow_law", "isothermal_glen");
73 m_config->set_number("flow_law.isothermal_Glen.ice_softness", pow(B0, -glen_n));
74
75 m_ssa = ssafactory(m_grid);
76 }
77
78 protected:
79 virtual void initializeSSACoefficients();
80
81 virtual void exactSolution(int i, int j,
82 double x, double y, double *u, double *v);
83
84
85 double H0; // Thickness
86 double L; // Half-width
87 double dhdx; // surface slope
88 double tauc0; // zero basal shear stress
89 double B0; // hardness
90 double glen_n;
91
92 bool dimensionless;
93
94 };
95
96 void SSATestCasePlug::initializeSSACoefficients() {
97
98 // The finite difference code uses the following flag to treat the non-periodic grid correctly.
99 m_config->set_flag("stress_balance.ssa.compute_surface_gradient_inward", true);
100 m_config->set_number("stress_balance.ssa.epsilon", 0.0);
101
102 // Ensure we never use the strength extension.
103 m_ssa->strength_extension->set_min_thickness(H0/2);
104
105 // Set constant coefficients.
106 m_geometry.ice_thickness.set(H0);
107 m_tauc.set(tauc0);
108
109 // Set boundary conditions (Dirichlet all the way around).
110 m_bc_mask.set(0.0);
111
112 IceModelVec::AccessList list{&m_bc_values, &m_bc_mask, &m_geometry.bed_elevation, &m_geometry.ice_surface_elevation};
113
114 for (Points p(*m_grid); p; p.next()) {
115 const int i = p.i(), j = p.j();
116
117 double myu, myv;
118 const double myx = m_grid->x(i), myy=m_grid->y(j);
119
120 m_geometry.bed_elevation(i,j) = -myx*(dhdx);
121 m_geometry.ice_surface_elevation(i,j) = m_geometry.bed_elevation(i,j) + H0;
122
123 bool edge = ((j == 0) || (j == (int)m_grid->My() - 1) ||
124 (i == 0) || (i == (int)m_grid->Mx() - 1));
125 if (edge) {
126 m_bc_mask(i,j) = 1;
127 exactSolution(i,j,myx,myy,&myu,&myv);
128 m_bc_values(i,j).u = myu;
129 m_bc_values(i,j).v = myv;
130 }
131 }
132
133 m_bc_values.update_ghosts();
134 m_bc_mask.update_ghosts();
135 m_geometry.bed_elevation.update_ghosts();
136 m_geometry.ice_surface_elevation.update_ghosts();
137 }
138
139 void SSATestCasePlug::exactSolution(int /*i*/, int /*j*/,
140 double /*x*/, double y,
141 double *u, double *v) {
142 double earth_grav = m_config->get_number("constants.standard_gravity"),
143 ice_rho = m_config->get_number("constants.ice.density");
144 double f = ice_rho * earth_grav * H0* dhdx;
145 double ynd = y/L;
146
147 *u = 0.5*pow(f,3)*pow(L,4)/pow(B0*H0,3)*(1-pow(ynd,4));
148 *v = 0;
149 }
150
151 } // end of namespace stressbalance
152 } // end of namespace pism
153
154 int main(int argc, char *argv[]) {
155
156 using namespace pism;
157 using namespace pism::stressbalance;
158
159 MPI_Comm com = MPI_COMM_WORLD; // won't be used except for rank,size
160 petsc::Initializer petsc(argc, argv, help);
161
162 com = PETSC_COMM_WORLD;
163
164 /* This explicit scoping forces destructors to be called before PetscFinalize() */
165 try {
166 Context::Ptr ctx = context_from_options(com, "ssa_test_plug");
167 Config::Ptr config = ctx->config();
168
169 std::string usage = "\n"
170 "usage of SSA_TEST_PLUG:\n"
171 " run ssa_test_plug -Mx <number> -My <number> -ssa_method <fd|fem>\n"
172 "\n";
173
174 bool stop = show_usage_check_req_opts(*ctx->log(), "ssa_test_plug", {}, usage);
175
176 if (stop) {
177 return 0;
178 }
179
180 // Parameters that can be overridden by command line options
181
182 unsigned int Mx = config->get_number("grid.Mx");
183 unsigned int My = config->get_number("grid.My");
184
185 auto method = config->get_string("stress_balance.ssa.method");
186 auto output_file = config->get_string("output.file_name");
187 auto glen_n = config->get_number("stress_balance.ssa.Glen_exponent");
188
189 // Determine the kind of solver to use.
190 SSAFactory ssafactory = NULL;
191 if (method == "fem") {
192 ssafactory = SSAFEMFactory;
193 } else if (method == "fd") {
194 ssafactory = SSAFDFactory;
195 } else {
196 /* can't happen */
197 }
198
199 SSATestCasePlug testcase(ctx, Mx, My, glen_n, ssafactory);
200 testcase.init();
201 testcase.run();
202 testcase.report("plug");
203 testcase.write(output_file);
204 }
205 catch (...) {
206 handle_fatal_errors(com);
207 return 1;
208 }
209
210 return 0;
211 }