PISM, A Parallel Ice Sheet Model  stable v2.0.5 committed by Constantine Khrulev on 2022-10-14 09:56:26 -0800
OptTillphiYieldStress.cc
Go to the documentation of this file.
1 // Copyright (C) 2004--2021 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 "OptTillphiYieldStress.hh"
20 
21 #include "pism/geometry/Geometry.hh"
22 #include "pism/util/Context.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/IceModelVec2CellType.hh"
25 #include "pism/util/MaxTimestep.hh"
26 #include "pism/util/Time.hh"
27 #include "pism/util/error_handling.hh"
28 #include "pism/util/io/File.hh"
29 #include "pism/util/pism_utilities.hh"
30 
31 namespace pism {
32 
33 /*! Optimization of till friction angle for given target surface elevation, analogous to
34  Pollard et al. (2012), TC 6(5), "A simple inverse method for the distribution of basal
35  sliding coefficients under ice sheets, applied to Antarctica"
36 */
38  : MohrCoulombYieldStress(grid),
39  m_mask(m_grid, "diff_mask", WITH_GHOSTS),
40  m_usurf_difference(m_grid, "usurf_difference", WITH_GHOSTS),
41  m_usurf_target(m_grid, "usurf", WITH_GHOSTS)
42 {
43  // In this model tillphi is NOT time-independent.
45 
46  m_name = "Iterative optimization of the till friction angle for the Mohr-Coulomb yield stress model";
47 
49  "target surface elevation for tillphi optimization",
50  "m", "m", "" /* no standard name */, 0);
52 
53  m_usurf_difference.set_attrs("diagnostic",
54  "difference between modeled and target"
55  " surface elevations",
56  "m", "m", "" /* no standard name */, 0);
58 
59  m_mask.set_attrs("diagnostic",
60  "one if the till friction angle was"
61  " updated by the last iteration, zero otherwise ",
62  "", "", // no units
63  "" /* no standard name */, 0);
64  m_mask.metadata()["flag_values"] = {0.0, 1.0};
65  m_mask.metadata()["flag_meanings"] = "no_update updated_during_last_iteration";
66 
67  {
68  // convergence threshold
69  m_dhdt_min = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.dhdt_min", "m / s");
70 
71  // scale used to compute tillphi adjustment using the surface elevation mismatch:
72  m_dphi_scale = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.dphi_scale", "degree / m");
73  // upper and lower bounds of the tillphi adjustment during an iteration:
74  m_dphi_max = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.dphi_max");
75  m_dphi_min = -2 * m_dphi_max;
76  // lower bound of tillphi:
77  m_phi0_min = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.phi0_min");
78  m_phi0_max = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.phi0_max");
79  m_topg_min = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.topg_min");
80  m_topg_max = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.topg_max");
81  // upper bound of tillphi:
82  m_phi_max = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.phi_max");
83 
84  if (m_phi0_min >= m_phi_max) {
86  "basal_yield_stress.mohr_coulomb.tillphi_opt: phi0_min >= phi_max");
87  }
88 
89  if (m_topg_min >= m_topg_max) {
91  "basal_yield_stress.mohr_coulomb.tillphi_opt: topg_min >= topg_max");
92  }
93  }
94 
95  {
96  m_time_name = m_config->get_string("time.dimension_name") + "_tillphi_opt";
97  m_t_last = m_grid->ctx()->time()->current();
98  m_update_interval = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.dt", "seconds");
99  m_t_eps = m_config->get_number("time_stepping.resolution", "seconds");
100  }
101 
102  m_log->message(2,
103  " Using iterative optimization of the till friction angle.\n"
104  " Lower bound phi0 of the till friction angle is a piecewise-linear function of bed elevation (b):\n"
105  " / %5.2f for b < %.f\n"
106  " phi0 = | %5.2f + (b - (%.f)) * (%.2f / %.f) for %.f < b < %.f\n"
107  " \\ %5.2f for %.f < b\n",
111 }
112 
113 /*!
114  * Initialize the last time tillphi was updated.
115  */
116 void OptTillphiYieldStress::init_t_last(const File &input_file) {
117  if (input_file.find_variable(m_time_name)) {
118  input_file.read_variable(m_time_name, {0}, {1}, &m_t_last);
119  } else {
120  m_t_last = m_grid->ctx()->time()->current();
121  }
122 }
123 
124 /*!
125  * Initialize the target ice surface elevation.
126  */
128  auto filename = m_config->get_string("basal_yield_stress.mohr_coulomb.tillphi_opt.file");
129 
130  if (not filename.empty()) {
131  m_usurf_target.regrid(filename, CRITICAL);
132  } else {
133  m_log->message(2, "* No file set to read target surface elevation from... using '%s'\n",
134  input_file.filename().c_str());
135 
136  m_usurf_target.regrid(input_file, CRITICAL);
137  }
138 
139  m_usurf_target.metadata().set_name("usurf_target");
140 }
141 
142 void OptTillphiYieldStress::restart_impl(const File &input_file, int record) {
143 
144  MohrCoulombYieldStress::restart_impl(input_file, record);
145 
146  init_t_last(input_file);
147 
148  init_usurf_target(input_file);
149 }
150 
151 //! Initialize the pseudo-plastic till mechanical model.
152 //! Target surface elevation and initial iteration time are set.
154  const YieldStressInputs &inputs) {
155 
156  MohrCoulombYieldStress::bootstrap_impl(input_file, inputs);
157 
158  init_t_last(input_file);
159 
160  init_usurf_target(input_file);
161 }
162 
164  (void) inputs;
166  "not implemented: till friction angle optimization "
167  "cannot be initialized without an input file");
168 }
169 
171  MaxTimestep dt_max;
172  {
173  if (t < m_t_last) {
175  "time %f is less than the previous time %f",
176  t, m_t_last);
177  }
178 
179  // Find the smallest time of the form m_t_last + k * m_update_interval that is greater
180  // than t
181  double k = ceil((t - m_t_last) / m_update_interval);
182 
183  double
184  t_next = m_t_last + k * m_update_interval,
185  dt = t_next - t;
186 
187  if (dt < m_t_eps) {
188  dt = m_update_interval;
189  }
190 
191  dt_max = MaxTimestep(dt, "tillphi_opt");
192  }
193 
194  auto dt_mohr_coulomb = MohrCoulombYieldStress::max_timestep_impl(t);
195 
196  return std::min(dt_max, dt_mohr_coulomb);
197 }
198 
200  double t, double dt) {
201 
202  double
203  t_next = m_t_last + m_update_interval,
204  t_final = t + dt;
205 
206  if (t_final < m_t_last) {
207  throw RuntimeError(PISM_ERROR_LOCATION, "cannot go back in time");
208  }
209 
210  if (std::abs(t_next - t_final) < m_t_eps) { // reached the next update time
212  inputs.geometry->bed_elevation,
213  inputs.geometry->cell_type);
214  m_t_last = t_final;
215  }
216 
218 }
219 
220 //! Perform an iteration to adjust the till friction angle according to the difference
221 //! between target and modeled surface elevations.
222 void OptTillphiYieldStress::update_tillphi(const IceModelVec2S &ice_surface_elevation,
223  const IceModelVec2S &bed_topography,
224  const IceModelVec2CellType &cell_type) {
225 
226  m_log->message(2, "* Updating till friction angle...\n");
227 
229  { &m_till_phi, &m_usurf_target, &m_usurf_difference, &m_mask, &ice_surface_elevation, &bed_topography, &cell_type };
230 
231  m_mask.set(0.0);
232 
233  double slope = (m_phi0_max - m_phi0_min) / (m_topg_max - m_topg_min);
234 
235  for (Points p(*m_grid); p; p.next()) {
236  const int i = p.i(), j = p.j();
237 
238  // Compute the lower bound of the till friction angle (default value corresponds
239  // to "bed_topography(i, j) > topg_max"):
240  double phi0 = m_phi0_max;
241  if (bed_topography(i, j) <= m_topg_min) {
242  phi0 = m_phi0_min;
243  } else if (bed_topography(i, j) <= m_topg_max) {
244  phi0 = m_phi0_min + (bed_topography(i, j) - m_topg_min) * slope;
245  }
246 
247  if (cell_type.grounded_ice(i, j)) {
248  double dh_previous = m_usurf_difference(i, j);
249  m_usurf_difference(i, j) = m_usurf_target(i, j) - ice_surface_elevation(i, j);
250  double dh_change = std::abs(m_usurf_difference(i, j) - dh_previous);
251 
252  if (dh_change / m_update_interval > m_dhdt_min) {
253  // Update tillphi if the rate of change of the surface elevation mismatch since
254  // the last iteration is greater than the convergence threshold m_dhdt_min.
255 
256  // Mark this location as "updated by the last iteration":
257  m_mask(i, j) = 1.0;
258 
259  // Compute (and clip) the tillphi adjustment
260  double dphi = m_usurf_difference(i, j) * m_dphi_scale;
261  dphi = pism::clip(dphi, m_dphi_min, m_dphi_max);
262 
263  // Update (and clip) the till friction angle:
264  m_till_phi(i, j) += dphi;
265  m_till_phi(i, j) = pism::clip(m_till_phi(i, j), phi0, m_phi_max);
266  }
267  } else if (cell_type.ocean(i, j)) {
268  // Floating and ice free ocean: use the bed-elevation-dependent lower bound of
269  // tillphi:
270  m_till_phi(i, j) = phi0;
271  }
272  } // end of the loop over grid points
273 }
274 
277 
278  if (not output.find_variable(m_time_name)) {
280 
281  output.write_attribute(m_time_name, "long_name",
282  "time of the last update of the till friction angle");
283  output.write_attribute(m_time_name, "calendar", m_grid->ctx()->time()->calendar());
284  output.write_attribute(m_time_name, "units", m_grid->ctx()->time()->units_string());
285  }
286 }
287 
290 
291  output.write_variable(m_time_name, {0}, {1}, &m_t_last);
292 }
293 
295 
296  return combine({{"tillphi", Diagnostic::wrap(m_till_phi)},
297  {"usurf_difference", Diagnostic::wrap(m_usurf_difference)},
298  {"usurf_target", Diagnostic::wrap(m_usurf_target)},
299  {"diff_mask", Diagnostic::wrap(m_mask)}},
301 }
302 
303 } // 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
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
static Ptr wrap(const T &input)
Definition: Diagnostic.hh:155
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:740
void define_variable(const std::string &name, IO_Type nctype, const std::vector< std::string > &dims) const
Define a variable.
Definition: File.cc:566
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:364
std::string filename() const
Definition: File.cc:310
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:753
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:631
High-level PISM I/O class.
Definition: File.hh:51
IceModelVec2S ice_surface_elevation
Definition: Geometry.hh:59
IceModelVec2CellType cell_type
Definition: Geometry.hh:57
IceModelVec2S bed_elevation
Definition: Geometry.hh:49
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
bool ocean(int i, int j) const
bool grounded_ice(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 set_time_independent(bool flag)
Set the time independent flag for all variables corresponding to this IceModelVec instance.
Definition: iceModelVec.cc:175
Combines the max. time step with the flag indicating if a restriction is active. Makes is possible to...
Definition: MaxTimestep.hh:31
void restart_impl(const File &input_file, int record)
void update_impl(const YieldStressInputs &inputs, double t, double dt)
void write_model_state_impl(const File &output) const
The default (empty implementation).
void define_model_state_impl(const File &output) const
MaxTimestep max_timestep_impl(double t) const
void bootstrap_impl(const File &input_file, const YieldStressInputs &inputs)
Initialize the pseudo-plastic till mechanical model.
PISM's default basal yield stress model which applies the Mohr-Coulomb model of deformable,...
void define_model_state_impl(const File &output) const
MaxTimestep max_timestep_impl(double t) const
void bootstrap_impl(const File &input_file, const YieldStressInputs &inputs)
void init_t_last(const File &input_file)
void init_impl(const YieldStressInputs &inputs)
OptTillphiYieldStress(IceGrid::ConstPtr g)
double m_update_interval
Update interval in seconds.
DiagnosticList diagnostics_impl() const
double m_t_last
time of the last till friction angle update
void update_tillphi(const IceModelVec2S &ice_surface_elevation, const IceModelVec2S &bed_topography, const IceModelVec2CellType &mask)
void restart_impl(const File &input_file, int record)
std::string m_time_name
Name of the variable used to store the last update time.
double m_t_eps
Temporal resolution to use when checking whether it's time to update.
void write_model_state_impl(const File &output) const
The default (empty implementation).
void update_impl(const YieldStressInputs &inputs, double t, double dt)
void init_usurf_target(const File &input_file)
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
void set_name(const std::string &name)
const Geometry * geometry
Definition: YieldStress.hh:33
std::string m_name
Definition: YieldStress.hh:77
DiagnosticList diagnostics_impl() const
Definition: YieldStress.cc:105
#define PISM_ERROR_LOCATION
@ PISM_DOUBLE
Definition: IO_Flags.hh:38
T clip(T x, T a, T b)
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:117
static const double k
Definition: exactTestP.cc:45
@ CRITICAL
Definition: IO_Flags.hh:70
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
T combine(const T &a, const T &b)
@ WITH_GHOSTS
Definition: iceModelVec.hh:49