PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
NullTransport.cc
Go to the documentation of this file.
1 // Copyright (C) 2012-2023 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 "pism/hydrology/NullTransport.hh"
20 #include "pism/geometry/Geometry.hh"
21 #include "pism/util/MaxTimestep.hh"
22 #include "pism/util/array/CellType.hh"
23 #include "pism/util/error_handling.hh"
24 #include "pism/util/pism_utilities.hh" // clip
25 
26 namespace pism {
27 namespace hydrology {
28 
29 NullTransport::NullTransport(std::shared_ptr<const Grid> g)
30  : Hydrology(g), m_Wtill_old(m_grid, "Wtill_old") {
31 
32  m_diffuse_tillwat = m_config->get_flag("hydrology.null_diffuse_till_water");
33  m_diffusion_time = m_config->get_number("hydrology.null_diffusion_time", "seconds");
34  m_diffusion_distance = m_config->get_number("hydrology.null_diffusion_distance", "meters");
35  m_tillwat_max = m_config->get_number("hydrology.tillwat_max", "meters");
36  m_tillwat_decay_rate = m_config->get_number("hydrology.tillwat_decay_rate", "m / second");
37 
38  if (m_tillwat_max < 0.0) {
40  "hydrology::NullTransport: hydrology.tillwat_max is negative.\n"
41  "This is not allowed.");
42  }
43 }
44 
46  m_log->message(2,
47  "* Initializing the null-transport (till only) subglacial hydrology model ...\n");
48 
49  if (m_diffuse_tillwat) {
50  m_log->message(
51  2, " [using lateral diffusion of stored till water as in Bueler and Brown, 2009]\n");
52  }
53 }
54 
55 void NullTransport::restart_impl(const File &input_file, int record) {
56  Hydrology::restart_impl(input_file, record);
57 }
58 
59 void NullTransport::bootstrap_impl(const File &input_file, const array::Scalar &ice_thickness) {
60  Hydrology::bootstrap_impl(input_file, ice_thickness);
61 }
62 
64  const array::Scalar &P) {
65  Hydrology::init_impl(W_till, W, P);
66 }
67 
69  (void)t;
70  if (m_diffuse_tillwat) {
71  const double dx2 = m_grid->dx() * m_grid->dx(), dy2 = m_grid->dy() * m_grid->dy(),
72  L = m_diffusion_distance, T = m_diffusion_time, K = L * L / (2.0 * T);
73 
74  return MaxTimestep(dx2 * dy2 / (2.0 * K * (dx2 + dy2)), "null-transport hydrology");
75  }
76 
77  return MaxTimestep("null-transport hydrology");
78 }
79 
80 //! Update the till water thickness by simply integrating the melt input.
81 /*!
82  Does a step of the trivial integration
83 
84  \f[ \frac{\partial W_{till}}{\partial t} = \frac{m}{\rho_w} - C\f]
85 
86  where \f$C=\f$`hydrology_tillwat_decay_rate`. Enforces bounds
87  \f$0 \le W_{till} \le W_{till}^{max}\f$ where the upper bound is
88  `hydrology_tillwat_max`. Here \f$m/\rho_w\f$ is `total_input`.
89 
90  Uses the current mass-continuity timestep `dt`. (Compare
91  hydrology::Routing::update_Wtill() which will generally be taking time steps
92  determined by the evolving transportable water layer in that model.)
93 
94  There is no attempt to report on conservation errors because this
95  hydrology::NullTransport model does not conserve water.
96 
97  There is no tranportable water thickness variable and no interaction with it.
98 */
99 void NullTransport::update_impl(double t, double dt, const Inputs& inputs) {
100  (void) t;
101 
102  bool add_surface_input = m_config->get_flag("hydrology.add_water_input_to_till_storage");
103 
104  // no transportable basal water
105  m_W.set(0.0);
106 
108  if (add_surface_input) {
110  }
111 
112  const double
113  water_density = m_config->get_number("constants.fresh_water.density"),
114  kg_per_m = m_grid->cell_area() * water_density; // kg m-1
115 
116  const auto &cell_type = inputs.geometry->cell_type;
117 
118  array::AccessScope list{&cell_type, &m_Wtill, &m_basal_melt_rate,
120 
121  if (add_surface_input) {
123  }
124 
125  for (auto p = m_grid->points(); p; p.next()) {
126  const int i = p.i(), j = p.j();
127 
128  double
129  W_old = m_Wtill(i, j),
130  dW_input = dt * m_basal_melt_rate(i, j);
131 
132  if (add_surface_input) {
133  dW_input += dt * m_surface_input_rate(i, j);
134  }
135 
136  if (W_old < 0.0) {
138  "negative subglacial water thickness of %f m at (%d, %d)",
139  W_old, i, j);
140  }
141 
142  // decay rate in areas under grounded ice
143  double dW_decay = dt * (- m_tillwat_decay_rate);
144 
145  if (not cell_type.grounded_ice(i, j)) {
146  dW_decay = 0.0;
147  }
148 
149  m_Wtill(i, j) = (W_old + dW_input) + dW_decay;
150 
151  // dW_decay is a "conservation error"
152  m_conservation_error_change(i, j) += dW_decay * kg_per_m;
153  }
154 
155  if (m_diffuse_tillwat) {
156  diffuse_till_water(dt);
157  }
158 
159  // remove water in ice-free areas and account for changes
161  inputs.no_model_mask,
162  m_tillwat_max, // global maximum till water thickness
163  m_tillwat_max, // till water thickness under the ocean
164  m_Wtill,
169 }
170 
172  // note: this call updates ghosts of m_Wtill_old
174 
175  const double
176  dx = m_grid->dx(),
177  dy = m_grid->dy(),
179  T = m_diffusion_time,
180  K = L * L / (2.0 * T),
181  Rx = K * dt / (dx * dx),
182  Ry = K * dt / (dy * dy);
183 
185  for (auto p = m_grid->points(); p; p.next()) {
186  const int i = p.i(), j = p.j();
187 
188  auto W = m_Wtill_old.star(i, j);
189 
190  double dWtill = (Rx * (W.w - 2.0 * W.c + W.e) + Ry * (W.s - 2.0 * W.c + W.n));
191 
192  m_Wtill(i, j) = W.c + dWtill;
193 
194  m_flow_change(i, j) = dWtill;
195  }
196 }
197 
198 } // end of namespace hydrology
199 } // end of namespace pism
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:158
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:162
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
High-level PISM I/O class.
Definition: File.hh:56
array::CellType2 cell_type
Definition: Geometry.hh:55
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
void copy_from(const Array2D< T > &source)
Definition: Array2D.hh:73
stencils::Star< T > star(int i, int j) const
Definition: Array2D.hh:79
void add(double alpha, const Array2D< T > &x)
Definition: Array2D.hh:65
void set(double c)
Result: v[j] <- c for all j.
Definition: Array.cc:707
array::Scalar m_flow_change
Definition: Hydrology.hh:196
virtual void bootstrap_impl(const File &input_file, const array::Scalar &ice_thickness)
Definition: Hydrology.cc:377
array::Scalar m_surface_input_rate
Definition: Hydrology.hh:179
array::Scalar m_grounding_line_change
Definition: Hydrology.hh:192
void enforce_bounds(const array::CellType &cell_type, const array::Scalar *no_model_mask, double max_thickness, double ocean_water_thickness, array::Scalar &water_thickness, array::Scalar &grounded_margin_change, array::Scalar &grounding_line_change, array::Scalar &conservation_error_change, array::Scalar &no_model_mask_change)
Correct the new water thickness based on boundary requirements.
Definition: Hydrology.cc:669
array::Scalar m_basal_melt_rate
Definition: Hydrology.hh:182
virtual void init_impl(const array::Scalar &W_till, const array::Scalar &W, const array::Scalar &P)
Definition: Hydrology.cc:387
array::Scalar1 m_W
effective thickness of transportable basal water
Definition: Hydrology.hh:173
virtual void restart_impl(const File &input_file, int record)
Definition: Hydrology.cc:370
array::Scalar m_input_change
Definition: Hydrology.hh:193
array::Scalar m_Wtill
effective thickness of basal water stored in till
Definition: Hydrology.hh:170
array::Scalar m_grounded_margin_change
Definition: Hydrology.hh:191
array::Scalar m_no_model_mask_change
Definition: Hydrology.hh:194
array::Scalar m_conservation_error_change
Definition: Hydrology.hh:190
The PISM subglacial hydrology model interface.
Definition: Hydrology.hh:109
const Geometry * geometry
Definition: Hydrology.hh:37
const array::Scalar1 * no_model_mask
Definition: Hydrology.hh:35
virtual void restart_impl(const File &input_file, int record)
virtual void init_impl(const array::Scalar &W_till, const array::Scalar &W, const array::Scalar &P)
virtual void bootstrap_impl(const File &input_file, const array::Scalar &ice_thickness)
virtual void update_impl(double t, double dt, const Inputs &inputs)
Solves an implicit step of a highly-simplified ODE.
virtual MaxTimestep max_timestep_impl(double t) const
NullTransport(std::shared_ptr< const Grid > g)
#define PISM_ERROR_LOCATION
static const double L
Definition: exactTestL.cc:40
static double K(double psi_x, double psi_y, double speed, double epsilon)
static const double g
Definition: exactTestP.cc:36
const double dy2
Definition: ssafd_code.cc:1
const double dx2
Definition: ssafd_code.cc:1