tsiafd_test.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
---
tsiafd_test.cc (14519B)
---
1 // Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 Ed Bueler and Constantine Khroulev
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 "\nSIAFD_TEST\n"
21 " Testing program for SIA, time-independent calculations separate from\n"
22 " IceModel. Uses verification test F. Also may be used in a PISM software"
23 "(regression) test.\n\n";
24
25 #include "SIAFD.hh"
26 #include "pism/basalstrength/basal_resistance.hh"
27 #include "pism/util/EnthalpyConverter.hh"
28 #include "pism/rheology/PatersonBuddCold.hh"
29 #include "pism/stressbalance/StressBalance.hh"
30 #include "pism/stressbalance/SSB_Modifier.hh"
31 #include "pism/stressbalance/ShallowStressBalance.hh"
32 #include "pism/util/IceGrid.hh"
33 #include "pism/util/Mask.hh"
34 #include "pism/util/Context.hh"
35 #include "pism/util/Time.hh"
36 #include "pism/util/VariableMetadata.hh"
37 #include "pism/util/error_handling.hh"
38 #include "pism/util/iceModelVec.hh"
39 #include "pism/util/io/File.hh"
40 #include "pism/util/petscwrappers/PetscInitializer.hh"
41 #include "pism/util/pism_utilities.hh"
42 #include "pism/util/pism_options.hh"
43 #include "pism/verification/tests/exactTestsFG.hh"
44 #include "pism/util/io/io_helpers.hh"
45 #include "pism/geometry/Geometry.hh"
46
47 namespace pism {
48
49 static void compute_strain_heating_errors(const IceModelVec3 &strain_heating,
50 const IceModelVec2S &thickness,
51 const IceGrid &grid,
52 double &gmax_strain_heating_err,
53 double &gav_strain_heating_err) {
54 double max_strain_heating_error = 0.0, av_strain_heating_error = 0.0, avcount = 0.0;
55 const int Mz = grid.Mz();
56
57 const double LforFG = 750000; // m
58
59 const double
60 ice_rho = grid.ctx()->config()->get_number("constants.ice.density"),
61 ice_c = grid.ctx()->config()->get_number("constants.ice.specific_heat_capacity");
62
63 IceModelVec::AccessList list{&thickness, &strain_heating};
64
65 ParallelSection loop(grid.com);
66 try {
67 for (Points p(grid); p; p.next()) {
68 const int i = p.i(), j = p.j();
69
70 double
71 xx = grid.x(i),
72 yy = grid.y(j),
73 r = sqrt(PetscSqr(xx) + PetscSqr(yy));
74
75 if ((r >= 1.0) && (r <= LforFG - 1.0)) {
76 // only evaluate error if inside sheet and not at central
77 // singularity
78 TestFGParameters F = exactFG(0.0, r, grid.z(), 0.0);
79
80 for (int k = 0; k < Mz; k++) {
81 F.Sig[k] *= ice_rho * ice_c; // scale exact strain_heating to J/(s m^3)
82 }
83 const int ks = grid.kBelowHeight(thickness(i, j));
84 const double *strain_heating_ij = strain_heating.get_column(i, j);
85 for (int k = 0; k < ks; k++) { // only eval error if below num surface
86 const double _strain_heating_error = fabs(strain_heating_ij[k] - F.Sig[k]);
87 max_strain_heating_error = std::max(max_strain_heating_error, _strain_heating_error);
88 avcount += 1.0;
89 av_strain_heating_error += _strain_heating_error;
90 }
91 }
92 }
93 } catch (...) {
94 loop.failed();
95 }
96 loop.check();
97
98 gmax_strain_heating_err = GlobalMax(grid.com, max_strain_heating_error);
99 gav_strain_heating_err = GlobalSum(grid.com, av_strain_heating_error);
100 double gavcount = GlobalSum(grid.com, avcount);
101 gav_strain_heating_err = gav_strain_heating_err/std::max(gavcount,1.0); // avoid div by zero
102 }
103
104
105 static void computeSurfaceVelocityErrors(const IceGrid &grid,
106 const IceModelVec2S &ice_thickness,
107 const IceModelVec3 &u3,
108 const IceModelVec3 &v3,
109 const IceModelVec3 &w3,
110 double &gmaxUerr,
111 double &gavUerr,
112 double &gmaxWerr,
113 double &gavWerr) {
114 double maxUerr = 0.0, maxWerr = 0.0, avUerr = 0.0, avWerr = 0.0;
115
116 const double LforFG = 750000; // m
117
118 IceModelVec::AccessList list{&ice_thickness, &u3, &v3, &w3};
119
120 for (Points p(grid); p; p.next()) {
121 const int i = p.i(), j = p.j();
122
123 double xx = grid.x(i), yy = grid.y(j),
124 r = sqrt(PetscSqr(xx) + PetscSqr(yy));
125 if ((r >= 1.0) && (r <= LforFG - 1.0)) { // only evaluate error if inside sheet
126 // and not at central singularity
127 const double H = ice_thickness(i, j);
128 std::vector<double> z(1, H);
129 TestFGParameters F = exactFG(0.0, r, z, 0.0);
130
131 const double uex = (xx/r) * F.U[0];
132 const double vex = (yy/r) * F.U[0];
133 // note that because getValZ does linear interpolation and H is not exactly at
134 // a grid point, this causes nonzero errors
135 const double Uerr = sqrt(PetscSqr(u3.getValZ(i,j,H) - uex) +
136 PetscSqr(v3.getValZ(i,j,H) - vex));
137 maxUerr = std::max(maxUerr,Uerr);
138 avUerr += Uerr;
139 const double Werr = fabs(w3.getValZ(i,j,H) - F.w[0]);
140 maxWerr = std::max(maxWerr,Werr);
141 avWerr += Werr;
142 }
143 }
144
145 gmaxUerr = GlobalMax(grid.com, maxUerr);
146 gmaxWerr = GlobalMax(grid.com, maxWerr);
147 gavUerr = GlobalSum(grid.com, avUerr);
148 gavUerr = gavUerr/(grid.Mx()*grid.My());
149 gavWerr = GlobalSum(grid.com, avWerr);
150 gavWerr = gavWerr/(grid.Mx()*grid.My());
151 }
152
153
154 static void enthalpy_from_temperature_cold(EnthalpyConverter &EC,
155 const IceGrid &grid,
156 const IceModelVec2S &thickness,
157 const IceModelVec3 &temperature,
158 IceModelVec3 &enthalpy) {
159
160 IceModelVec::AccessList list{&temperature, &enthalpy, &thickness};
161
162 for (Points p(grid); p; p.next()) {
163 const int i = p.i(), j = p.j();
164
165 const double *T_ij = temperature.get_column(i,j);
166 double *E_ij = enthalpy.get_column(i,j);
167
168 for (unsigned int k=0; k<grid.Mz(); ++k) {
169 double depth = thickness(i,j) - grid.z(k);
170 E_ij[k] = EC.enthalpy_permissive(T_ij[k], 0.0,
171 EC.pressure(depth));
172 }
173
174 }
175
176 enthalpy.update_ghosts();
177 }
178
179
180 //! \brief Set the test F initial state.
181 static void setInitStateF(IceGrid &grid,
182 EnthalpyConverter &EC,
183 IceModelVec2S &bed,
184 IceModelVec2Int &mask,
185 IceModelVec2S &surface,
186 IceModelVec2S &thickness,
187 IceModelVec3 &enthalpy) {
188
189 double
190 ST = 1.67e-5,
191 Tmin = 223.15, // K
192 LforFG = 750000; // m
193
194 bed.set(0.0);
195 mask.set(MASK_GROUNDED);
196
197 IceModelVec::AccessList list{&thickness, &enthalpy};
198
199 for (Points p(grid); p; p.next()) {
200 const int i = p.i(), j = p.j();
201
202 const double
203 r = std::max(radius(grid, i, j), 1.0), // avoid singularity at origin
204 Ts = Tmin + ST * r;
205
206 if (r > LforFG - 1.0) { // if (essentially) outside of sheet
207 thickness(i, j) = 0.0;
208 enthalpy.set_column(i, j, Ts);
209 } else {
210 TestFGParameters F = exactFG(0.0, r, grid.z(), 0.0);
211
212 thickness(i, j) = F.H;
213 enthalpy.set_column(i, j, &F.T[0]);
214 }
215 }
216
217
218 thickness.update_ghosts();
219
220 surface.copy_from(thickness);
221
222 enthalpy_from_temperature_cold(EC, grid, thickness, enthalpy,
223 enthalpy);
224 }
225
226 static void reportErrors(const IceGrid &grid,
227 units::System::Ptr unit_system,
228 const IceModelVec2S &thickness,
229 const IceModelVec3 &u_sia,
230 const IceModelVec3 &v_sia,
231 const IceModelVec3 &w_sia,
232 const IceModelVec3 &strain_heating) {
233
234 Logger::ConstPtr log = grid.ctx()->log();
235
236 // strain_heating errors if appropriate; reported in 10^6 J/(s m^3)
237 double max_strain_heating_error, av_strain_heating_error;
238 compute_strain_heating_errors(strain_heating, thickness,
239 grid,
240 max_strain_heating_error, av_strain_heating_error);
241
242 log->message(1,
243 "Sigma : maxSig avSig\n");
244 log->message(1, " %12.6f%12.6f\n",
245 max_strain_heating_error*1.0e6, av_strain_heating_error*1.0e6);
246
247 // surface velocity errors if exact values are available; reported in m year-1
248 double maxUerr, avUerr, maxWerr, avWerr;
249 computeSurfaceVelocityErrors(grid, thickness,
250 u_sia,
251 v_sia,
252 w_sia,
253 maxUerr, avUerr,
254 maxWerr, avWerr);
255
256 log->message(1,
257 "surf vels : maxUvec avUvec maxW avW\n");
258 log->message(1, " %12.6f%12.6f%12.6f%12.6f\n",
259 units::convert(unit_system, maxUerr, "m second-1", "m year-1"),
260 units::convert(unit_system, avUerr, "m second-1", "m year-1"),
261 units::convert(unit_system, maxWerr, "m second-1", "m year-1"),
262 units::convert(unit_system, avWerr, "m second-1", "m year-1"));
263 }
264
265 } // end of namespace pism
266
267 int main(int argc, char *argv[]) {
268
269 using namespace pism;
270 using namespace pism::stressbalance;
271
272 MPI_Comm com = MPI_COMM_WORLD;
273 petsc::Initializer petsc(argc, argv, help);
274
275 com = PETSC_COMM_WORLD;
276
277 /* This explicit scoping forces destructors to be called before PetscFinalize() */
278 try {
279 // set default verbosity
280 Context::Ptr ctx = context_from_options(com, "siafd_test");
281 Config::Ptr config = ctx->config();
282
283 config->set_flag("stress_balance.sia.grain_size_age_coupling", false);
284 config->set_string("stress_balance.sia.flow_law", "arr");
285
286 set_config_from_options(*config);
287
288 std::string usage = "\n"
289 "usage of SIAFD_TEST:\n"
290 " run siafd_test -Mx <number> -My <number> -Mz <number> -o foo.nc\n"
291 "\n";
292
293 bool stop = show_usage_check_req_opts(*ctx->log(), "siafd_test", {}, usage);
294
295 if (stop) {
296 return 0;
297 }
298
299 auto output_file = config->get_string("output.file_name");
300
301 GridParameters P(config);
302 P.Lx = 900e3;
303 P.Ly = P.Lx;
304 P.horizontal_size_from_options();
305
306 double Lz = 4000.0;
307 unsigned int Mz = config->get_number("grid.Mz");
308
309 P.z = IceGrid::compute_vertical_levels(Lz, Mz, EQUAL);
310 P.ownership_ranges_from_options(ctx->size());
311
312 // create grid and set defaults
313 IceGrid::Ptr grid(new IceGrid(ctx, P));
314 grid->report_parameters();
315
316 EnthalpyConverter::Ptr EC(new ColdEnthalpyConverter(*config));
317
318 const int WIDE_STENCIL = config->get_number("grid.max_stencil_width");
319
320 IceModelVec3
321 enthalpy(grid, "enthalpy", WITH_GHOSTS, WIDE_STENCIL),
322 age(grid, "age", WITHOUT_GHOSTS);
323
324 Geometry geometry(grid);
325 geometry.sea_level_elevation.set(0.0);
326
327 // age of the ice; is not used here
328 age.set_attrs("diagnostic", "age of the ice", "s", "s", "", 0);
329 age.set(0.0);
330
331 // enthalpy in the ice
332 enthalpy.set_attrs("model_state",
333 "ice enthalpy (includes sensible heat, latent heat, pressure)",
334 "J kg-1", "J kg-1", "", 0);
335 //
336 enthalpy.set(EC->enthalpy(263.15, 0.0, EC->pressure(1000.0)));
337
338 // Create the SIA solver object:
339
340 // We use SIA_Nonsliding and not SIAFD here because we need the z-component
341 // of the ice velocity, which is computed using incompressibility of ice in
342 // StressBalance::compute_vertical_velocity().
343 SIAFD *sia = new SIAFD(grid);
344 ZeroSliding *no_sliding = new ZeroSliding(grid);
345
346 // stress_balance will de-allocate no_sliding and sia.
347 StressBalance stress_balance(grid, no_sliding, sia);
348
349 // fill the fields:
350 setInitStateF(*grid, *EC,
351 geometry.bed_elevation,
352 geometry.cell_type,
353 geometry.ice_surface_elevation,
354 geometry.ice_thickness,
355 enthalpy);
356
357 geometry.ensure_consistency(config->get_number("geometry.ice_free_thickness_standard"));
358
359 // Initialize the SIA solver:
360 stress_balance.init();
361
362 IceModelVec2S melange_back_pressure;
363 melange_back_pressure.create(grid, "melange_back_pressure", WITHOUT_GHOSTS);
364 melange_back_pressure.set_attrs("boundary_condition",
365 "melange back pressure fraction", "", "", "", 0);
366 melange_back_pressure.set(0.0);
367
368 bool full_update = true;
369
370 stressbalance::Inputs inputs;
371 inputs.geometry = &geometry;
372 inputs.melange_back_pressure = &melange_back_pressure;
373 inputs.enthalpy = &enthalpy;
374 inputs.age = &age;
375
376 stress_balance.update(inputs, full_update);
377
378 // Report errors relative to the exact solution:
379 const IceModelVec3
380 &u3 = stress_balance.velocity_u(),
381 &v3 = stress_balance.velocity_v(),
382 &w3 = stress_balance.velocity_w();
383
384 const IceModelVec3 &sigma = stress_balance.volumetric_strain_heating();
385
386 reportErrors(*grid, ctx->unit_system(),
387 geometry.ice_thickness, u3, v3, w3, sigma);
388
389 // Write results to an output file:
390 File file(grid->com, output_file, PISM_NETCDF3, PISM_READWRITE_MOVE);
391 io::define_time(file, *ctx);
392 io::append_time(file, *ctx->config(), ctx->time()->current());
393
394 geometry.ice_surface_elevation.write(file);
395 geometry.ice_thickness.write(file);
396 geometry.cell_type.write(file);
397 geometry.bed_elevation.write(file);
398
399 sia->diffusivity().write(file);
400 u3.write(file);
401 v3.write(file);
402 w3.write(file);
403 sigma.write(file);
404 }
405 catch (...) {
406 handle_fatal_errors(com);
407 return 1;
408 }
409
410 return 0;
411 }