tTemperatureModel_Verification.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
---
tTemperatureModel_Verification.cc (4241B)
---
1 /* Copyright (C) 2016, 2017 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 "TemperatureModel_Verification.hh"
21
22 #include "pism/util/error_handling.hh"
23 #include "pism/verification/tests/exactTestsFG.hh"
24 #include "pism/verification/tests/exactTestK.h"
25 #include "pism/verification/tests/exactTestO.h"
26 #include "pism/energy/utilities.hh"
27 #include "pism/util/Time.hh"
28
29 namespace pism {
30 namespace energy {
31
32 static const double ApforG = 200; // m
33 static const double LforFG = 750000; // m
34 static const double ST = 1.67e-5;
35 static const double Tmin = 223.15; // K
36
37 TemperatureModel_Verification::TemperatureModel_Verification(IceGrid::ConstPtr grid,
38 stressbalance::StressBalance *stress_balance,
39 int testname, bool bedrock_is_ice)
40 : TemperatureModel(grid, stress_balance), m_testname(testname), m_bedrock_is_ice(bedrock_is_ice) {
41 // empty
42 }
43
44 void TemperatureModel_Verification::initialize_impl(const IceModelVec2S &basal_melt_rate,
45 const IceModelVec2S &ice_thickness,
46 const IceModelVec2S &surface_temperature,
47 const IceModelVec2S &climatic_mass_balance,
48 const IceModelVec2S &basal_heat_flux) {
49
50 // ignore provided basal melt rate
51 (void) basal_melt_rate;
52
53 m_basal_melt_rate.set(0.0);
54
55 switch (m_testname) {
56 case 'F':
57 case 'G':
58 initTestFG();
59 break;
60 case 'K':
61 case 'O':
62 initTestsKO();
63 break;
64 default:
65 TemperatureModel::initialize_impl(m_basal_melt_rate, ice_thickness, surface_temperature,
66 climatic_mass_balance, basal_heat_flux);
67 }
68
69 m_ice_temperature.update_ghosts();
70 m_basal_melt_rate.update_ghosts();
71
72 // this will update ghosts of m_ice_enthalpy
73 compute_enthalpy_cold(m_ice_temperature, ice_thickness, m_ice_enthalpy);
74 }
75
76 void TemperatureModel_Verification::initTestFG() {
77
78 IceModelVec::AccessList list{&m_ice_temperature};
79
80 const double time = m_testname == 'F' ? 0.0 : m_grid->ctx()->time()->current();
81 const double A = m_testname == 'F' ? 0.0 : ApforG;
82
83 for (Points p(*m_grid); p; p.next()) {
84 const int i = p.i(), j = p.j();
85
86 // avoid singularity at origin
87 const double r = std::max(radius(*m_grid, i, j), 1.0);
88
89 if (r > LforFG - 1.0) { // if (essentially) outside of sheet
90 m_ice_temperature.set_column(i, j, Tmin + ST * r);
91 } else {
92 TestFGParameters P = exactFG(time, r, m_grid->z(), A);
93 m_ice_temperature.set_column(i, j, &P.T[0]);
94 }
95 }
96 }
97
98 void TemperatureModel_Verification::initTestsKO() {
99
100 const unsigned int Mz = m_grid->Mz();
101
102 std::vector<double> T_column(Mz);
103
104 const double time = m_grid->ctx()->time()->current();
105
106 // evaluate exact solution in a column; all columns are the same
107 for (unsigned int k = 0; k < Mz; k++) {
108 if (m_testname == 'K') {
109 T_column[k] = exactK(time, m_grid->z(k), m_bedrock_is_ice).T;
110 } else {
111 T_column[k] = exactO(m_grid->z(k)).TT;
112 }
113 }
114
115 // fill m_ice_temperature
116 IceModelVec::AccessList list(m_ice_temperature);
117
118 ParallelSection loop(m_grid->com);
119 try {
120 for (Points p(*m_grid); p; p.next()) {
121 m_ice_temperature.set_column(p.i(), p.j(), &T_column[0]);
122 }
123 } catch (...) {
124 loop.failed();
125 }
126 loop.check();
127 }
128
129
130 } // end of namespace energy
131 } // end of namespace pism