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