24#include "pism/pism_config.hh"
26#include "pism/icemodel/IceModel.hh"
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"
52#include "pism/coupler/ocean/PyOceanModel.hh"
58 m_config(context->config()),
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),
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"),
71 m_ice_thickness_bc_mask(grid,
"thk_bc_mask"),
73 m_thickness_change(grid),
74 m_ts_times(new std::vector<
double>()) {
107 auto surface_input_file =
m_config->get_string(
"hydrology.surface_input.file");
108 if (not surface_input_file.empty()) {
110 int buffer_size =
static_cast<int>(
m_config->get_number(
"input.forcing.buffer_size"));
115 std::make_shared<array::Forcing>(
m_grid, file,
"water_input_rate",
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");
136 int year_increment =
static_cast<int>(
m_config->get_number(
"time_stepping.hit_multiples"));
138 if (year_increment > 0) {
142 int last_multiple = year - year % year_increment;
145 last_multiple -= year_increment *
static_cast<int>(year % year_increment < 0);
182 if (
m_config->get_flag(
"geometry.grounded_cell_fraction")) {
186 if (
m_config->get_flag(
"geometry.part_grid.enabled")) {
194 .
long_name(
"yield stress for basal till (plastic or pseudo-plastic model)")
201 .
long_name(
"temperature at the top surface of the bedrock thermal layer")
208 "ice basal melt rate from energy conservation and subshelf melt, in ice thickness per time")
218 .
long_name(
"Mask prescribing Dirichlet boundary locations for the sliding velocity")
228 double fill_value =
m_config->get_number(
"output.fill_value");
229 const double huge_value = 1e6;
232 "X-component of the SSA velocity boundary conditions");
234 "Y-component of the SSA velocity boundary conditions");
235 for (
int j : { 0, 1 }) {
245 .
long_name(
"Mask specifying locations where ice thickness is held constant")
305 for (
auto p =
m_grid->points(); p; p.next()) {
306 const int i = p.i(), j = p.j();
318 if (
m_config->get_flag(
"geometry.update.use_basal_melt_rate")) {
330 if (
m_config->get_flag(
"stress_balance.ssa.dirichlet_bc")) {
335 if (
m_config->get_flag(
"fracture_density.enabled")) {
375 const std::set<std::string> &additional_variables) {
376 std::string filename =
m_config->get_string(
"output.file");
378 if (filename.empty()) {
379 m_log->message(2,
"WARNING: output.file is empty. Using unnamed.nc instead.");
380 filename =
"unnamed.nc";
393 for (
const auto &v : additional_variables) {
413 double current_time =
m_time->current();
436 m_ocean->shelf_base_mass_flux(),
442 profiling.
begin(
"stress_balance");
444 profiling.
end(
"stress_balance");
448 e.
add_context(
"performing a time step. (Note: Model state was saved to '%s'.)",
449 output_file.c_str());
466 profiling.
begin(
"basal_yield_stress");
468 profiling.
end(
"basal_yield_stress");
485 profiling.
begin(
"age");
487 profiling.
end(
"age");
497 profiling.
begin(
"energy");
499 profiling.
end(
"energy");
506 if (
m_config->get_flag(
"fracture_density.enabled")) {
507 profiling.
begin(
"fracture_density");
509 profiling.
end(
"fracture_density");
515 if (do_mass_continuity) {
519 profiling.
begin(
"mass_transport");
545 {
"flux_staggered",
"flux_divergence"});
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());
556 profiling.
end(
"mass_transport");
559 profiling.
begin(
"front_retreat");
561 profiling.
end(
"front_retreat");
568 profiling.
begin(
"sea_level");
570 profiling.
end(
"sea_level");
572 profiling.
begin(
"ocean");
574 profiling.
end(
"ocean");
581 profiling.
begin(
"surface");
583 profiling.
end(
"surface");
586 if (do_mass_continuity) {
607 bool add_values =
true;
627 profiling.
begin(
"basal_hydrology");
629 profiling.
end(
"basal_hydrology");
634 int topg_state_counter =
m_beddef->bed_elevation().state_counter();
636 profiling.
begin(
"bed_deformation");
640 profiling.
end(
"bed_deformation");
647 if (
m_config->get_flag(
"time_stepping.assume_bed_elevation_changed")) {
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"
677 "-i %s -bootstrap -regrid_file %s -allow_extrapolation -Lz N [other options]\n"
679 m_grid->Lz(), o_file.c_str(), o_file.c_str(), o_file.c_str(),
m_grid->Lz());
705 }
else if (
m_config->get_flag(
"hydrology.surface_input_from_runoff")) {
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");
771 const double time =
m_time->current();
773 d.second->update(time, time);
777 m_log->message(2,
"running forward ...\n");
795 step(do_mass_conserve, do_skip);
801 bool tempAgeStep = updateAtDepth and (
m_age_model or do_energy);
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;
811 profiling.
begin(
"io");
817 if (stop_after_chekpoint) {
827 profiling.
stage_end(
"time-stepping loop");
829 return termination_reason;
844 profiling.
begin(
"initialization");
866 m_log->message(2,
"Computational domain and grid:\n");
867 m_grid->report_parameters();
874 profiling.
end(
"initialization");
931 const std::set<std::string> &vars,
932 const std::string &type,
933 const std::set<std::string> &available,
935 std::vector<std::string> missing;
936 for (
const auto &v : vars) {
937 if (available.find(v) == available.end()) {
938 missing.push_back(v);
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";
948 "%s variable%s %s %s not available!\n"
949 "Available variables:\n- %s",
952 join(missing,
",").c_str(),
954 set_join(available,
",\n- ").c_str());
958 "\nWARNING: %s variable%s %s %s not available!\n\n",
961 join(missing,
",").c_str(),
977 std::set<std::string> available;
979 available.insert(d.first);
982 auto extra_stop =
m_config->get_flag(
"output.extra.stop_missing");
986 auto requested =
set_split(
m_config->get_string(
"output.runtime.viewer.variables"),
',');
993 for (
const auto &v : available) {
994 if (requested.find(v) == requested.end()) {
1000 std::vector<std::string> missing;
1009 missing.push_back(v);
1016 if (not missing.empty()) {
1018 "requested scalar diagnostics %s are not available",
1019 join(missing,
",").c_str());
1032 d.second->update(
dt);
1035 const double time =
m_time->current();
1037 d.second->update(time -
dt, time);
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") {
1058 m_ocean = std::make_shared<ocean::PyOceanModelAdapter>(
m_grid, model);
High-level PISM I/O class.
array::Scalar1 sea_level_elevation
array::Scalar cell_grounded_fraction
void ensure_consistency(double ice_free_thickness_threshold)
array::Scalar2 ice_surface_elevation
array::Scalar1 ice_area_specific_volume
array::CellType2 cell_type
array::Scalar2 ice_thickness
array::Scalar2 bed_elevation
std::string m_adaptive_timestep_reason
std::map< std::string, const Component * > m_submodels
the list of sub-models, for writing model states and obtaining diagnostics
virtual energy::Inputs energy_model_inputs()
std::set< std::string > m_snapshot_vars
virtual int process_signals()
Catch signals -USR1, -USR2 and -TERM.
void enforce_consistency_of_geometry(ConsistencyFlag flag)
Update the surface elevation and the flow-type mask when the geometry has changed.
const Geometry & geometry() const
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.
std::shared_ptr< stressbalance::StressBalance > m_stress_balance
std::shared_ptr< Isochrones > m_isochrones
IceModelTerminationReason run()
std::shared_ptr< ocean::OceanModel > m_ocean
virtual void post_step_hook()
Virtual. Does nothing in IceModel. Derived classes can do more computation in each time step.
const ocean::OceanModel * ocean_model() const
std::map< std::string, TSDiagnostic::Ptr > m_ts_diagnostics
Requested scalar diagnostics.
virtual void energy_step()
Manage the solution of the energy equation, and related parallel communication.
std::shared_ptr< surface::SurfaceModel > m_surface
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
virtual stressbalance::Inputs stress_balance_inputs()
std::shared_ptr< FractureDensity > m_fracture
ThicknessChanges m_thickness_change
void init()
Manage the initialization of the IceModel object.
const Time::Ptr m_time
Time manager.
const GeometryEvolution & geometry_evolution() const
const stressbalance::StressBalance * stress_balance() const
virtual void reset_diagnostics()
static const int m_n_work2d
std::shared_ptr< YieldStress > m_basal_yield_stress_model
std::shared_ptr< calving::IcebergRemover > m_iceberg_remover
double m_timestep_hit_multiples_last_time
std::shared_ptr< Context > m_ctx
Execution context.
array::Scalar m_basal_melt_rate
rate of production of basal meltwater (ice-equivalent); no ghosts
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.
std::shared_ptr< array::Forcing > m_surface_input_for_hydrology
virtual void print_summary(bool tempAndAge)
virtual void update_diagnostics(double dt)
array::Scalar m_bedtoptemp
temperature at the top surface of the bedrock thermal layer
const energy::EnergyModel * energy_balance_model() const
std::string save_state_on_error(const std::string &suffix, const std::set< std::string > &additional_variables)
VariableMetadata m_output_global_attributes
stores global attributes saved in a PISM output file
const array::Scalar & frontal_melt() const
const array::Scalar & forced_retreat() const
std::vector< std::shared_ptr< array::Scalar2 > > m_work2d
virtual void process_options()
bool write_checkpoint()
Write a checkpoint (i.e. an intermediate result of a run).
virtual void allocate_storage()
Allocate all Arrays defined in IceModel.
std::set< std::string > m_checkpoint_vars
array::Scalar2 m_velocity_bc_mask
mask to determine Dirichlet boundary locations for the sliding velocity
std::string m_stdout_flags
Config::Ptr m_config
Configuration flags and parameters.
std::unique_ptr< GeometryEvolution > m_geometry_evolution
virtual void step(bool do_mass_continuity, bool do_skip)
The contents of the main PISM time-step.
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.
virtual void hydrology_step()
virtual void time_setup()
Initialize time from an input file or command-line options.
std::shared_ptr< energy::BedThermalUnit > m_btu
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
std::set< array::Array * > m_model_state
array::Scalar1 m_ice_thickness_bc_mask
Mask prescribing locations where ice thickness is held constant.
std::shared_ptr< AgeModel > m_age_model
const bed::BedDef * bed_deformation_model() const
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
array::Vector2 m_velocity_bc_values
Dirichlet boundary velocities.
virtual void prune_diagnostics()
std::string m_ts_filename
file to write scalar time-series to
virtual void update_fracture_density()
std::map< std::string, Diagnostic::Ptr > m_diagnostics
Requested spatially-variable diagnostics.
IceModel(std::shared_ptr< Grid > grid, const std::shared_ptr< Context > &context)
std::shared_ptr< ocean::sea_level::SeaLevel > m_sea_level
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...
virtual void update_viewers()
Update the runtime graphical viewers.
const array::Scalar & calving() const
virtual void pre_step_hook()
Virtual. Does nothing in IceModel. Derived classes can do more computation in each time step.
std::set< std::string > m_ts_vars
unsigned int m_skip_countdown
IceModelTerminationReason run_to(double run_end)
std::set< std::string > m_extra_vars
virtual YieldStressInputs yield_stress_inputs()
virtual void write_metadata(const File &file, MappingTreatment mapping_flag, HistoryTreatment history_flag) const
Write time-independent metadata to a file.
const YieldStress * basal_yield_stress_model() const
std::unique_ptr< hydrology::Hydrology > m_subglacial_hydrology
std::set< std::string > m_output_vars
const energy::BedThermalUnit * bedrock_thermal_model() const
void set_python_ocean_model(std::shared_ptr< ocean::PyOceanModel > model)
double dt_TempAge
enthalpy/temperature and age time-steps
double m_dt
mass continuity time step, s
std::shared_ptr< energy::EnergyModel > m_energy_model
const std::shared_ptr< Grid > m_grid
Computational grid.
std::shared_ptr< bed::BedDef > m_beddef
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
void message(int threshold, const char format[],...) const __attribute__((format(printf
Print a message to the log.
void begin(const char *name) const
void end(const char *name) const
void stage_end(const char *name) const
void stage_begin(const char *name) const
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
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.
void copy_from(const Array2D< T > &source)
void set_interpolation_type(InterpolationType type)
size_t size() const
Return the total number of elements in the owned part of an array.
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
void set(double c)
Result: v[j] <- c for all j.
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
bool next_to_ice(int i, int j) const
Ice-free margin (at least one of four neighbors has ice).
PISM bed deformation model (base class).
Given the temperature of the top of the bedrock, for the duration of one time-step,...
A very rudimentary PISM ocean model.
The class defining PISM's interface to the shallow stress balance code.
#define PISM_ERROR_LOCATION
@ PISM_READWRITE_MOVE
create a file for writing, move foo.nc to foo.nc~ if present
@ PISM_READONLY
open an existing file for reading only
double get_time(MPI_Comm comm)
io::Backend string_to_backend(const std::string &backend)
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
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)
std::string filename_add_suffix(const std::string &filename, const std::string &separator, const std::string &suffix)
Adds a suffix to a filename.
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
void pism_signal_handler(int sig)
array::Scalar forced_retreat
ThicknessChanges(const std::shared_ptr< const Grid > &grid)
array::Scalar frontal_melt