PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
LingleClark.cc
Go to the documentation of this file.
1 // Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2017, 2018, 2019, 2020, 2021, 2023 Constantine Khroulev
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/earth/LingleClark.hh"
20 
21 #include "pism/util/io/File.hh"
22 #include "pism/util/Time.hh"
23 #include "pism/util/Grid.hh"
24 #include "pism/util/ConfigInterface.hh"
25 #include "pism/util/error_handling.hh"
26 #include "pism/util/Vars.hh"
27 #include "pism/util/MaxTimestep.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/util/fftw_utilities.hh"
30 #include "pism/earth/LingleClarkSerial.hh"
31 #include "pism/util/Context.hh"
32 #include <memory>
33 
34 namespace pism {
35 namespace bed {
36 
37 LingleClark::LingleClark(std::shared_ptr<const Grid> grid)
38  : BedDef(grid),
39  m_total_displacement(m_grid, "bed_displacement"),
40  m_relief(m_grid, "bed_relief"),
41  m_load_thickness(grid, "load_thickness"),
42  m_elastic_displacement(grid, "elastic_bed_displacement") {
43 
44  m_time_name = m_config->get_string("time.dimension_name") + "_lingle_clark";
45  m_t_last = time().current();
46  m_update_interval = m_config->get_number("bed_deformation.lc.update_interval", "seconds");
47  m_t_eps = m_config->get_number("time_stepping.resolution", "seconds");
48 
49  if (m_update_interval < 1.0) {
51  "invalid bed_deformation.lc.update_interval = %f seconds",
53  }
54 
55  // A work vector. This storage is used to put thickness change on rank 0 and to get the plate
56  // displacement change back.
58  .long_name(
59  "total (viscous and elastic) displacement in the Lingle-Clark bed deformation model")
60  .units("meters");
61 
63 
65  .long_name("bed relief relative to the modeled bed displacement")
66  .units("meters");
67 
68  bool use_elastic_model = m_config->get_flag("bed_deformation.lc.elastic_model");
69 
71  .long_name(
72  "elastic part of the displacement in the Lingle-Clark bed deformation model; see :cite:`BLKfastearth`")
73  .units("meters");
75 
76  const int
77  Mx = m_grid->Mx(),
78  My = m_grid->My(),
79  Z = m_config->get_number("bed_deformation.lc.grid_size_factor"),
80  Nx = Z*(Mx - 1) + 1,
81  Ny = Z*(My - 1) + 1;
82 
83  const double
84  Lx = Z * (m_grid->x0() - m_grid->x(0)),
85  Ly = Z * (m_grid->y0() - m_grid->y(0));
86 
87  m_extended_grid = Grid::Shallow(m_grid->ctx(), Lx, Ly, m_grid->x0(), m_grid->y0(), Nx, Ny,
89 
91  std::make_shared<array::Scalar>(m_extended_grid, "viscous_bed_displacement");
92  m_viscous_displacement->metadata(0)
93  .long_name(
94  "bed displacement in the viscous half-space bed deformation model; see BuelerLingleBrown")
95  .units("meters");
96 
97  // coordinate variables of the extended grid should have different names
98  m_viscous_displacement->metadata().x().set_name("x_lc");
99  m_viscous_displacement->metadata().y().set_name("y_lc");
100 
101  // do not point to auxiliary coordinates "lon" and "lat".
102  m_viscous_displacement->metadata()["coordinates"] = "";
103 
104  m_viscous_displacement0 = m_viscous_displacement->allocate_proc0_copy();
105 
106  ParallelSection rank0(m_grid->com);
107  try {
108  if (m_grid->rank() == 0) {
109  m_serial_model.reset(new LingleClarkSerial(m_log, *m_config, use_elastic_model,
110  Mx, My,
111  m_grid->dx(), m_grid->dy(),
112  Nx, Ny));
113  }
114  } catch (...) {
115  rank0.failed();
116  }
117  rank0.check();
118 }
119 
121  // empty
122 }
123 
124 /*!
125  * Initialize the model by computing the viscous bed displacement using uplift and the elastic
126  * response using ice thickness.
127  *
128  * Then compute the bed relief as the difference between bed elevation and total bed displacement.
129  *
130  * This method has to initialize m_viscous_displacement, m_elastic_displacement,
131  * m_total_displacement, and m_relief.
132  */
133 void LingleClark::bootstrap_impl(const array::Scalar &bed_elevation,
134  const array::Scalar &bed_uplift,
135  const array::Scalar &ice_thickness,
136  const array::Scalar &sea_level_elevation) {
137  m_t_last = time().current();
138 
140 
141  compute_load(bed_elevation, ice_thickness, sea_level_elevation,
143 
144  std::shared_ptr<petsc::Vec> thickness0 = m_load_thickness.allocate_proc0_copy();
145 
146  // initialize the plate displacement
147  {
148  bed_uplift.put_on_proc0(*m_work0);
149  m_load_thickness.put_on_proc0(*thickness0);
150 
151  ParallelSection rank0(m_grid->com);
152  try {
153  if (m_grid->rank() == 0) {
154  PetscErrorCode ierr = 0;
155 
156  m_serial_model->bootstrap(*thickness0, *m_work0);
157 
158  ierr = VecCopy(m_serial_model->total_displacement(), *m_work0);
159  PISM_CHK(ierr, "VecCopy");
160 
161  ierr = VecCopy(m_serial_model->viscous_displacement(), *m_viscous_displacement0);
162  PISM_CHK(ierr, "VecCopy");
163 
164  ierr = VecCopy(m_serial_model->elastic_displacement(), *m_elastic_displacement0);
165  PISM_CHK(ierr, "VecCopy");
166  }
167  } catch (...) {
168  rank0.failed();
169  }
170  rank0.check();
171  }
172 
174 
176 
178 
179  // compute bed relief
181 }
182 
183 /*!
184  * Return the load response matrix for the elastic response.
185  *
186  * This method is used for testing only.
187  */
188 std::shared_ptr<array::Scalar> LingleClark::elastic_load_response_matrix() const {
189  std::shared_ptr<array::Scalar> result(new array::Scalar(m_extended_grid, "lrm"));
190 
191  int
192  Nx = m_extended_grid->Mx(),
193  Ny = m_extended_grid->My();
194 
195  auto lrm0 = result->allocate_proc0_copy();
196 
197  {
198  ParallelSection rank0(m_grid->com);
199  try {
200  if (m_grid->rank() == 0) {
201  std::vector<std::complex<double> > array(Nx * Ny);
202 
203  m_serial_model->compute_load_response_matrix((fftw_complex*)array.data());
204 
205  get_real_part((fftw_complex*)array.data(), 1.0, Nx, Ny, Nx, Ny, 0, 0, *lrm0);
206  }
207  } catch (...) {
208  rank0.failed();
209  }
210  rank0.check();
211  }
212 
213  result->get_from_proc0(*lrm0);
214 
215  return result;
216 }
217 
218 /*! Initialize the Lingle-Clark bed deformation model.
219  *
220  * Inputs:
221  *
222  * - bed topography,
223  * - ice thickness,
224  * - plate displacement (either read from a file or bootstrapped using uplift) and
225  * possibly re-gridded.
226  */
227 void LingleClark::init_impl(const InputOptions &opts, const array::Scalar &ice_thickness,
228  const array::Scalar &sea_level_elevation) {
229  m_log->message(2, "* Initializing the Lingle-Clark bed deformation model...\n");
230 
231  if (opts.type == INIT_RESTART or opts.type == INIT_BOOTSTRAP) {
232  File input_file(m_grid->com, opts.filename, io::PISM_NETCDF3, io::PISM_READONLY);
233 
234  if (input_file.find_variable(m_time_name)) {
235  input_file.read_variable(m_time_name, {0}, {1}, &m_t_last);
236  } else {
237  m_t_last = time().current();
238  }
239  } else {
240  m_t_last = time().current();
241  }
242 
243  // Initialize bed topography and uplift maps.
244  BedDef::init_impl(opts, ice_thickness, sea_level_elevation);
245 
247 
248  if (opts.type == INIT_RESTART) {
249  // Set viscous displacement by reading from the input file.
250  m_viscous_displacement->read(opts.filename, opts.record);
251  // Set elastic displacement by reading from the input file.
253  } else if (opts.type == INIT_BOOTSTRAP) {
254  this->bootstrap(m_topg, m_uplift, ice_thickness, sea_level_elevation);
255  } else {
256  // do nothing
257  }
258 
259  // Try re-gridding plate_displacement.
260  regrid("Lingle-Clark bed deformation model",
262  regrid("Lingle-Clark bed deformation model",
264 
265  compute_load(m_topg, ice_thickness, sea_level_elevation,
267 
268  // Now that viscous displacement and elastic displacement are finally initialized,
269  // put them on rank 0 and initialize the serial model itself.
270  {
273 
274  ParallelSection rank0(m_grid->com);
275  try {
276  if (m_grid->rank() == 0) { // only processor zero does the work
277  PetscErrorCode ierr = 0;
278 
280 
281  ierr = VecCopy(m_serial_model->total_displacement(), *m_work0);
282  PISM_CHK(ierr, "VecCopy");
283  }
284  } catch (...) {
285  rank0.failed();
286  }
287  rank0.check();
288  }
289 
291 
292  // compute bed relief
294 }
295 
297 
298  if (t < m_t_last) {
300  "time %f is less than the previous time %f",
301  t, m_t_last);
302  }
303 
304  // Find the smallest time of the form m_t_last + k * m_update_interval that is greater
305  // than t
306  double k = ceil((t - m_t_last) / m_update_interval);
307 
308  double
309  t_next = m_t_last + k * m_update_interval,
310  dt_max = t_next - t;
311 
312  if (dt_max < m_t_eps) {
313  dt_max = m_update_interval;
314  }
315 
316  return MaxTimestep(dt_max, "bed_def lc");
317 }
318 
319 /*!
320  * Get total bed displacement on the PISM grid.
321  */
323  return m_total_displacement;
324 }
325 
327  return *m_viscous_displacement;
328 }
329 
331  return m_elastic_displacement;
332 }
333 
335  return m_relief;
336 }
337 
338 void LingleClark::step(const array::Scalar &ice_thickness,
339  const array::Scalar &sea_level_elevation,
340  double dt) {
341 
342  compute_load(m_topg, ice_thickness, sea_level_elevation,
344 
346 
347  ParallelSection rank0(m_grid->com);
348  try {
349  if (m_grid->rank() == 0) { // only processor zero does the step
350  PetscErrorCode ierr = 0;
351 
352  m_serial_model->step(dt, *m_work0);
353 
354  ierr = VecCopy(m_serial_model->total_displacement(), *m_work0);
355  PISM_CHK(ierr, "VecCopy");
356 
357  ierr = VecCopy(m_serial_model->viscous_displacement(), *m_viscous_displacement0);
358  PISM_CHK(ierr, "VecCopy");
359 
360  ierr = VecCopy(m_serial_model->elastic_displacement(), *m_elastic_displacement0);
361  PISM_CHK(ierr, "VecCopy");
362  }
363  } catch (...) {
364  rank0.failed();
365  }
366  rank0.check();
367 
369 
371 
373 
374  // Update bed elevation using bed displacement and relief.
375  {
377  // Increment the topg state counter. SIAFD relies on this!
379  }
380 
381  //! Finally, we need to update bed uplift and topg_last.
384 }
385 
386 //! Update the Lingle-Clark bed deformation model.
387 void LingleClark::update_impl(const array::Scalar &ice_thickness,
388  const array::Scalar &sea_level_elevation,
389  double t, double dt) {
390 
391  double
392  t_next = m_t_last + m_update_interval,
393  t_final = t + dt;
394 
395  if (t_final < m_t_last) {
396  throw RuntimeError(PISM_ERROR_LOCATION, "cannot go back in time");
397  }
398 
399  if (std::abs(t_next - t_final) < m_t_eps) { // reached the next update time
400  double dt_beddef = t_final - m_t_last;
401  step(ice_thickness, sea_level_elevation, dt_beddef);
402  m_t_last = t_final;
403  }
404 }
405 
406 void LingleClark::define_model_state_impl(const File &output) const {
408  m_viscous_displacement->define(output, io::PISM_DOUBLE);
410 
411  if (not output.find_variable(m_time_name)) {
413 
414  output.write_attribute(m_time_name, "long_name",
415  "time of the last update of the Lingle-Clark bed deformation model");
416  output.write_attribute(m_time_name, "calendar", time().calendar());
417  output.write_attribute(m_time_name, "units", time().units_string());
418  }
419 }
420 
421 void LingleClark::write_model_state_impl(const File &output) const {
423 
424  m_viscous_displacement->write(output);
426 
427  output.write_variable(m_time_name, {0}, {1}, &m_t_last);
428 }
429 
431  DiagnosticList result = {
432  {"viscous_bed_displacement", Diagnostic::wrap(*m_viscous_displacement)},
433  {"elastic_bed_displacement", Diagnostic::wrap(m_elastic_displacement)},
434  };
435  return combine(result, BedDef::diagnostics_impl());
436 }
437 
438 } // end of namespace bed
439 } // end of namespace pism
const Time & time() const
Definition: Component.cc:109
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
static Ptr wrap(const T &input)
Definition: Diagnostic.hh:160
void read_variable(const std::string &variable_name, const std::vector< unsigned int > &start, const std::vector< unsigned int > &count, double *ip) const
Definition: File.cc:747
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
void write_variable(const std::string &variable_name, const std::vector< unsigned int > &start, const std::vector< unsigned int > &count, const double *op) const
Definition: File.cc:760
void define_variable(const std::string &name, io::Type nctype, const std::vector< std::string > &dims) const
Define a variable.
Definition: File.cc:573
void write_attribute(const std::string &var_name, const std::string &att_name, io::Type nctype, const std::vector< double > &values) const
Write a multiple-valued double attribute.
Definition: File.cc:638
High-level PISM I/O class.
Definition: File.hh:56
static std::shared_ptr< Grid > Shallow(std::shared_ptr< const Context > ctx, double Lx, double Ly, double x0, double y0, unsigned int Mx, unsigned int My, grid::Registration r, grid::Periodicity p)
Initialize a uniform, shallow (3 z-levels) grid with half-widths (Lx,Ly) and Mx by My nodes.
Definition: Grid.cc:134
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
double current() const
Current time, in seconds.
Definition: Time.cc:465
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
void copy_from(const Array2D< T > &source)
Definition: Array2D.hh:73
void add(double alpha, const Array2D< T > &x)
Definition: Array2D.hh:65
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 get_from_proc0(petsc::Vec &onp0)
Gets a local Array2 from processor 0.
Definition: Array.cc:1095
void write(const std::string &filename) const
Definition: Array.cc:800
void inc_state_counter()
Increment the object state counter.
Definition: Array.cc:151
std::shared_ptr< petsc::Vec > allocate_proc0_copy() const
Definition: Array.cc:963
void put_on_proc0(petsc::Vec &onp0) const
Puts a local array::Scalar on processor 0.
Definition: Array.cc:1052
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: Array.cc:553
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: BedDef.cc:64
virtual DiagnosticList diagnostics_impl() const
Definition: BedDef.cc:69
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: BedDef.cc:59
array::Scalar2 m_topg_last
bed elevation at the time of the last update
Definition: BedDef.hh:83
void bootstrap(const array::Scalar &bed_elevation, const array::Scalar &bed_uplift, const array::Scalar &ice_thickness, const array::Scalar &sea_level_elevation)
Initialize using provided bed elevation and uplift.
Definition: BedDef.cc:82
array::Scalar2 m_topg
current bed elevation
Definition: BedDef.hh:80
virtual void init_impl(const InputOptions &opts, const array::Scalar &ice_thickness, const array::Scalar &sea_level_elevation)
Initialize from the context (input file and the "variables" database).
Definition: BedDef.cc:105
array::Scalar m_uplift
bed uplift rate
Definition: BedDef.hh:86
const array::Scalar & bed_elevation() const
Definition: BedDef.cc:51
void compute_uplift(const array::Scalar &bed, const array::Scalar &bed_last, double dt, array::Scalar &result)
Compute bed uplift (dt is in seconds).
Definition: BedDef.cc:167
PISM bed deformation model (base class).
Definition: BedDef.hh:39
Class implementing the bed deformation model described in [BLKfastearth].
void step(const array::Scalar &ice_thickness, const array::Scalar &sea_level_elevation, double dt)
Definition: LingleClark.cc:338
std::shared_ptr< petsc::Vec > m_elastic_displacement0
rank 0 storage for the elastic displacement
Definition: LingleClark.hh:95
MaxTimestep max_timestep_impl(double t) const
Definition: LingleClark.cc:296
void bootstrap_impl(const array::Scalar &bed_elevation, const array::Scalar &bed_uplift, const array::Scalar &ice_thickness, const array::Scalar &sea_level_elevation)
Definition: LingleClark.cc:133
void init_impl(const InputOptions &opts, const array::Scalar &ice_thickness, const array::Scalar &sea_level_elevation)
Definition: LingleClark.cc:227
double m_t_eps
Temporal resolution to use when checking whether it's time to update.
Definition: LingleClark.hh:102
std::shared_ptr< petsc::Vec > m_work0
Definition: LingleClark.hh:73
array::Scalar m_elastic_displacement
Elastic bed displacement (part of the model state)
Definition: LingleClark.hh:93
const array::Scalar & elastic_displacement() const
Definition: LingleClark.cc:330
std::shared_ptr< Grid > m_extended_grid
extended grid for the viscous plate displacement
Definition: LingleClark.hh:85
std::string m_time_name
Name of the variable used to store the last update time.
Definition: LingleClark.hh:104
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: LingleClark.cc:421
std::shared_ptr< array::Scalar > m_viscous_displacement
Viscous displacement on the extended grid (part of the model state).
Definition: LingleClark.hh:88
double m_t_last
time of the last bed deformation update
Definition: LingleClark.hh:98
array::Scalar m_relief
Bed relief relative to the bed displacement.
Definition: LingleClark.hh:76
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: LingleClark.cc:406
std::shared_ptr< petsc::Vec > m_viscous_displacement0
rank 0 storage using the extended grid
Definition: LingleClark.hh:90
const array::Scalar & relief() const
Definition: LingleClark.cc:334
LingleClark(std::shared_ptr< const Grid > g)
Definition: LingleClark.cc:37
DiagnosticList diagnostics_impl() const
Definition: LingleClark.cc:430
std::unique_ptr< LingleClarkSerial > m_serial_model
Serial viscoelastic bed deformation model.
Definition: LingleClark.hh:82
const array::Scalar & total_displacement() const
Definition: LingleClark.cc:322
array::Scalar m_total_displacement
Total (viscous and elastic) bed displacement.
Definition: LingleClark.hh:69
const array::Scalar & viscous_displacement() const
Definition: LingleClark.cc:326
void update_impl(const array::Scalar &ice_thickness, const array::Scalar &sea_level_elevation, double t, double dt)
Update the Lingle-Clark bed deformation model.
Definition: LingleClark.cc:387
double m_update_interval
Update interval in seconds.
Definition: LingleClark.hh:100
std::shared_ptr< array::Scalar > elastic_load_response_matrix() const
Definition: LingleClark.cc:188
array::Scalar m_load_thickness
Ice-equivalent load thickness.
Definition: LingleClark.hh:79
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
double compute_load(double bed, double ice_thickness, double sea_level, double ice_density, double ocean_density)
Definition: BedDef.cc:174
@ NOT_PERIODIC
Definition: Grid.hh:51
@ CELL_CORNER
Definition: Grid.hh:53
@ PISM_NETCDF3
Definition: IO_Flags.hh:57
@ PISM_READONLY
open an existing file for reading only
Definition: IO_Flags.hh:72
@ PISM_DOUBLE
Definition: IO_Flags.hh:52
@ INIT_BOOTSTRAP
Definition: Component.hh:56
@ INIT_RESTART
Definition: Component.hh:56
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:125
static const double k
Definition: exactTestP.cc:42
static std::string calendar(const File *input_file, const Config &config, const Logger &log)
Definition: Time.cc:146
T combine(const T &a, const T &b)
void get_real_part(fftw_complex *input, double normalization, int Mx, int My, int Nx, int Ny, int i0, int j0, petsc::Vec &output)
Get the real part of input and put it in output.
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