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