Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.2-d6b3a29ca committed by Constantine Khrulev on 2025-03-28
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
EnergyModel.cc
Go to the documentation of this file.
1/* Copyright (C) 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024 PISM Authors
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
20#include "pism/energy/EnergyModel.hh"
21#include "pism/energy/utilities.hh"
22#include "pism/stressbalance/StressBalance.hh"
23#include "pism/util/MaxTimestep.hh"
24#include "pism/util/Profiling.hh"
25#include "pism/util/Vars.hh"
26#include "pism/util/array/CellType.hh"
27#include "pism/util/error_handling.hh"
28#include "pism/util/io/File.hh"
29#include "pism/util/pism_utilities.hh"
30
31namespace pism {
32namespace energy {
33
34static void check_input(const array::Array *ptr, const char *name) {
35 if (ptr == NULL) {
36 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "energy balance model input %s was not provided", name);
37 }
38}
39
42 basal_heat_flux = NULL;
43 cell_type = NULL;
44 ice_thickness = NULL;
46 shelf_base_temp = NULL;
47 surface_temp = NULL;
49
51 u3 = NULL;
52 v3 = NULL;
53 w3 = NULL;
54
55 no_model_mask = NULL;
56}
57
58void Inputs::check() const {
59 check_input(cell_type, "cell_type");
60 check_input(basal_frictional_heating, "basal_frictional_heating");
61 check_input(basal_heat_flux, "basal_heat_flux");
62 check_input(ice_thickness, "ice_thickness");
63 check_input(surface_liquid_fraction, "surface_liquid_fraction");
64 check_input(shelf_base_temp, "shelf_base_temp");
65 check_input(surface_temp, "surface_temp");
66 check_input(till_water_thickness, "till_water_thickness");
67
68 check_input(volumetric_heating_rate, "volumetric_heating_rate");
69 check_input(u3, "u3");
70 check_input(v3, "v3");
71 check_input(w3, "w3");
72}
73
80
88
89
90bool marginal(const array::Scalar1 &thickness, int i, int j, double threshold) {
91 auto H = thickness.box(i, j);
92
93 return ((H.e < threshold) or (H.ne < threshold) or (H.n < threshold) or (H.nw < threshold) or
94 (H.w < threshold) or (H.sw < threshold) or (H.s < threshold) or (H.se < threshold));
95}
96
97
104
105
106EnergyModel::EnergyModel(std::shared_ptr<const Grid> grid,
107 std::shared_ptr<const stressbalance::StressBalance> stress_balance)
108 : Component(grid),
109 m_ice_enthalpy(m_grid, "enthalpy", array::WITH_GHOSTS, m_grid->z(),
110 m_config->get_number("grid.max_stencil_width")),
111 m_work(m_grid, "work_vector", array::WITHOUT_GHOSTS, m_grid->z()),
112 m_basal_melt_rate(m_grid, "basal_melt_rate_grounded"),
113 m_stress_balance(stress_balance) {
114
115 // POSSIBLE standard name = land_ice_enthalpy
117 .long_name("ice enthalpy (includes sensible heat, latent heat, pressure)")
118 .units("J kg^-1");
119
120 {
121 // ghosted to allow the "redundant" computation of tauc
123 .long_name(
124 "ice basal melt rate from energy conservation, in ice thickness per time (valid in grounded areas)")
125 .units("m s^-1");
126 // We could use land_ice_basal_melt_rate, but that way both basal_melt_rate_grounded and bmelt
127 // have this standard name.
128 m_basal_melt_rate.metadata()["comment"] = "positive basal melt rate corresponds to ice loss";
129 }
130
131 // a 3d work vector
132 m_work.metadata(0).long_name("usually new values of temperature or enthalpy during time step");
133}
134
135void EnergyModel::init_enthalpy(const File &input_file, bool do_regrid, int record) {
136
137 if (input_file.variable_exists("enthalpy")) {
138 if (do_regrid) {
140 } else {
141 m_ice_enthalpy.read(input_file, record);
142 }
143 } else if (input_file.variable_exists("temp")) {
144 array::Array3D &temp = m_work, &liqfrac = m_ice_enthalpy;
145
146 {
147 temp.set_name("temp");
148 temp.metadata(0).set_name("temp");
149 temp.metadata(0)
150 .long_name("ice temperature")
151 .units("kelvin")
152 .standard_name("land_ice_temperature");
153
154 if (do_regrid) {
155 temp.regrid(input_file, io::Default::Nil());
156 } else {
157 temp.read(input_file, record);
158 }
159 }
160
161 const array::Scalar &ice_thickness = *m_grid->variables().get_2d_scalar("land_ice_thickness");
162
163 if (input_file.variable_exists("liqfrac")) {
164 SpatialVariableMetadata enthalpy_metadata = m_ice_enthalpy.metadata();
165
166 liqfrac.set_name("liqfrac");
167 liqfrac.metadata(0).set_name("liqfrac");
168 liqfrac.metadata(0).long_name("ice liquid water fraction").units("1");
169
170 if (do_regrid) {
171 liqfrac.regrid(input_file, io::Default::Nil());
172 } else {
173 liqfrac.read(input_file, record);
174 }
175
176 m_ice_enthalpy.set_name(enthalpy_metadata.get_name());
177 m_ice_enthalpy.metadata() = enthalpy_metadata;
178
179 m_log->message(2,
180 " - Computing enthalpy using ice temperature,"
181 " liquid water fraction and thickness...\n");
182 compute_enthalpy(temp, liqfrac, ice_thickness, m_ice_enthalpy);
183 } else {
184 m_log->message(2,
185 " - Computing enthalpy using ice temperature and thickness "
186 "and assuming zero liquid water fraction...\n");
187 compute_enthalpy_cold(temp, ice_thickness, m_ice_enthalpy);
188 }
189 } else {
191 "neither enthalpy nor temperature was found in '%s'.\n",
192 input_file.name().c_str());
193 }
194}
195
196/*!
197 * The `-regrid_file` may contain enthalpy, temperature, or *both* temperature and liquid water
198 * fraction.
199 */
201
202 auto regrid_filename = m_config->get_string("input.regrid.file");
203 auto regrid_vars = set_split(m_config->get_string("input.regrid.vars"), ',');
204
205
206 if (regrid_filename.empty()) {
207 return;
208 }
209
210 std::string enthalpy_name = m_ice_enthalpy.metadata().get_name();
211
212 if (regrid_vars.empty() or member(enthalpy_name, regrid_vars)) {
213 File regrid_file(m_grid->com, regrid_filename, io::PISM_GUESS, io::PISM_READONLY);
214 init_enthalpy(regrid_file, true, 0);
215 }
216}
217
218
219void EnergyModel::restart(const File &input_file, int record) {
220 this->restart_impl(input_file, record);
221}
222
223void EnergyModel::bootstrap(const File &input_file,
224 const array::Scalar &ice_thickness,
225 const array::Scalar &surface_temperature,
226 const array::Scalar &climatic_mass_balance,
227 const array::Scalar &basal_heat_flux) {
228 this->bootstrap_impl(input_file,
229 ice_thickness, surface_temperature,
230 climatic_mass_balance, basal_heat_flux);
231}
232
233void EnergyModel::initialize(const array::Scalar &basal_melt_rate,
234 const array::Scalar &ice_thickness,
235 const array::Scalar &surface_temperature,
236 const array::Scalar &climatic_mass_balance,
237 const array::Scalar &basal_heat_flux) {
238 this->initialize_impl(basal_melt_rate,
239 ice_thickness,
240 surface_temperature,
241 climatic_mass_balance,
242 basal_heat_flux);
243}
244
245void EnergyModel::update(double t, double dt, const Inputs &inputs) {
246 // reset standard out flags at the beginning of every time step
247 m_stdout_flags = "";
249
250 profiling().begin("ice_energy");
251 {
252 // this call should fill m_work with new values of enthalpy
253 this->update_impl(t, dt, inputs);
254
256 }
257 profiling().end("ice_energy");
258
259 // globalize m_stats and update m_stdout_flags
260 {
261 m_stats.sum(m_grid->com);
262
263 if (m_stats.reduced_accuracy_counter > 0.0) { // count of when BOMBPROOF switches to lower accuracy
264 const double reduced_accuracy_percentage = 100.0 * (m_stats.reduced_accuracy_counter / (m_grid->Mx() * m_grid->My()));
265 const double reporting_threshold = 5.0; // only report if above 5%
266
267 if (reduced_accuracy_percentage > reporting_threshold and m_log->get_threshold() > 2) {
268 m_stdout_flags = (pism::printf(" [BPsacr=%.4f%%] ", reduced_accuracy_percentage) +
270 }
271 }
272
273 if (m_stats.bulge_counter > 0) {
274 // count of when advection bulges are limited; frequently it is identically zero
277 }
278 }
279}
280
282 // silence a compiler warning
283 (void) t;
284
285 if (m_stress_balance == NULL) {
287 "EnergyModel: no stress balance provided."
288 " Cannot compute max. time step.");
289 }
290
291 return MaxTimestep(m_stress_balance->max_timestep_cfl_3d().dt_max.value(), "energy");
292}
293
294const std::string& EnergyModel::stdout_flags() const {
295 return m_stdout_flags;
296}
297
299 return m_stats;
300}
301
303 return m_ice_enthalpy;
304}
305
306/*! @brief Basal melt rate in grounded areas. (It is set to zero elsewhere.) */
310
311/*! @brief Ice loss "flux" due to ice liquefaction. */
312class LiquifiedIceFlux : public TSDiag<TSFluxDiagnostic,EnergyModel> {
313public:
315 : TSDiag<TSFluxDiagnostic, EnergyModel>(m, "liquified_ice_flux") {
316
317 set_units("m^3 / second", "m^3 / year");
318 m_variable["long_name"] =
319 "rate of ice loss due to liquefaction, averaged over the reporting interval";
320 m_variable["comment"] = "positive means ice loss";
321 m_variable["cell_methods"] = "time: mean";
322 }
323protected:
324 double compute() {
325 // liquified ice volume during the last time step
326 return model->stats().liquified_ice_volume;
327 }
328};
329
331 DiagnosticList result;
332 result = {
333 {"enthalpy", Diagnostic::wrap(m_ice_enthalpy)},
334 {"basal_melt_rate_grounded", Diagnostic::wrap(m_basal_melt_rate)}
335 };
336 return result;
337}
338
340 return {
341 {"liquified_ice_flux", TSDiagnostic::Ptr(new LiquifiedIceFlux(this))}
342 };
343}
344
345} // end of namespace energy
346
347} // end of namespace pism
const Config::ConstPtr m_config
configuration database used by this component
Definition Component.hh:158
const Logger::ConstPtr m_log
logger (for easy access)
Definition Component.hh:162
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition Component.hh:156
const Profiling & profiling() const
Definition Component.cc:113
A class defining a common interface for most PISM sub-models.
Definition Component.hh:118
static Ptr wrap(const T &input)
bool variable_exists(const std::string &short_name) const
Checks if a variable exists.
Definition File.cc:378
std::string name() const
Definition File.cc:274
High-level PISM I/O class.
Definition File.hh:55
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
void begin(const char *name) const
Definition Profiling.cc:75
void end(const char *name) const
Definition Profiling.cc:91
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
Spatial NetCDF variable (corresponding to a 2D or 3D scalar field).
void set_units(const std::string &units, const std::string &output_units)
VariableMetadata m_variable
std::shared_ptr< TSDiagnostic > Ptr
Scalar diagnostic reporting a "flux".
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
VariableMetadata & set_name(const std::string &name)
VariableMetadata & standard_name(const std::string &input)
std::string get_name() const
stencils::Box< T > box(int i, int j) const
Definition Array2D.hh:93
void copy_from(const Array3D &input)
Definition Array3D.cc:209
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition Array3D.hh:33
void read(const std::string &filename, unsigned int time)
Definition Array.cc:731
void set_name(const std::string &name)
Sets the variable name to name.
Definition Array.cc:342
void regrid(const std::string &filename, io::Default default_value)
Definition Array.cc:736
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition Array.cc:476
Abstract class for reading, writing, allocating, and accessing a DA-based PETSc Vec (2D and 3D fields...
Definition Array.hh:207
EnergyModelStats & operator+=(const EnergyModelStats &other)
void init_enthalpy(const File &input_file, bool regrid, int record)
Initialize enthalpy by reading it from a file, or by reading temperature and liquid water fraction,...
const EnergyModelStats & stats() const
void update(double t, double dt, const Inputs &inputs)
const array::Array3D & enthalpy() const
void bootstrap(const File &input_file, const array::Scalar &ice_thickness, const array::Scalar &surface_temperature, const array::Scalar &climatic_mass_balance, const array::Scalar &basal_heat_flux)
Bootstrapping using heuristics.
void initialize(const array::Scalar &basal_melt_rate, const array::Scalar &ice_thickness, const array::Scalar &surface_temperature, const array::Scalar &climatic_mass_balance, const array::Scalar &basal_heat_flux)
Initialize using formulas (for runs using synthetic data).
virtual MaxTimestep max_timestep_impl(double t) const
EnergyModelStats m_stats
std::shared_ptr< const stressbalance::StressBalance > m_stress_balance
const std::string & stdout_flags() const
virtual void initialize_impl(const array::Scalar &basal_melt_rate, const array::Scalar &ice_thickness, const array::Scalar &surface_temperature, const array::Scalar &climatic_mass_balance, const array::Scalar &basal_heat_flux)=0
void restart(const File &input_file, int record)
virtual void restart_impl(const File &input_file, int record)=0
virtual DiagnosticList diagnostics_impl() const
const array::Scalar & basal_melt_rate() const
Basal melt rate in grounded areas. (It is set to zero elsewhere.)
array::Array3D m_ice_enthalpy
virtual void bootstrap_impl(const File &input_file, const array::Scalar &ice_thickness, const array::Scalar &surface_temperature, const array::Scalar &climatic_mass_balance, const array::Scalar &basal_heat_flux)=0
EnergyModel(std::shared_ptr< const Grid > grid, std::shared_ptr< const stressbalance::StressBalance > stress_balance)
virtual void update_impl(double t, double dt, const Inputs &inputs)=0
virtual TSDiagnosticList ts_diagnostics_impl() const
array::Scalar m_basal_melt_rate
void regrid_enthalpy()
Regrid enthalpy from the -regrid_file.
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::Scalar * no_model_mask
const array::Array3D * u3
const array::Scalar * surface_temp
const array::CellType * cell_type
const array::Array3D * volumetric_heating_rate
LiquifiedIceFlux(const EnergyModel *m)
Ice loss "flux" due to ice liquefaction.
static Default Nil()
Definition IO_Flags.hh:93
#define PISM_ERROR_LOCATION
static void check_input(const array::Array *ptr, const char *name)
bool marginal(const array::Scalar1 &thickness, int i, int j, double threshold)
void compute_enthalpy_cold(const array::Array3D &temperature, const array::Scalar &ice_thickness, array::Array3D &result)
Compute ice enthalpy from temperature temperature by assuming the ice has zero liquid fraction.
Definition utilities.cc:48
void compute_enthalpy(const array::Array3D &temperature, const array::Array3D &liquid_water_fraction, const array::Scalar &ice_thickness, array::Array3D &result)
Compute result (enthalpy) from temperature and liquid fraction.
Definition utilities.cc:105
@ PISM_GUESS
Definition IO_Flags.hh:56
@ PISM_READONLY
open an existing file for reading only
Definition IO_Flags.hh:68
std::map< std::string, TSDiagnostic::Ptr > TSDiagnosticList
std::string printf(const char *format,...)
std::map< std::string, Diagnostic::Ptr > DiagnosticList
std::set< std::string > set_split(const std::string &input, char separator)
Transform a separator-separated list (a string) into a set of strings.
bool member(const std::string &string, const std::set< std::string > &set)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)