PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
OptTillphiYieldStress.cc
Go to the documentation of this file.
1 // Copyright (C) 2004--2023 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 "pism/basalstrength/OptTillphiYieldStress.hh"
20 
21 #include "pism/geometry/Geometry.hh"
22 #include "pism/util/Context.hh"
23 #include "pism/util/Grid.hh"
24 #include "pism/util/array/CellType.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 */
37 OptTillphiYieldStress::OptTillphiYieldStress(std::shared_ptr<const Grid> grid)
38  : MohrCoulombYieldStress(grid),
39  m_mask(m_grid, "diff_mask"),
40  m_usurf_difference(m_grid, "usurf_difference"),
41  m_usurf_target(m_grid, "usurf")
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  .long_name("target surface elevation for tillphi optimization")
50  .units("m")
51  .set_time_independent(true);
52 
54  .long_name("difference between modeled and target surface elevations")
55  .units("m");
56 
58 
60  .long_name(
61  "one if the till friction angle was updated by the last iteration, zero otherwise ");
62 
63  m_mask.metadata()["flag_values"] = {0.0, 1.0};
64  m_mask.metadata()["flag_meanings"] = "no_update updated_during_last_iteration";
65 
66  {
67  // convergence threshold
68  m_dhdt_min = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.dhdt_min", "m / s");
69 
70  // scale used to compute tillphi adjustment using the surface elevation mismatch:
71  m_dphi_scale = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.dphi_scale", "degree / m");
72  // upper and lower bounds of the tillphi adjustment during an iteration:
73  m_dphi_max = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.dphi_max");
74  m_dphi_min = -2 * m_dphi_max;
75  // lower bound of tillphi:
76  m_phi0_min = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.phi0_min");
77  m_phi0_max = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.phi0_max");
78  m_topg_min = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.topg_min");
79  m_topg_max = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.topg_max");
80  // upper bound of tillphi:
81  m_phi_max = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.phi_max");
82 
83  if (m_phi0_min >= m_phi_max) {
85  "basal_yield_stress.mohr_coulomb.tillphi_opt: phi0_min >= phi_max");
86  }
87 
88  if (m_topg_min >= m_topg_max) {
90  "basal_yield_stress.mohr_coulomb.tillphi_opt: topg_min >= topg_max");
91  }
92  }
93 
94  {
95  m_time_name = m_config->get_string("time.dimension_name") + "_tillphi_opt";
96  m_t_last = time().current();
97  m_update_interval = m_config->get_number("basal_yield_stress.mohr_coulomb.tillphi_opt.dt", "seconds");
98  m_t_eps = m_config->get_number("time_stepping.resolution", "seconds");
99  }
100 
101  m_log->message(2,
102  " Using iterative optimization of the till friction angle.\n"
103  " Lower bound phi0 of the till friction angle is a piecewise-linear function of bed elevation (b):\n"
104  " / %5.2f for b < %.f\n"
105  " phi0 = | %5.2f + (b - (%.f)) * (%.2f / %.f) for %.f < b < %.f\n"
106  " \\ %5.2f for %.f < b\n",
110 }
111 
112 /*!
113  * Initialize the last time tillphi was updated.
114  */
115 void OptTillphiYieldStress::init_t_last(const File &input_file) {
116  if (input_file.find_variable(m_time_name)) {
117  input_file.read_variable(m_time_name, {0}, {1}, &m_t_last);
118  } else {
119  m_t_last = time().current();
120  }
121 }
122 
123 /*!
124  * Initialize the target ice surface elevation.
125  */
127  auto filename = m_config->get_string("basal_yield_stress.mohr_coulomb.tillphi_opt.file");
128 
129  if (not filename.empty()) {
131  } else {
132  m_log->message(2, "* No file set to read target surface elevation from... using '%s'\n",
133  input_file.filename().c_str());
134 
135  m_usurf_target.regrid(input_file, io::Default::Nil());
136  }
137 
138  m_usurf_target.metadata().set_name("usurf_target");
139 }
140 
141 void OptTillphiYieldStress::restart_impl(const File &input_file, int record) {
142 
143  MohrCoulombYieldStress::restart_impl(input_file, record);
144 
145  init_t_last(input_file);
146 
147  init_usurf_target(input_file);
148 }
149 
150 //! Initialize the pseudo-plastic till mechanical model.
151 //! Target surface elevation and initial iteration time are set.
153  const YieldStressInputs &inputs) {
154 
155  MohrCoulombYieldStress::bootstrap_impl(input_file, inputs);
156 
157  init_t_last(input_file);
158 
159  init_usurf_target(input_file);
160 }
161 
163  (void) inputs;
165  "not implemented: till friction angle optimization "
166  "cannot be initialized without an input file");
167 }
168 
170  MaxTimestep dt_max;
171  {
172  if (t < m_t_last) {
174  "time %f is less than the previous time %f",
175  t, m_t_last);
176  }
177 
178  // Find the smallest time of the form m_t_last + k * m_update_interval that is greater
179  // than t
180  double k = ceil((t - m_t_last) / m_update_interval);
181 
182  double
183  t_next = m_t_last + k * m_update_interval,
184  dt = t_next - t;
185 
186  if (dt < m_t_eps) {
187  dt = m_update_interval;
188  }
189 
190  dt_max = MaxTimestep(dt, "tillphi_opt");
191  }
192 
193  auto dt_mohr_coulomb = MohrCoulombYieldStress::max_timestep_impl(t);
194 
195  return std::min(dt_max, dt_mohr_coulomb);
196 }
197 
199  double t, double dt) {
200 
201  double
202  t_next = m_t_last + m_update_interval,
203  t_final = t + dt;
204 
205  if (t_final < m_t_last) {
206  throw RuntimeError(PISM_ERROR_LOCATION, "cannot go back in time");
207  }
208 
209  if (std::abs(t_next - t_final) < m_t_eps) { // reached the next update time
211  inputs.geometry->bed_elevation,
212  inputs.geometry->cell_type);
213  m_t_last = t_final;
214  }
215 
217 }
218 
219 //! Perform an iteration to adjust the till friction angle according to the difference
220 //! between target and modeled surface elevations.
221 void OptTillphiYieldStress::update_tillphi(const array::Scalar &ice_surface_elevation,
222  const array::Scalar &bed_topography,
223  const array::CellType &cell_type) {
224 
225  m_log->message(2, "* Updating till friction angle...\n");
226 
227  array::AccessScope list
228  { &m_till_phi, &m_usurf_target, &m_usurf_difference, &m_mask, &ice_surface_elevation, &bed_topography, &cell_type };
229 
230  m_mask.set(0.0);
231 
232  double slope = (m_phi0_max - m_phi0_min) / (m_topg_max - m_topg_min);
233 
234  for (auto p = m_grid->points(); p; p.next()) {
235  const int i = p.i(), j = p.j();
236 
237  // Compute the lower bound of the till friction angle (default value corresponds
238  // to "bed_topography(i, j) > topg_max"):
239  double phi0 = m_phi0_max;
240  if (bed_topography(i, j) <= m_topg_min) {
241  phi0 = m_phi0_min;
242  } else if (bed_topography(i, j) <= m_topg_max) {
243  phi0 = m_phi0_min + (bed_topography(i, j) - m_topg_min) * slope;
244  }
245 
246  if (cell_type.grounded_ice(i, j)) {
247  double dh_previous = m_usurf_difference(i, j);
248  m_usurf_difference(i, j) = m_usurf_target(i, j) - ice_surface_elevation(i, j);
249  double dh_change = std::abs(m_usurf_difference(i, j) - dh_previous);
250 
251  if (dh_change / m_update_interval > m_dhdt_min) {
252  // Update tillphi if the rate of change of the surface elevation mismatch since
253  // the last iteration is greater than the convergence threshold m_dhdt_min.
254 
255  // Mark this location as "updated by the last iteration":
256  m_mask(i, j) = 1.0;
257 
258  // Compute (and clip) the tillphi adjustment
259  double dphi = m_usurf_difference(i, j) * m_dphi_scale;
260  dphi = pism::clip(dphi, m_dphi_min, m_dphi_max);
261 
262  // Update (and clip) the till friction angle:
263  m_till_phi(i, j) += dphi;
264  m_till_phi(i, j) = pism::clip(m_till_phi(i, j), phi0, m_phi_max);
265  }
266  } else if (cell_type.ocean(i, j)) {
267  // Floating and ice free ocean: use the bed-elevation-dependent lower bound of
268  // tillphi:
269  m_till_phi(i, j) = phi0;
270  }
271  } // end of the loop over grid points
272 }
273 
276 
277  if (not output.find_variable(m_time_name)) {
279 
280  output.write_attribute(m_time_name, "long_name",
281  "time of the last update of the till friction angle");
282  output.write_attribute(m_time_name, "calendar", time().calendar());
283  output.write_attribute(m_time_name, "units", time().units_string());
284  }
285 }
286 
289 
290  output.write_variable(m_time_name, {0}, {1}, &m_t_last);
291 }
292 
294 
295  return combine({{"tillphi", Diagnostic::wrap(m_till_phi)},
296  {"usurf_difference", Diagnostic::wrap(m_usurf_difference)},
297  {"usurf_target", Diagnostic::wrap(m_usurf_target)},
298  {"diff_mask", Diagnostic::wrap(m_mask)}},
300 }
301 
302 } // 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
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
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
std::string filename() const
Definition: File.cc:307
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
array::Scalar2 ice_surface_elevation
Definition: Geometry.hh:57
array::CellType2 cell_type
Definition: Geometry.hh:55
array::Scalar2 bed_elevation
Definition: Geometry.hh:47
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)
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 array::Scalar &ice_surface_elevation, const array::Scalar &bed_topography, const array::CellType &mask)
OptTillphiYieldStress(std::shared_ptr< const Grid > g)
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
double current() const
Current time, in seconds.
Definition: Time.cc:465
VariableMetadata & set_name(const std::string &name)
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
VariableMetadata & set_time_independent(bool flag)
const Geometry * geometry
Definition: YieldStress.hh:32
std::string m_name
Definition: YieldStress.hh:76
DiagnosticList diagnostics_impl() const
Definition: YieldStress.cc:104
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
void set(double c)
Result: v[j] <- c for all j.
Definition: Array.cc:707
void regrid(const std::string &filename, io::Default default_value)
Definition: Array.cc:814
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: Array.cc:553
bool grounded_ice(int i, int j) const
Definition: CellType.hh:46
bool ocean(int i, int j) const
Definition: CellType.hh:34
"Cell type" mask. Adds convenience methods to array::Scalar.
Definition: CellType.hh:30
static Default Nil()
Definition: IO_Flags.hh:97
#define PISM_ERROR_LOCATION
double min(const array::Scalar &input)
Finds minimum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:193
@ PISM_DOUBLE
Definition: IO_Flags.hh:52
T clip(T x, T a, T b)
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)