PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
PointwiseIsostasy.cc
Go to the documentation of this file.
1 // Copyright (C) 2010, 2011, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2022, 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/BedDef.hh"
20 #include "pism/util/Grid.hh"
21 #include "pism/util/Time.hh"
22 #include "pism/util/ConfigInterface.hh"
23 #include "pism/util/Vars.hh"
24 #include "pism/util/MaxTimestep.hh"
25 
26 namespace pism {
27 namespace bed {
28 
29 PointwiseIsostasy::PointwiseIsostasy(std::shared_ptr<const Grid> grid)
30  : BedDef(grid), m_load_last(m_grid, "load_last") {
31  // empty
32 }
33 
34 void PointwiseIsostasy::init_impl(const InputOptions &opts, const array::Scalar &ice_thickness,
35  const array::Scalar &sea_level_elevation) {
36 
37  m_log->message(2,
38  "* Initializing the pointwise isostasy bed deformation model...\n");
39 
40  BedDef::init_impl(opts, ice_thickness, sea_level_elevation);
41 
42  // store the initial load
43  compute_load(m_topg, ice_thickness, sea_level_elevation, m_load_last);
44 }
45 
46 
48  const array::Scalar &bed_uplift,
49  const array::Scalar &ice_thickness,
50  const array::Scalar &sea_level_elevation) {
51  BedDef::bootstrap_impl(bed_elevation, bed_uplift, ice_thickness, sea_level_elevation);
52 
53  // store initial load and bed elevation
54  compute_load(bed_elevation, ice_thickness, sea_level_elevation, m_load_last);
56 }
57 
59  (void) t;
60  return MaxTimestep("bed_def iso");
61 }
62 
63 //! Updates the pointwise isostasy model.
65  const array::Scalar &sea_level_elevation,
66  double t, double dt) {
67  (void) t;
68 
69  const double
70  mantle_density = m_config->get_number("bed_deformation.mantle_density"),
71  load_density = m_config->get_number("constants.ice.density"),
72  ocean_density = m_config->get_number("constants.sea_water.density"),
73  ice_density = m_config->get_number("constants.ice.density"),
74  f = load_density / mantle_density;
75 
76  //! Our goal: topg = topg_last - f*(load - load_last)
77 
79  &ice_thickness, &sea_level_elevation, &m_load_last};
80 
81  ParallelSection loop(m_grid->com);
82  try {
83  for (auto p = m_grid->points(); p; p.next()) {
84  const int i = p.i(), j = p.j();
85 
86  double load = compute_load(m_topg(i, j),
87  ice_thickness(i, j),
88  sea_level_elevation(i, j),
89  ice_density, ocean_density);
90 
91  m_topg(i, j) = m_topg_last(i, j) - f * (load - m_load_last(i, j));
92  m_load_last(i, j) = load;
93  }
94  } catch (...) {
95  loop.failed();
96  }
97  loop.check();
98 
99  //! Finally, we need to update bed uplift, topg_last and load_last.
101 
103 
104  //! Increment the topg state counter. SIAFD relies on this!
106 }
107 
108 } // end of namespace bed
109 } // 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
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.
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
void inc_state_counter()
Increment the object state counter.
Definition: Array.cc:151
virtual 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: BedDef.cc:88
array::Scalar2 m_topg_last
bed elevation at the time of the last update
Definition: BedDef.hh:83
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
void bootstrap_impl(const array::Scalar &bed_elevation, const array::Scalar &bed_uplift, const array::Scalar &ice_thickness, const array::Scalar &sea_level_elevation)
array::Scalar m_load_last
last ice load (ice-equivalent thickness)
Definition: BedDef.hh:121
MaxTimestep max_timestep_impl(double t) const
PointwiseIsostasy(std::shared_ptr< const Grid > g)
void update_impl(const array::Scalar &ice_thickness, const array::Scalar &sea_level_elevation, double t, double dt)
Updates the pointwise isostasy model.
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).
double compute_load(double bed, double ice_thickness, double sea_level, double ice_density, double ocean_density)
Definition: BedDef.cc:174