tPSVerification.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
---
tPSVerification.cc (6939B)
---
1 /* Copyright (C) 2014, 2015, 2016, 2017, 2018 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 <algorithm>
21 #include <vector>
22
23 #include "PSVerification.hh"
24 #include "pism/coupler/AtmosphereModel.hh"
25 #include "pism/rheology/PatersonBuddCold.hh"
26 #include "pism/util/EnthalpyConverter.hh"
27 #include "pism/util/Time.hh"
28 #include "pism/util/ConfigInterface.hh"
29
30 #include "tests/exactTestsABCD.h"
31 #include "tests/exactTestsFG.hh"
32 #include "tests/exactTestH.h"
33 #include "tests/exactTestL.hh"
34
35 #include "pism/util/error_handling.hh"
36 #include "pism/util/IceGrid.hh"
37 #include "pism/util/MaxTimestep.hh"
38
39 namespace pism {
40 namespace surface {
41
42 const double Verification::ablationRateOutside = 0.02; // m year-1
43 const double Verification::secpera = 3.15569259747e7;
44
45 const double Verification::ST = 1.67e-5;
46 const double Verification::Tmin = 223.15; // K
47 const double Verification::LforFG = 750000; // m
48 const double Verification::ApforG = 200; // m
49
50 Verification::Verification(IceGrid::ConstPtr g,
51 EnthalpyConverter::Ptr EC, int test)
52 : PSFormulas(g), m_testname(test), m_EC(EC) {
53 // empty
54 }
55
56 Verification::~Verification() {
57 // empty
58 }
59
60 void Verification::init_impl(const Geometry &geometry) {
61 // Make sure that ice surface temperature and climatic mass balance
62 // get initialized at the beginning of the run (as far as I can tell
63 // this affects zero-length runs only).
64 update(geometry, m_grid->ctx()->time()->current(), 0);
65 }
66
67 void Verification::define_model_state_impl(const File &output) const {
68 m_mass_flux->define(output);
69 m_temperature->define(output);
70 }
71
72 void Verification::write_model_state_impl(const File &output) const {
73 m_mass_flux->write(output);
74 m_temperature->write(output);
75 }
76
77 MaxTimestep Verification::max_timestep_impl(double t) const {
78 (void) t;
79 return MaxTimestep("verification surface model");
80 }
81
82 /** Initialize climate inputs of tests K and O.
83 *
84 * @return 0 on success
85 */
86 void Verification::update_KO() {
87
88 m_mass_flux->set(0.0);
89 m_temperature->set(223.15);
90 }
91
92 /** Update the test L climate input (once).
93 *
94 * Unlike other `update_...()` methods, this one uses [kg m-2 s-1]
95 * as units of the climatic_mass_balance.
96 *
97 * @return 0 on success
98 */
99 void Verification::update_L() {
100 double A0, T0;
101
102 rheology::PatersonBuddCold tgaIce("stress_balance.sia.", *m_config, m_EC);
103
104 // compute T so that A0 = A(T) = Acold exp(-Qcold/(R T)) (i.e. for PatersonBuddCold);
105 // set all temps to this constant
106 A0 = 1.0e-16/secpera; // = 3.17e-24 1/(Pa^3 s); (EISMINT value) flow law parameter
107 T0 = tgaIce.tempFromSoftness(A0);
108
109 m_temperature->set(T0);
110
111 const double
112 ice_density = m_config->get_number("constants.ice.density"),
113 a0 = units::convert(m_sys, 0.3, "m year-1", "m second-1"),
114 L = 750e3,
115 Lsqr = L * L;
116
117 IceModelVec::AccessList list(*m_mass_flux);
118 for (Points p(*m_grid); p; p.next()) {
119 const int i = p.i(), j = p.j();
120
121 double r = radius(*m_grid, i, j);
122 (*m_mass_flux)(i, j) = a0 * (1.0 - (2.0 * r * r / Lsqr));
123
124 (*m_mass_flux)(i, j) *= ice_density; // convert to [kg m-2 s-1]
125 }
126 }
127
128 void Verification::update_V() {
129
130 // initialize temperature; the value used does not matter
131 m_temperature->set(273.15);
132
133 // initialize mass balance:
134 m_mass_flux->set(0.0);
135 }
136
137 void Verification::update_impl(const Geometry &geometry, double t, double dt) {
138 (void) geometry;
139 (void) dt;
140
141 switch (m_testname) {
142 case 'A':
143 case 'B':
144 case 'C':
145 case 'D':
146 case 'H':
147 update_ABCDH(t);
148 break;
149 case 'F':
150 case 'G':
151 update_FG(t);
152 break;
153 case 'K':
154 case 'O':
155 update_KO();
156 break;
157 case 'L':
158 {
159 update_L();
160 // return here; note update_L() uses correct units
161 return;
162 }
163 case 'V':
164 update_V();
165 break;
166 default:
167 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Test %c is not implemented.", m_testname);
168 }
169
170 // convert from [m second-1] to [kg m-2 s-1]
171 m_mass_flux->scale(m_config->get_number("constants.ice.density"));
172 }
173
174 /** Update climate inputs for tests A, B, C, D, E, H.
175 *
176 * @return 0 on success
177 */
178 void Verification::update_ABCDH(double time) {
179 double A0, T0, accum;
180
181 double f = m_config->get_number("constants.ice.density") / m_config->get_number("bed_deformation.mantle_density");
182
183 rheology::PatersonBuddCold tgaIce("stress_balance.sia.", *m_config, m_EC);
184
185 // compute T so that A0 = A(T) = Acold exp(-Qcold/(R T)) (i.e. for PatersonBuddCold);
186 // set all temps to this constant
187 A0 = 1.0e-16/secpera; // = 3.17e-24 1/(Pa^3 s); (EISMINT value) flow law parameter
188 T0 = tgaIce.tempFromSoftness(A0);
189
190 m_temperature->set(T0);
191
192 IceModelVec::AccessList list(*m_mass_flux);
193 ParallelSection loop(m_grid->com);
194 try {
195 for (Points p(*m_grid); p; p.next()) {
196 const int i = p.i(), j = p.j();
197
198 const double r = radius(*m_grid, i, j);
199 switch (m_testname) {
200 case 'A':
201 accum = exactA(r).M;
202 break;
203 case 'B':
204 accum = exactB(time, r).M;
205 break;
206 case 'C':
207 accum = exactC(time, r).M;
208 break;
209 case 'D':
210 accum = exactD(time, r).M;
211 break;
212 case 'H':
213 accum = exactH(f, time, r).M;
214 break;
215 default:
216 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "test must be A, B, C, D, or H, got %c",
217 m_testname);
218 }
219 (*m_mass_flux)(i, j) = accum;
220 }
221 } catch (...) {
222 loop.failed();
223 }
224 loop.check();
225 }
226
227 void Verification::update_FG(double time) {
228
229 const double t = m_testname == 'F' ? 0.0 : time;
230 const double A = m_testname == 'F' ? 0.0 : ApforG;
231
232 IceModelVec::AccessList list{m_mass_flux.get(), m_temperature.get()};
233
234 for (Points p(*m_grid); p; p.next()) {
235 const int i = p.i(), j = p.j();
236
237 // avoid singularity at origin
238 const double r = std::max(radius(*m_grid, i, j), 1.0);
239
240 (*m_temperature)(i, j) = Tmin + ST * r;
241
242 if (r > LforFG - 1.0) {
243 // if (essentially) outside of sheet
244 (*m_mass_flux)(i, j) = - ablationRateOutside / secpera;
245 } else {
246 (*m_mass_flux)(i, j) = exactFG(t, r, m_grid->z(), A).M;
247 }
248 }
249 }
250
251 } // end of namespace surface
252 } // end of namespace pism