tBTU_Full.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
---
tBTU_Full.cc (8458B)
---
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 #include "BTU_Full.hh"
20
21 #include "pism/util/pism_options.hh"
22 #include "pism/util/io/File.hh"
23 #include "pism/util/error_handling.hh"
24 #include "pism/util/pism_utilities.hh"
25 #include "pism/util/MaxTimestep.hh"
26 #include "BedrockColumn.hh"
27
28 namespace pism {
29 namespace energy {
30
31
32 BTU_Full::BTU_Full(IceGrid::ConstPtr g, const BTUGrid &grid)
33 : BedThermalUnit(g),
34 m_bootstrapping_needed(false) {
35
36 m_k = m_config->get_number("energy.bedrock_thermal.conductivity");
37
38 const double
39 rho = m_config->get_number("energy.bedrock_thermal.density"),
40 c = m_config->get_number("energy.bedrock_thermal.specific_heat_capacity");
41 // build constant diffusivity for heat equation
42 m_D = m_k / (rho * c);
43
44 // validate Lbz
45 if (grid.Lbz <= 0.0) {
46 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Invalid bedrock thermal layer depth: %f m",
47 grid.Lbz);
48 }
49
50 // validate Mbz
51 if (grid.Mbz < 2) {
52 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "Invalid number of layers of the bedrock thermal layer: %d",
53 grid.Mbz);
54 }
55
56 {
57 m_Mbz = grid.Mbz;
58 m_Lbz = grid.Lbz;
59
60 std::map<std::string, std::string> attrs;
61 attrs["units"] = "m";
62 attrs["long_name"] = "Z-coordinate in bedrock";
63 attrs["axis"] = "Z";
64 attrs["positive"] = "up";
65
66 std::vector<double> z(m_Mbz);
67 double dz = m_Lbz / (m_Mbz - 1);
68 for (unsigned int k = 0; k < m_Mbz; ++k) {
69 z[k] = -m_Lbz + k * dz;
70 }
71 z.back() = 0.0;
72
73 m_temp.reset(new IceModelVec3Custom(m_grid, "litho_temp", "zb", z, attrs));
74
75 m_temp->set_attrs("model_state",
76 "lithosphere (bedrock) temperature, in BTU_Full",
77 "K", "K", "", 0);
78 m_temp->metadata().set_number("valid_min", 0.0);
79 }
80
81 m_column.reset(new BedrockColumn("bedrock_column", *m_config, vertical_spacing(), Mz()));
82 }
83
84 BTU_Full::~BTU_Full() {
85 // empty
86 }
87
88
89 //! \brief Initialize the bedrock thermal unit.
90 void BTU_Full::init_impl(const InputOptions &opts) {
91
92 m_log->message(2, "* Initializing the bedrock thermal unit...\n");
93
94 // 2D initialization. Takes care of the flux through the bottom surface of the thermal layer.
95 BedThermalUnit::init_impl(opts);
96
97 // Initialize the temperature field.
98 {
99 // store the current "revision number" of the temperature field
100 const int temp_revision = m_temp->state_counter();
101
102 if (opts.type == INIT_RESTART) {
103 File input_file(m_grid->com, opts.filename, PISM_GUESS, PISM_READONLY);
104
105 if (input_file.find_variable("litho_temp")) {
106 m_temp->read(input_file, opts.record);
107 }
108 // otherwise the bedrock temperature is either interpolated from a -regrid_file or filled
109 // using bootstrapping (below)
110 }
111
112 regrid("bedrock thermal layer", *m_temp, REGRID_WITHOUT_REGRID_VARS);
113
114 if (m_temp->state_counter() == temp_revision) {
115 m_bootstrapping_needed = true;
116 } else {
117 m_bootstrapping_needed = false;
118 }
119 }
120
121 update_flux_through_top_surface();
122 }
123
124
125 /** Returns the vertical spacing used by the bedrock grid.
126 */
127 double BTU_Full::vertical_spacing_impl() const {
128 return m_Lbz / (m_Mbz - 1.0);
129 }
130
131 unsigned int BTU_Full::Mz_impl() const {
132 return m_Mbz;
133 }
134
135
136 double BTU_Full::depth_impl() const {
137 return m_Lbz;
138 }
139
140 void BTU_Full::define_model_state_impl(const File &output) const {
141 m_bottom_surface_flux.define(output);
142 m_temp->define(output);
143 }
144
145 void BTU_Full::write_model_state_impl(const File &output) const {
146 m_bottom_surface_flux.write(output);
147 m_temp->write(output);
148 }
149
150 MaxTimestep BTU_Full::max_timestep_impl(double t) const {
151 (void) t;
152 // No time step restriction: we are using an unconditionally stable method.
153 return MaxTimestep("bedrock thermal layer");
154 }
155
156
157 /** Perform a step of the bedrock thermal model.
158 */
159 void BTU_Full::update_impl(const IceModelVec2S &bedrock_top_temperature,
160 double t, double dt) {
161 (void) t;
162
163 if (m_bootstrapping_needed) {
164 bootstrap(bedrock_top_temperature);
165 m_bootstrapping_needed = false;
166 }
167
168 if (dt < 0) {
169 throw RuntimeError(PISM_ERROR_LOCATION, "dt < 0 is not allowed");
170 }
171
172 IceModelVec::AccessList list{m_temp.get(), &m_bottom_surface_flux, &bedrock_top_temperature};
173
174 ParallelSection loop(m_grid->com);
175 try {
176 for (Points p(*m_grid); p; p.next()) {
177 const int i = p.i(), j = p.j();
178
179 double *T = m_temp->get_column(i, j);
180
181 m_column->solve(dt, m_bottom_surface_flux(i, j), bedrock_top_temperature(i, j),
182 T, // input
183 T); // output
184
185 // Check that T is positive:
186 for (unsigned int k = 0; k < m_Mbz; ++k) {
187 if (T[k] <= 0.0) {
188 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
189 "invalid bedrock temperature: %f Kelvin at %d,%d,%d",
190 T[k], i, j, k);
191 }
192 }
193 }
194 } catch (...) {
195 loop.failed();
196 }
197 loop.check();
198
199 update_flux_through_top_surface();
200 }
201
202 /*! Computes the heat flux from the bedrock thermal layer upward into the
203 ice/bedrock interface:
204 \f[G_0 = -k_b \frac{\partial T_b}{\partial z}\big|_{z=0}.\f]
205 Uses the second-order finite difference expression
206 \f[\frac{\partial T_b}{\partial z}\big|_{z=0} \approx \frac{3 T_b(0) - 4 T_b(-\Delta z) + T_b(-2\Delta z)}{2 \Delta z}\f]
207 where \f$\Delta z\f$ is the equal spacing in the bedrock.
208
209 The above expression only makes sense when `Mbz` = `temp.n_levels` >= 3.
210 When `Mbz` = 2 we use first-order differencing. When temp was not created,
211 the `Mbz` <= 1 cases, we return the stored geothermal flux.
212 */
213 void BTU_Full::update_flux_through_top_surface() {
214
215 if (m_bootstrapping_needed) {
216 m_top_surface_flux.copy_from(m_bottom_surface_flux);
217 return;
218 }
219
220 double dz = this->vertical_spacing();
221 const int k0 = m_Mbz - 1; // Tb[k0] = ice/bed interface temp, at z=0
222
223 IceModelVec::AccessList list{m_temp.get(), &m_top_surface_flux};
224
225 if (m_Mbz >= 3) {
226
227 for (Points p(*m_grid); p; p.next()) {
228 const int i = p.i(), j = p.j();
229
230 const double *Tb = m_temp->get_column(i, j);
231 m_top_surface_flux(i, j) = - m_k * (3 * Tb[k0] - 4 * Tb[k0-1] + Tb[k0-2]) / (2 * dz);
232 }
233
234 } else {
235
236 for (Points p(*m_grid); p; p.next()) {
237 const int i = p.i(), j = p.j();
238
239 const double *Tb = m_temp->get_column(i, j);
240 m_top_surface_flux(i, j) = - m_k * (Tb[k0] - Tb[k0-1]) / dz;
241 }
242
243 }
244 }
245
246 const IceModelVec3Custom& BTU_Full::temperature() const {
247 if (m_bootstrapping_needed) {
248 throw RuntimeError(PISM_ERROR_LOCATION, "bedrock temperature is not available (bootstrapping is needed)");
249 }
250
251 return *m_temp;
252 }
253
254 void BTU_Full::bootstrap(const IceModelVec2S &bedrock_top_temperature) {
255
256 m_log->message(2,
257 " bootstrapping to fill lithosphere temperatures in the bedrock thermal layer\n"
258 " using temperature at the top bedrock surface and geothermal flux\n"
259 " (bedrock temperature is linear in depth)...\n");
260
261 double dz = this->vertical_spacing();
262 const int k0 = m_Mbz - 1; // Tb[k0] = ice/bedrock interface temp
263
264 IceModelVec::AccessList list{&bedrock_top_temperature, &m_bottom_surface_flux, m_temp.get()};
265 for (Points p(*m_grid); p; p.next()) {
266 const int i = p.i(), j = p.j();
267
268 double *Tb = m_temp->get_column(i, j); // Tb points into temp memory
269
270 Tb[k0] = bedrock_top_temperature(i, j);
271 for (int k = k0-1; k >= 0; k--) {
272 Tb[k] = Tb[k+1] + dz * m_bottom_surface_flux(i, j) / m_k;
273 }
274 }
275
276 m_temp->inc_state_counter(); // mark as modified
277 }
278
279 } // end of namespace energy
280 } // end of namespace pism