tIceModel.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tIceModel.hh (15136B)
---
1 // Copyright (C) 2004-2020 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 #ifndef __iceModel_hh
20 #define __iceModel_hh
21
22 //! \file IceModel.hh Definition of class IceModel.
23 /*! \file IceModel.hh
24 IceModel is a big class which is an ice flow model. It contains all parts that
25 are not well-defined, separated components. Such components are better places
26 to put sub-models that have a clear, general interface to the rest of an ice
27 sheet model.
28
29 IceModel has pointers to well-defined components, when they exist.
30
31 IceModel generally interprets user options, and initializes components based on
32 such options. It manages the initialization sequences (%e.g. a restart from a
33 file containing a complete model state, versus bootstrapping).
34 */
35
36 #include <map>
37 #include <set>
38 #include <string>
39 #include <vector>
40 #include <memory>
41
42 // IceModel owns a bunch of fields, so we have to include this.
43 #include "pism/util/iceModelVec.hh"
44 #include "pism/util/IceModelVec2CellType.hh"
45 #include "pism/util/ConfigInterface.hh"
46 #include "pism/util/Context.hh"
47 #include "pism/util/Logger.hh"
48 #include "pism/util/Time.hh"
49 #include "pism/util/Diagnostic.hh"
50 #include "pism/util/MaxTimestep.hh"
51 #include "pism/geometry/Geometry.hh"
52 #include "pism/geometry/GeometryEvolution.hh"
53 #include "pism/stressbalance/StressBalance.hh"
54 #include "pism/basalstrength/YieldStress.hh"
55
56 namespace pism {
57
58 namespace ocean {
59 class OceanModel;
60 namespace sea_level {
61 class SeaLevel;
62 }
63 }
64
65 namespace surface {
66 class SurfaceModel;
67 }
68
69 namespace hydrology {
70 class Hydrology;
71 }
72
73 namespace calving {
74 class EigenCalving;
75 class vonMisesCalving;
76 class FloatKill;
77 class HayhurstCalving;
78 class CalvingAtThickness;
79 class IcebergRemover;
80 }
81
82 class FractureDensity;
83
84 namespace energy {
85 class BedThermalUnit;
86 class Inputs;
87 class EnergyModelStats;
88 class EnergyModel;
89 }
90
91 namespace frontalmelt {
92 class FrontalMelt;
93 }
94
95 namespace bed {
96 class BedDef;
97 }
98
99 class IceGrid;
100 class AgeModel;
101 class IceModelVec2CellType;
102 class IceModelVec2T;
103 class Component;
104 class FrontRetreat;
105 class PrescribedRetreat;
106
107 //! The base class for PISM. Contains all essential variables, parameters, and flags for modelling
108 //! an ice sheet.
109 class IceModel {
110 public:
111 IceModel(IceGrid::Ptr g, Context::Ptr context);
112
113 // the destructor must be virtual merely because some members are virtual
114 virtual ~IceModel();
115
116 IceGrid::Ptr grid() const;
117 Context::Ptr ctx() const;
118
119 void init();
120
121 /** Run PISM in the "standalone" mode. */
122 virtual void run();
123
124 /** Advance the current PISM run to a specific time */
125 virtual void run_to(double time);
126
127 virtual void save_results();
128
129 void list_diagnostics() const;
130 void list_diagnostics_json() const;
131 std::map<std::string, std::vector<VariableMetadata>> describe_diagnostics() const;
132 std::map<std::string, std::vector<VariableMetadata>> describe_ts_diagnostics() const;
133
134 const IceModelVec2S &calving() const;
135 const IceModelVec2S &frontal_melt() const;
136 const IceModelVec2S &forced_retreat() const;
137
138 double ice_volume_temperate(double thickness_threshold) const;
139 double ice_volume_cold(double thickness_threshold) const;
140 double ice_area_temperate(double thickness_threshold) const;
141 double ice_area_cold(double thickness_threshold) const;
142
143 const stressbalance::StressBalance* stress_balance() const;
144 const ocean::OceanModel* ocean_model() const;
145 const frontalmelt::FrontalMelt* frontalmelt_model() const;
146 const bed::BedDef* bed_model() const;
147 const energy::BedThermalUnit* bedrock_thermal_model() const;
148 const energy::EnergyModel* energy_balance_model() const;
149
150 const Geometry& geometry() const;
151 const GeometryEvolution& geometry_evolution() const;
152
153 double dt() const;
154
155 protected:
156 virtual void allocate_submodels();
157 virtual void allocate_stressbalance();
158 virtual void allocate_age_model();
159 virtual void allocate_bed_deformation();
160 virtual void allocate_bedrock_thermal_unit();
161 virtual void allocate_energy_model();
162 virtual void allocate_subglacial_hydrology();
163 virtual void allocate_basal_yield_stress();
164 virtual void allocate_couplers();
165 virtual void allocate_geometry_evolution();
166 virtual void allocate_iceberg_remover();
167
168 virtual stressbalance::Inputs stress_balance_inputs();
169
170 virtual energy::Inputs energy_model_inputs();
171
172 virtual YieldStressInputs yield_stress_inputs();
173
174 virtual void time_setup();
175 virtual void model_state_setup();
176 virtual void misc_setup();
177 virtual void init_diagnostics();
178 virtual void init_calving();
179 virtual void init_frontal_melt();
180 virtual void init_front_retreat();
181 virtual void prune_diagnostics();
182 virtual void update_diagnostics(double dt);
183 virtual void reset_diagnostics();
184
185 virtual void step(bool do_mass_continuity, bool do_skip);
186 virtual void pre_step_hook();
187 virtual void post_step_hook();
188
189 void reset_counters();
190
191 // see iMbootstrap.cc
192 virtual void bootstrap_2d(const File &input_file);
193
194 // see iMoptions.cc
195 virtual void process_options();
196 virtual std::set<std::string> output_variables(const std::string &keyword);
197
198 virtual void compute_lat_lon();
199
200 // see iMIO.cc
201 virtual void restart_2d(const File &input_file, unsigned int record);
202 virtual void initialize_2d();
203
204 enum OutputKind {INCLUDE_MODEL_STATE = 0, JUST_DIAGNOSTICS};
205 virtual void save_variables(const File &file,
206 OutputKind kind,
207 const std::set<std::string> &variables,
208 double time,
209 IO_Type default_diagnostics_type = PISM_FLOAT);
210
211 virtual void define_model_state(const File &file);
212 virtual void write_model_state(const File &file);
213
214 enum HistoryTreatment {OVERWRITE_HISTORY = 0, PREPEND_HISTORY};
215 enum MappingTreatment {WRITE_MAPPING = 0, SKIP_MAPPING};
216 virtual void write_metadata(const File &file, MappingTreatment mapping_flag,
217 HistoryTreatment history_flag);
218
219 virtual void write_mapping(const File &file);
220 virtual void write_run_stats(const File &file);
221
222
223 virtual void define_diagnostics(const File &file,
224 const std::set<std::string> &variables,
225 IO_Type default_type);
226 virtual void write_diagnostics(const File &file,
227 const std::set<std::string> &variables);
228
229 //! Computational grid
230 const IceGrid::Ptr m_grid;
231 //! Configuration flags and parameters
232 const Config::Ptr m_config;
233 //! Execution context
234 const Context::Ptr m_ctx;
235 //! Unit system
236 const units::System::Ptr m_sys;
237 //! Logger
238 const Logger::Ptr m_log;
239 //! Time manager
240 const Time::Ptr m_time;
241
242 //! stores global attributes saved in a PISM output file
243 VariableMetadata m_output_global_attributes;
244
245 //! run statistics
246 VariableMetadata m_run_stats;
247
248 //! the list of sub-models, for writing model states and obtaining diagnostics
249 std::map<std::string,const Component*> m_submodels;
250
251 std::unique_ptr<hydrology::Hydrology> m_subglacial_hydrology;
252 std::shared_ptr<YieldStress> m_basal_yield_stress_model;
253
254 std::shared_ptr<IceModelVec2T> m_surface_input_for_hydrology;
255
256 energy::BedThermalUnit *m_btu;
257 energy::EnergyModel *m_energy_model;
258
259 std::shared_ptr<AgeModel> m_age_model;
260
261 std::shared_ptr<calving::IcebergRemover> m_iceberg_remover;
262 std::shared_ptr<calving::FloatKill> m_float_kill_calving;
263 std::shared_ptr<calving::CalvingAtThickness> m_thickness_threshold_calving;
264 std::shared_ptr<calving::EigenCalving> m_eigen_calving;
265 std::shared_ptr<calving::HayhurstCalving> m_hayhurst_calving;
266 std::shared_ptr<calving::vonMisesCalving> m_vonmises_calving;
267 std::shared_ptr<PrescribedRetreat> m_prescribed_retreat;
268
269 std::shared_ptr<FrontRetreat> m_front_retreat;
270
271 std::shared_ptr<surface::SurfaceModel> m_surface;
272 std::shared_ptr<ocean::OceanModel> m_ocean;
273 std::shared_ptr<frontalmelt::FrontalMelt> m_frontal_melt;
274 std::shared_ptr<ocean::sea_level::SeaLevel> m_sea_level;
275
276 bed::BedDef *m_beddef;
277
278 // state variables and some diagnostics/internals
279
280 Geometry m_geometry;
281 std::unique_ptr<GeometryEvolution> m_geometry_evolution;
282 bool m_new_bed_elevation;
283
284 //! ghosted
285 IceModelVec2S m_basal_yield_stress;
286 //! rate of production of basal meltwater (ice-equivalent); no ghosts
287 IceModelVec2S m_basal_melt_rate;
288 //! temperature at the top surface of the bedrock thermal layer
289 IceModelVec2S m_bedtoptemp;
290
291 std::shared_ptr<FractureDensity> m_fracture;
292
293 //! mask to determine Dirichlet boundary locations
294 IceModelVec2Int m_ssa_dirichlet_bc_mask;
295 //! Dirichlet boundary velocities
296 IceModelVec2V m_ssa_dirichlet_bc_values;
297
298 // parameters
299 //! mass continuity time step, s
300 double m_dt;
301 //! time of last update for enthalpy/temperature
302 double t_TempAge;
303 //! enthalpy/temperature and age time-steps
304 double dt_TempAge;
305
306 unsigned int m_skip_countdown;
307
308 std::string m_adaptive_timestep_reason;
309
310 std::string m_stdout_flags;
311
312 // see iceModel.cc
313 virtual void allocate_storage();
314
315 virtual MaxTimestep max_timestep_diffusivity();
316 virtual void max_timestep(double &dt_result, unsigned int &skip_counter);
317 virtual unsigned int skip_counter(double input_dt, double input_dt_diffusivity);
318
319 // see energy.cc
320 virtual void bedrock_thermal_model_step();
321 virtual void energy_step();
322
323 virtual void hydrology_step();
324
325 virtual void combine_basal_melt_rate(const Geometry &geometry,
326 const IceModelVec2S &shelf_base_mass_flux,
327 const IceModelVec2S &grounded_basal_melt_rate,
328 IceModelVec2S &result);
329
330 enum ConsistencyFlag {REMOVE_ICEBERGS, DONT_REMOVE_ICEBERGS};
331 void enforce_consistency_of_geometry(ConsistencyFlag flag);
332
333 virtual void front_retreat_step();
334
335 void compute_geometry_change(const IceModelVec2S &thickness,
336 const IceModelVec2S &Href,
337 const IceModelVec2S &thickness_old,
338 const IceModelVec2S &Href_old,
339 bool add_values,
340 IceModelVec2S &output);
341
342 // see iMIO.cc
343 virtual void regrid();
344
345 // see iMfractures.cc
346 virtual void update_fracture_density();
347
348 // see iMreport.cc
349 virtual double compute_temperate_base_fraction(double ice_area);
350 virtual double compute_original_ice_fraction(double ice_volume);
351 virtual void print_summary(bool tempAndAge);
352 virtual void print_summary_line(bool printPrototype, bool tempAndAge,
353 double delta_t,
354 double volume, double area,
355 double meltfrac, double max_diffusivity);
356
357
358 // see iMutil.cc
359 virtual int process_signals();
360 virtual void prepend_history(const std::string &string);
361 virtual void update_run_stats();
362
363 // working space (a convenience)
364 static const int m_n_work2d = 3;
365 mutable IceModelVec2S m_work2d[m_n_work2d];
366
367 std::shared_ptr<stressbalance::StressBalance> m_stress_balance;
368
369 struct ThicknessChanges {
370 ThicknessChanges(IceGrid::ConstPtr grid);
371
372 // calving during the last time step
373 IceModelVec2S calving;
374
375 // frontal melt during the last time step
376 IceModelVec2S frontal_melt;
377
378 // parameterized retreat
379 IceModelVec2S forced_retreat;
380 };
381
382 ThicknessChanges m_thickness_change;
383
384 /*!
385 * The set of variables that the "state" of IceModel consists of.
386 */
387 std::set<IceModelVec*> m_model_state;
388 //! Requested spatially-variable diagnostics.
389 std::map<std::string,Diagnostic::Ptr> m_diagnostics;
390 //! Requested scalar diagnostics.
391 std::map<std::string,TSDiagnostic::Ptr> m_ts_diagnostics;
392
393 // Set of variables to put in the output file:
394 std::set<std::string> m_output_vars;
395
396 // This is related to the snapshot saving feature
397 std::string m_snapshots_filename;
398 bool m_save_snapshots, m_snapshots_file_is_ready, m_split_snapshots;
399 std::vector<double> m_snapshot_times;
400 std::set<std::string> m_snapshot_vars;
401 unsigned int m_current_snapshot;
402 void init_snapshots();
403 void write_snapshot();
404 MaxTimestep save_max_timestep(double my_t);
405
406 //! file to write scalar time-series to
407 std::string m_ts_filename;
408 //! requested times for scalar time-series
409 std::shared_ptr<std::vector<double>> m_ts_times;
410 std::set<std::string> m_ts_vars;
411 void init_timeseries();
412 void flush_timeseries();
413 MaxTimestep ts_max_timestep(double my_t);
414
415 // spatially-varying time-series
416 bool m_save_extra, m_extra_file_is_ready, m_split_extra;
417 std::string m_extra_filename;
418 std::vector<double> m_extra_times;
419 unsigned int m_next_extra;
420 double m_last_extra;
421 std::set<std::string> m_extra_vars;
422 TimeBoundsMetadata m_extra_bounds;
423 std::unique_ptr<File> m_extra_file;
424 void init_extras();
425 void write_extras();
426 MaxTimestep extras_max_timestep(double my_t);
427
428 // automatic backups
429 std::string m_backup_filename;
430 double m_last_backup_time;
431 std::set<std::string> m_backup_vars;
432 void init_backups();
433 void write_backup();
434
435 // last time at which PISM hit a multiple of X years, see the configuration parameter
436 // time_stepping.hit_multiples
437 double m_timestep_hit_multiples_last_time;
438
439 // diagnostic viewers; see iMviewers.cc
440 virtual void update_viewers();
441 virtual void view_field(const IceModelVec *field);
442 std::map<std::string,petsc::Viewer::Ptr> m_viewers;
443
444 private:
445 TimeseriesMetadata m_timestamp;
446 double m_start_time; // this is used in the wall-clock-time backup code
447 };
448
449 MaxTimestep reporting_max_timestep(const std::vector<double> ×, double t,
450 const std::string &description);
451
452 void check_minimum_ice_thickness(const IceModelVec2S &ice_thickness);
453 bool check_maximum_ice_thickness(const IceModelVec2S &ice_thickness);
454
455 void bedrock_surface_temperature(const IceModelVec2S &sea_level,
456 const IceModelVec2CellType &cell_type,
457 const IceModelVec2S &bed_topography,
458 const IceModelVec2S &ice_thickness,
459 const IceModelVec2S &basal_enthalpy,
460 const IceModelVec2S &ice_surface_temperature,
461 IceModelVec2S &result);
462
463 } // end of namespace pism
464
465 #endif /* __iceModel_hh */