tEnergyModel.cc - 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
---
tEnergyModel.cc (13772B)
---
1 /* Copyright (C) 2016, 2017, 2018, 2019 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 "EnergyModel.hh"
21 #include "pism/util/MaxTimestep.hh"
22 #include "pism/stressbalance/StressBalance.hh"
23 #include "pism/util/io/File.hh"
24 #include "pism/util/Vars.hh"
25 #include "utilities.hh"
26 #include "pism/util/EnthalpyConverter.hh"
27 #include "bootstrapping.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/util/error_handling.hh"
30 #include "pism/util/IceModelVec2CellType.hh"
31 #include "pism/util/pism_options.hh"
32 #include "pism/util/Profiling.hh"
33
34 namespace pism {
35 namespace energy {
36
37 static void check_input(const IceModelVec *ptr, const char *name) {
38 if (ptr == NULL) {
39 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "energy balance model input %s was not provided", name);
40 }
41 }
42
43 Inputs::Inputs() {
44 basal_frictional_heating = NULL;
45 basal_heat_flux = NULL;
46 cell_type = NULL;
47 ice_thickness = NULL;
48 surface_liquid_fraction = NULL;
49 shelf_base_temp = NULL;
50 surface_temp = NULL;
51 till_water_thickness = NULL;
52
53 volumetric_heating_rate = NULL;
54 u3 = NULL;
55 v3 = NULL;
56 w3 = NULL;
57
58 no_model_mask = NULL;
59 }
60
61 void Inputs::check() const {
62 check_input(cell_type, "cell_type");
63 check_input(basal_frictional_heating, "basal_frictional_heating");
64 check_input(basal_heat_flux, "basal_heat_flux");
65 check_input(ice_thickness, "ice_thickness");
66 check_input(surface_liquid_fraction, "surface_liquid_fraction");
67 check_input(shelf_base_temp, "shelf_base_temp");
68 check_input(surface_temp, "surface_temp");
69 check_input(till_water_thickness, "till_water_thickness");
70
71 check_input(volumetric_heating_rate, "volumetric_heating_rate");
72 check_input(u3, "u3");
73 check_input(v3, "v3");
74 check_input(w3, "w3");
75 }
76
77 EnergyModelStats::EnergyModelStats() {
78 bulge_counter = 0;
79 reduced_accuracy_counter = 0;
80 low_temperature_counter = 0;
81 liquified_ice_volume = 0.0;
82 }
83
84 EnergyModelStats& EnergyModelStats::operator+=(const EnergyModelStats &other) {
85 bulge_counter += other.bulge_counter;
86 reduced_accuracy_counter += other.reduced_accuracy_counter;
87 low_temperature_counter += other.low_temperature_counter;
88 liquified_ice_volume += other.liquified_ice_volume;
89 return *this;
90 }
91
92
93 bool marginal(const IceModelVec2S &thickness, int i, int j, double threshold) {
94 int
95 n = j + 1,
96 e = i + 1,
97 s = j - 1,
98 w = i - 1;
99
100 const double
101 N = thickness(i, n),
102 E = thickness(e, j),
103 S = thickness(i, s),
104 W = thickness(w, j),
105 NW = thickness(w, n),
106 SW = thickness(w, s),
107 NE = thickness(e, n),
108 SE = thickness(e, s);
109
110 return ((E < threshold) or
111 (NE < threshold) or
112 (N < threshold) or
113 (NW < threshold) or
114 (W < threshold) or
115 (SW < threshold) or
116 (S < threshold) or
117 (SE < threshold));
118 }
119
120
121 void EnergyModelStats::sum(MPI_Comm com) {
122 bulge_counter = GlobalSum(com, bulge_counter);
123 reduced_accuracy_counter = GlobalSum(com, reduced_accuracy_counter);
124 low_temperature_counter = GlobalSum(com, low_temperature_counter);
125 liquified_ice_volume = GlobalSum(com, liquified_ice_volume);
126 }
127
128
129 EnergyModel::EnergyModel(IceGrid::ConstPtr grid,
130 stressbalance::StressBalance *stress_balance)
131 : Component(grid), m_stress_balance(stress_balance) {
132
133 const unsigned int WIDE_STENCIL = m_config->get_number("grid.max_stencil_width");
134
135 {
136 m_ice_enthalpy.create(m_grid, "enthalpy", WITH_GHOSTS, WIDE_STENCIL);
137 // POSSIBLE standard name = land_ice_enthalpy
138 m_ice_enthalpy.set_attrs("model_state",
139 "ice enthalpy (includes sensible heat, latent heat, pressure)",
140 "J kg-1", "J kg-1", "", 0);
141 }
142
143 {
144 m_basal_melt_rate.create(m_grid, "basal_melt_rate_grounded", WITHOUT_GHOSTS);
145 // ghosted to allow the "redundant" computation of tauc
146 m_basal_melt_rate.set_attrs("model_state",
147 "ice basal melt rate from energy conservation, in ice thickness per time (valid in grounded areas)",
148 "m s-1", "m year-1", "", 0);
149 // We could use land_ice_basal_melt_rate, but that way both basal_melt_rate_grounded and bmelt
150 // have this standard name.
151 m_basal_melt_rate.metadata().set_string("comment", "positive basal melt rate corresponds to ice loss");
152 }
153
154 // a 3d work vector
155 {
156 m_work.create(m_grid, "work_vector", WITHOUT_GHOSTS);
157 m_work.set_attrs("internal",
158 "usually new values of temperature or enthalpy during time step",
159 "", "", "", 0);
160 }
161 }
162
163 void EnergyModel::init_enthalpy(const File &input_file, bool do_regrid, int record) {
164
165 if (input_file.find_variable("enthalpy")) {
166 if (do_regrid) {
167 m_ice_enthalpy.regrid(input_file, CRITICAL);
168 } else {
169 m_ice_enthalpy.read(input_file, record);
170 }
171 } else if (input_file.find_variable("temp")) {
172 IceModelVec3
173 &temp = m_work,
174 &liqfrac = m_ice_enthalpy;
175
176 {
177 temp.set_name("temp");
178 temp.metadata(0).set_name("temp");
179 temp.set_attrs("temporary", "ice temperature",
180 "Kelvin", "Kelvin", "land_ice_temperature", 0);
181
182 if (do_regrid) {
183 temp.regrid(input_file, CRITICAL);
184 } else {
185 temp.read(input_file, record);
186 }
187 }
188
189 const IceModelVec2S & ice_thickness = *m_grid->variables().get_2d_scalar("land_ice_thickness");
190
191 if (input_file.find_variable("liqfrac")) {
192 SpatialVariableMetadata enthalpy_metadata = m_ice_enthalpy.metadata();
193
194 liqfrac.set_name("liqfrac");
195 liqfrac.metadata(0).set_name("liqfrac");
196 liqfrac.set_attrs("temporary", "ice liquid water fraction",
197 "1", "1", "", 0);
198
199 if (do_regrid) {
200 liqfrac.regrid(input_file, CRITICAL);
201 } else {
202 liqfrac.read(input_file, record);
203 }
204
205 m_ice_enthalpy.set_name(enthalpy_metadata.get_name());
206 m_ice_enthalpy.metadata() = enthalpy_metadata;
207
208 m_log->message(2,
209 " - Computing enthalpy using ice temperature,"
210 " liquid water fraction and thickness...\n");
211 compute_enthalpy(temp, liqfrac, ice_thickness, m_ice_enthalpy);
212 } else {
213 m_log->message(2,
214 " - Computing enthalpy using ice temperature and thickness "
215 "and assuming zero liquid water fraction...\n");
216 compute_enthalpy_cold(temp, ice_thickness, m_ice_enthalpy);
217 }
218 } else {
219 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
220 "neither enthalpy nor temperature was found in '%s'.\n",
221 input_file.filename().c_str());
222 }
223 }
224
225 /*!
226 * The `-regrid_file` may contain enthalpy, temperature, or *both* temperature and liquid water
227 * fraction.
228 */
229 void EnergyModel::regrid_enthalpy() {
230
231 auto regrid_filename = m_config->get_string("input.regrid.file");
232 auto regrid_vars = set_split(m_config->get_string("input.regrid.vars"), ',');
233
234
235 if (regrid_filename.empty()) {
236 return;
237 }
238
239 std::string enthalpy_name = m_ice_enthalpy.metadata().get_name();
240
241 if (regrid_vars.empty() or member(enthalpy_name, regrid_vars)) {
242 File regrid_file(m_grid->com, regrid_filename, PISM_GUESS, PISM_READONLY);
243 init_enthalpy(regrid_file, true, 0);
244 }
245 }
246
247
248 void EnergyModel::restart(const File &input_file, int record) {
249 this->restart_impl(input_file, record);
250 }
251
252 void EnergyModel::bootstrap(const File &input_file,
253 const IceModelVec2S &ice_thickness,
254 const IceModelVec2S &surface_temperature,
255 const IceModelVec2S &climatic_mass_balance,
256 const IceModelVec2S &basal_heat_flux) {
257 this->bootstrap_impl(input_file,
258 ice_thickness, surface_temperature,
259 climatic_mass_balance, basal_heat_flux);
260 }
261
262 void EnergyModel::initialize(const IceModelVec2S &basal_melt_rate,
263 const IceModelVec2S &ice_thickness,
264 const IceModelVec2S &surface_temperature,
265 const IceModelVec2S &climatic_mass_balance,
266 const IceModelVec2S &basal_heat_flux) {
267 this->initialize_impl(basal_melt_rate,
268 ice_thickness,
269 surface_temperature,
270 climatic_mass_balance,
271 basal_heat_flux);
272 }
273
274 void EnergyModel::update(double t, double dt, const Inputs &inputs) {
275 // reset standard out flags at the beginning of every time step
276 m_stdout_flags = "";
277 m_stats = EnergyModelStats();
278
279 const Profiling &profiling = m_grid->ctx()->profiling();
280
281 profiling.begin("ice_energy");
282 {
283 // this call should fill m_work with new values of enthalpy
284 this->update_impl(t, dt, inputs);
285
286 m_work.update_ghosts(m_ice_enthalpy);
287 }
288 profiling.end("ice_energy");
289
290 // globalize m_stats and update m_stdout_flags
291 {
292 char buffer[50] = "";
293
294 m_stats.sum(m_grid->com);
295
296 if (m_stats.reduced_accuracy_counter > 0.0) { // count of when BOMBPROOF switches to lower accuracy
297 const double reduced_accuracy_percentage = 100.0 * (m_stats.reduced_accuracy_counter / (m_grid->Mx() * m_grid->My()));
298 const double reporting_threshold = 5.0; // only report if above 5%
299
300 if (reduced_accuracy_percentage > reporting_threshold and m_log->get_threshold() > 2) {
301 snprintf(buffer, 50, " [BPsacr=%.4f%%] ", reduced_accuracy_percentage);
302 m_stdout_flags = buffer + m_stdout_flags;
303 }
304 }
305
306 if (m_stats.bulge_counter > 0) {
307 // count of when advection bulges are limited; frequently it is identically zero
308 snprintf(buffer, 50, " BULGE=%d ", m_stats.bulge_counter);
309 m_stdout_flags = buffer + m_stdout_flags;
310 }
311 }
312 }
313
314 MaxTimestep EnergyModel::max_timestep_impl(double t) const {
315 // silence a compiler warning
316 (void) t;
317
318 if (m_stress_balance == NULL) {
319 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
320 "EnergyModel: no stress balance provided."
321 " Cannot compute max. time step.");
322 }
323
324 return MaxTimestep(m_stress_balance->max_timestep_cfl_3d().dt_max.value(), "energy");
325 }
326
327 const std::string& EnergyModel::stdout_flags() const {
328 return m_stdout_flags;
329 }
330
331 const EnergyModelStats& EnergyModel::stats() const {
332 return m_stats;
333 }
334
335 const IceModelVec3 & EnergyModel::enthalpy() const {
336 return m_ice_enthalpy;
337 }
338
339 /*! @brief Basal melt rate in grounded areas. (It is set to zero elsewhere.) */
340 const IceModelVec2S & EnergyModel::basal_melt_rate() const {
341 return m_basal_melt_rate;
342 }
343
344 /*! @brief Ice loss "flux" due to ice liquefaction. */
345 class LiquifiedIceFlux : public TSDiag<TSFluxDiagnostic,EnergyModel> {
346 public:
347 LiquifiedIceFlux(const EnergyModel *m)
348 : TSDiag<TSFluxDiagnostic, EnergyModel>(m, "liquified_ice_flux") {
349
350 set_units("m3 / second", "m3 / year");
351 m_ts.variable().set_string("long_name",
352 "rate of ice loss due to liquefaction,"
353 " averaged over the reporting interval");
354 m_ts.variable().set_string("comment", "positive means ice loss");
355 m_ts.variable().set_string("cell_methods", "time: mean");
356 }
357 protected:
358 double compute() {
359 // liquified ice volume during the last time step
360 return model->stats().liquified_ice_volume;
361 }
362 };
363
364 namespace diagnostics {
365 /*! @brief Report ice enthalpy. */
366 class Enthalpy : public Diag<EnergyModel>
367 {
368 public:
369 Enthalpy(const EnergyModel *m)
370 : Diag<EnergyModel>(m) {
371 m_vars = {model->enthalpy().metadata()};
372 }
373
374 protected:
375 IceModelVec::Ptr compute_impl() const {
376
377 IceModelVec3::Ptr result(new IceModelVec3(m_grid, "enthalpy", WITHOUT_GHOSTS));
378 result->metadata(0) = m_vars[0];
379
380 const IceModelVec3 &input = model->enthalpy();
381
382 // FIXME: implement IceModelVec3::copy_from()
383
384 IceModelVec::AccessList list {result.get(), &input};
385 ParallelSection loop(m_grid->com);
386 try {
387 for (Points p(*m_grid); p; p.next()) {
388 const int i = p.i(), j = p.j();
389
390 result->set_column(i, j, input.get_column(i, j));
391 }
392 } catch (...) {
393 loop.failed();
394 }
395 loop.check();
396
397
398 return result;
399 }
400 };
401
402 } // end of namespace diagnostics
403
404 DiagnosticList EnergyModel::diagnostics_impl() const {
405 DiagnosticList result;
406 result = {
407 {"enthalpy", Diagnostic::Ptr(new diagnostics::Enthalpy(this))},
408 {"basal_melt_rate_grounded", Diagnostic::wrap(m_basal_melt_rate)}
409 };
410 return result;
411 }
412
413 TSDiagnosticList EnergyModel::ts_diagnostics_impl() const {
414 return {
415 {"liquified_ice_flux", TSDiagnostic::Ptr(new LiquifiedIceFlux(this))}
416 };
417 }
418
419 } // end of namespace energy
420
421 } // end of namespace pism