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"
40 m_usurf_difference(m_grid,
"usurf_difference",
WITH_GHOSTS),
46 m_name =
"Iterative optimization of the till friction angle for the Mohr-Coulomb yield stress model";
49 "target surface elevation for tillphi optimization",
54 "difference between modeled and target"
55 " surface elevations",
60 "one if the till friction angle was"
61 " updated by the last iteration, zero otherwise ",
65 m_mask.
metadata()[
"flag_meanings"] =
"no_update updated_during_last_iteration";
69 m_dhdt_min =
m_config->get_number(
"basal_yield_stress.mohr_coulomb.tillphi_opt.dhdt_min",
"m / s");
72 m_dphi_scale =
m_config->get_number(
"basal_yield_stress.mohr_coulomb.tillphi_opt.dphi_scale",
"degree / m");
74 m_dphi_max =
m_config->get_number(
"basal_yield_stress.mohr_coulomb.tillphi_opt.dphi_max");
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");
82 m_phi_max =
m_config->get_number(
"basal_yield_stress.mohr_coulomb.tillphi_opt.phi_max");
86 "basal_yield_stress.mohr_coulomb.tillphi_opt: phi0_min >= phi_max");
91 "basal_yield_stress.mohr_coulomb.tillphi_opt: topg_min >= topg_max");
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",
128 auto filename =
m_config->get_string(
"basal_yield_stress.mohr_coulomb.tillphi_opt.file");
130 if (not filename.empty()) {
133 m_log->message(2,
"* No file set to read target surface elevation from... using '%s'\n",
166 "not implemented: till friction angle optimization "
167 "cannot be initialized without an input file");
175 "time %f is less than the previous time %f",
196 return std::min(dt_max, dt_mohr_coulomb);
200 double t,
double dt) {
210 if (std::abs(t_next - t_final) <
m_t_eps) {
226 m_log->message(2,
"* Updating till friction angle...\n");
236 const int i = p.i(), j = p.j();
243 }
else if (bed_topography(i, j) <=
m_topg_max) {
267 }
else if (cell_type.
ocean(i, j)) {
282 "time of the last update of the till friction angle");
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
const Config::ConstPtr m_config
configuration database used by this component
const Logger::ConstPtr m_log
logger (for easy access)
const IceGrid::ConstPtr m_grid
grid used by this component
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
void define_variable(const std::string &name, IO_Type nctype, const std::vector< std::string > &dims) const
Define a variable.
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.
std::string filename() const
void write_variable(const std::string &variable_name, const std::vector< unsigned int > &start, const std::vector< unsigned int > &count, const double *op) const
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.
High-level PISM I/O class.
IceModelVec2S ice_surface_elevation
IceModelVec2CellType cell_type
IceModelVec2S bed_elevation
std::shared_ptr< const IceGrid > ConstPtr
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.
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.
void regrid(const std::string &filename, RegriddingFlag flag, double default_value=0.0)
void set(double c)
Result: v[j] <- c for all j.
void set_time_independent(bool flag)
Set the time independent flag for all variables corresponding to this IceModelVec instance.
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)
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
IceModelVec2S m_usurf_target
void update_tillphi(const IceModelVec2S &ice_surface_elevation, const IceModelVec2S &bed_topography, const IceModelVec2CellType &mask)
IceModelVec2S m_usurf_difference
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
DiagnosticList diagnostics_impl() const
#define PISM_ERROR_LOCATION
std::map< std::string, Diagnostic::Ptr > DiagnosticList
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)