PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
BTU_Full.cc
Go to the documentation of this file.
1 /* Copyright (C) 2016, 2017, 2018, 2019, 2020, 2021 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 
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::vector<double> z(m_Mbz);
61  double dz = m_Lbz / (m_Mbz - 1);
62  for (unsigned int k = 0; k < m_Mbz; ++k) {
63  z[k] = -m_Lbz + k * dz;
64  }
65  z.back() = 0.0;
66 
67  m_temp.reset(new IceModelVec3(m_grid, "litho_temp", WITHOUT_GHOSTS, z));
68  m_temp->metadata(0).z().set_name("zb");
69 
70  std::map<std::string, std::string> z_attrs =
71  {{"units", "m"},
72  {"long_name", "Z-coordinate in bedrock"},
73  {"axis", "Z"},
74  {"positive", "up"}};
75  for (const auto &z_attr : z_attrs) {
76  m_temp->metadata(0).z().set_string(z_attr.first, z_attr.second);
77  }
78 
79  m_temp->set_attrs("model_state",
80  "lithosphere (bedrock) temperature, in BTU_Full",
81  "K", "K", "", 0);
82  m_temp->metadata()["valid_min"] = {0.0};
83  }
84 
85  m_column.reset(new BedrockColumn("bedrock_column", *m_config, vertical_spacing(), Mz()));
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.
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 
122 }
123 
124 
125 /** Returns the vertical spacing used by the bedrock grid.
126  */
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 {
142  m_temp->define(output);
143 }
144 
145 void BTU_Full::write_model_state_impl(const File &output) const {
147  m_temp->write(output);
148 }
149 
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 
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) {
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 
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 */
214 
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 
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 
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
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
IceGrid::ConstPtr grid() const
Definition: Component.cc:105
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:142
@ REGRID_WITHOUT_REGRID_VARS
Definition: Component.hh:131
virtual void regrid(const std::string &module_name, IceModelVec &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS)
Definition: Component.cc:151
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
VariableLookupData find_variable(const std::string &short_name, const std::string &std_name) const
Find a variable using its standard name and/or short name.
Definition: File.cc:364
High-level PISM I/O class.
Definition: File.hh:51
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
void copy_from(const IceModelVec2S &source)
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: iceModelVec.hh:404
void define(const File &file, IO_Type default_type=PISM_DOUBLE) const
Define variables corresponding to an IceModelVec in a file opened using file.
Definition: iceModelVec.cc:523
void write(const std::string &filename) const
Definition: iceModelVec.cc:822
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
void failed()
Indicates a failure of a parallel section.
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
virtual MaxTimestep max_timestep_impl(double my_t) const
Definition: BTU_Full.cc:150
unsigned int m_Mbz
number of vertical levels within the bedrock
Definition: BTU_Full.hh:119
void update_flux_through_top_surface()
Definition: BTU_Full.cc:213
virtual unsigned int Mz_impl() const
Definition: BTU_Full.cc:131
double m_Lbz
thickness of the bedrock layer, in meters
Definition: BTU_Full.hh:121
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: BTU_Full.cc:145
std::shared_ptr< BedrockColumn > m_column
Definition: BTU_Full.hh:128
virtual void update_impl(const IceModelVec2S &bedrock_top_temperature, double t, double dt)=0
bool m_bootstrapping_needed
true if the model needs to "bootstrap" the temperature field during the first time step
Definition: BTU_Full.hh:124
virtual double depth_impl() const
Definition: BTU_Full.cc:136
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: BTU_Full.cc:140
IceModelVec3::Ptr m_temp
Definition: BTU_Full.hh:111
virtual double vertical_spacing_impl() const
Definition: BTU_Full.cc:127
virtual void bootstrap(const IceModelVec2S &bedrock_top_temperature)
Definition: BTU_Full.cc:254
BTU_Full(IceGrid::ConstPtr g, const BTUGrid &vertical_grid)
Definition: BTU_Full.cc:32
double m_D
diffusivity of the heat flow within the bedrock layer
Definition: BTU_Full.hh:116
virtual void init_impl(const InputOptions &opts)
Initialize the bedrock thermal unit.
Definition: BTU_Full.cc:90
const IceModelVec3 & temperature() const
Bedrock thermal layer temperature field.
Definition: BTU_Full.cc:246
double m_k
bedrock thermal conductivity
Definition: BTU_Full.hh:114
IceModelVec2S m_bottom_surface_flux
upward heat flux through the bottom surface of the bed thermal layer
virtual void init_impl(const InputOptions &opts)
Initialize the bedrock thermal unit.
IceModelVec2S m_top_surface_flux
upward heat flux through the top surface of the bed thermal layer
Given the temperature of the top of the bedrock, for the duration of one time-step,...
#define PISM_ERROR_LOCATION
#define rho
Definition: exactTestM.c:35
static const double g
Definition: exactTestP.cc:39
@ INIT_RESTART
Definition: Component.hh:39
static const double k
Definition: exactTestP.cc:45
@ PISM_GUESS
Definition: IO_Flags.hh:41
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:49
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
InitializationType type
initialization type
Definition: Component.hh:44
std::string filename
name of the input file (if applicable)
Definition: Component.hh:46
unsigned int record
index of the record to re-start from
Definition: Component.hh:48