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