PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
BedDef.cc
Go to the documentation of this file.
1 // Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 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 "BedDef.hh"
20 #include "pism/util/pism_utilities.hh"
21 #include "pism/util/Time.hh"
22 #include "pism/util/Vars.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/ConfigInterface.hh"
25 #include "pism/util/Context.hh"
26 
27 namespace pism {
28 namespace bed {
29 
31  : Component(grid),
32  m_wide_stencil(m_config->get_number("grid.max_stencil_width")),
33  m_topg(m_grid, "topg", WITH_GHOSTS, m_wide_stencil),
34  m_topg_last(m_grid, "topg", WITH_GHOSTS, m_wide_stencil),
35  m_uplift(m_grid, "dbdt", WITHOUT_GHOSTS)
36 {
37 
38  m_topg.set_attrs("model_state", "bedrock surface elevation",
39  "m", "m", "bedrock_altitude", 0);
40 
41  m_topg_last.set_attrs("model_state", "bedrock surface elevation",
42  "m", "m", "bedrock_altitude", 0);
43 
44  m_uplift.set_attrs("model_state", "bedrock uplift rate",
45  "m s-1", "mm year-1", "tendency_of_bedrock_altitude", 0);
46 }
47 
49  return m_topg;
50 }
51 
52 const IceModelVec2S& BedDef::uplift() const {
53  return m_uplift;
54 }
55 
56 void BedDef::define_model_state_impl(const File &output) const {
57  m_uplift.define(output);
58  m_topg.define(output);
59 }
60 
61 void BedDef::write_model_state_impl(const File &output) const {
62  m_uplift.write(output);
63  m_topg.write(output);
64 }
65 
67  DiagnosticList result;
68  result = {
69  {"dbdt", Diagnostic::wrap(m_uplift)},
70  {"topg", Diagnostic::wrap(m_topg)}
71  };
72 
73  return result;
74 }
75 
76 void BedDef::init(const InputOptions &opts, const IceModelVec2S &ice_thickness,
77  const IceModelVec2S &sea_level_elevation) {
78  this->init_impl(opts, ice_thickness, sea_level_elevation);
79 }
80 
81 //! Initialize using provided bed elevation and uplift.
82 void BedDef::bootstrap(const IceModelVec2S &bed_elevation,
83  const IceModelVec2S &bed_uplift,
84  const IceModelVec2S &ice_thickness,
85  const IceModelVec2S &sea_level_elevation) {
86  this->bootstrap_impl(bed_elevation, bed_uplift, ice_thickness, sea_level_elevation);
87 }
88 
89 void BedDef::bootstrap_impl(const IceModelVec2S &bed_elevation,
90  const IceModelVec2S &bed_uplift,
91  const IceModelVec2S &ice_thickness,
92  const IceModelVec2S &sea_level_elevation) {
94  m_uplift.copy_from(bed_uplift);
95 
96  // suppress a compiler warning:
97  (void) ice_thickness;
98  (void) sea_level_elevation;
99 }
100 
101 void BedDef::update(const IceModelVec2S &ice_thickness,
102  const IceModelVec2S &sea_level_elevation,
103  double t, double dt) {
104  this->update_impl(ice_thickness, sea_level_elevation, t, dt);
105 }
106 
107 //! Initialize from the context (input file and the "variables" database).
108 void BedDef::init_impl(const InputOptions &opts, const IceModelVec2S &ice_thickness,
109  const IceModelVec2S &sea_level_elevation) {
110  (void) ice_thickness;
111  (void) sea_level_elevation;
112 
113  switch (opts.type) {
114  case INIT_RESTART:
115  // read bed elevation and uplift rate from file
116  m_log->message(2,
117  " reading bed topography and uplift from %s ... \n",
118  opts.filename.c_str());
119  // re-starting
120  m_topg.read(opts.filename, opts.record); // fails if not found!
121  m_uplift.read(opts.filename, opts.record); // fails if not found!
122  break;
123  case INIT_BOOTSTRAP:
124  // bootstrapping
126  m_config->get_number("bootstrapping.defaults.bed"));
128  m_config->get_number("bootstrapping.defaults.uplift"));
129  break;
130  case INIT_OTHER:
131  default:
132  {
133  // do nothing
134  }
135  }
136 
137  // process -regrid_file and -regrid_vars
138  regrid("bed deformation", m_topg);
139  // uplift is not a part of the model state, but the user may want to take it from a -regrid_file
140  // during bootstrapping
141  regrid("bed deformation", m_uplift);
142 
143  std::string uplift_file = m_config->get_string("bed_deformation.bed_uplift_file");
144  if (not uplift_file.empty()) {
145  m_log->message(2,
146  " reading bed uplift from %s ... \n",
147  uplift_file.c_str());
148  m_uplift.regrid(uplift_file, CRITICAL);
149  }
150 
151  std::string correction_file = m_config->get_string("bed_deformation.bed_topography_delta_file");
152  if (not correction_file.empty()) {
153  apply_topg_offset(correction_file);
154  }
155 
156  // this should be the last thing we do here
158 }
159 
160 /*!
161  * Apply a correction to the bed topography by reading topg_delta from filename.
162  */
163 void BedDef::apply_topg_offset(const std::string &filename) {
164  m_log->message(2, " Adding a bed topography correction read in from %s...\n",
165  filename.c_str());
166 
167  IceModelVec2S topg_delta(m_grid, "topg_delta", WITHOUT_GHOSTS);
168  topg_delta.set_attrs("internal", "bed topography correction",
169  "meters", "meters", "", 0);
170 
171  topg_delta.regrid(filename, CRITICAL);
172 
173  m_topg.add(1.0, topg_delta);
174 }
175 
176 //! Compute bed uplift (dt is in seconds).
177 void BedDef::compute_uplift(const IceModelVec2S &bed, const IceModelVec2S &bed_last,
178  double dt, IceModelVec2S &result) {
179  bed.add(-1, bed_last, result);
180  //! uplift = (topg - topg_last) / dt
181  result.scale(1.0 / dt);
182 }
183 
184 double compute_load(double bed, double ice_thickness, double sea_level,
185  double ice_density, double ocean_density) {
186 
187  double
188  ice_load = ice_thickness,
189  ocean_depth = std::max(sea_level - bed, 0.0),
190  ocean_load = (ocean_density / ice_density) * ocean_depth;
191 
192  // this excludes the load of ice shelves
193  return ice_load > ocean_load ? ice_load : 0.0;
194 }
195 
196 /*! Compute the load on the bedrock in units of ice-equivalent thickness.
197  *
198  */
199 void compute_load(const IceModelVec2S &bed_elevation,
200  const IceModelVec2S &ice_thickness,
201  const IceModelVec2S &sea_level_elevation,
202  IceModelVec2S &result) {
203 
204  Config::ConstPtr config = result.grid()->ctx()->config();
205 
206  const double
207  ice_density = config->get_number("constants.ice.density"),
208  ocean_density = config->get_number("constants.sea_water.density");
209 
210  IceModelVec::AccessList list{&bed_elevation, &ice_thickness, &sea_level_elevation, &result};
211 
212  for (Points p(*result.grid()); p; p.next()) {
213  const int i = p.i(), j = p.j();
214 
215  result(i, j) = compute_load(bed_elevation(i, j),
216  ice_thickness(i, j),
217  sea_level_elevation(i, j),
218  ice_density, ocean_density);
219  }
220 }
221 
222 } // end of namespace bed
223 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
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
virtual void regrid(const std::string &module_name, IceModelVec &variable, RegriddingFlag flag=NO_REGRID_WITHOUT_REGRID_VARS)
Definition: Component.cc:151
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
A class defining a common interface for most PISM sub-models.
Definition: Component.hh:101
std::shared_ptr< const Config > ConstPtr
static Ptr wrap(const T &input)
Definition: Diagnostic.hh:155
High-level PISM I/O class.
Definition: File.hh:51
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
void add(double alpha, const IceModelVec2S &x)
void copy_from(const IceModelVec2S &source)
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 scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
Definition: iceModelVec.cc:252
void read(const std::string &filename, unsigned int time)
Definition: iceModelVec.cc:833
IceGrid::ConstPtr grid() const
Definition: iceModelVec.cc:128
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
virtual void init_impl(const InputOptions &opts, const IceModelVec2S &ice_thickness, const IceModelVec2S &sea_level_elevation)
Initialize from the context (input file and the "variables" database).
Definition: BedDef.cc:108
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
Definition: BedDef.cc:61
virtual DiagnosticList diagnostics_impl() const
Definition: BedDef.cc:66
void update(const IceModelVec2S &ice_thickness, const IceModelVec2S &sea_level_elevation, double t, double dt)
Definition: BedDef.cc:101
const IceModelVec2S & uplift() const
Definition: BedDef.cc:52
virtual void update_impl(const IceModelVec2S &ice_thickness, const IceModelVec2S &sea_level_elevation, double t, double dt)=0
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
Definition: BedDef.cc:56
IceModelVec2S m_topg
current bed elevation
Definition: BedDef.hh:81
void init(const InputOptions &opts, const IceModelVec2S &ice_thickness, const IceModelVec2S &sea_level_elevation)
Definition: BedDef.cc:76
BedDef(IceGrid::ConstPtr g)
Definition: BedDef.cc:30
virtual void bootstrap_impl(const IceModelVec2S &bed_elevation, const IceModelVec2S &bed_uplift, const IceModelVec2S &ice_thickness, const IceModelVec2S &sea_level_elevation)
Definition: BedDef.cc:89
void compute_uplift(const IceModelVec2S &bed, const IceModelVec2S &bed_last, double dt, IceModelVec2S &result)
Compute bed uplift (dt is in seconds).
Definition: BedDef.cc:177
IceModelVec2S m_uplift
bed uplift rate
Definition: BedDef.hh:87
IceModelVec2S m_topg_last
bed elevation at the time of the last update
Definition: BedDef.hh:84
virtual void apply_topg_offset(const std::string &filename)
Definition: BedDef.cc:163
const IceModelVec2S & bed_elevation() const
Definition: BedDef.cc:48
void bootstrap(const IceModelVec2S &bed_elevation, const IceModelVec2S &bed_uplift, const IceModelVec2S &ice_thickness, const IceModelVec2S &sea_level_elevation)
Initialize using provided bed elevation and uplift.
Definition: BedDef.cc:82
double compute_load(double bed, double ice_thickness, double sea_level, double ice_density, double ocean_density)
Definition: BedDef.cc:184
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
@ INIT_BOOTSTRAP
Definition: Component.hh:39
@ INIT_OTHER
Definition: Component.hh:39
@ INIT_RESTART
Definition: Component.hh:39
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:117
@ OPTIONAL
Definition: IO_Flags.hh:70
@ CRITICAL
Definition: IO_Flags.hh:70
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
@ WITH_GHOSTS
Definition: iceModelVec.hh:49
InitializationType type
initialization type
Definition: Component.hh:44
std::string filename
name of the input file (if applicable)
Definition: Component.hh:46
unsigned int record
index of the record to re-start from
Definition: Component.hh:48