PISM, A Parallel Ice Sheet Model  stable v2.0.6 committed by Constantine Khrulev on 2023-01-23 15:14:38 -0900
FlowLaw.hh
Go to the documentation of this file.
1 // Copyright (C) 2004-2018, 2020, 2021 Jed Brown, Ed Bueler, and Constantine Khroulev
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 #ifndef __flowlaws_hh
20 #define __flowlaws_hh
21 
22 #include <string>
23 
24 #include "pism/util/EnthalpyConverter.hh"
25 #include "pism/util/Vector2.hh"
26 
27 namespace pism {
28 
29 class IceModelVec2S;
30 class IceModelVec3;
31 
32 class Config;
33 
34 /*!
35  * This uses the definition of squared second invariant from Hutter and several others, namely the
36  * output is @f$ D^2 = \frac 1 2 D_{ij} D_{ij} @f$ where incompressibility is used to compute
37  * @f$ D_{zz}. @f$
38  *
39  * This is the approximation of the full second invariant corresponding to the shallow shelf
40  * approximation. In particular, we assume that @f$ u @f$ and @f$ v @f$ are depth-independent (@f$
41  * u_z = v_z = 0 @f$) and neglect horizontal derivatives of the vertical velocity (@f$ w_x = w_y = 0
42  * @f$).
43  */
44 static inline double secondInvariant_2D(const Vector2 &U_x, const Vector2 &U_y) {
45  const double
46  u_x = U_x.u,
47  u_y = U_y.u,
48  v_x = U_x.v,
49  v_y = U_y.v,
50  w_z = -(u_x + v_y); // w_z is computed assuming incompressibility of ice
51  return 0.5 * (u_x * u_x + v_y * v_y + w_z * w_z + 0.5*(u_y + v_x)*(u_y + v_x));
52 }
53 
54 //! Ice flow laws.
55 namespace rheology {
56 
57 //! Abstract class containing the constitutive relation for the flow of ice (of
58 //! the Paterson-Budd type).
59 /*!
60  This is the interface which most of PISM uses for rheology.
61 
62  The current implementation of stress-balance computations in PISM restrict
63  possible choices of rheologies to ones that
64 
65  - are power laws
66 
67  - allow factoring out a temperature- (or enthalpy-) dependent ice hardness
68  factor
69 
70  - can be represented in the viscosity form
71 
72  @note FlowLaw derived classes should implement hardness() in
73  terms of softness(). That way in many cases we only need to
74  re-implement softness... to turn one flow law into another.
75 */
76 class FlowLaw {
77 public:
78  FlowLaw(const std::string &prefix, const Config &config,
80  virtual ~FlowLaw() = default;
81 
82  void effective_viscosity(double hardness, double gamma,
83  double *nu, double *dnu) const;
84 
85  void effective_viscosity(double hardness, double gamma, double eps,
86  double *nu, double *dnu) const;
87 
88  std::string name() const;
89  double exponent() const;
90 
91  EnthalpyConverter::Ptr EC() const;
92 
93  double hardness(double E, double p) const;
94  void hardness_n(const double *enthalpy, const double *pressure,
95  unsigned int n, double *result) const;
96 
97  double softness(double E, double p) const;
98 
99  double flow(double stress, double E, double pressure, double grainsize) const;
100  void flow_n(const double *stress, const double *E,
101  const double *pressure, const double *grainsize,
102  unsigned int n, double *result) const;
103 
104 protected:
105  virtual double flow_impl(double stress, double E,
106  double pressure, double grainsize) const;
107  virtual void flow_n_impl(const double *stress, const double *E,
108  const double *pressure, const double *grainsize,
109  unsigned int n, double *result) const;
110  virtual double hardness_impl(double E, double p) const;
111  virtual void hardness_n_impl(const double *enthalpy, const double *pressure,
112  unsigned int n, double *result) const;
113  virtual double softness_impl(double E, double p) const = 0;
114 
115 protected:
116  std::string m_name;
117 
118  //! ice density
119  double m_rho;
120  //! Clausius-Clapeyron gradient
122  //! melting point temperature (for water, 273.15 K)
124 
126 
127  double softness_paterson_budd(double T_pa) const;
128 
129  //! regularization parameter for @f$ \gamma @f$
130  double m_schoofReg;
131 
132  //! @f$ (1 - n) / (2n) @f$; used to compute viscosity
134  //! @f$ - 1 / n @f$; used to compute hardness
136 
137  //! Paterson-Budd softness, cold case
138  double m_A_cold;
139  //! Paterson-Budd softness, warm case
140  double m_A_warm;
141  //! Activation energy, cold case
142  double m_Q_cold;
143  //! Activation energy, warm case
144  double m_Q_warm;
145  //! critical temperature (cold -- warm transition)
146  double m_crit_temp;
147 
148  //! acceleration due to gravity
150  //! ideal gas constant
152  //! power law exponent
153  double m_n;
154 };
155 
156 double averaged_hardness(const FlowLaw &ice,
157  double ice_thickness,
158  int kbelowH,
159  const double *zlevels,
160  const double *enthalpy);
161 
162 void averaged_hardness_vec(const FlowLaw &ice,
163  const IceModelVec2S &ice_thickness,
164  const IceModelVec3 &enthalpy,
165  IceModelVec2S &result);
166 
167 bool FlowLawUsesGrainSize(const FlowLaw &flow_law);
168 
169 } // end of namespace rheology
170 } // end of namespace pism
171 
172 #endif // __flowlaws_hh
A class for storing and accessing PISM configuration flags and parameters.
std::shared_ptr< EnthalpyConverter > Ptr
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: iceModelVec.hh:404
double v
Definition: Vector2.hh:103
double u
Definition: Vector2.hh:103
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2.hh:29
double m_crit_temp
critical temperature (cold – warm transition)
Definition: FlowLaw.hh:146
double softness_paterson_budd(double T_pa) const
Return the softness parameter A(T) for a given temperature T.
Definition: FlowLaw.cc:81
double m_schoofReg
regularization parameter for
Definition: FlowLaw.hh:130
virtual double hardness_impl(double E, double p) const
Definition: FlowLaw.cc:134
double m_A_warm
Paterson-Budd softness, warm case.
Definition: FlowLaw.hh:140
double m_A_cold
Paterson-Budd softness, cold case.
Definition: FlowLaw.hh:138
double m_Q_cold
Activation energy, cold case.
Definition: FlowLaw.hh:142
double m_standard_gravity
acceleration due to gravity
Definition: FlowLaw.hh:149
double m_viscosity_power
; used to compute viscosity
Definition: FlowLaw.hh:133
virtual ~FlowLaw()=default
void flow_n(const double *stress, const double *E, const double *pressure, const double *grainsize, unsigned int n, double *result) const
Definition: FlowLaw.cc:99
virtual void hardness_n_impl(const double *enthalpy, const double *pressure, unsigned int n, double *result) const
Definition: FlowLaw.cc:127
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 ...
Definition: FlowLaw.cc:163
void hardness_n(const double *enthalpy, const double *pressure, unsigned int n, double *result) const
Definition: FlowLaw.cc:122
std::string name() const
Definition: FlowLaw.cc:67
double exponent() const
Definition: FlowLaw.cc:75
FlowLaw(const std::string &prefix, const Config &config, EnthalpyConverter::Ptr EC)
Definition: FlowLaw.cc:36
std::string m_name
Definition: FlowLaw.hh:116
double flow(double stress, double E, double pressure, double grainsize) const
The flow law itself.
Definition: FlowLaw.cc:89
double softness(double E, double p) const
Definition: FlowLaw.cc:114
virtual double softness_impl(double E, double p) const =0
double m_rho
ice density
Definition: FlowLaw.hh:119
double m_hardness_power
; used to compute hardness
Definition: FlowLaw.hh:135
double m_melting_point_temp
melting point temperature (for water, 273.15 K)
Definition: FlowLaw.hh:123
EnthalpyConverter::Ptr EC() const
Definition: FlowLaw.cc:71
double m_beta_CC_grad
Clausius-Clapeyron gradient.
Definition: FlowLaw.hh:121
double m_ideal_gas_constant
ideal gas constant
Definition: FlowLaw.hh:151
virtual double flow_impl(double stress, double E, double pressure, double grainsize) const
Definition: FlowLaw.cc:94
double m_Q_warm
Activation energy, warm case.
Definition: FlowLaw.hh:144
double m_n
power law exponent
Definition: FlowLaw.hh:153
double hardness(double E, double p) const
Definition: FlowLaw.cc:118
EnthalpyConverter::Ptr m_EC
Definition: FlowLaw.hh:125
virtual void flow_n_impl(const double *stress, const double *E, const double *pressure, const double *grainsize, unsigned int n, double *result) const
Definition: FlowLaw.cc:105
#define n
Definition: exactTestM.c:37
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 .
Definition: FlowLaw.cc:214
void averaged_hardness_vec(const FlowLaw &ice, const IceModelVec2S &thickness, const IceModelVec3 &enthalpy, IceModelVec2S &result)
Definition: FlowLaw.cc:182
bool FlowLawUsesGrainSize(const FlowLaw &flow_law)
Definition: FlowLaw.cc:263
static double secondInvariant_2D(const Vector2 &U_x, const Vector2 &U_y)
Definition: FlowLaw.hh:44