19 #include "pism/basalstrength/OptTillphiYieldStress.hh"
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"
39 m_mask(m_grid,
"diff_mask"),
40 m_usurf_difference(m_grid,
"usurf_difference"),
41 m_usurf_target(m_grid,
"usurf")
46 m_name =
"Iterative optimization of the till friction angle for the Mohr-Coulomb yield stress model";
49 .
long_name(
"target surface elevation for tillphi optimization")
54 .
long_name(
"difference between modeled and target surface elevations")
61 "one if the till friction angle was updated by the last iteration, zero otherwise ");
64 m_mask.
metadata()[
"flag_meanings"] =
"no_update updated_during_last_iteration";
68 m_dhdt_min =
m_config->get_number(
"basal_yield_stress.mohr_coulomb.tillphi_opt.dhdt_min",
"m / s");
71 m_dphi_scale =
m_config->get_number(
"basal_yield_stress.mohr_coulomb.tillphi_opt.dphi_scale",
"degree / m");
73 m_dphi_max =
m_config->get_number(
"basal_yield_stress.mohr_coulomb.tillphi_opt.dphi_max");
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");
81 m_phi_max =
m_config->get_number(
"basal_yield_stress.mohr_coulomb.tillphi_opt.phi_max");
85 "basal_yield_stress.mohr_coulomb.tillphi_opt: phi0_min >= phi_max");
90 "basal_yield_stress.mohr_coulomb.tillphi_opt: topg_min >= topg_max");
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",
127 auto filename =
m_config->get_string(
"basal_yield_stress.mohr_coulomb.tillphi_opt.file");
129 if (not filename.empty()) {
132 m_log->message(2,
"* No file set to read target surface elevation from... using '%s'\n",
165 "not implemented: till friction angle optimization "
166 "cannot be initialized without an input file");
174 "time %f is less than the previous time %f",
195 return std::min(dt_max, dt_mohr_coulomb);
199 double t,
double dt) {
209 if (std::abs(t_next - t_final) <
m_t_eps) {
225 m_log->message(2,
"* Updating till friction angle...\n");
234 for (
auto p =
m_grid->points(); p; p.next()) {
235 const int i = p.i(), j = p.j();
242 }
else if (bed_topography(i, j) <=
m_topg_max) {
266 }
else if (cell_type.
ocean(i, j)) {
281 "time of the last update of the till friction angle");
const Time & time() const
const Config::ConstPtr m_config
configuration database used by this component
const Logger::ConstPtr m_log
logger (for easy access)
const std::shared_ptr< const Grid > 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
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 define_variable(const std::string &name, io::Type nctype, const std::vector< std::string > &dims) const
Define a variable.
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.
array::Scalar2 ice_surface_elevation
array::CellType2 cell_type
array::Scalar2 bed_elevation
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)
array::Scalar1 m_usurf_target
void init_impl(const YieldStressInputs &inputs)
double m_update_interval
Update interval in seconds.
array::Scalar1 m_usurf_difference
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.
DiagnosticList diagnostics_impl() const
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
void set(double c)
Result: v[j] <- c for all j.
void regrid(const std::string &filename, io::Default default_value)
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
bool grounded_ice(int i, int j) const
bool ocean(int i, int j) const
"Cell type" mask. Adds convenience methods to array::Scalar.
#define PISM_ERROR_LOCATION
double min(const array::Scalar &input)
Finds minimum over all the values in an array::Scalar object. Ignores ghosts.
std::map< std::string, Diagnostic::Ptr > DiagnosticList
static std::string calendar(const File *input_file, const Config &config, const Logger &log)
T combine(const T &a, const T &b)