23 #include "pism/util/pism_utilities.hh"
24 #include "pism/util/EnthalpyConverter.hh"
25 #include "pism/util/pism_options.hh"
26 #include "pism/util/iceModelVec.hh"
28 #include "pism/util/ConfigInterface.hh"
29 #include "pism/util/IceGrid.hh"
31 #include "pism/util/error_handling.hh"
61 schoofLen = config.
get_number(
"flow_law.Schoof_regularizing_length",
"m"),
62 schoofVel = config.
get_number(
"flow_law.Schoof_regularizing_velocity",
"m second-1");
90 double pressure,
double gs)
const {
91 return this->
flow_impl(stress, enthalpy, pressure, gs);
95 double pressure,
double )
const {
96 return softness(enthalpy, pressure) * pow(stress,
m_n-1);
100 const double *pressure,
const double *grainsize,
101 unsigned int n,
double *result)
const {
102 this->
flow_n_impl(stress, enthalpy, pressure, grainsize,
n, result);
106 const double *pressure,
const double *grainsize,
107 unsigned int n,
double *result)
const {
108 for (
unsigned int k = 0;
k <
n; ++
k) {
109 result[
k] = this->
flow(stress[
k], enthalpy[
k], pressure[
k], grainsize[
k]);
123 unsigned int n,
double *result)
const {
128 unsigned int n,
double *result)
const {
129 for (
unsigned int k = 0;
k <
n; ++
k) {
130 result[
k] = this->
hardness(enthalpy[
k], pressure[
k]);
164 double *nu,
double *dnu)
const {
169 double *nu,
double *dnu)
const {
173 if (PetscLikely(nu != NULL)) {
177 if (PetscLikely(dnu != NULL)) {
194 const int i = p.i(), j = p.j();
197 double H = thickness(i,j);
198 const double *enthColumn = enthalpy.
get_column(i, j);
200 &(grid.
z()[0]), enthColumn);
215 double thickness,
int kbelowH,
216 const double *zlevels,
217 const double *enthalpy) {
229 for (
int i = 1; i <= kbelowH; ++i) {
231 p1 = EC.
pressure(thickness - zlevels[i]),
236 B += (zlevels[i] - zlevels[i-1]) * (
h0 + h1);
248 depth = thickness - zlevels[kbelowH],
251 B += depth * ice.
hardness(enthalpy[kbelowH], p);
264 static const double gs[] = {1e-4, 1e-3, 1e-2, 1}, s=1e4, E=400000, p=1e6;
265 double ref = flow_law.
flow(s, E, p, gs[0]);
266 for (
int i=1; i<4; i++) {
267 if (flow_law.
flow(s, E, p, gs[i]) != ref) {
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
double get_number(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
A class for storing and accessing PISM configuration flags and parameters.
std::shared_ptr< EnthalpyConverter > Ptr
double pressure(double depth) const
Get pressure in ice from depth below surface using the hydrostatic assumption.
Converts between specific enthalpy and temperature or liquid content.
unsigned int kBelowHeight(double height) const
Return the index k into zlevels[] so that zlevels[k] <= height < zlevels[k+1] and k < Mz.
const std::vector< double > & z() const
Z-coordinates within the ice.
Describes the PISM grid and the distribution of data across processors.
double * get_column(int i, int j)
A virtual class collecting methods common to ice and bedrock 3D fields.
void update_ghosts()
Updates ghost points.
IceGrid::ConstPtr grid() const
void failed()
Indicates a failure of a parallel section.
double m_crit_temp
critical temperature (cold – warm transition)
double softness_paterson_budd(double T_pa) const
Return the softness parameter A(T) for a given temperature T.
double m_schoofReg
regularization parameter for
virtual double hardness_impl(double E, double p) const
double m_A_warm
Paterson-Budd softness, warm case.
double m_A_cold
Paterson-Budd softness, cold case.
double m_Q_cold
Activation energy, cold case.
double m_standard_gravity
acceleration due to gravity
double m_viscosity_power
; used to compute viscosity
void flow_n(const double *stress, const double *E, const double *pressure, const double *grainsize, unsigned int n, double *result) const
virtual void hardness_n_impl(const double *enthalpy, const double *pressure, unsigned int n, double *result) const
void effective_viscosity(double hardness, double gamma, double *nu, double *dnu) const
Computes the regularized effective viscosity and its derivative with respect to the second invariant ...
void hardness_n(const double *enthalpy, const double *pressure, unsigned int n, double *result) const
FlowLaw(const std::string &prefix, const Config &config, EnthalpyConverter::Ptr EC)
double flow(double stress, double E, double pressure, double grainsize) const
The flow law itself.
double softness(double E, double p) const
virtual double softness_impl(double E, double p) const =0
double m_hardness_power
; used to compute hardness
double m_melting_point_temp
melting point temperature (for water, 273.15 K)
EnthalpyConverter::Ptr EC() const
double m_beta_CC_grad
Clausius-Clapeyron gradient.
double m_ideal_gas_constant
ideal gas constant
virtual double flow_impl(double stress, double E, double pressure, double grainsize) const
double m_Q_warm
Activation energy, warm case.
double m_n
power law exponent
double hardness(double E, double p) const
EnthalpyConverter::Ptr m_EC
virtual void flow_n_impl(const double *stress, const double *E, const double *pressure, const double *grainsize, unsigned int n, double *result) const
#define PISM_ERROR_LOCATION
double averaged_hardness(const FlowLaw &ice, double thickness, int kbelowH, const double *zlevels, const double *enthalpy)
Computes vertical average of B(E, p) ice hardness, namely .
void averaged_hardness_vec(const FlowLaw &ice, const IceModelVec2S &thickness, const IceModelVec3 &enthalpy, IceModelVec2S &result)
bool FlowLawUsesGrainSize(const FlowLaw &flow_law)