tiCMthermo.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
---
tiCMthermo.cc (18118B)
---
1 // Copyright (C) 2004-2018 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
21 #include "tests/exactTestsFG.hh"
22 #include "tests/exactTestK.h"
23 #include "tests/exactTestO.h"
24 #include "iceCompModel.hh"
25
26 #include "pism/stressbalance/StressBalance.hh"
27 #include "pism/util/Time.hh"
28 #include "pism/util/IceGrid.hh"
29 #include "pism/util/pism_options.hh"
30 #include "pism/util/error_handling.hh"
31 #include "pism/earth/BedDef.hh"
32 #include "pism/util/ConfigInterface.hh"
33 #include "pism/util/pism_utilities.hh"
34 #include "BTU_Verification.hh"
35 #include "pism/energy/TemperatureModel.hh"
36 #include "pism/coupler/SurfaceModel.hh"
37 #include "pism/coupler/OceanModel.hh"
38 #include "pism/hydrology/Hydrology.hh"
39
40 namespace pism {
41
42 // boundary conditions for tests F, G (same as EISMINT II Experiment F)
43 const double IceCompModel::m_ST = 1.67e-5;
44 const double IceCompModel::m_Tmin = 223.15; // K
45 const double IceCompModel::m_LforFG = 750000; // m
46 const double IceCompModel::m_ApforG = 200; // m
47
48 void IceCompModel::energy_step() {
49
50 energy::EnergyModelStats stats;
51
52 IceModelVec2S &bedtoptemp = m_work2d[1];
53 IceModelVec2S &basal_enthalpy = m_work2d[2];
54 m_energy_model->enthalpy().getHorSlice(basal_enthalpy, 0.0);
55
56 bedrock_surface_temperature(m_geometry.sea_level_elevation,
57 m_geometry.cell_type,
58 m_geometry.bed_elevation,
59 m_geometry.ice_thickness,
60 basal_enthalpy,
61 m_surface->temperature(),
62 bedtoptemp);
63
64 m_btu->update(bedtoptemp, t_TempAge, dt_TempAge);
65
66 energy::Inputs inputs;
67 {
68 inputs.basal_frictional_heating = &m_stress_balance->basal_frictional_heating();
69 inputs.basal_heat_flux = &m_btu->flux_through_top_surface(); // bedrock thermal layer
70 inputs.cell_type = &m_geometry.cell_type; // geometry
71 inputs.ice_thickness = &m_geometry.ice_thickness; // geometry
72 inputs.shelf_base_temp = &m_ocean->shelf_base_temperature(); // ocean model
73 inputs.surface_liquid_fraction = &m_surface->liquid_water_fraction(); // surface model
74 inputs.surface_temp = &m_surface->temperature(); // surface model
75 inputs.till_water_thickness = &m_subglacial_hydrology->till_water_thickness();
76
77 inputs.volumetric_heating_rate = &m_stress_balance->volumetric_strain_heating();
78 inputs.u3 = &m_stress_balance->velocity_u();
79 inputs.v3 = &m_stress_balance->velocity_v();
80 inputs.w3 = &m_stress_balance->velocity_w();
81
82 inputs.check(); // make sure all data members were set
83 }
84
85 if ((m_testname == 'F') || (m_testname == 'G')) {
86 // Compute compensatory strain heating (fills strain_heating3_comp).
87 getCompSourcesTestFG();
88
89 // Add computed strain heating to the compensatory part.
90 m_strain_heating3_comp.add(1.0, *inputs.volumetric_heating_rate);
91
92 // Use the result.
93 inputs.volumetric_heating_rate = &m_strain_heating3_comp;
94 }
95
96 m_energy_model->update(t_TempAge, dt_TempAge, inputs);
97
98 m_stdout_flags = m_energy_model->stdout_flags() + m_stdout_flags;
99 }
100
101 void IceCompModel::initTestFG() {
102
103 IceModelVec::AccessList list{&m_geometry.ice_thickness};
104
105 const double time = m_testname == 'F' ? 0.0 : m_time->current();
106 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
107
108 for (Points p(*m_grid); p; p.next()) {
109 const int i = p.i(), j = p.j();
110
111 // avoid singularity at origin
112 const double r = std::max(radius(*m_grid, i, j), 1.0);
113
114 if (r > m_LforFG - 1.0) { // if (essentially) outside of sheet
115 m_geometry.ice_thickness(i, j) = 0.0;
116 } else {
117 m_geometry.ice_thickness(i, j) = exactFG(time, r, m_grid->z(), A).H;
118 }
119 }
120
121 m_geometry.ice_thickness.update_ghosts();
122
123 {
124 IceModelVec2S bed_topography(m_grid, "topg", WITHOUT_GHOSTS);
125 bed_topography.set(0.0);
126
127 IceModelVec2S bed_uplift(m_grid, "uplift", WITHOUT_GHOSTS);
128 bed_uplift.set(0.0);
129
130 IceModelVec2S sea_level(m_grid, "sea_level", WITHOUT_GHOSTS);
131 sea_level.set(0.0);
132
133 m_beddef->bootstrap(bed_topography, bed_uplift, m_geometry.ice_thickness,
134 sea_level);
135 }
136 }
137
138 void IceCompModel::initTestsKO() {
139
140 IceModelVec2S bed_topography(m_grid, "topg", WITHOUT_GHOSTS);
141 bed_topography.set(0.0);
142
143 IceModelVec2S bed_uplift(m_grid, "uplift", WITHOUT_GHOSTS);
144 bed_uplift.set(0.0);
145
146 IceModelVec2S sea_level(m_grid, "sea_level", WITHOUT_GHOSTS);
147 sea_level.set(0.0);
148
149 m_geometry.ice_thickness.set(3000.0);
150
151 m_beddef->bootstrap(bed_topography, bed_uplift, m_geometry.ice_thickness,
152 sea_level);
153 }
154
155 void IceCompModel::getCompSourcesTestFG() {
156
157 const double
158 ice_rho = m_config->get_number("constants.ice.density"),
159 ice_c = m_config->get_number("constants.ice.specific_heat_capacity");
160
161 // before temperature and flow step, set strain_heating_c from exact values
162
163 IceModelVec::AccessList list{&m_strain_heating3_comp};
164
165 const double time = m_testname == 'F' ? 0.0 : m_time->current();
166 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
167
168 for (Points p(*m_grid); p; p.next()) {
169 const int i = p.i(), j = p.j();
170
171 double r = std::max(radius(*m_grid, i, j), 1.0); // avoid singularity at origin
172
173 if (r > m_LforFG - 1.0) { // outside of sheet
174 m_strain_heating3_comp.set_column(i, j, 0.0);
175 } else {
176 TestFGParameters P = exactFG(time, r, m_grid->z(), A);
177
178 m_strain_heating3_comp.set_column(i, j, &P.Sigc[0]);
179 }
180 }
181
182 // scale strain_heating to J/(s m^3)
183 m_strain_heating3_comp.scale(ice_rho * ice_c);
184 }
185
186 void IceCompModel::computeTemperatureErrors(double &gmaxTerr,
187 double &gavTerr) {
188 double maxTerr = 0.0, avTerr = 0.0, avcount = 0.0;
189
190 if (m_testname != 'F' and m_testname != 'G') {
191 throw RuntimeError(PISM_ERROR_LOCATION, "temperature errors only computable for tests F and G");
192 }
193
194 const double time = m_testname == 'F' ? 0.0 : m_time->current();
195 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
196
197 energy::TemperatureModel *m = dynamic_cast<energy::TemperatureModel*>(m_energy_model);
198 const IceModelVec3 &ice_temperature = m->temperature();
199
200 IceModelVec::AccessList list{&m_geometry.ice_thickness, &ice_temperature};
201
202 ParallelSection loop(m_grid->com);
203 try {
204 for (Points p(*m_grid); p; p.next()) {
205 const int i = p.i(), j = p.j();
206
207 const double r = radius(*m_grid, i, j);
208 const double *T = ice_temperature.get_column(i, j);
209
210 // only evaluate error if inside sheet and not at central
211 // singularity
212 if ((r >= 1.0) and (r <= m_LforFG - 1.0)) {
213 TestFGParameters P = exactFG(time, r, m_grid->z(), A);
214
215 // only evaluate error if below ice surface
216 const int ks = m_grid->kBelowHeight(m_geometry.ice_thickness(i, j));
217 for (int k = 0; k < ks; k++) {
218 const double Terr = fabs(T[k] - P.T[k]);
219 maxTerr = std::max(maxTerr, Terr);
220 avcount += 1.0;
221 avTerr += Terr;
222 }
223 }
224 }
225 } catch (...) {
226 loop.failed();
227 }
228 loop.check();
229
230 gmaxTerr = GlobalMax(m_grid->com, maxTerr);
231 gavTerr = GlobalSum(m_grid->com, avTerr);
232 double gavcount = GlobalSum(m_grid->com, avcount);
233 gavTerr = gavTerr / std::max(gavcount, 1.0); // avoid division by zero
234 }
235
236
237 void IceCompModel::computeIceBedrockTemperatureErrors(double &gmaxTerr, double &gavTerr,
238 double &gmaxTberr, double &gavTberr) {
239
240 if (m_testname != 'K' and m_testname != 'O') {
241 throw RuntimeError(PISM_ERROR_LOCATION, "ice and bedrock temperature errors only computable for tests K and O");
242 }
243
244 double maxTerr = 0.0, avTerr = 0.0, avcount = 0.0;
245 double maxTberr = 0.0, avTberr = 0.0, avbcount = 0.0;
246
247 const double *Tb, *T;
248 std::vector<double> Tex(m_grid->Mz());
249
250 energy::BTU_Verification *my_btu = dynamic_cast<energy::BTU_Verification*>(m_btu);
251 if (my_btu == NULL) {
252 throw RuntimeError(PISM_ERROR_LOCATION, "BTU_Verification is required");
253 }
254 const IceModelVec3Custom &bedrock_temp = my_btu->temperature();
255
256 std::vector<double> zblevels = bedrock_temp.levels();
257 unsigned int Mbz = (unsigned int)zblevels.size();
258 std::vector<double> Tbex(Mbz);
259
260 switch (m_testname) {
261 case 'K':
262 for (unsigned int k = 0; k < m_grid->Mz(); k++) {
263 TestKParameters K = exactK(m_time->current(), m_grid->z(k), m_bedrock_is_ice_forK);
264 Tex[k] = K.T;
265 }
266 for (unsigned int k = 0; k < Mbz; k++) {
267 TestKParameters K = exactK(m_time->current(), zblevels[k], m_bedrock_is_ice_forK);
268 Tbex[k] = K.T;
269 }
270 break;
271 case 'O':
272 for (unsigned int k = 0; k < m_grid->Mz(); k++) {
273 Tex[k] = exactO(m_grid->z(k)).TT;
274 }
275 for (unsigned int k = 0; k < Mbz; k++) {
276 Tbex[k] = exactO(zblevels[k]).TT;
277 }
278 break;
279 default:
280 throw RuntimeError(PISM_ERROR_LOCATION, "ice and bedrock temperature errors only for tests K and O");
281 }
282
283 energy::TemperatureModel *m = dynamic_cast<energy::TemperatureModel*>(m_energy_model);
284 const IceModelVec3 &ice_temperature = m->temperature();
285
286 IceModelVec::AccessList list{&ice_temperature, &bedrock_temp};
287
288 for (Points p(*m_grid); p; p.next()) {
289 const int i = p.i(), j = p.j();
290
291 Tb = bedrock_temp.get_column(i, j);
292 for (unsigned int kb = 0; kb < Mbz; kb++) {
293 const double Tberr = fabs(Tb[kb] - Tbex[kb]);
294 maxTberr = std::max(maxTberr, Tberr);
295 avbcount += 1.0;
296 avTberr += Tberr;
297 }
298 T = ice_temperature.get_column(i, j);
299 for (unsigned int k = 0; k < m_grid->Mz(); k++) {
300 const double Terr = fabs(T[k] - Tex[k]);
301 maxTerr = std::max(maxTerr, Terr);
302 avcount += 1.0;
303 avTerr += Terr;
304 }
305 }
306
307 gmaxTerr = GlobalMax(m_grid->com, maxTerr);
308 gavTerr = GlobalSum(m_grid->com, avTerr);
309 double gavcount = GlobalSum(m_grid->com, avcount);
310 gavTerr = gavTerr/std::max(gavcount, 1.0); // avoid division by zero
311
312 gmaxTberr = GlobalMax(m_grid->com, maxTberr);
313 gavTberr = GlobalSum(m_grid->com, avTberr);
314 double gavbcount = GlobalSum(m_grid->com, avbcount);
315 gavTberr = gavTberr/std::max(gavbcount, 1.0); // avoid division by zero
316 }
317
318
319 void IceCompModel::computeBasalTemperatureErrors(double &gmaxTerr, double &gavTerr, double ¢erTerr) {
320
321 if (m_testname != 'F' and m_testname != 'G') {
322 throw RuntimeError(PISM_ERROR_LOCATION, "temperature errors only computable for tests F and G");
323 }
324
325 double
326 Texact = 0.0,
327 domeT = 0.0,
328 domeTexact = 0.0,
329 Terr = 0.0,
330 avTerr = 0.0;
331
332 const double time = m_testname == 'F' ? 0.0 : m_time->current();
333 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
334 std::vector<double> z(1, 0.0);
335
336 energy::TemperatureModel *m = dynamic_cast<energy::TemperatureModel*>(m_energy_model);
337 const IceModelVec3 &ice_temperature = m->temperature();
338
339 IceModelVec::AccessList list(ice_temperature);
340
341 ParallelSection loop(m_grid->com);
342 try {
343 for (Points p(*m_grid); p; p.next()) {
344 const int i = p.i(), j = p.j();
345
346 double r = std::max(radius(*m_grid, i, j), 1.0);
347
348 if (r > m_LforFG - 1.0) { // outside of sheet
349 Texact = m_Tmin + m_ST * r; // = Ts
350 } else {
351 Texact = exactFG(time, r, z, A).T[0];
352 }
353
354 const double Tbase = ice_temperature.get_column(i,j)[0];
355 if (i == ((int)m_grid->Mx() - 1) / 2 and
356 j == ((int)m_grid->My() - 1) / 2) {
357 domeT = Tbase;
358 domeTexact = Texact;
359 }
360 // compute maximum errors
361 Terr = std::max(Terr, fabs(Tbase - Texact));
362 // add to sums for average errors
363 avTerr += fabs(Tbase - Texact);
364 }
365 } catch (...) {
366 loop.failed();
367 }
368 loop.check();
369
370 double gdomeT, gdomeTexact;
371
372 gmaxTerr = GlobalMax(m_grid->com, Terr);
373 gavTerr = GlobalSum(m_grid->com, avTerr);
374 gavTerr = gavTerr/(m_grid->Mx()*m_grid->My());
375 gdomeT = GlobalMax(m_grid->com, domeT);
376 gdomeTexact = GlobalMax(m_grid->com, domeTexact);
377 centerTerr = fabs(gdomeT - gdomeTexact);
378 }
379
380
381 void IceCompModel::compute_strain_heating_errors(double &gmax_strain_heating_err, double &gav_strain_heating_err) {
382 double max_strain_heating_err = 0.0, av_strain_heating_err = 0.0, avcount = 0.0;
383
384 if (m_testname != 'F' and m_testname != 'G') {
385 throw RuntimeError(PISM_ERROR_LOCATION, "strain-heating (strain_heating) errors only computable for tests F and G");
386 }
387
388 const double
389 ice_rho = m_config->get_number("constants.ice.density"),
390 ice_c = m_config->get_number("constants.ice.specific_heat_capacity");
391
392 const IceModelVec3 &strain_heating3 = m_stress_balance->volumetric_strain_heating();
393
394 IceModelVec::AccessList list{&m_geometry.ice_thickness, &strain_heating3};
395
396 const double time = m_testname == 'F' ? 0.0 : m_time->current();
397 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
398
399 ParallelSection loop(m_grid->com);
400 try {
401 for (Points p(*m_grid); p; p.next()) {
402 const int i = p.i(), j = p.j();
403
404 double r = radius(*m_grid, i, j);
405 if ((r >= 1.0) && (r <= m_LforFG - 1.0)) {
406 // only evaluate error if inside sheet and not at central singularity
407
408 TestFGParameters P = exactFG(time, r, m_grid->z(), A);
409
410 for (unsigned int k = 0; k < m_grid->Mz(); k++) {
411 // scale exact strain_heating to J/(s m^3)
412 P.Sig[k] *= ice_rho * ice_c;
413 }
414
415 const unsigned int ks = m_grid->kBelowHeight(m_geometry.ice_thickness(i, j));
416 const double *strain_heating = strain_heating3.get_column(i, j);
417
418 for (unsigned int k = 0; k < ks; k++) { // only evaluate error if below ice surface
419 const double strain_heating_err = fabs(strain_heating[k] - P.Sig[k]);
420 max_strain_heating_err = std::max(max_strain_heating_err, strain_heating_err);
421 avcount += 1.0;
422 av_strain_heating_err += strain_heating_err;
423 }
424 }
425 }
426 } catch (...) {
427 loop.failed();
428 }
429 loop.check();
430
431 gmax_strain_heating_err = GlobalMax(m_grid->com, max_strain_heating_err);
432 gav_strain_heating_err = GlobalSum(m_grid->com, av_strain_heating_err);
433 double gavcount = GlobalSum(m_grid->com, avcount);
434 gav_strain_heating_err = gav_strain_heating_err/std::max(gavcount, 1.0); // avoid div by zero
435 }
436
437
438 void IceCompModel::computeSurfaceVelocityErrors(double &gmaxUerr, double &gavUerr,
439 double &gmaxWerr, double &gavWerr) {
440 double
441 maxUerr = 0.0,
442 maxWerr = 0.0,
443 avUerr = 0.0,
444 avWerr = 0.0;
445
446 const IceModelVec3
447 &u3 = m_stress_balance->velocity_u(),
448 &v3 = m_stress_balance->velocity_v(),
449 &w3 = m_stress_balance->velocity_w();
450
451 IceModelVec::AccessList list{&m_geometry.ice_thickness, &u3, &v3, &w3};
452
453 const double time = m_testname == 'F' ? 0.0 : m_time->current();
454 const double A = m_testname == 'F' ? 0.0 : m_ApforG;
455
456 ParallelSection loop(m_grid->com);
457 try {
458 for (Points p(*m_grid); p; p.next()) {
459 const int i = p.i(), j = p.j();
460
461 const double H = m_geometry.ice_thickness(i, j);
462 std::vector<double> z(1, H);
463
464 const double
465 x = m_grid->x(i),
466 y = m_grid->y(j),
467 r = radius(*m_grid, i, j);
468
469 if ((r >= 1.0) and (r <= m_LforFG - 1.0)) {
470 // only evaluate error if inside sheet and not at central singularity
471
472 TestFGParameters P = exactFG(time, r, z, A);
473
474 const double
475 uex = (x/r) * P.U[0],
476 vex = (y/r) * P.U[0];
477 // note that because getValZ does linear interpolation and H(i, j) is not exactly at
478 // a grid point, this causes nonzero errors
479 const double Uerr = sqrt(PetscSqr(u3.getValZ(i, j, H) - uex) +
480 PetscSqr(v3.getValZ(i, j, H) - vex));
481 maxUerr = std::max(maxUerr, Uerr);
482 avUerr += Uerr;
483 const double Werr = fabs(w3.getValZ(i, j, H) - P.w[0]);
484 maxWerr = std::max(maxWerr, Werr);
485 avWerr += Werr;
486 }
487 }
488 } catch (...) {
489 loop.failed();
490 }
491 loop.check();
492
493 gmaxUerr = GlobalMax(m_grid->com, maxUerr);
494 gmaxWerr = GlobalMax(m_grid->com, maxWerr);
495 gavUerr = GlobalSum(m_grid->com, avUerr);
496 gavUerr = gavUerr/(m_grid->Mx()*m_grid->My());
497 gavWerr = GlobalSum(m_grid->com, avWerr);
498 gavWerr = gavWerr/(m_grid->Mx()*m_grid->My());
499 }
500
501
502 void IceCompModel::computeBasalMeltRateErrors(double &gmaxbmelterr, double &gminbmelterr) {
503 double
504 maxbmelterr = -9.99e40,
505 minbmelterr = 9.99e40;
506
507 if (m_testname != 'O') {
508 throw RuntimeError(PISM_ERROR_LOCATION, "basal melt rate errors are only computable for test O");
509 }
510
511 // we just need one constant from exact solution:
512 double bmelt = exactO(0.0).bmelt;
513
514 const IceModelVec2S &basal_melt_rate = m_energy_model->basal_melt_rate();
515
516 IceModelVec::AccessList list(basal_melt_rate);
517
518 for (Points p(*m_grid); p; p.next()) {
519 const int i = p.i(), j = p.j();
520
521 double err = fabs(basal_melt_rate(i, j) - bmelt);
522 maxbmelterr = std::max(maxbmelterr, err);
523 minbmelterr = std::min(minbmelterr, err);
524 }
525
526 gmaxbmelterr = GlobalMax(m_grid->com, maxbmelterr);
527 gminbmelterr = GlobalMin(m_grid->com, minbmelterr);
528 }
529
530 } // end of namespace pism