tAgeModel.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tAgeModel.cc (8050B)
---
1 /* Copyright (C) 2016, 2017, 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 "AgeModel.hh"
21
22 #include "pism/age/AgeColumnSystem.hh"
23 #include "pism/util/error_handling.hh"
24 #include "pism/util/Vars.hh"
25 #include "pism/util/io/File.hh"
26
27 namespace pism {
28
29 AgeModelInputs::AgeModelInputs(const IceModelVec2S *thickness,
30 const IceModelVec3 *u,
31 const IceModelVec3 *v,
32 const IceModelVec3 *w)
33 : ice_thickness(thickness), u3(u), v3(v), w3(w) {
34 // empty
35 }
36
37 AgeModelInputs::AgeModelInputs() {
38 ice_thickness = NULL;
39 u3 = NULL;
40 v3 = NULL;
41 w3 = NULL;
42 }
43
44 static void check_input(const IceModelVec *ptr, const char *name) {
45 if (ptr == NULL) {
46 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
47 "ice age model input %s was not provided", name);
48 }
49 }
50
51 void AgeModelInputs::check() const {
52 check_input(ice_thickness, "ice_thickness");
53 check_input(u3, "u3");
54 check_input(v3, "v3");
55 check_input(w3, "w3");
56 }
57
58 AgeModel::AgeModel(IceGrid::ConstPtr grid, stressbalance::StressBalance *stress_balance)
59 : Component(grid),
60 // FIXME: should be able to use width=1...
61 m_ice_age(m_grid, "age", WITH_GHOSTS, m_config->get_number("grid.max_stencil_width")),
62 m_work(m_grid, "work_vector", WITHOUT_GHOSTS),
63 m_stress_balance(stress_balance) {
64
65 m_ice_age.set_attrs("model_state", "age of ice",
66 "s", "years", "" /* no standard name*/, 0);
67
68 m_ice_age.metadata().set_number("valid_min", 0.0);
69
70 m_work.set_attrs("internal", "new values of age during time step",
71 "s", "s", "", 0);
72 }
73
74 /*!
75 Let \f$\tau(t,x,y,z)\f$ be the age of the ice. Denote the three-dimensional
76 velocity field within the ice fluid as \f$(u,v,w)\f$. The age equation
77 is \f$d\tau/dt = 1\f$, that is, ice may move but it gets one year older in one
78 year. Thus
79 \f[ \frac{\partial \tau}{\partial t} + u \frac{\partial \tau}{\partial x}
80 + v \frac{\partial \tau}{\partial y} + w \frac{\partial \tau}{\partial z} = 1 \f]
81 This equation is purely advective and hyperbolic. The right-hand side is "1" as
82 long as age \f$\tau\f$ and time \f$t\f$ are measured in the same units.
83 Because the velocity field is incompressible, \f$\nabla \cdot (u,v,w) = 0\f$,
84 we can rewrite the equation as
85 \f[ \frac{\partial \tau}{\partial t} + \nabla \left( (u,v,w) \tau \right) = 1 \f]
86 There is a conservative first-order numerical method; see AgeColumnSystem::solveThisColumn().
87
88 The boundary condition is that when the ice falls as snow it has age zero.
89 That is, \f$\tau(t,x,y,h(t,x,y)) = 0\f$ in accumulation areas. There is no
90 boundary condition elsewhere on the ice upper surface, as the characteristics
91 go outward in the ablation zone. If the velocity in the bottom cell of ice
92 is upward (\f$w>0\f$) then we also apply a zero age boundary condition,
93 \f$\tau(t,x,y,0) = 0\f$. This is the case where ice freezes on at the base,
94 either grounded basal ice freezing on stored water in till, or marine basal ice.
95 (Note that the water that is frozen-on as ice might be quite "old" in the sense
96 that its most recent time in the atmosphere was long ago; this comment is
97 relevant to any analysis which relates isotope ratios to modeled age.)
98
99 The numerical method is a conservative form of first-order upwinding, but the
100 vertical advection term is computed implicitly. Thus there is no CFL-type
101 stability condition from the vertical velocity; CFL is only for the horizontal
102 velocity. We use a finely-spaced, equally-spaced vertical grid in the
103 calculation. Note that the columnSystemCtx methods coarse_to_fine() and
104 fine_to_coarse() interpolate back and forth between this fine grid and
105 the storage grid. The storage grid may or may not be equally-spaced. See
106 AgeColumnSystem::solve() for the actual method.
107 */
108 void AgeModel::update(double t, double dt, const AgeModelInputs &inputs) {
109
110 // fix a compiler warning
111 (void) t;
112
113 inputs.check();
114
115 const IceModelVec2S &ice_thickness = *inputs.ice_thickness;
116
117 const IceModelVec3
118 &u3 = *inputs.u3,
119 &v3 = *inputs.v3,
120 &w3 = *inputs.w3;
121
122 AgeColumnSystem system(m_grid->z(), "age",
123 m_grid->dx(), m_grid->dy(), dt,
124 m_ice_age, u3, v3, w3); // linear system to solve in each column
125
126 size_t Mz_fine = system.z().size();
127 std::vector<double> x(Mz_fine); // space for solution
128
129 IceModelVec::AccessList list{&ice_thickness, &u3, &v3, &w3, &m_ice_age, &m_work};
130
131 unsigned int Mz = m_grid->Mz();
132
133 ParallelSection loop(m_grid->com);
134 try {
135 for (Points p(*m_grid); p; p.next()) {
136 const int i = p.i(), j = p.j();
137
138 system.init(i, j, ice_thickness(i, j));
139
140 if (system.ks() == 0) {
141 // if no ice, set the entire column to zero age
142 m_work.set_column(i, j, 0.0);
143 } else {
144 // general case: solve advection PDE
145
146 // solve the system for this column; call checks that params set
147 system.solve(x);
148
149 // put solution in IceModelVec3
150 system.fine_to_coarse(x, i, j, m_work);
151
152 // Ensure that the age of the ice is non-negative.
153 //
154 // FIXME: this is a kludge. We need to ensure that our numerical method has the maximum
155 // principle instead. (We may still need this for correctness, though.)
156 double *column = m_work.get_column(i, j);
157 for (unsigned int k = 0; k < Mz; ++k) {
158 if (column[k] < 0.0) {
159 column[k] = 0.0;
160 }
161 }
162 }
163 }
164 } catch (...) {
165 loop.failed();
166 }
167 loop.check();
168
169 m_work.update_ghosts(m_ice_age);
170 }
171
172 const IceModelVec3 & AgeModel::age() const {
173 return m_ice_age;
174 }
175
176 MaxTimestep AgeModel::max_timestep_impl(double t) const {
177 // fix a compiler warning
178 (void) t;
179
180 if (m_stress_balance == NULL) {
181 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
182 "AgeModel: no stress balance provided."
183 " Cannot compute max. time step.");
184 }
185
186 return MaxTimestep(m_stress_balance->max_timestep_cfl_3d().dt_max.value(), "age model");
187 }
188
189 void AgeModel::init(const InputOptions &opts) {
190
191 m_log->message(2, "* Initializing the age model...\n");
192
193
194 double initial_age_years = m_config->get_number("age.initial_value", "years");
195
196 if (opts.type == INIT_RESTART) {
197 File input_file(m_grid->com, opts.filename, PISM_GUESS, PISM_READONLY);
198
199 if (input_file.find_variable("age")) {
200 m_ice_age.read(input_file, opts.record);
201 } else {
202 m_log->message(2,
203 "PISM WARNING: input file '%s' does not have the 'age' variable.\n"
204 " Setting it to %f years...\n",
205 opts.filename.c_str(), initial_age_years);
206 m_ice_age.set(m_config->get_number("age.initial_value", "seconds"));
207 }
208 } else {
209 m_log->message(2, " - setting initial age to %.4f years\n", initial_age_years);
210 m_ice_age.set(m_config->get_number("age.initial_value", "seconds"));
211 }
212
213 regrid("Age Model", m_ice_age, REGRID_WITHOUT_REGRID_VARS);
214 }
215
216 void AgeModel::define_model_state_impl(const File &output) const {
217 m_ice_age.define(output);
218 }
219
220 void AgeModel::write_model_state_impl(const File &output) const {
221 m_ice_age.write(output);
222 }
223
224 } // end of namespace pism