tssa_test_linear.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_linear.cc (6486B)
---
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: linear flow. The rheology is
20 linear (i.e. n=1 in the Glen flow law) and the basal shear stress is also
21 linear viscous flow. The geometry consists of a constant surface slope in
22 the positive x-direction, and dirichlet conditions leading to an exponential
23 solution are imposed along the entire boundary.
24 */
25
26
27 static char help[] =
28 "\nSSA_TEST_EXP\n"
29 " Testing program for the finite element implementation of the SSA.\n"
30 " Does a time-independent calculation. Does not run IceModel or a derived\n"
31 " class thereof.Also may be used in a PISM\n"
32 " software (regression) test.\n\n";
33
34 #include <cmath>
35
36 #include "pism/basalstrength/basal_resistance.hh" // IceBasalResistancePlasticLaw
37 #include "pism/stressbalance/ssa/SSAFD.hh"
38 #include "pism/stressbalance/ssa/SSAFEM.hh"
39 #include "pism/stressbalance/ssa/SSATestCase.hh"
40 #include "pism/util/Context.hh"
41 #include "pism/util/VariableMetadata.hh"
42 #include "pism/util/error_handling.hh"
43 #include "pism/util/iceModelVec.hh"
44 #include "pism/util/io/File.hh"
45 #include "pism/util/petscwrappers/PetscInitializer.hh"
46 #include "pism/util/pism_utilities.hh"
47 #include "pism/util/pism_options.hh"
48 #include "pism/verification/tests/exactTestsIJ.h"
49
50 namespace pism {
51 namespace stressbalance {
52
53 class SSATestCaseExp: public SSATestCase
54 {
55 public:
56 SSATestCaseExp(Context::Ptr ctx, int Mx, int My, SSAFactory ssafactory)
57 : SSATestCase(ctx, Mx, My, 50e3, 50e3, CELL_CORNER, NOT_PERIODIC) {
58 L = units::convert(ctx->unit_system(), 50, "km", "m"); // 50km half-width
59 H0 = 500; // meters
60 dhdx = 0.005; // pure number
61 nu0 = units::convert(ctx->unit_system(), 30.0, "MPa year", "Pa s");
62 tauc0 = 1.e4; // 1kPa
63
64 // Use a pseudo-plastic law with linear till
65 m_config->set_flag("basal_resistance.pseudo_plastic.enabled", true);
66 m_config->set_number("basal_resistance.pseudo_plastic.q", 1.0);
67
68 // The following is irrelevant because we will force linear rheology later.
69 m_enthalpyconverter = EnthalpyConverter::Ptr(new EnthalpyConverter(*m_config));
70
71 m_ssa = ssafactory(m_grid);
72 }
73
74 protected:
75 virtual void initializeSSACoefficients();
76
77 virtual void exactSolution(int i, int j,
78 double x, double y, double *u, double *v);
79
80 double L, H0, dhdx, nu0, tauc0;
81 };
82
83 void SSATestCaseExp::initializeSSACoefficients() {
84
85 // Force linear rheology
86 m_ssa->strength_extension->set_notional_strength(nu0 * H0);
87 m_ssa->strength_extension->set_min_thickness(4000*10);
88
89 // The finite difference code uses the following flag to treat the non-periodic grid correctly.
90 m_config->set_flag("stress_balance.ssa.compute_surface_gradient_inward", true);
91
92 // Set constants for most coefficients.
93 m_geometry.ice_thickness.set(H0);
94 m_geometry.ice_surface_elevation.set(H0);
95 m_geometry.bed_elevation.set(0.0);
96 m_tauc.set(tauc0);
97
98
99 // Set boundary conditions (Dirichlet all the way around).
100 m_bc_mask.set(0.0);
101
102 IceModelVec::AccessList list{&m_bc_values, &m_bc_mask};
103
104 for (Points p(*m_grid); p; p.next()) {
105 const int i = p.i(), j = p.j();
106
107 double myu, myv;
108 const double myx = m_grid->x(i), myy=m_grid->y(j);
109
110 bool edge = ((j == 0) || (j == (int)m_grid->My() - 1) ||
111 (i == 0) || (i == (int)m_grid->Mx() - 1));
112 if (edge) {
113 m_bc_mask(i,j) = 1;
114 exactSolution(i,j,myx,myy,&myu,&myv);
115 m_bc_values(i,j).u = myu;
116 m_bc_values(i,j).v = myv;
117 }
118 }
119
120 m_bc_values.update_ghosts();
121 m_bc_mask.update_ghosts();
122 }
123
124
125 void SSATestCaseExp::exactSolution(int /*i*/, int /*j*/, double x, double /*y*/,
126 double *u, double *v) {
127 double tauc_threshold_velocity = m_config->get_number("basal_resistance.pseudo_plastic.u_threshold",
128 "m second-1");
129 double v0 = units::convert(m_sys, 100.0, "m year-1", "m second-1");
130 // double alpha=log(2.)/(2*L);
131 double alpha = sqrt((tauc0/tauc_threshold_velocity) / (4*nu0*H0));
132 *u = v0*exp(-alpha*(x-L));
133 *v = 0;
134 }
135
136 } // end of namespace stressbalance
137 } // end of namespace pism
138
139 int main(int argc, char *argv[]) {
140
141 using namespace pism;
142 using namespace pism::stressbalance;
143
144 MPI_Comm com = MPI_COMM_WORLD; // won't be used except for rank,size
145 petsc::Initializer petsc(argc, argv, help);
146
147 com = PETSC_COMM_WORLD;
148
149 /* This explicit scoping forces destructors to be called before PetscFinalize() */
150 try {
151 Context::Ptr ctx = context_from_options(com, "ssa_test_linear");
152 Config::Ptr config = ctx->config();
153
154 std::string usage = "\n"
155 "usage:\n"
156 " run ssa_test_linear -Mx <number> -My <number> -ssa_method <fd|fem>\n"
157 "\n";
158
159 bool stop = show_usage_check_req_opts(*ctx->log(), "ssa_test_linear", {}, usage);
160
161 if (stop) {
162 return 0;
163 }
164
165 // Parameters that can be overridden by command line options
166 unsigned int Mx = config->get_number("grid.Mx");
167 unsigned int My = config->get_number("grid.My");
168
169 auto method = config->get_string("stress_balance.ssa.method");
170 auto output_file = config->get_string("output.file_name");
171
172 // Determine the kind of solver to use.
173 SSAFactory ssafactory = NULL;
174 if (method == "fem") {
175 ssafactory = SSAFEMFactory;
176 } else if (method == "fd") {
177 ssafactory = SSAFDFactory;
178 } else {
179 /* can't happen */
180 }
181
182 SSATestCaseExp testcase(ctx, Mx, My, ssafactory);
183 testcase.init();
184 testcase.run();
185 testcase.report("linear");
186 testcase.write(output_file);
187 }
188 catch (...) {
189 handle_fatal_errors(com);
190 return 1;
191 }
192
193 return 0;
194 }