tssa_test_cfbc.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_cfbc.cc (7038B)
---
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_TESTCFBC\n"
21 " Testing program for PISM's implementations of the SSA.\n"
22 " Does a time-independent calculation. Does not run IceModel or a derived\n"
23 " class thereof. Uses the van der Veen flow-line shelf geometry. Also may be\n"
24 " used in a PISM 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/SSAFD_diagnostics.hh"
29 #include "pism/stressbalance/ssa/SSATestCase.hh"
30 #include "pism/stressbalance/ssa/SSAFEM.hh"
31 #include "pism/util/Mask.hh"
32 #include "pism/util/Context.hh"
33 #include "pism/util/VariableMetadata.hh"
34 #include "pism/util/error_handling.hh"
35 #include "pism/util/iceModelVec.hh"
36 #include "pism/util/io/File.hh"
37 #include "pism/util/petscwrappers/PetscInitializer.hh"
38 #include "pism/util/pism_utilities.hh"
39 #include "pism/util/pism_options.hh"
40
41 namespace pism {
42 namespace stressbalance {
43
44 // thickness profile in the van der Veen solution
45 static double H_exact(double V0, double H0, double C, double x) {
46 const double Q0 = V0*H0;
47 return pow(4 * C / Q0 * x + 1/pow(H0, 4), -0.25);
48 }
49
50 // velocity profile; corresponds to constant flux
51 static double u_exact(double V0, double H0, double C, double x) {
52 const double Q0 = V0*H0;
53 return Q0 / H_exact(V0, H0, C, x);
54 }
55
56 class SSATestCaseCFBC: public SSATestCase {
57 public:
58 SSATestCaseCFBC(Context::Ptr ctx, int Mx, int My, SSAFactory ssafactory)
59 : SSATestCase(ctx, Mx, My, 250e3, 250e3, CELL_CENTER, Y_PERIODIC) {
60 V0 = units::convert(ctx->unit_system(), 300.0, "m year-1", "m second-1");
61 H0 = 600.0; // meters
62 C = 2.45e-18;
63
64 m_config->set_number("flow_law.isothermal_Glen.ice_softness",
65 pow(1.9e8, -m_config->get_number("stress_balance.ssa.Glen_exponent")));
66 m_config->set_flag("stress_balance.ssa.compute_surface_gradient_inward", false);
67 m_config->set_flag("stress_balance.calving_front_stress_bc", true);
68 m_config->set_string("stress_balance.ssa.flow_law", "isothermal_glen");
69
70 m_enthalpyconverter = EnthalpyConverter::Ptr(new EnthalpyConverter(*m_config));
71
72 m_ssa = ssafactory(m_grid);
73 }
74
75 virtual void write_nuH(const std::string &filename);
76
77 protected:
78 virtual void initializeSSACoefficients();
79
80 virtual void exactSolution(int i, int j,
81 double x, double y, double *u, double *v);
82
83 double V0, //!< grounding line vertically-averaged velocity
84 H0, //!< grounding line thickness (meters)
85 C; //!< "typical constant ice parameter"
86 };
87
88 void SSATestCaseCFBC::write_nuH(const std::string &filename) {
89
90 SSAFD *ssafd = dynamic_cast<SSAFD*>(m_ssa);
91 if (ssafd != NULL) {
92 SSAFD_nuH(ssafd).compute()->write(filename);
93 }
94 }
95
96 void SSATestCaseCFBC::initializeSSACoefficients() {
97
98 m_tauc.set(0.0); // irrelevant
99 m_geometry.bed_elevation.set(-1000.0); // assures shelf is floating
100
101 double enth0 = m_enthalpyconverter->enthalpy(273.15, 0.01, 0.0); // 0.01 water fraction
102 m_ice_enthalpy.set(enth0);
103
104 IceModelVec::AccessList list{&m_geometry.ice_thickness,
105 &m_geometry.ice_surface_elevation, &m_bc_mask, &m_bc_values, &m_geometry.cell_type};
106
107 double ocean_rho = m_config->get_number("constants.sea_water.density"),
108 ice_rho = m_config->get_number("constants.ice.density");
109
110 const double x_min = m_grid->x(0);
111
112 for (Points p(*m_grid); p; p.next()) {
113 const int i = p.i(), j = p.j();
114
115 const double x = m_grid->x(i);
116
117 if (i != (int)m_grid->Mx() - 1) {
118 m_geometry.ice_thickness(i, j) = H_exact(V0, H0, C, x - x_min);
119 m_geometry.cell_type(i, j) = MASK_FLOATING;
120 } else {
121 m_geometry.ice_thickness(i, j) = 0;
122 m_geometry.cell_type(i, j) = MASK_ICE_FREE_OCEAN;
123 }
124
125 m_geometry.ice_surface_elevation(i,j) = (1.0 - ice_rho / ocean_rho) * m_geometry.ice_thickness(i, j);
126
127 if (i == 0) {
128 m_bc_mask(i, j) = 1;
129 m_bc_values(i, j).u = V0;
130 m_bc_values(i, j).v = 0;
131 } else {
132 m_bc_mask(i, j) = 0;
133 m_bc_values(i, j).u = 0;
134 m_bc_values(i, j).v = 0;
135 }
136 }
137
138
139 // communicate what we have set
140 m_geometry.ice_surface_elevation.update_ghosts();
141
142 m_geometry.ice_thickness.update_ghosts();
143
144 m_bc_mask.update_ghosts();
145
146 m_geometry.cell_type.update_ghosts();
147
148 m_bc_values.update_ghosts();
149 }
150
151 void SSATestCaseCFBC::exactSolution(int i, int /*j*/,
152 double x, double /*y*/,
153 double *u, double *v) {
154 const double x_min = m_grid->x(0);
155
156 if (i != (int)m_grid->Mx() - 1) {
157 *u = u_exact(V0, H0, C, x - x_min);
158 } else {
159 *u = 0;
160 }
161
162 *v = 0;
163 }
164
165 } // end of namespace stressbalance
166 } // end of namespace pism
167
168 int main(int argc, char *argv[]) {
169
170 using namespace pism;
171 using namespace pism::stressbalance;
172
173 MPI_Comm com = MPI_COMM_WORLD;
174 petsc::Initializer petsc(argc, argv, help);
175
176 com = PETSC_COMM_WORLD;
177
178 /* This explicit scoping forces destructors to be called before PetscFinalize() */
179 try {
180 Context::Ptr ctx = context_from_options(com, "ssa_test_cfbc");
181 Config::Ptr config = ctx->config();
182
183 std::string usage = "\n"
184 "usage of SSA_TEST_CFBC:\n"
185 " run ssa_test_cfbc -Mx <number> -My <number>\n"
186 "\n";
187
188 bool stop = show_usage_check_req_opts(*ctx->log(), "ssa_test_cfbc", {}, usage);
189
190 if (stop) {
191 return 0;
192 }
193
194 // Parameters that can be overridden by command line options
195 unsigned int Mx = config->get_number("grid.Mx");
196 unsigned int My = config->get_number("grid.My");
197
198 auto method = config->get_string("stress_balance.ssa.method");
199 auto output_file = config->get_string("output.file_name");
200
201 // Determine the kind of solver to use.
202 SSAFactory ssafactory = NULL;
203 if (method == "fem") {
204 ssafactory = SSAFEMFactory;
205 } else if (method == "fd") {
206 ssafactory = SSAFDFactory;
207 } else {
208 /* can't happen */
209 }
210
211 SSATestCaseCFBC testcase(ctx, Mx, My, ssafactory);
212 testcase.init();
213 testcase.run();
214 testcase.report("V");
215 testcase.write(output_file);
216 testcase.write_nuH(output_file);
217 }
218 catch (...) {
219 handle_fatal_errors(com);
220 return 1;
221 }
222
223 return 0;
224 }