Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
IceModel.cc
Go to the documentation of this file.
1// Copyright (C) 2004-2024 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#include <algorithm>
22#include <petscsys.h>
23
24#include "pism/pism_config.hh"
25
26#include "pism/icemodel/IceModel.hh"
27
28#include "pism/basalstrength/YieldStress.hh"
29#include "pism/frontretreat/util/IcebergRemover.hh"
30#include "pism/energy/BedThermalUnit.hh"
31#include "pism/hydrology/Hydrology.hh"
32#include "pism/stressbalance/StressBalance.hh"
33#include "pism/util/Grid.hh"
34#include "pism/util/ConfigInterface.hh"
35#include "pism/util/Diagnostic.hh"
36#include "pism/util/error_handling.hh"
37#include "pism/coupler/SeaLevel.hh"
38#include "pism/coupler/OceanModel.hh"
39#include "pism/coupler/SurfaceModel.hh"
40#include "pism/earth/BedDef.hh"
41#include "pism/util/pism_signal.h"
42#include "pism/util/Vars.hh"
43#include "pism/util/Profiling.hh"
44#include "pism/util/pism_utilities.hh"
45#include "pism/age/AgeModel.hh"
46#include "pism/age/Isochrones.hh"
47#include "pism/energy/EnergyModel.hh"
48#include "pism/util/io/File.hh"
49#include "pism/util/array/Forcing.hh"
50#include "pism/fracturedensity/FractureDensity.hh"
51#include "pism/coupler/util/options.hh" // ForcingOptions
52#include "pism/coupler/ocean/PyOceanModel.hh"
53
54namespace pism {
55
56IceModel::IceModel(std::shared_ptr<Grid> grid, const std::shared_ptr<Context> &context)
57 : m_grid(grid),
58 m_config(context->config()),
59 m_ctx(context),
60 m_sys(context->unit_system()),
61 m_log(context->log()),
62 m_time(context->time()),
63 m_output_global_attributes("PISM_GLOBAL", m_sys),
64 m_geometry(m_grid),
65 m_new_bed_elevation(true),
66 m_basal_yield_stress(m_grid, "tauc"),
67 m_basal_melt_rate(m_grid, "bmelt"),
68 m_bedtoptemp(m_grid, "bedtoptemp"),
69 m_velocity_bc_mask(m_grid, "vel_bc_mask"),
70 m_velocity_bc_values(m_grid, "_bc"), // u_bc and v_bc
71 m_ice_thickness_bc_mask(grid, "thk_bc_mask"),
72 m_step_counter(0),
73 m_thickness_change(grid),
74 m_ts_times(new std::vector<double>()) {
75
78
79 pism_signal = 0;
80 signal(SIGTERM, pism_signal_handler);
81 signal(SIGUSR1, pism_signal_handler);
82 signal(SIGUSR2, pism_signal_handler);
83
84 m_surface = nullptr;
85 m_ocean = nullptr;
86 m_beddef = nullptr;
87
88 m_btu = nullptr;
89 m_energy_model = nullptr;
90
91 m_output_global_attributes["Conventions"] = "CF-1.6";
93
94 m_fracture = nullptr;
95
97
98 // allocate temporary storage
99 {
100 // 2d work vectors
101 for (int j = 0; j < m_n_work2d; j++) {
102 m_work2d.push_back(
103 std::make_shared<array::Scalar2>(m_grid, pism::printf("work_vector_%d", j)));
104 }
105 }
106
107 auto surface_input_file = m_config->get_string("hydrology.surface_input.file");
108 if (not surface_input_file.empty()) {
109 ForcingOptions surface_input(*m_ctx, "hydrology.surface_input");
110 int buffer_size = static_cast<int>(m_config->get_number("input.forcing.buffer_size"));
111
112 File file(m_grid->com, surface_input.filename, io::PISM_NETCDF3, io::PISM_READONLY);
113
115 std::make_shared<array::Forcing>(m_grid, file, "water_input_rate",
116 "", // no standard name
117 buffer_size, surface_input.periodic);
119 .long_name("water input rate for the subglacial hydrology model")
120 .units("kg m^-2 s^-1")
121 .output_units("kg m^-2 year^-1");
122 m_surface_input_for_hydrology->metadata()["valid_min"] = { 0.0 };
123 }
124}
125
126double IceModel::dt() const {
127 return m_dt;
128}
129
131 dt_TempAge = 0.0;
132 m_dt = 0.0;
134
135 {
136 int year_increment = static_cast<int>(m_config->get_number("time_stepping.hit_multiples"));
137
138 if (year_increment > 0) {
139 // start year:
140 auto year = m_time->units().date(m_time->current(), m_time->calendar()).year;
141
142 int last_multiple = year - year % year_increment;
143 // correct last_multiple if 'year' is negative
144 // and not a multiple of year_increment:
145 last_multiple -= year_increment * static_cast<int>(year % year_increment < 0);
146
147 units::DateTime last_date{ last_multiple, 1, 1, 0, 0, 0.0 };
148
149 m_timestep_hit_multiples_last_time = m_time->units().time(last_date, m_time->calendar());
150 } else {
152 }
153 }
154}
155
156
158 // empty; defined here to be able to use more forward-declared classes in IceModel.hh
159}
160
161
162//! Allocate all Arrays defined in IceModel.
163/*!
164 This procedure allocates the memory used to store model state, diagnostic and
165 work vectors and sets metadata.
166
167 Default values should not be set here; please use set_vars_from_options().
168
169 All the memory allocated here is freed by Arrays' destructors.
170*/
172
173 // FIXME: this should do for now, but we should pass a const reference to Geometry to sub-models
174 // as a function argument.
175 m_grid->variables().add(m_geometry.ice_surface_elevation);
176 m_grid->variables().add(m_geometry.ice_thickness);
177 m_grid->variables().add(m_geometry.cell_type);
178 m_grid->variables().add(m_geometry.sea_level_elevation);
179 m_grid->variables().add(m_geometry.longitude);
180 m_grid->variables().add(m_geometry.latitude);
181
182 if (m_config->get_flag("geometry.grounded_cell_fraction")) {
183 m_grid->variables().add(m_geometry.cell_grounded_fraction);
184 }
185
186 if (m_config->get_flag("geometry.part_grid.enabled")) {
188 }
189
190 // yield stress for basal till (plastic or pseudo-plastic model)
191 {
192 // PROPOSED standard_name = land_ice_basal_material_yield_stress
194 .long_name("yield stress for basal till (plastic or pseudo-plastic model)")
195 .units("Pa");
196 m_grid->variables().add(m_basal_yield_stress);
197 }
198
199 {
201 .long_name("temperature at the top surface of the bedrock thermal layer")
202 .units("kelvin");
203 }
204
205 // basal melt rate
207 .long_name(
208 "ice basal melt rate from energy conservation and subshelf melt, in ice thickness per time")
209 .units("m s^-1")
210 .output_units("m year^-1")
211 .standard_name("land_ice_basal_melt_rate");
212 m_basal_melt_rate.metadata()["comment"] = "positive basal melt rate corresponds to ice loss";
213 m_grid->variables().add(m_basal_melt_rate);
214
215 // Sliding velocity (usually SSA) Dirichlet B.C. locations and values
216 {
218 .long_name("Mask prescribing Dirichlet boundary locations for the sliding velocity")
221 m_velocity_bc_mask.metadata()["flag_values"] = { 0, 1 };
222 m_velocity_bc_mask.metadata()["flag_meanings"] = "no_data boundary_condition";
223
225 }
226 // SSA Dirichlet B.C. values
227 {
228 double fill_value = m_config->get_number("output.fill_value");
229 const double huge_value = 1e6;
230 // vel_bc
232 "X-component of the SSA velocity boundary conditions");
234 "Y-component of the SSA velocity boundary conditions");
235 for (int j : { 0, 1 }) {
236 m_velocity_bc_values.metadata(j)["valid_range"] = { -huge_value, huge_value };
237 m_velocity_bc_values.metadata(j)["_FillValue"] = { fill_value };
239 }
240 }
241
242 // Ice thickness BC mask
243 {
245 .long_name("Mask specifying locations where ice thickness is held constant")
248 m_ice_thickness_bc_mask.metadata()["flag_values"] = {0, 1};
249 m_ice_thickness_bc_mask.metadata()["flag_meanings"] = "no_data boundary_condition";
250
252 }
253
254 // Add some variables to the list of "model state" fields.
257
259
264}
265
266//! Update the surface elevation and the flow-type mask when the geometry has changed.
267/*!
268 This procedure should be called whenever necessary to maintain consistency of geometry.
269
270 For instance, it should be called when either ice thickness or bed elevation change.
271 In particular we always want \f$h = H + b\f$ to apply at grounded points, and, on the
272 other hand, we want the mask to reflect that the ice is floating if the flotation
273 criterion applies at a point.
274
275 If `flag == REMOVE_ICEBERGS`, also calls the code which removes icebergs, to avoid
276 stress balance solver problems caused by ice that is not attached to the grounded ice
277 sheet.
278*/
280
281 m_geometry.bed_elevation.copy_from(m_beddef->bed_elevation());
283
284 if (m_iceberg_remover and flag == REMOVE_ICEBERGS) {
285 // The iceberg remover has to use the same mask as the stress balance code, hence the
286 // stress-balance-related threshold here.
287 m_geometry.ensure_consistency(m_config->get_number("stress_balance.ice_free_thickness_standard"));
288
292 // The call above modifies ice thickness and updates the mask accordingly, but we re-compute the
293 // mask (we need to use a different threshold).
294 }
295
296 // This will ensure that ice area specific volume is zero if ice thickness is greater
297 // than zero, then compute new surface elevation and mask.
298 m_geometry.ensure_consistency(m_config->get_number("geometry.ice_free_thickness_standard"));
299
300 if (flag == REMOVE_ICEBERGS) {
301 // clean up partially-filled cells that are not next to ice
304
305 for (auto p = m_grid->points(); p; p.next()) {
306 const int i = p.i(), j = p.j();
307
308 if (m_geometry.ice_area_specific_volume(i, j) > 0.0 and
309 not m_geometry.cell_type.next_to_ice(i, j)) {
311 }
312 }
313 }
314}
315
318 if (m_config->get_flag("geometry.update.use_basal_melt_rate")) {
320 }
321
323 result.geometry = &m_geometry;
325 result.enthalpy = &m_energy_model->enthalpy();
326 result.age = m_age_model ? &m_age_model->age() : nullptr;
327
328 result.water_column_pressure = &m_ocean->average_water_column_pressure();
329
330 if (m_config->get_flag("stress_balance.ssa.dirichlet_bc")) {
331 result.bc_mask = &m_velocity_bc_mask;
333 }
334
335 if (m_config->get_flag("fracture_density.enabled")) {
336 result.fracture_density = &m_fracture->density();
337 }
338
339 return result;
340}
341
343 energy::Inputs result;
344
345 result.basal_frictional_heating = &m_stress_balance->basal_frictional_heating();
346 result.basal_heat_flux = &m_btu->flux_through_top_surface(); // bedrock thermal layer
347 result.cell_type = &m_geometry.cell_type; // geometry
348 result.ice_thickness = &m_geometry.ice_thickness; // geometry
349 result.shelf_base_temp = &m_ocean->shelf_base_temperature(); // ocean model
350 result.till_water_thickness = &m_subglacial_hydrology->till_water_thickness();
351 result.surface_liquid_fraction = &m_surface->liquid_water_fraction(); // surface model
352 result.surface_temp = &m_surface->temperature(); // surface model
353
354 result.volumetric_heating_rate = &m_stress_balance->volumetric_strain_heating();
355 result.u3 = &m_stress_balance->velocity_u();
356 result.v3 = &m_stress_balance->velocity_v();
357 result.w3 = &m_stress_balance->velocity_w();
358
359 result.check(); // make sure all data members were set
360
361 return result;
362}
363
365 YieldStressInputs result;
366
367 result.geometry = &m_geometry;
368 result.till_water_thickness = &m_subglacial_hydrology->till_water_thickness();
369 result.subglacial_water_thickness = &m_subglacial_hydrology->subglacial_water_thickness();
370
371 return result;
372}
373
374std::string IceModel::save_state_on_error(const std::string &suffix,
375 const std::set<std::string> &additional_variables) {
376 std::string filename = m_config->get_string("output.file");
377
378 if (filename.empty()) {
379 m_log->message(2, "WARNING: output.file is empty. Using unnamed.nc instead.");
380 filename = "unnamed.nc";
381 }
382
383 filename = filename_add_suffix(filename, suffix, "");
384
385 File file(m_grid->com,
386 filename,
387 string_to_backend(m_config->get_string("output.format")),
389
391
392 auto variables = output_variables("small");
393 for (const auto &v : additional_variables) {
394 variables.insert(v);
395 }
396
397 save_variables(file, INCLUDE_MODEL_STATE, variables, m_time->current());
398
399 return filename;
400}
401
402//! The contents of the main PISM time-step.
403/*!
404During the time-step we perform the following actions:
405 */
406void IceModel::step(bool do_mass_continuity,
407 bool do_skip) {
408
410
411 const Profiling &profiling = m_ctx->profiling();
412
413 double current_time = m_time->current();
414
415 //! \li call pre_step_hook() to let derived classes do more
417
418 //! \li update the velocity field; in some cases the whole three-dimensional
419 //! field is updated and in some cases just the vertically-averaged
420 //! horizontal velocity is updated
421
422 // always "update" ice velocity (possibly trivially); only update
423 // SSA and only update velocities at depth if suggested by temp and age
424 // stability criterion; note *lots* of communication is avoided by skipping
425 // SSA (and temp/age)
426
427 const bool updateAtDepth = (m_skip_countdown == 0);
428
429 // Combine basal melt rate in grounded (computed during the energy
430 // step) and floating (provided by an ocean model) areas.
431 //
432 // Basal melt rate may be used by a stress balance model to compute vertical velocity of
433 // ice.
434 {
436 m_ocean->shelf_base_mass_flux(),
437 m_energy_model->basal_melt_rate(),
439 }
440
441 try {
442 profiling.begin("stress_balance");
443 m_stress_balance->update(stress_balance_inputs(), updateAtDepth);
444 profiling.end("stress_balance");
445 } catch (RuntimeError &e) {
446 std::string output_file = save_state_on_error("_stressbalance_failed", {});
447
448 e.add_context("performing a time step. (Note: Model state was saved to '%s'.)",
449 output_file.c_str());
450 throw;
451 }
452
453
454 m_stdout_flags += m_stress_balance->stdout_report();
455
456 m_stdout_flags += (updateAtDepth ? "v" : "V");
457
458 //! \li determine the time step according to a variety of stability criteria
459 auto dt_info = max_timestep(m_skip_countdown);
460 m_dt = dt_info.dt;
461 m_adaptive_timestep_reason = dt_info.reason;
462 m_skip_countdown = dt_info.skip_counter;
463
464 //! \li update the yield stress for the plastic till model (if appropriate)
466 profiling.begin("basal_yield_stress");
467 m_basal_yield_stress_model->update(yield_stress_inputs(), current_time, m_dt);
468 profiling.end("basal_yield_stress");
469 m_basal_yield_stress.copy_from(m_basal_yield_stress_model->basal_material_yield_stress());
470 m_stdout_flags += "y";
471 } else {
472 m_stdout_flags += "$";
473 }
474
475 dt_TempAge += m_dt;
476
477 //! \li update the age of the ice (if appropriate)
478 if (m_age_model and updateAtDepth) {
479 AgeModelInputs inputs;
481 inputs.u3 = &m_stress_balance->velocity_u();
482 inputs.v3 = &m_stress_balance->velocity_v();
483 inputs.w3 = &m_stress_balance->velocity_w();
484
485 profiling.begin("age");
486 m_age_model->update(current_time, dt_TempAge, inputs);
487 profiling.end("age");
488 m_stdout_flags += "a";
489 } else {
490 m_stdout_flags += "$";
491 }
492
493 //! \li update the enthalpy (or temperature) field according to the conservation of
494 //! energy model based (especially) on the new velocity field; see
495 //! energy_step()
496 if (updateAtDepth) { // do the energy step
497 profiling.begin("energy");
498 energy_step();
499 profiling.end("energy");
500 m_stdout_flags += "E";
501 } else {
502 m_stdout_flags += "$";
503 }
504
505 //! \li update the fracture density field; see update_fracture_density()
506 if (m_config->get_flag("fracture_density.enabled")) {
507 profiling.begin("fracture_density");
509 profiling.end("fracture_density");
510 }
511
512 //! \li update the thickness of the ice according to the mass conservation model and calving
513 //! parameterizations
514
515 if (do_mass_continuity) {
516 // reset the conservation error field:
517 m_geometry_evolution->reset();
518
519 profiling.begin("mass_transport");
520 {
521 // Note that there are three adaptive time-stepping criteria. Two of them (using max.
522 // diffusion and 2D CFL) are limiting the mass-continuity time-step and the third (3D
523 // CFL) limits the energy and age time-steps.
524
525 // The mass-continuity time-step is usually smaller, and the skipping mechanism lets us
526 // do several mass-continuity steps for each energy step.
527
528 // When -no_mass is set, mass-continuity-related time-step restrictions are disabled,
529 // making "skipping" unnecessary.
530
531 // This is why the following two lines appear here and are executed only if
532 // do_mass_continuity is true.
533 if (do_skip and m_skip_countdown > 0) {
535 }
536
537 try {
539 m_dt,
540 m_stress_balance->advective_velocity(),
541 m_stress_balance->diffusive_flux(),
543 } catch (RuntimeError &e) {
544 std::string output_file = save_state_on_error("_mass_transport_failed",
545 {"flux_staggered", "flux_divergence"});
546
547 e.add_context("performing a mass transport time step (dt=%f s). (Note: Model state was saved to '%s'.)",
548 m_dt, output_file.c_str());
549 throw;
550 }
551
552 m_geometry_evolution->apply_flux_divergence(m_geometry);
553
555 }
556 profiling.end("mass_transport");
557
558 // calving, frontal melt, and discharge accounting
559 profiling.begin("front_retreat");
561 profiling.end("front_retreat");
562
563 m_stdout_flags += "h";
564 } else {
565 m_stdout_flags += "$";
566 }
567
568 profiling.begin("sea_level");
569 m_sea_level->update(m_geometry, current_time, m_dt);
570 profiling.end("sea_level");
571
572 profiling.begin("ocean");
573 m_ocean->update(m_geometry, current_time, m_dt);
574 profiling.end("ocean");
575
576 // The sea level elevation might have changed, so we need to update the mask, etc. Note
577 // that THIS MAY PRODUCE ICEBERGS, but we assume that the surface model does not care.
579
580 //! \li Update surface and ocean models.
581 profiling.begin("surface");
582 m_surface->update(m_geometry, current_time, m_dt);
583 profiling.end("surface");
584
585
586 if (do_mass_continuity) {
587 // compute and apply effective surface and basal mass balance
588
589 m_geometry_evolution->source_term_step(m_geometry, m_dt,
591 m_surface->mass_flux(),
593 m_geometry_evolution->apply_mass_fluxes(m_geometry);
594
595 // add removed icebergs to discharge due to calving
596 {
597 auto &old_H = *m_work2d[0], &old_Href = *m_work2d[1];
598
599 {
600 old_H.copy_from(m_geometry.ice_thickness);
601 old_Href.copy_from(m_geometry.ice_area_specific_volume);
602 }
603
604 // the last call has to remove icebergs
606
607 bool add_values = true;
610 old_H, old_Href,
611 add_values,
613 }
614 }
615
616 if (m_isochrones) {
617 m_isochrones->update(current_time, m_dt,
618 m_stress_balance->velocity_u(),
619 m_stress_balance->velocity_v(),
621 m_geometry_evolution->top_surface_mass_balance(),
622 m_geometry_evolution->bottom_surface_mass_balance());
623 }
624
625 //! \li update the state variables in the subglacial hydrology model (typically
626 //! water thickness and sometimes pressure)
627 profiling.begin("basal_hydrology");
629 profiling.end("basal_hydrology");
630
631 //! \li compute the bed deformation, which depends on current thickness, bed elevation,
632 //! and sea level
633 if (m_beddef) {
634 int topg_state_counter = m_beddef->bed_elevation().state_counter();
635
636 profiling.begin("bed_deformation");
639 current_time, m_dt);
640 profiling.end("bed_deformation");
641
642 m_new_bed_elevation = m_beddef->bed_elevation().state_counter() != topg_state_counter;
643 } else {
644 m_new_bed_elevation = false;
645 }
646
647 if (m_config->get_flag("time_stepping.assume_bed_elevation_changed")) {
648 m_new_bed_elevation = true;
649 }
650
653 m_stdout_flags += "b";
654 } else {
655 m_stdout_flags += " ";
656 }
657
658 //! \li call post_step_hook() to let derived classes do more
660
661 // Done with the step; now adopt the new time.
662 m_time->step(m_dt);
663
664 if (updateAtDepth) {
665 t_TempAge = m_time->current();
666 dt_TempAge = 0.0;
667 }
668
669 // Check if the ice thickness exceeded the height of the computational box and stop if it did.
670 if (max(m_geometry.ice_thickness) > m_grid->Lz()) {
671 auto o_file = save_state_on_error("_max_thickness", {});
672
674 "Ice thickness exceeds the height of the computational box (%7.4f m).\n"
675 "The model state was saved to '%s'. To continue this simulation,\n"
676 "run with\n"
677 "-i %s -bootstrap -regrid_file %s -allow_extrapolation -Lz N [other options]\n"
678 "where N > %7.4f.",
679 m_grid->Lz(), o_file.c_str(), o_file.c_str(), o_file.c_str(), m_grid->Lz());
680 }
681
682 // end the flag line
684}
685
686/*!
687 * Note: don't forget to update IceRegionalModel::hydrology_step() if necessary.
688 */
690 hydrology::Inputs inputs;
691
692 array::Scalar &sliding_speed = *m_work2d[0];
693 compute_magnitude(m_stress_balance->advective_velocity(), sliding_speed);
694
695 inputs.no_model_mask = nullptr;
696 inputs.geometry = &m_geometry;
697 inputs.surface_input_rate = nullptr;
699 inputs.ice_sliding_speed = &sliding_speed;
700
702 m_surface_input_for_hydrology->update(m_time->current(), m_dt);
703 m_surface_input_for_hydrology->average(m_time->current(), m_dt);
705 } else if (m_config->get_flag("hydrology.surface_input_from_runoff")) {
706 // convert [kg m-2] to [kg m-2 s-1]
707 array::Scalar &surface_input_rate = *m_work2d[1];
708 surface_input_rate.copy_from(m_surface->runoff());
709 surface_input_rate.scale(1.0 / m_dt);
710 inputs.surface_input_rate = &surface_input_rate;
711 }
712
713 m_subglacial_hydrology->update(m_time->current(), m_dt, inputs);
714}
715
716//! Virtual. Does nothing in `IceModel`. Derived classes can do more computation in each time step.
718 // empty
719}
720
721//! Virtual. Does nothing in `IceModel`. Derived classes can do more computation in each time step.
723 // empty
724}
725
726/**
727 * Run the time-stepping loop from the current model time to `time`.
728 *
729 * This should be called by the coupler controlling PISM when it is
730 * running alongside a GCM.
731 *
732 * @param run_end model time (in seconds) to run to
733 *
734 * @return 0 on success
735 */
737
738 m_time->set_end(run_end);
739
740 return run();
741}
742
743
744/**
745 * Run the time-stepping loop from the current time until the time
746 * specified by the IceModel::grid::time object.
747 *
748 * This is the method used by PISM in the "standalone" mode.
749 */
751 const Profiling &profiling = m_ctx->profiling();
752
753 bool do_mass_conserve = m_config->get_flag("geometry.update.enabled");
754 bool do_energy = member(m_config->get_string("energy.model"), {"cold", "enthalpy"});
755 bool do_skip = m_config->get_flag("time_stepping.skip.enabled");
756
757 // de-allocate diagnostics that are not needed
759
760 // Enforce consistency *and* remove icebergs. During time-stepping we remove icebergs at
761 // the end of the time step, so we need to ensure that ice geometry is "OK" before the
762 // first step.
764
765 // Update spatially-variable diagnostics at the beginning of the run.
766 write_extras();
767
768 // Update scalar time series to remember the state at the beginning of the run.
769 // This is needed to compute rates of change of the ice mass, volume, etc.
770 {
771 const double time = m_time->current();
772 for (const auto &d : m_ts_diagnostics) {
773 d.second->update(time, time);
774 }
775 }
776
777 m_log->message(2, "running forward ...\n");
778
779 m_stdout_flags.erase(); // clear it out
780 print_summary_line(true, do_energy, 0.0, 0.0, 0.0, 0.0, 0.0);
781 print_summary(do_energy); // report starting state
782
783 t_TempAge = m_time->current();
784 dt_TempAge = 0.0;
785
786 IceModelTerminationReason termination_reason = PISM_DONE;
787 // main loop for time evolution
788 // IceModel::step calls Time::step(dt), ensuring that this while loop
789 // will terminate
790 profiling.stage_begin("time-stepping loop");
791 while (m_time->current() < m_time->end()) {
792
793 m_stdout_flags.erase(); // clear it out
794
795 step(do_mass_conserve, do_skip);
796
798
799 // report a summary for major steps or the last one
800 bool updateAtDepth = m_skip_countdown == 0;
801 bool tempAgeStep = updateAtDepth and (m_age_model or do_energy);
802
803 double time_resolution = m_config->get_number("time_stepping.resolution");
804 const bool show_step = tempAgeStep or fabs(m_time->current() - m_time->end()) < time_resolution;
805 print_summary(show_step);
806
807 // update viewers before writing extras because writing extras resets diagnostics
809
810 // writing these fields here ensures that we do it after the last time-step
811 profiling.begin("io");
813 write_extras();
814 bool stop_after_chekpoint = write_checkpoint();
815 profiling.end("io");
816
817 if (stop_after_chekpoint) {
818 termination_reason = PISM_CHEKPOINT;
819 break;
820 }
821
822 if (process_signals() != 0) {
823 termination_reason = PISM_SIGNAL;
824 break;
825 }
826 } // end of the time-stepping loop
827 profiling.stage_end("time-stepping loop");
828
829 return termination_reason;
830}
831
832//! Manage the initialization of the IceModel object.
833/*!
834Please see the documenting comments of the functions called below to find
835explanations of their intended uses.
836 */
838 // Get the start time in seconds and ensure that it is consistent
839 // across all processors.
841
842 const Profiling &profiling = m_ctx->profiling();
843
844 profiling.begin("initialization");
845
846 //! The IceModel initialization sequence is this:
847
848 //! 1) Initialize model time:
849 time_setup();
850
851 //! 2) Process the options:
853
854 //! 3) Memory allocation:
856
857 //! 4) Allocate PISM components modeling some physical processes.
859
860 //! 6) Initialize coupler models and fill the model state variables
861 //! (from a PISM output file, from a bootstrapping file using some
862 //! modeling choices or using formulas). Calls IceModel::regrid()
864
865 //! 7) Report grid parameters:
866 m_log->message(2, "Computational domain and grid:\n");
867 m_grid->report_parameters();
868
869 //! 8) Miscellaneous stuff: set up the bed deformation model, initialize the
870 //! basal till model, initialize snapshots. This has to happen *after*
871 //! regridding.
872 misc_setup();
873
874 profiling.end("initialization");
875}
876
878 return m_geometry;
879}
880
884
886 return this->m_stress_balance.get();
887}
888
890 return m_ocean.get();
891}
892
894 return m_btu.get();
895}
896
900
904
906 return m_beddef.get();
907}
908
909/*!
910 * Return thickness change due to calving (over the last time step).
911 */
915
916/*!
917 * Return thickness change due to frontal melt (over the last time step).
918 */
922
923/*!
924 * Return thickness change due to forced retreat (over the last time step).
925 */
929
931 const std::set<std::string> &vars,
932 const std::string &type,
933 const std::set<std::string> &available,
934 bool stop) {
935 std::vector<std::string> missing;
936 for (const auto &v : vars) {
937 if (available.find(v) == available.end()) {
938 missing.push_back(v);
939 }
940 }
941
942 if (not missing.empty()) {
943 size_t N = missing.size();
944 const char *ending = N > 1 ? "s" : "";
945 const char *verb = N > 1 ? "are" : "is";
946 if (stop) {
948 "%s variable%s %s %s not available!\n"
949 "Available variables:\n- %s",
950 type.c_str(),
951 ending,
952 join(missing, ",").c_str(),
953 verb,
954 set_join(available, ",\n- ").c_str());
955 }
956
957 log.message(2,
958 "\nWARNING: %s variable%s %s %s not available!\n\n",
959 type.c_str(),
960 ending,
961 join(missing, ",").c_str(),
962 verb);
963 }
964}
965
966/*!
967 * De-allocate diagnostics that were not requested.
968 *
969 * Checks viewers, -extra_vars, -checkpoint, -save_vars, and regular output.
970 *
971 * FIXME: I need to make sure that these reporting mechanisms are active. It is possible that
972 * variables are on a list, but that list is not actually used.
973 */
975
976 // get the list of available diagnostics
977 std::set<std::string> available;
978 for (const auto &d : m_diagnostics) {
979 available.insert(d.first);
980 }
981
982 auto extra_stop = m_config->get_flag("output.extra.stop_missing");
983 warn_about_missing(*m_log, m_extra_vars, "diagnostic", available, extra_stop);
984
985 // get the list of requested diagnostics
986 auto requested = set_split(m_config->get_string("output.runtime.viewer.variables"), ',');
987 requested = combine(requested, m_output_vars);
988 requested = combine(requested, m_snapshot_vars);
989 requested = combine(requested, m_extra_vars);
990 requested = combine(requested, m_checkpoint_vars);
991
992 // de-allocate diagnostics that were not requested
993 for (const auto &v : available) {
994 if (requested.find(v) == requested.end()) {
995 m_diagnostics.erase(v);
996 }
997 }
998
999 // scalar time series
1000 std::vector<std::string> missing;
1001 if (not m_ts_filename.empty() and m_ts_vars.empty()) {
1002 // use all diagnostics
1003 } else {
1004 TSDiagnosticList diagnostics;
1005 for (const auto &v : m_ts_vars) {
1006 if (m_ts_diagnostics.find(v) != m_ts_diagnostics.end()) {
1007 diagnostics[v] = m_ts_diagnostics[v];
1008 } else {
1009 missing.push_back(v);
1010 }
1011 }
1012 // replace m_ts_diagnostics with requested diagnostics, de-allocating the rest
1013 m_ts_diagnostics = diagnostics;
1014 }
1015
1016 if (not missing.empty()) {
1018 "requested scalar diagnostics %s are not available",
1019 join(missing, ",").c_str());
1020 }
1021}
1022
1023/*!
1024 * Update diagnostics.
1025 *
1026 * This usually involves accumulating data needed to computed time-averaged quantities.
1027 *
1028 * Call this after prune_diagnostics() to avoid unnecessary work.
1029 */
1031 for (const auto &d : m_diagnostics) {
1032 d.second->update(dt);
1033 }
1034
1035 const double time = m_time->current();
1036 for (const auto &d : m_ts_diagnostics) {
1037 d.second->update(time - dt, time);
1038 }
1039}
1040
1041/*!
1042 * Reset accumulators in diagnostics that compute time-averaged quantities.
1043 */
1045 for (auto &d : m_diagnostics) {
1046 d.second->reset();
1047 }
1048}
1049
1050IceModel::ThicknessChanges::ThicknessChanges(const std::shared_ptr<const Grid> &grid)
1051 : calving(grid, "thickness_change_due_to_calving"),
1052 frontal_melt(grid, "thickness_change_due_to_frontal_melt"),
1053 forced_retreat(grid, "thickness_change_due_to_forced_retreat") {
1054 // empty
1055}
1056
1057void IceModel::set_python_ocean_model(std::shared_ptr<ocean::PyOceanModel> model) {
1058 m_ocean = std::make_shared<ocean::PyOceanModelAdapter>(m_grid, model);
1059 m_submodels["ocean model"] = m_ocean.get();
1060}
1061
1062} // end of namespace pism
const array::Array3D * w3
Definition AgeModel.hh:41
const array::Array3D * v3
Definition AgeModel.hh:40
const array::Array3D * u3
Definition AgeModel.hh:39
const array::Scalar * ice_thickness
Definition AgeModel.hh:38
High-level PISM I/O class.
Definition File.hh:55
array::Scalar1 sea_level_elevation
Definition Geometry.hh:48
array::Scalar cell_grounded_fraction
Definition Geometry.hh:56
void ensure_consistency(double ice_free_thickness_threshold)
Definition Geometry.cc:114
array::Scalar2 ice_surface_elevation
Definition Geometry.hh:57
array::Scalar1 ice_area_specific_volume
Definition Geometry.hh:52
array::CellType2 cell_type
Definition Geometry.hh:55
array::Scalar2 ice_thickness
Definition Geometry.hh:51
array::Scalar longitude
Definition Geometry.hh:42
array::Scalar2 bed_elevation
Definition Geometry.hh:47
array::Scalar latitude
Definition Geometry.hh:41
std::string m_adaptive_timestep_reason
Definition IceModel.hh:319
std::map< std::string, const Component * > m_submodels
the list of sub-models, for writing model states and obtaining diagnostics
Definition IceModel.hh:252
virtual energy::Inputs energy_model_inputs()
Definition IceModel.cc:342
std::set< std::string > m_snapshot_vars
Definition IceModel.hh:429
virtual int process_signals()
Catch signals -USR1, -USR2 and -TERM.
Definition utilities.cc:47
void enforce_consistency_of_geometry(ConsistencyFlag flag)
Update the surface elevation and the flow-type mask when the geometry has changed.
Definition IceModel.cc:279
const Geometry & geometry() const
Definition IceModel.cc:877
void compute_geometry_change(const array::Scalar &thickness, const array::Scalar &Href, const array::Scalar &thickness_old, const array::Scalar &Href_old, bool add_values, array::Scalar &output)
virtual void model_state_setup()
Sets the starting values of model state variables.
double m_start_time
Definition IceModel.hh:474
std::shared_ptr< stressbalance::StressBalance > m_stress_balance
Definition IceModel.hh:395
std::shared_ptr< Isochrones > m_isochrones
Definition IceModel.hh:263
IceModelTerminationReason run()
Definition IceModel.cc:750
std::shared_ptr< ocean::OceanModel > m_ocean
Definition IceModel.hh:280
virtual void post_step_hook()
Virtual. Does nothing in IceModel. Derived classes can do more computation in each time step.
Definition IceModel.cc:722
const ocean::OceanModel * ocean_model() const
Definition IceModel.cc:889
std::map< std::string, TSDiagnostic::Ptr > m_ts_diagnostics
Requested scalar diagnostics.
Definition IceModel.hh:419
virtual void energy_step()
Manage the solution of the energy equation, and related parallel communication.
Definition energy.cc:60
std::shared_ptr< surface::SurfaceModel > m_surface
Definition IceModel.hh:279
void write_extras()
Write spatially-variable diagnostic quantities.
void write_snapshot()
Writes a snapshot of the model state (if necessary)
unsigned int m_step_counter
Definition IceModel.hh:323
virtual stressbalance::Inputs stress_balance_inputs()
Definition IceModel.cc:316
std::shared_ptr< FractureDensity > m_fracture
Definition IceModel.hh:299
ThicknessChanges m_thickness_change
Definition IceModel.hh:410
virtual ~IceModel()
Definition IceModel.cc:157
void init()
Manage the initialization of the IceModel object.
Definition IceModel.cc:837
const Time::Ptr m_time
Time manager.
Definition IceModel.hh:246
Geometry m_geometry
Definition IceModel.hh:288
const GeometryEvolution & geometry_evolution() const
Definition IceModel.cc:881
const stressbalance::StressBalance * stress_balance() const
Definition IceModel.cc:885
bool m_new_bed_elevation
Definition IceModel.hh:290
virtual void reset_diagnostics()
Definition IceModel.cc:1044
static const int m_n_work2d
Definition IceModel.hh:392
std::shared_ptr< YieldStress > m_basal_yield_stress_model
Definition IceModel.hh:255
std::shared_ptr< calving::IcebergRemover > m_iceberg_remover
Definition IceModel.hh:265
double m_timestep_hit_multiples_last_time
Definition IceModel.hh:465
std::shared_ptr< Context > m_ctx
Execution context.
Definition IceModel.hh:240
array::Scalar m_basal_melt_rate
rate of production of basal meltwater (ice-equivalent); no ghosts
Definition IceModel.hh:295
virtual void front_retreat_step()
virtual void misc_setup()
Miscellaneous initialization tasks plus tasks that need the fields that can come from regridding.
const Logger::Ptr m_log
Logger.
Definition IceModel.hh:244
std::shared_ptr< array::Forcing > m_surface_input_for_hydrology
Definition IceModel.hh:257
virtual void print_summary(bool tempAndAge)
Definition printout.cc:89
virtual void update_diagnostics(double dt)
Definition IceModel.cc:1030
array::Scalar m_bedtoptemp
temperature at the top surface of the bedrock thermal layer
Definition IceModel.hh:297
const energy::EnergyModel * energy_balance_model() const
Definition IceModel.cc:897
std::string save_state_on_error(const std::string &suffix, const std::set< std::string > &additional_variables)
Definition IceModel.cc:374
VariableMetadata m_output_global_attributes
stores global attributes saved in a PISM output file
Definition IceModel.hh:249
const array::Scalar & frontal_melt() const
Definition IceModel.cc:919
const array::Scalar & forced_retreat() const
Definition IceModel.cc:926
std::vector< std::shared_ptr< array::Scalar2 > > m_work2d
Definition IceModel.hh:393
virtual void process_options()
double dt() const
Definition IceModel.cc:126
bool write_checkpoint()
Write a checkpoint (i.e. an intermediate result of a run).
virtual void allocate_storage()
Allocate all Arrays defined in IceModel.
Definition IceModel.cc:171
std::set< std::string > m_checkpoint_vars
Definition IceModel.hh:459
array::Scalar2 m_velocity_bc_mask
mask to determine Dirichlet boundary locations for the sliding velocity
Definition IceModel.hh:302
void reset_counters()
Definition IceModel.cc:130
std::string m_stdout_flags
Definition IceModel.hh:321
Config::Ptr m_config
Configuration flags and parameters.
Definition IceModel.hh:238
std::unique_ptr< GeometryEvolution > m_geometry_evolution
Definition IceModel.hh:289
virtual void step(bool do_mass_continuity, bool do_skip)
The contents of the main PISM time-step.
Definition IceModel.cc:406
virtual void combine_basal_melt_rate(const Geometry &geometry, const array::Scalar &shelf_base_mass_flux, const array::Scalar &grounded_basal_melt_rate, array::Scalar &result)
Combine basal melt rate in grounded and floating areas.
Definition energy.cc:85
virtual void hydrology_step()
Definition IceModel.cc:689
virtual void time_setup()
Initialize time from an input file or command-line options.
std::shared_ptr< energy::BedThermalUnit > m_btu
Definition IceModel.hh:259
virtual void save_variables(const File &file, OutputKind kind, const std::set< std::string > &variables, double time, io::Type default_diagnostics_type=io::PISM_FLOAT) const
Definition output.cc:150
std::set< array::Array * > m_model_state
Definition IceModel.hh:415
array::Scalar1 m_ice_thickness_bc_mask
Mask prescribing locations where ice thickness is held constant.
Definition IceModel.hh:307
std::shared_ptr< AgeModel > m_age_model
Definition IceModel.hh:262
const bed::BedDef * bed_deformation_model() const
Definition IceModel.cc:905
virtual std::set< std::string > output_variables(const std::string &keyword)
Assembles a list of diagnostics corresponding to an output file size.
double t_TempAge
time of last update for enthalpy/temperature
Definition IceModel.hh:313
array::Vector2 m_velocity_bc_values
Dirichlet boundary velocities.
Definition IceModel.hh:304
virtual void prune_diagnostics()
Definition IceModel.cc:974
std::string m_ts_filename
file to write scalar time-series to
Definition IceModel.hh:436
virtual void update_fracture_density()
std::map< std::string, Diagnostic::Ptr > m_diagnostics
Requested spatially-variable diagnostics.
Definition IceModel.hh:417
IceModel(std::shared_ptr< Grid > grid, const std::shared_ptr< Context > &context)
Definition IceModel.cc:56
std::shared_ptr< ocean::sea_level::SeaLevel > m_sea_level
Definition IceModel.hh:282
virtual void print_summary_line(bool printPrototype, bool tempAndAge, double delta_t, double volume, double area, double meltfrac, double max_diffusivity)
Print a line to stdout which summarizes the state of the modeled ice sheet at the end of the time ste...
Definition printout.cc:164
virtual void update_viewers()
Update the runtime graphical viewers.
Definition viewers.cc:58
const array::Scalar & calving() const
Definition IceModel.cc:912
virtual void pre_step_hook()
Virtual. Does nothing in IceModel. Derived classes can do more computation in each time step.
Definition IceModel.cc:717
std::set< std::string > m_ts_vars
Definition IceModel.hh:439
unsigned int m_skip_countdown
Definition IceModel.hh:317
IceModelTerminationReason run_to(double run_end)
Definition IceModel.cc:736
std::set< std::string > m_extra_vars
Definition IceModel.hh:450
virtual YieldStressInputs yield_stress_inputs()
Definition IceModel.cc:364
virtual void write_metadata(const File &file, MappingTreatment mapping_flag, HistoryTreatment history_flag) const
Write time-independent metadata to a file.
Definition output.cc:72
const YieldStress * basal_yield_stress_model() const
Definition IceModel.cc:901
std::unique_ptr< hydrology::Hydrology > m_subglacial_hydrology
Definition IceModel.hh:254
std::set< std::string > m_output_vars
Definition IceModel.hh:422
const energy::BedThermalUnit * bedrock_thermal_model() const
Definition IceModel.cc:893
void set_python_ocean_model(std::shared_ptr< ocean::PyOceanModel > model)
Definition IceModel.cc:1057
double dt_TempAge
enthalpy/temperature and age time-steps
Definition IceModel.hh:315
double m_dt
mass continuity time step, s
Definition IceModel.hh:311
std::shared_ptr< energy::EnergyModel > m_energy_model
Definition IceModel.hh:260
const std::shared_ptr< Grid > m_grid
Computational grid.
Definition IceModel.hh:236
std::shared_ptr< bed::BedDef > m_beddef
Definition IceModel.hh:284
virtual TimesteppingInfo max_timestep(unsigned int counter)
Use various stability criteria to determine the time step for an evolution run.
virtual void allocate_submodels()
Allocate PISM's sub-models implementing some physical processes.
array::Scalar2 m_basal_yield_stress
ghosted
Definition IceModel.hh:293
void message(int threshold, const char format[],...) const __attribute__((format(printf
Print a message to the log.
Definition Logger.cc:49
A basic logging class.
Definition Logger.hh:40
void begin(const char *name) const
Definition Profiling.cc:75
void end(const char *name) const
Definition Profiling.cc:91
void stage_end(const char *name) const
Definition Profiling.cc:119
void stage_begin(const char *name) const
Definition Profiling.cc:103
void add_context(const std::string &message)
Add a message providing some context. This way we can (sort of) get a stack trace even though C++ exc...
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
VariableMetadata & standard_name(const std::string &input)
VariableMetadata & output_units(const std::string &input)
VariableMetadata & set_output_type(io::Type type)
VariableMetadata & set_time_independent(bool flag)
const array::Scalar * till_water_thickness
const array::Scalar * subglacial_water_thickness
const Geometry * geometry
The PISM basal yield stress model interface (virtual base class)
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition Array.hh:64
void copy_from(const Array2D< T > &source)
Definition Array2D.hh:73
void set_interpolation_type(InterpolationType type)
Definition Array.cc:178
size_t size() const
Return the total number of elements in the owned part of an array.
Definition Array.cc:909
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
Definition Array.cc:224
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition Array.cc:476
bool next_to_ice(int i, int j) const
Ice-free margin (at least one of four neighbors has ice).
Definition CellType.hh:87
PISM bed deformation model (base class).
Definition BedDef.hh:37
Given the temperature of the top of the bedrock, for the duration of one time-step,...
const array::Scalar * surface_liquid_fraction
const array::Scalar1 * ice_thickness
const array::Array3D * w3
const array::Array3D * v3
const array::Scalar * shelf_base_temp
const array::Scalar * basal_heat_flux
const array::Scalar * basal_frictional_heating
const array::Scalar * till_water_thickness
const array::Array3D * u3
const array::Scalar * surface_temp
const array::CellType * cell_type
const array::Array3D * volumetric_heating_rate
const Geometry * geometry
Definition Hydrology.hh:37
const array::Scalar * ice_sliding_speed
Definition Hydrology.hh:41
const array::Scalar * surface_input_rate
Definition Hydrology.hh:39
const array::Scalar1 * no_model_mask
Definition Hydrology.hh:35
const array::Scalar * basal_melt_rate
Definition Hydrology.hh:40
A very rudimentary PISM ocean model.
Definition OceanModel.hh:33
const array::Scalar1 * fracture_density
const array::Scalar * basal_yield_stress
const array::Scalar * bc_mask
const array::Array3D * age
const array::Scalar * water_column_pressure
const array::Array3D * enthalpy
const array::Scalar * basal_melt_rate
const array::Vector * bc_values
The class defining PISM's interface to the shallow stress balance code.
#define PISM_ERROR_LOCATION
@ 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:74
@ PISM_READONLY
open an existing file for reading only
Definition IO_Flags.hh:68
double get_time(MPI_Comm comm)
io::Backend string_to_backend(const std::string &backend)
Definition File.cc:56
std::string set_join(const std::set< std::string > &input, const std::string &separator)
std::map< std::string, TSDiagnostic::Ptr > TSDiagnosticList
std::string printf(const char *format,...)
IceModelTerminationReason
Definition IceModel.hh:114
@ PISM_SIGNAL
Definition IceModel.hh:114
@ PISM_DONE
Definition IceModel.hh:114
@ PISM_CHEKPOINT
Definition IceModel.hh:114
std::set< std::string > set_split(const std::string &input, char separator)
Transform a separator-separated list (a string) into a set of strings.
void warn_about_missing(const Logger &log, const std::set< std::string > &vars, const std::string &type, const std::set< std::string > &available, bool stop)
Definition IceModel.cc:930
std::string filename_add_suffix(const std::string &filename, const std::string &separator, const std::string &suffix)
Adds a suffix to a filename.
std::string version()
bool member(const std::string &string, const std::set< std::string > &set)
T combine(const T &a, const T &b)
std::string join(const std::vector< std::string > &strings, const std::string &separator)
Concatenate strings, inserting separator between elements.
volatile sig_atomic_t pism_signal
Definition pism_signal.c:23
void pism_signal_handler(int sig)
Definition pism_signal.c:25
std::string filename
Definition options.hh:33
ThicknessChanges(const std::shared_ptr< const Grid > &grid)
Definition IceModel.cc:1050