Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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
31namespace 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*/
37OptTillphiYieldStress::OptTillphiYieldStress(std::shared_ptr<const Grid> 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")
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");
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 */
116 if (input_file.variable_exists(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.name().c_str());
134
136 }
137
138 m_usurf_target.metadata().set_name("usurf_target");
139}
140
141void 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) {
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.
221void 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
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.variable_exists(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
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)
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:699
bool variable_exists(const std::string &short_name) const
Checks if a variable exists.
Definition File.cc:378
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:712
void define_variable(const std::string &name, io::Type nctype, const std::vector< std::string > &dims) const
Define a variable.
Definition File.cc:543
std::string name() const
Definition File.cc:274
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:590
High-level PISM I/O class.
Definition File.hh:55
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...
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:461
VariableMetadata & long_name(const std::string &input)
VariableMetadata & units(const std::string &input)
VariableMetadata & set_name(const std::string &name)
VariableMetadata & set_time_independent(bool flag)
const Geometry * geometry
std::string m_name
DiagnosticList diagnostics_impl() const
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition Array.hh:64
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
void regrid(const std::string &filename, io::Default default_value)
Definition Array.cc:736
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition Array.cc:476
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:93
#define PISM_ERROR_LOCATION
@ PISM_DOUBLE
Definition IO_Flags.hh:52
T clip(T x, T a, T b)
std::map< std::string, Diagnostic::Ptr > DiagnosticList
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)