PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
ForceThickness.cc
Go to the documentation of this file.
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 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 "ForceThickness.hh"
20 #include "pism/util/IceGrid.hh"
21 #include "pism/util/Vars.hh"
22 
23 #include "pism/util/ConfigInterface.hh"
24 #include "pism/util/Mask.hh"
25 #include "pism/util/pism_options.hh"
26 #include "pism/util/error_handling.hh"
27 #include "pism/util/io/io_helpers.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/util/IceModelVec2CellType.hh"
30 #include "pism/util/MaxTimestep.hh"
31 #include "pism/util/io/File.hh"
32 #include "pism/geometry/Geometry.hh"
33 
34 namespace pism {
35 namespace surface {
36 
37 ///// "Force-to-thickness" mechanism
38 ForceThickness::ForceThickness(IceGrid::ConstPtr g, std::shared_ptr<SurfaceModel> input)
39  : SurfaceModel(g, input),
40  m_target_thickness(m_grid, "thk", WITHOUT_GHOSTS),
41  m_ftt_mask(m_grid, "ftt_mask", WITHOUT_GHOSTS)
42 {
43 
44  m_alpha = m_config->get_number("surface.force_to_thickness.alpha", "s-1");
45  m_alpha_ice_free_factor = m_config->get_number("surface.force_to_thickness.ice_free_alpha_factor");
46  m_ice_free_thickness_threshold = m_config->get_number("surface.force_to_thickness.ice_free_thickness_threshold");
47  m_start_time = m_config->get_number("surface.force_to_thickness.start_time", "seconds");
48 
49  m_ftt_mask.set_attrs("diagnostic",
50  "mask specifying where to apply the force-to-thickness mechanism",
51  "", "", "", 0); // no units and no standard name
52  m_ftt_mask.set(1.0); // default: applied in whole domain
55 
57 
61 }
62 
63 void ForceThickness::init_impl(const Geometry &geometry) {
64 
65  m_input_model->init(geometry);
66 
67  m_log->message(2,
68  "* Initializing force-to-thickness mass-balance modifier...\n");
69 
70  std::string input_file = m_config->get_string("surface.force_to_thickness_file");
71 
72  if (input_file.empty()) {
73  throw RuntimeError(PISM_ERROR_LOCATION, "surface.force_to_thickness_file cannot be empty");
74  }
75 
76  m_log->message(2,
77  " alpha = %.6f year-1 for -force_to_thickness mechanism\n"
78  " alpha = %.6f year-1 in areas with target ice thickness of less than %.3f meters\n",
79  units::convert(m_sys, m_alpha, "s-1", "year-1"),
82 
83  // check of the input file is really there and regrid the target thickness
84  File file(m_grid->com, input_file, PISM_GUESS, PISM_READONLY);
85 
86  m_log->message(2,
87  " reading target thickness 'thk' from %s ...\n"
88  " (this field will appear in output file as 'ftt_target_thk')\n",
89  input_file.c_str());
90  {
91  m_target_thickness.metadata(0).set_name("thk"); // name to read by
92  // set attributes for the read stage; see below for reset
93  m_target_thickness.set_attrs("diagnostic",
94  "target thickness for force-to-thickness mechanism (hit this at end of run)",
95  "m", "m",
96  "land_ice_thickness", 0); // standard_name *to read by*
97 
98  m_target_thickness.regrid(input_file, CRITICAL);
99 
100  // reset name to avoid confusion; set attributes again to overwrite "read by" choices above
101  m_target_thickness.metadata(0).set_name("ftt_target_thk");
102  m_target_thickness.set_attrs("diagnostic",
103  "target thickness for force-to-thickness mechanism (wants to hit this at end of run)",
104  "m", "m",
105  "", 0); // no CF standard_name, to put it mildly
106  }
107 
108  {
109  m_log->message(2,
110  " reading force-to-thickness mask 'ftt_mask' from %s ...\n",
111  input_file.c_str());
112  m_ftt_mask.regrid(input_file, CRITICAL);
113  }
114 }
115 
116 /*!
117 If `-force_to_thickness_file` `foo.nc` is in use then vthktarget will have a target ice thickness
118 map. Let \f$H_{\text{tar}}\f$ be this target thickness,
119 and let \f$H\f$ be the current model thickness. Recall that the mass continuity
120 equation solved by GeometryEvolution is
121  \f[ \frac{\partial H}{\partial t} = M - S - \nabla\cdot \mathbf{q} \f]
122 and that this procedure is supposed to produce \f$M\f$.
123 In this context, the semantics of `-force_to_thickness` are that \f$M\f$ is modified
124 by a multiple of the difference between the target thickness and the current thickness.
125 In particular, the \f$\Delta M\f$ that is produced here is
126  \f[\Delta M = \alpha (H_{\text{tar}} - H)\f]
127 where \f$\alpha>0\f$ is determined below. Note \f$\Delta M\f$ is positive in
128 areas where \f$H_{\text{tar}} > H\f$, so we are adding mass there, and we are ablating
129 in the other case.
130 
131 Let \f$t_s\f$ be the start time and \f$t_e\f$ the end time for the run.
132 Without flow or basal mass balance, or any surface mass balance other than the
133 \f$\Delta M\f$ computed here, we are solving
134  \f[ \frac{\partial H}{\partial t} = \alpha (H_{\text{tar}} - H) \f]
135 Let's assume \f$H(t_s)=H_0\f$. This initial value problem has solution
136 \f$H(t) = H_{\text{tar}} + (H_0 - H_{\text{tar}}) e^{-\alpha (t-t_s)}\f$
137 and so
138  \f[ H(t_e) = H_{\text{tar}} + (H_0 - H_{\text{tar}}) e^{-\alpha (t_e-t_s)} \f]
139 
140 The constant \f$\alpha\f$ has a default value `pism_config:surface.force_to_thickness.alpha`.
141 
142 The next example uses files generated from the EISMINT-Greenland experiment;
143 see the corresponding chapter of the User's Manual.
144 
145 Suppose we regard the SSL2 run as a spin-up to reach a better temperature field.
146 It is a spinup in which the surface was allowed to evolve. Assume the
147 early file `green20km_y1.nc` has the target thickness, because it essentially
148 has the input thickness. This script adds a 500 a run, to finalize the spinup,
149 in which the ice sheet geometry goes from the the thickness values in
150 `green_ssl2_110ka.nc` to values very close to those in `green20km_y1.nc`:
151 \code
152 #!/bin/bash
153 
154 NN=8 # default number of processors
155 if [ $# -gt 0 ] ; then # if user says "test_ftt.sh 8" then NN = 8
156  NN="$1"
157 fi
158 
159 # set MPIDO if using different MPI execution command, for example:
160 # $ export PISM_MPIDO="aprun -n "
161 if [ -n "${PISM_MPIDO:+1}" ] ; then # check if env var is already set
162  echo "$SCRIPTNAME PISM_MPIDO = $PISM_MPIDO (already set)"
163 else
164  PISM_MPIDO="mpiexec -n "
165  echo "$SCRIPTNAME PISM_MPIDO = $PISM_MPIDO"
166 fi
167 
168 # check if env var PISM_DO was set (i.e. PISM_DO=echo for a 'dry' run)
169 if [ -n "${PISM_DO:+1}" ] ; then # check if env var DO is already set
170  echo "$SCRIPTNAME PISM_DO = $PISM_DO (already set)"
171 else
172  PISM_DO=""
173 fi
174 
175 # prefix to pism (not to executables)
176 if [ -n "${PISM_PREFIX:+1}" ] ; then # check if env var is already set
177  echo "$SCRIPTNAME PISM_PREFIX = $PISM_PREFIX (already set)"
178 else
179  PISM_PREFIX="" # just a guess
180  echo "$SCRIPTNAME PISM_PREFIX = $PISM_PREFIX"
181 fi
182 
183 # set PISM_EXEC if using different executables, for example:
184 # $ export PISM_EXEC="pismr -energy cold"
185 if [ -n "${PISM_EXEC:+1}" ] ; then # check if env var is already set
186  echo "$SCRIPTNAME PISM_EXEC = $PISM_EXEC (already set)"
187 else
188  PISM_EXEC="pismr"
189  echo "$SCRIPTNAME PISM_EXEC = $PISM_EXEC"
190 fi
191 
192 
193 PISM="${PISM_PREFIX}${PISM_EXEC}"
194 
195 cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
196  -surface pdd -pdd_fausto \
197  -o no_force.nc -ts_file ts_no_force.nc -ts_times -1000:yearly:0"
198 $PISM_DO $cmd
199 
200 echo
201 
202 cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
203  -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc \
204  -o default_force.nc -ts_file ts_default_force.nc -ts_times -1000:yearly:0"
205 $PISM_DO $cmd
206 
207 echo
208 
209 cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
210  -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc -force_to_thickness_alpha 0.005 \
211  -o weak_force.nc -ts_file ts_weak_force.nc -ts_times -1000:yearly:0"
212 $PISM_DO $cmd
213 
214 
215 cmd="$PISM_MPIDO $NN $PISM -ys -1000.0 -ye 0 -skip 5 -i green_ssl2_110ka.nc -atmosphere searise_greenland \
216  -surface pdd,forcing -pdd_fausto -force_to_thickness_file green20km_y1.nc -force_to_thickness_alpha 0.05 \
217  -o strong_force.nc -ts_file ts_strong_force.nc -ts_times -1000:yearly:0"
218 $PISM_DO $cmd
219 
220 \endcode
221 The script also has a run with no forcing, one with forcing at a lower alpha value,
222 a factor of five smaller than the default, and one with a forcing at a higher alpha value, a factor of five higher.
223  */
225  const IceModelVec2S &ice_thickness,
226  const IceModelVec2CellType &cell_type,
227  IceModelVec2S &result) const {
228 
229  if (time < m_start_time) {
230  return;
231  }
232 
233  m_log->message(5,
234  " updating surface mass balance using -force_to_thickness mechanism ...");
235 
236  double ice_density = m_config->get_number("constants.ice.density");
237 
238  IceModelVec::AccessList list{&cell_type, &ice_thickness,
239  &m_target_thickness, &m_ftt_mask, &result};
240 
241  for (Points p(*m_grid); p; p.next()) {
242  const int i = p.i(), j = p.j();
243 
244  if (m_ftt_mask(i,j) > 0.5 and cell_type.grounded(i, j)) {
246  result(i,j) += ice_density * m_alpha * (m_target_thickness(i,j) - ice_thickness(i,j));
247  } else {
248  result(i,j) += ice_density * m_alpha * m_alpha_ice_free_factor * (m_target_thickness(i,j) - ice_thickness(i,j));
249  }
250  }
251  }
252  // no communication needed
253 }
254 
255 void ForceThickness::update_impl(const Geometry &geometry, double t, double dt) {
256  m_input_model->update(geometry, t, dt);
257 
258  m_mass_flux->copy_from(m_input_model->mass_flux());
259 
261  geometry.ice_thickness,
262  geometry.cell_type,
263  *m_mass_flux);
264 
268 }
269 
271  return *m_mass_flux;
272 }
273 
275  return *m_accumulation;
276 }
277 
279  return *m_melt;
280 }
281 
283  return *m_runoff;
284 }
285 
286 /*!
287 The timestep restriction is, by direct analogy, the same as for
288  \f[\frac{dy}{dt} = - \alpha y\f]
289 with explicit (forward) Euler. If \f$\Delta t\f$ is the time step then Euler is
290 \f$y_{n+1} = (1-\alpha \Delta t) y_n\f$. We require for stability that
291 \f$|y_{n+1}|\le |y_n|\f$, which is to say \f$|1-\alpha \Delta t|\le 1\f$.
292 Equivalently (since \f$\alpha \Delta t>0\f$),
293  \f[\alpha \Delta t\le 2\f]
294 Therefore we set here
295  \f[\Delta t = \frac{2}{\alpha}.\f]
296  */
298  double max_dt = 2.0 / m_alpha;
299  MaxTimestep input_max_dt = m_input_model->max_timestep(my_t);
300 
301  return std::min(input_max_dt, MaxTimestep(max_dt, "surface forcing"));
302 }
303 
305  m_ftt_mask.define(output);
306  m_target_thickness.define(output);
307 
308  if (m_input_model != NULL) {
309  m_input_model->define_model_state(output);
310  }
311 }
312 
313 void ForceThickness::write_model_state_impl(const File &output) const {
314  m_ftt_mask.write(output);
315  m_target_thickness.write(output);
316 
317  if (m_input_model != NULL) {
318  m_input_model->write_model_state(output);
319  }
320 }
321 
322 } // end of namespace surface
323 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
const units::System::Ptr m_sys
unit system used by this component
Definition: Component.hh:140
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
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
High-level PISM I/O class.
Definition: File.hh:51
IceModelVec2CellType cell_type
Definition: Geometry.hh:57
IceModelVec2S ice_thickness
Definition: Geometry.hh:53
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
bool grounded(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
void set_attrs(const std::string &pism_intent, const std::string &long_name, const std::string &units, const std::string &glaciological_units, const std::string &standard_name, unsigned int component)
Sets NetCDF attributes of an IceModelVec object.
Definition: iceModelVec.cc:399
void regrid(const std::string &filename, RegriddingFlag flag, double default_value=0.0)
Definition: iceModelVec.cc:838
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
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 set_time_independent(bool flag)
void set_name(const std::string &name)
void set_output_type(IO_Type type)
ForceThickness(IceGrid::ConstPtr g, std::shared_ptr< SurfaceModel > input)
const IceModelVec2S & runoff_impl() const
void write_model_state_impl(const File &output) const
The default (empty implementation).
MaxTimestep max_timestep_impl(double t) const
void define_model_state_impl(const File &output) const
The default (empty implementation).
void init_impl(const Geometry &geometry)
const IceModelVec2S & mass_flux_impl() const
void adjust_mass_flux(double time, const IceModelVec2S &ice_thickness, const IceModelVec2CellType &cell_type, IceModelVec2S &result) const
const IceModelVec2S & accumulation_impl() const
const IceModelVec2S & melt_impl() const
void update_impl(const Geometry &geometry, double t, double dt)
IceModelVec2S::Ptr m_mass_flux
static IceModelVec2S::Ptr allocate_mass_flux(IceGrid::ConstPtr grid)
Definition: SurfaceModel.cc:76
IceModelVec2S::Ptr m_melt
static IceModelVec2S::Ptr allocate_melt(IceGrid::ConstPtr grid)
IceModelVec2S::Ptr m_accumulation
std::shared_ptr< SurfaceModel > m_input_model
void dummy_accumulation(const IceModelVec2S &smb, IceModelVec2S &result)
static IceModelVec2S::Ptr allocate_accumulation(IceGrid::ConstPtr grid)
static IceModelVec2S::Ptr allocate_runoff(IceGrid::ConstPtr grid)
IceModelVec2S::Ptr m_runoff
void dummy_runoff(const IceModelVec2S &smb, IceModelVec2S &result)
void dummy_melt(const IceModelVec2S &smb, IceModelVec2S &result)
The interface of PISM's surface models.
Definition: SurfaceModel.hh:45
#define PISM_ERROR_LOCATION
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition: Units.cc:72
@ PISM_INT
Definition: IO_Flags.hh:36
static const double g
Definition: exactTestP.cc:39
@ CRITICAL
Definition: IO_Flags.hh:70
@ PISM_GUESS
Definition: IO_Flags.hh:41
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:49
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49