ticeCompModel.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
---
ticeCompModel.cc (32417B)
---
1 // Copyright (C) 2004-2019 Jed Brown, 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 #include <cmath>
20 #include <cstring>
21
22 #include <vector> // STL vector container; sortable; used in test L
23 #include <algorithm> // required by sort(...) in test L
24
25 #include "tests/exactTestsABCD.h"
26 #include "tests/exactTestsFG.hh"
27 #include "tests/exactTestH.h"
28 #include "tests/exactTestL.hh"
29
30 #include "iceCompModel.hh"
31 #include "pism/stressbalance/sia/SIAFD.hh"
32 #include "pism/stressbalance/ShallowStressBalance.hh"
33 #include "pism/rheology/PatersonBuddCold.hh"
34 #include "pism/stressbalance/StressBalance.hh"
35 #include "pism/util/EnthalpyConverter.hh"
36 #include "pism/util/io/File.hh"
37 #include "pism/util/pism_options.hh"
38 #include "pism/coupler/ocean/Constant.hh"
39 #include "pism/coupler/SeaLevel.hh"
40 #include "PSVerification.hh"
41 #include "pism/util/Mask.hh"
42 #include "pism/util/error_handling.hh"
43 #include "pism/earth/BedDef.hh"
44 #include "pism/util/IceGrid.hh"
45 #include "pism/util/Time.hh"
46 #include "pism/util/ConfigInterface.hh"
47 #include "pism/util/Context.hh"
48 #include "pism/util/io/io_helpers.hh"
49 #include "pism/util/Logger.hh"
50 #include "pism/util/pism_utilities.hh"
51 #include "BTU_Verification.hh"
52 #include "pism/energy/BTU_Minimal.hh"
53 #include "TemperatureModel_Verification.hh"
54
55 namespace pism {
56
57 using units::convert;
58
59 IceCompModel::IceCompModel(IceGrid::Ptr g, Context::Ptr context, int mytest)
60 : IceModel(g, context), m_testname(mytest), m_bedrock_is_ice_forK(false) {
61
62 m_log->message(2, "starting Test %c ...\n", m_testname);
63
64 // Override some defaults from parent class
65 m_config->set_number("stress_balance.sia.enhancement_factor", 1.0);
66 // none use bed smoothing & bed roughness parameterization
67 m_config->set_number("stress_balance.sia.bed_smoother.range", 0.0);
68
69 // set values of flags in run()
70 m_config->set_flag("geometry.update.enabled", true);
71 m_config->set_flag("geometry.update.use_basal_melt_rate", false);
72
73 // flow law settings
74 switch (m_testname) {
75 case 'A':
76 case 'B':
77 case 'C':
78 case 'D':
79 case 'H':
80 case 'L':
81 {
82 m_config->set_string("stress_balance.sia.flow_law", "isothermal_glen");
83 const double year = convert(m_sys, 1.0, "year", "seconds");
84 m_config->set_number("flow_law.isothermal_Glen.ice_softness", 1.0e-16 / year);
85 break;
86 }
87 case 'V':
88 {
89 m_config->set_string("stress_balance.ssa.flow_law", "isothermal_glen");
90 const double
91 hardness = 1.9e8,
92 softness = pow(hardness,
93 -m_config->get_number("stress_balance.ssa.Glen_exponent"));
94 m_config->set_number("flow_law.isothermal_Glen.ice_softness", softness);
95 break;
96 }
97 case 'F':
98 case 'G':
99 case 'K':
100 case 'O':
101 default:
102 {
103 m_config->set_string("stress_balance.sia.flow_law", "arr");
104 break;
105 }
106 }
107
108 if (m_testname == 'H') {
109 m_config->set_string("bed_deformation.model", "iso");
110 } else {
111 m_config->set_string("bed_deformation.model", "none");
112 }
113
114 if ((m_testname == 'F') || (m_testname == 'G') || (m_testname == 'K') || (m_testname == 'O')) {
115 m_config->set_flag("energy.enabled", true);
116 // essentially turn off run-time reporting of extremely low computed
117 // temperatures; *they will be reported as errors* anyway
118 m_config->set_number("energy.minimum_allowed_temperature", 0.0);
119 m_config->set_number("energy.max_low_temperature_count", 1000000);
120 } else {
121 m_config->set_flag("energy.enabled", false);
122 }
123
124 m_config->set_flag("ocean.always_grounded", true);
125
126 // special considerations for K and O wrt thermal bedrock and pressure-melting
127 if ((m_testname == 'K') || (m_testname == 'O')) {
128 m_config->set_flag("energy.allow_temperature_above_melting", false);
129 } else {
130 // note temps are generally allowed to go above pressure melting in verify
131 m_config->set_flag("energy.allow_temperature_above_melting", true);
132 }
133
134 if (m_testname == 'V') {
135 // no sub-shelf melting
136 m_config->set_flag("geometry.update.use_basal_melt_rate", false);
137
138 // this test is isothermal
139 m_config->set_flag("energy.enabled", false);
140
141 // use the SSA solver
142 m_config->set_string("stress_balance_model", "ssa");
143
144 // this certainly is not a "dry simulation"
145 m_config->set_flag("ocean.always_grounded", false);
146
147 m_config->set_flag("stress_balance.ssa.dirichlet_bc", true);
148 }
149
150 m_config->set_flag("energy.temperature_based", true);
151 }
152
153 void IceCompModel::allocate_storage() {
154
155 IceModel::allocate_storage();
156
157 m_HexactL.create(m_grid, "HexactL", WITH_GHOSTS, 2);
158
159 m_strain_heating3_comp.create(m_grid,"strain_heating_comp", WITHOUT_GHOSTS);
160 m_strain_heating3_comp.set_attrs("internal","rate of compensatory strain heating in ice",
161 "W m-3", "W m-3", "", 0);
162 }
163
164 void IceCompModel::allocate_bedrock_thermal_unit() {
165
166 if (m_btu != NULL) {
167 return;
168 }
169
170 // this switch changes Test K to make material properties for bedrock the same as for ice
171 bool biiSet = options::Bool("-bedrock_is_ice", "set bedrock properties to those of ice");
172 if (biiSet) {
173 if (m_testname == 'K') {
174 m_log->message(1,
175 "setting material properties of bedrock to those of ice in Test K\n");
176 m_config->set_number("energy.bedrock_thermal.density", m_config->get_number("constants.ice.density"));
177 m_config->set_number("energy.bedrock_thermal.conductivity", m_config->get_number("constants.ice.thermal_conductivity"));
178 m_config->set_number("energy.bedrock_thermal.specific_heat_capacity", m_config->get_number("constants.ice.specific_heat_capacity"));
179 m_bedrock_is_ice_forK = true;
180 } else {
181 m_log->message(1,
182 "IceCompModel WARNING: option -bedrock_is_ice ignored; only applies to Test K\n");
183 }
184 }
185
186 if (m_testname != 'K') {
187 // now make bedrock have same material properties as ice
188 // (note Mbz=1 also, by default, but want ice/rock interface to see
189 // pure ice from the point of view of applying geothermal boundary
190 // condition, especially in tests F and G)
191 m_config->set_number("energy.bedrock_thermal.density", m_config->get_number("constants.ice.density"));
192 m_config->set_number("energy.bedrock_thermal.conductivity", m_config->get_number("constants.ice.thermal_conductivity"));
193 m_config->set_number("energy.bedrock_thermal.specific_heat_capacity", m_config->get_number("constants.ice.specific_heat_capacity"));
194 }
195
196 energy::BTUGrid bed_vertical_grid = energy::BTUGrid::FromOptions(m_grid->ctx());
197
198 if (bed_vertical_grid.Mbz > 1) {
199 m_btu = new energy::BTU_Verification(m_grid, bed_vertical_grid,
200 m_testname, m_bedrock_is_ice_forK);
201 } else {
202 m_btu = new energy::BTU_Minimal(m_grid);
203 }
204
205 m_submodels["bedrock thermal layer"] = m_btu;
206 }
207
208 void IceCompModel::allocate_energy_model() {
209
210 if (m_energy_model != NULL) {
211 return;
212 }
213
214 m_log->message(2, "# Allocating an energy balance model...\n");
215
216 // this switch changes Test K to make material properties for bedrock the same as for ice
217 bool bedrock_is_ice = options::Bool("-bedrock_is_ice", "set bedrock properties to those of ice");
218 m_energy_model = new energy::TemperatureModel_Verification(m_grid, m_stress_balance.get(), m_testname, bedrock_is_ice);
219
220 m_submodels["energy balance model"] = m_energy_model;
221 }
222
223
224 void IceCompModel::allocate_bed_deformation() {
225
226 IceModel::allocate_bed_deformation();
227
228 // for simple isostasy
229 m_f = m_config->get_number("constants.ice.density") / m_config->get_number("bed_deformation.mantle_density");
230
231 std::string bed_def_model = m_config->get_string("bed_deformation.model");
232
233 if ((m_testname == 'H') && bed_def_model != "iso") {
234 m_log->message(1,
235 "IceCompModel WARNING: Test H should be run with option\n"
236 " '-bed_def iso' for the reported errors to be correct.\n");
237 }
238 }
239
240 void IceCompModel::allocate_couplers() {
241 EnthalpyConverter::Ptr EC = m_ctx->enthalpy_converter();
242
243 // Climate will always come from verification test formulas.
244 m_surface.reset(new surface::Verification(m_grid, EC, m_testname));
245 m_submodels["surface process model"] = m_surface.get();
246
247 m_ocean.reset(new ocean::Constant(m_grid));
248 m_submodels["ocean model"] = m_ocean.get();
249
250 m_sea_level.reset(new ocean::sea_level::SeaLevel(m_grid));
251 m_submodels["sea level forcing"] = m_sea_level.get();
252 }
253
254 void IceCompModel::bootstrap_2d(const File &input_file) {
255 (void) input_file;
256 throw RuntimeError(PISM_ERROR_LOCATION, "pismv (IceCompModel) does not support bootstrapping.");
257 }
258
259 void IceCompModel::initialize_2d() {
260 m_log->message(3, "initializing Test %c from formulas ...\n",m_testname);
261
262 m_geometry.bed_elevation.set(0.0);
263 m_geometry.sea_level_elevation.set(0.0);
264
265 IceModelVec2S uplift(m_grid, "uplift", WITHOUT_GHOSTS);
266 uplift.set(0.0);
267
268 m_beddef->bootstrap(m_geometry.bed_elevation,
269 uplift,
270 m_geometry.ice_thickness,
271 m_geometry.sea_level_elevation);
272
273 // Test-specific initialization:
274 switch (m_testname) {
275 case 'A':
276 case 'B':
277 case 'C':
278 case 'D':
279 case 'H':
280 initTestABCDH();
281 break;
282 case 'F':
283 case 'G':
284 initTestFG(); // see iCMthermo.cc
285 break;
286 case 'K':
287 case 'O':
288 initTestsKO(); // see iCMthermo.cc
289 break;
290 case 'L':
291 initTestL();
292 break;
293 case 'V':
294 test_V_init();
295 break;
296 default:
297 throw RuntimeError(PISM_ERROR_LOCATION, "Desired test not implemented by IceCompModel.");
298 }
299 }
300
301 void IceCompModel::initTestABCDH() {
302
303 const double time = m_time->current();
304
305 m_geometry.cell_type.set(MASK_GROUNDED);
306
307 IceModelVec::AccessList list(m_geometry.ice_thickness);
308
309 ParallelSection loop(m_grid->com);
310 try {
311 for (Points p(*m_grid); p; p.next()) {
312 const int i = p.i(), j = p.j();
313
314 const double r = radius(*m_grid, i, j);
315 switch (m_testname) {
316 case 'A':
317 m_geometry.ice_thickness(i, j) = exactA(r).H;
318 break;
319 case 'B':
320 m_geometry.ice_thickness(i, j) = exactB(time, r).H;
321 break;
322 case 'C':
323 m_geometry.ice_thickness(i, j) = exactC(time, r).H;
324 break;
325 case 'D':
326 m_geometry.ice_thickness(i, j) = exactD(time, r).H;
327 break;
328 case 'H':
329 m_geometry.ice_thickness(i, j) = exactH(m_f, time, r).H;
330 break;
331 default:
332 throw RuntimeError(PISM_ERROR_LOCATION, "test must be A, B, C, D, or H");
333 }
334 }
335 } catch (...) {
336 loop.failed();
337 }
338 loop.check();
339
340 m_geometry.ice_thickness.update_ghosts();
341
342 {
343 IceModelVec2S bed_uplift(m_grid, "uplift", WITHOUT_GHOSTS);
344 bed_uplift.set(0.0);
345
346 if (m_testname == 'H') {
347 m_geometry.bed_elevation.copy_from(m_geometry.ice_thickness);
348 m_geometry.bed_elevation.scale(-m_f);
349 } else { // flat bed case otherwise
350 m_geometry.bed_elevation.set(0.0);
351 }
352 m_geometry.sea_level_elevation.set(0.0);
353 m_beddef->bootstrap(m_geometry.bed_elevation, bed_uplift, m_geometry.ice_thickness,
354 m_geometry.sea_level_elevation);
355 }
356 }
357
358 //! Class used initTestL() in generating sorted list for ODE solver.
359 class rgrid {
360 public:
361 double r;
362 int i,j;
363 };
364
365 //! Comparison used initTestL() in generating sorted list for ODE solver.
366 struct rgridReverseSort {
367 bool operator()(rgrid a, rgrid b) {
368 return (a.r > b.r);
369 }
370 };
371
372 void IceCompModel::initTestL() {
373
374 m_log->message(2, "* Initializing ice thickness and bed topography for test L...\n");
375
376 // setup to evaluate test L; requires solving an ODE numerically
377 // using sorted list of radii, sorted in decreasing radius order
378 const int MM = m_grid->xm() * m_grid->ym();
379
380 std::vector<rgrid> rrv(MM);
381 int k = 0;
382 for (Points p(*m_grid); p; p.next()) {
383 const int i = p.i(), j = p.j();
384
385 rrv[k].i = i;
386 rrv[k].j = j;
387 rrv[k].r = radius(*m_grid, i,j);
388
389 k += 1;
390 }
391
392 std::sort(rrv.begin(), rrv.end(), rgridReverseSort()); // so rrv[k].r > rrv[k+1].r
393
394 // get soln to test L at these radii; solves ODE only once (on each processor)
395 std::vector<double> rr(MM), HH(MM), bb(MM), aa(MM);
396
397 for (k = 0; k < MM; k++) {
398 rr[k] = rrv[k].r;
399 }
400
401 ExactLParameters L = exactL(rr);
402
403 {
404 IceModelVec2S bed_uplift(m_grid, "uplift", WITHOUT_GHOSTS);
405
406 IceModelVec::AccessList list{&m_geometry.ice_thickness, &m_geometry.bed_elevation};
407
408 for (k = 0; k < MM; k++) {
409 m_geometry.ice_thickness(rrv[k].i, rrv[k].j) = L.H[k];
410 m_geometry.bed_elevation(rrv[k].i, rrv[k].j) = L.b[k];
411 }
412
413 m_geometry.ice_thickness.update_ghosts();
414
415 bed_uplift.set(0.0);
416
417 m_beddef->bootstrap(m_geometry.bed_elevation, bed_uplift, m_geometry.ice_thickness,
418 m_geometry.sea_level_elevation);
419 }
420
421 // store copy of ice_thickness for evaluating geometry errors
422 m_HexactL.copy_from(m_geometry.ice_thickness);
423 }
424
425 //! \brief Tests A and E have a thickness B.C. (ice_thickness == 0 outside a circle of radius 750km).
426 void IceCompModel::reset_thickness_test_A() {
427 const double LforAE = 750e3; // m
428
429 IceModelVec::AccessList list(m_geometry.ice_thickness);
430
431 for (Points p(*m_grid); p; p.next()) {
432 const int i = p.i(), j = p.j();
433
434 if (radius(*m_grid, i, j) > LforAE) {
435 m_geometry.ice_thickness(i, j) = 0;
436 }
437 }
438
439 m_geometry.ice_thickness.update_ghosts();
440 }
441
442 void IceCompModel::computeGeometryErrors(double &gvolexact, double &gareaexact,
443 double &gdomeHexact, double &volerr,
444 double &areaerr, double &gmaxHerr,
445 double &gavHerr, double &gmaxetaerr,
446 double ¢erHerr) {
447 // compute errors in thickness, eta=thickness^{(2n+2)/n}, volume, area
448
449 const double time = m_time->current();
450 double
451 Hexact = 0.0,
452 vol = 0.0,
453 area = 0.0,
454 domeH = 0.0,
455 volexact = 0.0,
456 areaexact = 0.0,
457 domeHexact = 0.0;
458 double
459 Herr = 0.0,
460 avHerr = 0.0,
461 etaerr = 0.0;
462
463 IceModelVec::AccessList list(m_geometry.ice_thickness);
464 if (m_testname == 'L') {
465 list.add(m_HexactL);
466 }
467
468 double
469 seawater_density = m_config->get_number("constants.sea_water.density"),
470 ice_density = m_config->get_number("constants.ice.density"),
471 Glen_n = m_config->get_number("stress_balance.sia.Glen_exponent"),
472 standard_gravity = m_config->get_number("constants.standard_gravity");
473
474 // area of grid square in square km:
475 const double a = m_grid->dx() * m_grid->dy() * 1e-3 * 1e-3;
476 const double m = (2.0 * Glen_n + 2.0) / Glen_n;
477
478 ParallelSection loop(m_grid->com);
479 try {
480 for (Points p(*m_grid); p; p.next()) {
481 const int i = p.i(), j = p.j();
482
483 if (m_geometry.ice_thickness(i,j) > 0) {
484 area += a;
485 vol += a * m_geometry.ice_thickness(i,j) * 1e-3;
486 }
487 double xx = m_grid->x(i), r = radius(*m_grid, i,j);
488 switch (m_testname) {
489 case 'A':
490 Hexact = exactA(r).H;
491 break;
492 case 'B':
493 Hexact = exactB(time, r).H;
494 break;
495 case 'C':
496 Hexact = exactC(time, r).H;
497 break;
498 case 'D':
499 Hexact = exactD(time, r).H;
500 break;
501 case 'F':
502 if (r > m_LforFG - 1.0) { // outside of sheet
503 Hexact=0.0;
504 } else {
505 r=std::max(r,1.0);
506 std::vector<double> z(1, 0.0);
507 Hexact = exactFG(0.0, r, z, 0.0).H;
508 }
509 break;
510 case 'G':
511 if (r > m_LforFG -1.0) { // outside of sheet
512 Hexact=0.0;
513 } else {
514 r=std::max(r,1.0);
515 std::vector<double> z(1, 0.0);
516 Hexact = exactFG(time, r, z, m_ApforG).H;
517 }
518 break;
519 case 'H':
520 Hexact = exactH(m_f, time, r).H;
521 break;
522 case 'K':
523 case 'O':
524 Hexact = 3000.0;
525 break;
526 case 'L':
527 Hexact = m_HexactL(i,j);
528 break;
529 case 'V':
530 {
531 double
532 H0 = 600.0,
533 v0 = convert(m_sys, 300.0, "m year-1", "m second-1"),
534 Q0 = H0 * v0,
535 B0 = m_stress_balance->shallow()->flow_law()->hardness(0, 0),
536 C = pow(ice_density * standard_gravity * (1.0 - ice_density/seawater_density) / (4 * B0), 3);
537
538 Hexact = pow(4 * C / Q0 * xx + 1/pow(H0, 4), -0.25);
539 }
540 break;
541 default:
542 throw RuntimeError(PISM_ERROR_LOCATION, "test must be A, B, C, D, F, G, H, K, L, or O");
543 }
544
545 if (Hexact > 0) {
546 areaexact += a;
547 volexact += a * Hexact * 1e-3;
548 }
549 if (i == ((int)m_grid->Mx() - 1)/2 and
550 j == ((int)m_grid->My() - 1)/2) {
551 domeH = m_geometry.ice_thickness(i,j);
552 domeHexact = Hexact;
553 }
554 // compute maximum errors
555 Herr = std::max(Herr,fabs(m_geometry.ice_thickness(i,j) - Hexact));
556 etaerr = std::max(etaerr,fabs(pow(m_geometry.ice_thickness(i,j),m) - pow(Hexact,m)));
557 // add to sums for average errors
558 avHerr += fabs(m_geometry.ice_thickness(i,j) - Hexact);
559 }
560 } catch (...) {
561 loop.failed();
562 }
563 loop.check();
564
565 // globalize (find errors over all processors)
566 double gvol, garea, gdomeH;
567 gvolexact = GlobalSum(m_grid->com, volexact);
568 gdomeHexact = GlobalMax(m_grid->com, domeHexact);
569 gareaexact = GlobalSum(m_grid->com, areaexact);
570
571 gvol = GlobalSum(m_grid->com, vol);
572 garea = GlobalSum(m_grid->com, area);
573 volerr = fabs(gvol - gvolexact);
574 areaerr = fabs(garea - gareaexact);
575
576 gmaxHerr = GlobalMax(m_grid->com, Herr);
577 gavHerr = GlobalSum(m_grid->com, avHerr);
578 gavHerr = gavHerr/(m_grid->Mx()*m_grid->My());
579 gmaxetaerr = GlobalMax(m_grid->com, etaerr);
580
581 gdomeH = GlobalMax(m_grid->com, domeH);
582 centerHerr = fabs(gdomeH - gdomeHexact);
583 }
584
585 void IceCompModel::post_step_hook() {
586 if (m_testname == 'A') {
587 reset_thickness_test_A();
588 }
589 }
590
591
592 void IceCompModel::print_summary(bool /* tempAndAge */) {
593 // we always show a summary at every step
594 IceModel::print_summary(true);
595 }
596
597
598 void IceCompModel::reportErrors() {
599 // geometry errors to report (for all tests except K and O):
600 // -- max thickness error
601 // -- average (at each grid point on whole grid) thickness error
602 // -- max (thickness)^(2n+2)/n error
603 // -- volume error
604 // -- area error
605 // and temperature errors (for tests F & G & K & O):
606 // -- max T error over 3D domain of ice
607 // -- av T error over 3D domain of ice
608 // and basal temperature errors (for tests F & G):
609 // -- max basal temp error
610 // -- average (at each grid point on whole grid) basal temp error
611 // and bedrock temperature errors (for tests K & O):
612 // -- max Tb error over 3D domain of bedrock
613 // -- av Tb error over 3D domain of bedrock
614 // and strain-heating (Sigma) errors (for tests F & G):
615 // -- max Sigma error over 3D domain of ice (in 10^-3 K a^-1)
616 // -- av Sigma error over 3D domain of ice (in 10^-3 K a^-1)
617 // and basal melt rate error (for test O):
618 // -- max bmelt error over base of ice
619 // and surface velocity errors (for tests F & G):
620 // -- max |<us,vs> - <usex,vsex>| error
621 // -- av |<us,vs> - <usex,vsex>| error
622 // -- max ws error
623 // -- av ws error
624 // and basal sliding errors (for test E):
625 // -- max ub error
626 // -- max vb error
627 // -- max |<ub,vb> - <ubexact,vbexact>| error
628 // -- av |<ub,vb> - <ubexact,vbexact>| error
629
630 bool no_report = options::Bool("-no_report", "Don't report numerical errors");
631
632 if (no_report) {
633 return;
634 }
635
636 double maxbmelterr, minbmelterr, volexact, areaexact, domeHexact,
637 volerr, areaerr, maxHerr, avHerr, maxetaerr, centerHerr;
638 double maxTerr, avTerr, basemaxTerr, baseavTerr, basecenterTerr, maxTberr, avTberr;
639 double max_strain_heating_error, av_strain_heating_error;
640 double maxUerr, avUerr, maxWerr, avWerr;
641
642 const rheology::FlowLaw &flow_law = *m_stress_balance->modifier()->flow_law();
643 const double m = (2.0 * flow_law.exponent() + 2.0) / flow_law.exponent();
644
645 EnthalpyConverter::Ptr EC = m_ctx->enthalpy_converter();
646
647 if ((m_testname == 'F' or m_testname == 'G') and m_testname != 'V' and
648 not FlowLawIsPatersonBuddCold(flow_law, *m_config, EC)) {
649 m_log->message(1,
650 "pismv WARNING: flow law must be cold part of Paterson-Budd ('-siafd_flow_law arr')\n"
651 " for reported errors in test %c to be meaningful!\n",
652 m_testname);
653 }
654
655 m_log->message(1,
656 "NUMERICAL ERRORS evaluated at final time (relative to exact solution):\n");
657
658
659 // geometry (thickness, vol) errors if appropriate; reported in m except for relmaxETA
660 if ((m_testname != 'K') && (m_testname != 'O')) {
661 computeGeometryErrors(volexact,areaexact,domeHexact,
662 volerr,areaerr,maxHerr,avHerr,maxetaerr,centerHerr);
663 m_log->message(1,
664 "geometry : prcntVOL maxH avH relmaxETA\n"); // no longer reporting centerHerr
665 m_log->message(1, " %12.6f%12.6f%12.6f%12.6f\n",
666 100*volerr/volexact, maxHerr, avHerr,
667 maxetaerr/pow(domeHexact,m));
668
669 }
670
671 // temperature errors for F and G
672 if ((m_testname == 'F') || (m_testname == 'G')) {
673 computeTemperatureErrors(maxTerr, avTerr);
674 computeBasalTemperatureErrors(basemaxTerr, baseavTerr, basecenterTerr);
675 m_log->message(1,
676 "temp : maxT avT basemaxT baseavT\n"); // no longer reporting basecenterT
677 m_log->message(1, " %12.6f%12.6f%12.6f%12.6f\n",
678 maxTerr, avTerr, basemaxTerr, baseavTerr);
679
680 } else if ((m_testname == 'K') || (m_testname == 'O')) {
681 computeIceBedrockTemperatureErrors(maxTerr, avTerr, maxTberr, avTberr);
682 m_log->message(1,
683 "temp : maxT avT maxTb avTb\n");
684 m_log->message(1, " %12.6f%12.6f%12.6f%12.6f\n",
685 maxTerr, avTerr, maxTberr, avTberr);
686
687 }
688
689 // strain_heating errors if appropriate; reported in 10^6 J/(s m^3)
690 if ((m_testname == 'F') || (m_testname == 'G')) {
691 compute_strain_heating_errors(max_strain_heating_error, av_strain_heating_error);
692 m_log->message(1,
693 "Sigma : maxSig avSig\n");
694 m_log->message(1, " %12.6f%12.6f\n",
695 max_strain_heating_error*1.0e6, av_strain_heating_error*1.0e6);
696 }
697
698 // surface velocity errors if exact values are available; reported in m year-1
699 if ((m_testname == 'F') || (m_testname == 'G')) {
700 computeSurfaceVelocityErrors(maxUerr, avUerr, maxWerr, avWerr);
701 m_log->message(1,
702 "surf vels : maxUvec avUvec maxW avW\n");
703 m_log->message(1,
704 " %12.6f%12.6f%12.6f%12.6f\n",
705 convert(m_sys, maxUerr, "m second-1", "m year-1"),
706 convert(m_sys, avUerr, "m second-1", "m year-1"),
707 convert(m_sys, maxWerr, "m second-1", "m year-1"),
708 convert(m_sys, avWerr, "m second-1", "m year-1"));
709 }
710
711 // basal melt rate errors if appropriate; reported in m year-1
712 if (m_testname == 'O') {
713 computeBasalMeltRateErrors(maxbmelterr, minbmelterr);
714 if (maxbmelterr != minbmelterr) {
715 m_log->message(1,
716 "IceCompModel WARNING: unexpected Test O situation: max and min of bmelt error\n"
717 " are different: maxbmelterr = %f, minbmelterr = %f\n",
718 convert(m_sys, maxbmelterr, "m second-1", "m year-1"),
719 convert(m_sys, minbmelterr, "m second-1", "m year-1"));
720 }
721 m_log->message(1,
722 "basal melt: max\n");
723 m_log->message(1, " %11.5f\n",
724 convert(m_sys, maxbmelterr, "m second-1", "m year-1"));
725
726 }
727
728 m_log->message(1, "NUM ERRORS DONE\n");
729
730 options::String report_file("-report_file", "NetCDF error report file");
731 bool append = options::Bool("-append", "Append the NetCDF error report");
732
733 IO_Mode mode = append ? PISM_READWRITE : PISM_READWRITE_MOVE;
734
735 if (report_file.is_set()) {
736 unsigned int start;
737 TimeseriesMetadata err("N", "N", m_sys);
738
739 err.set_string("units", "1");
740
741 m_log->message(2, "Also writing errors to '%s'...\n", report_file->c_str());
742
743 // Find the number of records in this file:
744 File file(m_grid->com, report_file, PISM_NETCDF3, mode); // OK to use netcdf3
745
746 start = file.dimension_length("N");
747
748 io::write_attributes(file, m_output_global_attributes, PISM_DOUBLE);
749
750 // Write the dimension variable:
751 io::write_timeseries(file, err, (size_t)start, (double)(start + 1), PISM_INT);
752
753 // Always write grid parameters:
754 err.set_name("dx");
755 err.set_string("units", "meters");
756 io::write_timeseries(file, err, (size_t)start, m_grid->dx());
757 err.set_name("dy");
758 io::write_timeseries(file, err, (size_t)start, m_grid->dy());
759 err.set_name("dz");
760 io::write_timeseries(file, err, (size_t)start, m_grid->dz_max());
761
762 // Always write the test name:
763 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
764 err.set_name("test");
765 io::write_timeseries(file, err, (size_t)start, (double)m_testname, PISM_BYTE);
766
767 if ((m_testname != 'K') && (m_testname != 'O')) {
768 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
769 err.set_name("relative_volume");
770 err.set_string("units", "percent");
771 err.set_string("long_name", "relative ice volume error");
772 io::write_timeseries(file, err, (size_t)start, 100*volerr/volexact);
773
774 err.set_name("relative_max_eta");
775 err.set_string("units", "1");
776 err.set_string("long_name", "relative $\\eta$ error");
777 io::write_timeseries(file, err, (size_t)start, maxetaerr/pow(domeHexact,m));
778
779 err.set_name("maximum_thickness");
780 err.set_string("units", "meters");
781 err.set_string("long_name", "maximum ice thickness error");
782 io::write_timeseries(file, err, (size_t)start, maxHerr);
783
784 err.set_name("average_thickness");
785 err.set_string("units", "meters");
786 err.set_string("long_name", "average ice thickness error");
787 io::write_timeseries(file, err, (size_t)start, avHerr);
788 }
789
790 if ((m_testname == 'F') || (m_testname == 'G')) {
791 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
792 err.set_name("maximum_temperature");
793 err.set_string("units", "Kelvin");
794 err.set_string("long_name", "maximum ice temperature error");
795 io::write_timeseries(file, err, (size_t)start, maxTerr);
796
797 err.set_name("average_temperature");
798 err.set_string("long_name", "average ice temperature error");
799 io::write_timeseries(file, err, (size_t)start, avTerr);
800
801 err.set_name("maximum_basal_temperature");
802 err.set_string("long_name", "maximum basal temperature error");
803 io::write_timeseries(file, err, (size_t)start, basemaxTerr);
804 err.set_name("average_basal_temperature");
805 err.set_string("long_name", "average basal temperature error");
806 io::write_timeseries(file, err, (size_t)start, baseavTerr);
807
808 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
809 err.set_name("maximum_sigma");
810 err.set_string("units", "J s-1 m-3");
811 err.set_string("glaciological_units", "1e6 J s-1 m-3");
812 err.set_string("long_name", "maximum strain heating error");
813 io::write_timeseries(file, err, (size_t)start, max_strain_heating_error);
814
815 err.set_name("average_sigma");
816 err.set_string("long_name", "average strain heating error");
817 io::write_timeseries(file, err, (size_t)start, av_strain_heating_error);
818
819 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
820 err.set_name("maximum_surface_velocity");
821 err.set_string("long_name", "maximum ice surface horizontal velocity error");
822 err.set_string("units", "m second-1");
823 err.set_string("glaciological_units", "meters year-1");
824 io::write_timeseries(file, err, (size_t)start, maxUerr);
825
826 err.set_name("average_surface_velocity");
827 err.set_string("long_name", "average ice surface horizontal velocity error");
828 io::write_timeseries(file, err, (size_t)start, avUerr);
829
830 err.set_name("maximum_surface_w");
831 err.set_string("long_name", "maximum ice surface vertical velocity error");
832 io::write_timeseries(file, err, (size_t)start, maxWerr);
833
834 err.set_name("average_surface_w");
835 err.set_string("long_name", "average ice surface vertical velocity error");
836 io::write_timeseries(file, err, (size_t)start, avWerr);
837 }
838
839 if ((m_testname == 'K') || (m_testname == 'O')) {
840 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
841 err.set_name("maximum_temperature");
842 err.set_string("units", "Kelvin");
843 err.set_string("long_name", "maximum ice temperature error");
844 io::write_timeseries(file, err, (size_t)start, maxTerr);
845
846 err.set_name("average_temperature");
847 err.set_string("long_name", "average ice temperature error");
848 io::write_timeseries(file, err, (size_t)start, avTerr);
849
850 err.set_name("maximum_bedrock_temperature");
851 err.set_string("long_name", "maximum bedrock temperature error");
852 io::write_timeseries(file, err, (size_t)start, maxTberr);
853
854 err.set_name("average_bedrock_temperature");
855 err.set_string("long_name", "average bedrock temperature error");
856 io::write_timeseries(file, err, (size_t)start, avTberr);
857 }
858
859 if (m_testname == 'O') {
860 err.clear_all_strings(); err.clear_all_doubles(); err.set_string("units", "1");
861 err.set_name("maximum_basal_melt_rate");
862 err.set_string("units", "m second-1");
863 err.set_string("glaciological_units", "meters year-1");
864 io::write_timeseries(file, err, (size_t)start, maxbmelterr);
865 }
866 }
867
868 }
869
870 //! \brief Initialize test V.
871 /*
872 Try
873
874 pismv -test V -y 1000 -part_grid -ssa_method fd -cfbc -o fig4-blue.nc
875 pismv -test V -y 1000 -part_grid -ssa_method fd -o fig4-green.nc
876
877 to try to reproduce Figure 4.
878
879 Try
880
881 pismv -test V -y 3000 -ssa_method fd -cfbc -o fig5.nc -thickness_calving_threshold 250 -part_grid
882
883 with -Mx 51, -Mx 101, -Mx 201 for figure 5,
884
885 pismv -test V -y 300 -ssa_method fd -o fig6-ab.nc
886
887 for 6a and 6b,
888
889 pismv -test V -y 300 -ssa_method fd -cfbc -part_grid -o fig6-cd.nc
890
891 for 6c and 6d,
892
893 pismv -test V -y 300 -ssa_method fd -cfbc -part_grid -o fig6-ef.nc
894
895 for 6e and 6f.
896
897 */
898 void IceCompModel::test_V_init() {
899
900 {
901 IceModelVec2S bed_uplift(m_grid, "uplift", WITHOUT_GHOSTS);
902 bed_uplift.set(0.0);
903 m_geometry.bed_elevation.set(-1000);
904
905 m_beddef->bootstrap(m_geometry.bed_elevation, bed_uplift, m_geometry.ice_thickness,
906 m_geometry.sea_level_elevation);
907 }
908
909 // set SSA boundary conditions:
910 double upstream_velocity = convert(m_sys, 300.0, "m year-1", "m second-1"),
911 upstream_thk = 600.0;
912
913 IceModelVec::AccessList list{&m_geometry.ice_thickness, &m_ssa_dirichlet_bc_mask,
914 &m_ssa_dirichlet_bc_values};
915
916 for (Points p(*m_grid); p; p.next()) {
917 const int i = p.i(), j = p.j();
918
919 if (i <= 2) {
920 m_ssa_dirichlet_bc_mask(i,j) = 1;
921 m_ssa_dirichlet_bc_values(i,j) = Vector2(upstream_velocity, 0.0);
922 m_geometry.ice_thickness(i, j) = upstream_thk;
923 } else {
924 m_ssa_dirichlet_bc_mask(i,j) = 0;
925 m_ssa_dirichlet_bc_values(i,j) = Vector2(0.0, 0.0);
926 m_geometry.ice_thickness(i, j) = 0;
927 }
928 }
929
930 m_ssa_dirichlet_bc_mask.update_ghosts();
931
932 m_ssa_dirichlet_bc_values.update_ghosts();
933
934 m_geometry.ice_thickness.update_ghosts();
935 }
936
937 } // end of namespace pism