PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
AgeModel.cc
Go to the documentation of this file.
1/* Copyright (C) 2016, 2017, 2019, 2020, 2022, 2023, 2024 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
26namespace 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
43static 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
51 check_input(ice_thickness, "ice_thickness");
52 check_input(u3, "u3");
53 check_input(v3, "v3");
54 check_input(w3, "w3");
55}
56
57AgeModel::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/*!
75Let \f$\tau(t,x,y,z)\f$ be the age of the ice. Denote the three-dimensional
76velocity field within the ice fluid as \f$(u,v,w)\f$. The age equation
77is \f$d\tau/dt = 1\f$, that is, ice may move but it gets one year older in one
78year. 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]
81This equation is purely advective and hyperbolic. The right-hand side is "1" as
82long as age \f$\tau\f$ and time \f$t\f$ are measured in the same units.
83Because the velocity field is incompressible, \f$\nabla \cdot (u,v,w) = 0\f$,
84we can rewrite the equation as
85 \f[ \frac{\partial \tau}{\partial t} + \nabla \left( (u,v,w) \tau \right) = 1 \f]
86There is a conservative first-order numerical method; see AgeColumnSystem::solveThisColumn().
87
88The boundary condition is that when the ice falls as snow it has age zero.
89That is, \f$\tau(t,x,y,h(t,x,y)) = 0\f$ in accumulation areas. There is no
90boundary condition elsewhere on the ice upper surface, as the characteristics
91go outward in the ablation zone. If the velocity in the bottom cell of ice
92is 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,
94either 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
96that its most recent time in the atmosphere was long ago; this comment is
97relevant to any analysis which relates isotope ratios to modeled age.)
98
99The numerical method is a conservative form of first-order upwinding, but the
100vertical advection term is computed implicitly. Thus there is no CFL-type
101stability condition from the vertical velocity; CFL is only for the horizontal
102velocity. We use a finely-spaced, equally-spaced vertical grid in the
103calculation. Note that the columnSystemCtx methods coarse_to_fine() and
104fine_to_coarse() interpolate back and forth between this fine grid and
105the storage grid. The storage grid may or may not be equally-spaced. See
106AgeColumnSystem::solve() for the actual method.
107 */
108void 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
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
187void 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.variable_exists("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
214void AgeModel::define_model_state_impl(const File &output) const {
216}
217
218void 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
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
bool variable_exists(const std::string &short_name) const
Checks if a variable exists.
Definition File.cc:378
High-level PISM I/O class.
Definition File.hh:55
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
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:64
void copy_from(const Array3D &input)
Definition Array3D.cc:209
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:50
double * get_column(int i, int j)
Definition Array3D.cc:121
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:731
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:463
void write(const std::string &filename) const
Definition Array.cc:722
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition Array.cc:476
Abstract class for reading, writing, allocating, and accessing a DA-based PETSc Vec (2D and 3D fields...
Definition Array.hh:207
const std::vector< double > & z() const
void fine_to_coarse(const std::vector< double > &input, int i, int j, array::Array3D &output) const
unsigned int ks() const
#define PISM_ERROR_LOCATION
@ PISM_GUESS
Definition IO_Flags.hh:56
@ PISM_READONLY
open an existing file for reading only
Definition IO_Flags.hh:68
@ 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