PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
AgeModel.cc
Go to the documentation of this file.
1 /* Copyright (C) 2016, 2017, 2019, 2020, 2022, 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 
20 #include "pism/age/AgeModel.hh"
21 #include "pism/age/AgeColumnSystem.hh"
22 #include "pism/util/error_handling.hh"
23 #include "pism/util/io/File.hh"
24 #include <memory>
25 
26 namespace pism {
27 
29  const array::Array3D *u,
30  const array::Array3D *v,
31  const array::Array3D *w)
32  : ice_thickness(thickness), u3(u), v3(v), w3(w) {
33  // empty
34 }
35 
37  ice_thickness = NULL;
38  u3 = NULL;
39  v3 = NULL;
40  w3 = NULL;
41 }
42 
43 static void check_input(const array::Array *ptr, const char *name) {
44  if (ptr == NULL) {
46  "ice age model input %s was not provided", name);
47  }
48 }
49 
50 void AgeModelInputs::check() const {
51  check_input(ice_thickness, "ice_thickness");
52  check_input(u3, "u3");
53  check_input(v3, "v3");
54  check_input(w3, "w3");
55 }
56 
57 AgeModel::AgeModel(std::shared_ptr<const Grid> grid,
58  std::shared_ptr<const stressbalance::StressBalance> stress_balance)
59  : Component(grid),
60  // FIXME: should be able to use width=1...
61  m_ice_age(m_grid, "age", array::WITH_GHOSTS, m_grid->z(), m_config->get_number("grid.max_stencil_width")),
62  m_work(m_grid, "work_vector", array::WITHOUT_GHOSTS, m_grid->z()),
63  m_stress_balance(stress_balance) {
64 
66  .long_name("age of ice")
67  .units("s");
68 
69  m_ice_age.metadata()["valid_min"] = {0.0};
70 
71  m_work.metadata().units("s");
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 array::Scalar &ice_thickness = *inputs.ice_thickness;
116 
117  const array::Array3D
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  array::AccessScope 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 (auto p = m_grid->points(); 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 array::Array3D
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 
170 }
171 
172 const array::Array3D & AgeModel::age() const {
173  return m_ice_age;
174 }
175 
177 
178  if (m_stress_balance == nullptr) {
180  "AgeModel: no stress balance provided."
181  " Cannot compute max. time step.");
182  }
183 
184  return MaxTimestep(m_stress_balance->max_timestep_cfl_3d().dt_max.value(), "age model");
185 }
186 
187 void AgeModel::init(const InputOptions &opts) {
188 
189  m_log->message(2, "* Initializing the age model...\n");
190 
191 
192  double initial_age_years = m_config->get_number("age.initial_value", "years");
193 
194  if (opts.type == INIT_RESTART) {
195  File input_file(m_grid->com, opts.filename, io::PISM_GUESS, io::PISM_READONLY);
196 
197  if (input_file.find_variable("age")) {
198  m_ice_age.read(input_file, opts.record);
199  } else {
200  m_log->message(2,
201  "PISM WARNING: input file '%s' does not have the 'age' variable.\n"
202  " Setting it to %f years...\n",
203  opts.filename.c_str(), initial_age_years);
204  m_ice_age.set(m_config->get_number("age.initial_value", "seconds"));
205  }
206  } else {
207  m_log->message(2, " - setting initial age to %.4f years\n", initial_age_years);
208  m_ice_age.set(m_config->get_number("age.initial_value", "seconds"));
209  }
210 
212 }
213 
214 void AgeModel::define_model_state_impl(const File &output) const {
216 }
217 
218 void AgeModel::write_model_state_impl(const File &output) const {
219  m_ice_age.write(output);
220 }
221 
222 } // end of namespace pism
void init(int i, int j, double thickness)
void solve(std::vector< double > &x)
First-order upwind scheme with implicit in the vertical: one column solve.
Tridiagonal linear system for vertical column of age (pure advection) problem.
const array::Array3D * w3
Definition: AgeModel.hh:41
void check() const
Definition: AgeModel.cc:50
const array::Array3D * v3
Definition: AgeModel.hh:40
const array::Array3D * u3
Definition: AgeModel.hh:39
const array::Scalar * ice_thickness
Definition: AgeModel.hh:38
void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: AgeModel.cc:214
array::Array3D m_ice_age
Definition: AgeModel.hh:59
AgeModel(std::shared_ptr< const Grid > grid, std::shared_ptr< const stressbalance::StressBalance > stress_balance)
Definition: AgeModel.cc:57
std::shared_ptr< const stressbalance::StressBalance > m_stress_balance
Definition: AgeModel.hh:61
array::Array3D m_work
Definition: AgeModel.hh:60
MaxTimestep max_timestep_impl(double t) const
Definition: AgeModel.cc:176
const array::Array3D & age() const
Definition: AgeModel.cc:172
void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: AgeModel.cc:218
void init(const InputOptions &opts)
Definition: AgeModel.cc:187
void update(double t, double dt, const AgeModelInputs &inputs)
Definition: AgeModel.cc:108
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
@ REGRID_WITHOUT_REGRID_VARS
Definition: Component.hh:151
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
void regrid(const std::string &module_name, array::Array &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS)
Definition: Component.cc:159
A class defining a common interface for most PISM sub-models.
Definition: Component.hh:118
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:361
High-level PISM I/O class.
Definition: File.hh:56
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
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
void copy_from(const Array3D &input)
Definition: Array3D.cc:208
void set_column(int i, int j, double c)
Set all values of scalar quantity to given a single value in a particular column.
Definition: Array3D.cc:49
double * get_column(int i, int j)
Definition: Array3D.cc:120
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: Array3D.hh:33
void read(const std::string &filename, unsigned int time)
Definition: Array.cc:809
void define(const File &file, io::Type default_type) const
Define variables corresponding to an Array in a file opened using file.
Definition: Array.cc:540
void write(const std::string &filename) const
Definition: Array.cc:800
void set(double c)
Result: v[j] <- c for all j.
Definition: Array.cc:707
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: Array.cc:553
Abstract class for reading, writing, allocating, and accessing a DA-based PETSc Vec (2D and 3D fields...
Definition: Array.hh:208
const std::vector< double > & z() const
unsigned int ks() const
void fine_to_coarse(const std::vector< double > &fine, int i, int j, array::Array3D &coarse) const
#define PISM_ERROR_LOCATION
@ WITH_GHOSTS
Definition: Array.hh:62
@ WITHOUT_GHOSTS
Definition: Array.hh:62
static Vector3 column(const double A[3][3], size_t k)
Definition: Element.cc:51
@ PISM_GUESS
Definition: IO_Flags.hh:56
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:72
@ PISM_DOUBLE
Definition: IO_Flags.hh:52
static void check_input(const array::Array *ptr, const char *name)
Definition: AgeModel.cc:43
@ INIT_RESTART
Definition: Component.hh:56
static const double k
Definition: exactTestP.cc:42
InitializationType type
initialization type
Definition: Component.hh:61
std::string filename
name of the input file (if applicable)
Definition: Component.hh:63
unsigned int record
index of the record to re-start from
Definition: Component.hh:65