tFlowLaw.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tFlowLaw.hh (5873B)
---
1 // Copyright (C) 2004-2018 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,
79 EnthalpyConverter::Ptr EC);
80 virtual ~FlowLaw();
81
82 void effective_viscosity(double hardness, double gamma,
83 double *nu, double *dnu) const;
84
85 std::string name() const;
86 double exponent() const;
87 double enhancement_factor() const;
88 double enhancement_factor_interglacial() const;
89
90 EnthalpyConverter::Ptr EC() const;
91
92 double hardness(double E, double p) const;
93 void hardness_n(const double *enthalpy, const double *pressure,
94 unsigned int n, double *result) const;
95
96 double softness(double E, double p) const;
97
98 double flow(double stress, double E, double pressure, double grainsize) const;
99 void flow_n(const double *stress, const double *E,
100 const double *pressure, const double *grainsize,
101 unsigned int n, double *result) const;
102
103 protected:
104 virtual double flow_impl(double stress, double E,
105 double pressure, double grainsize) const;
106 virtual void flow_n_impl(const double *stress, const double *E,
107 const double *pressure, const double *grainsize,
108 unsigned int n, double *result) const;
109 virtual double hardness_impl(double E, double p) const;
110 virtual void hardness_n_impl(const double *enthalpy, const double *pressure,
111 unsigned int n, double *result) const;
112 virtual double softness_impl(double E, double p) const = 0;
113
114 protected:
115 std::string m_name;
116
117 double m_rho, //!< ice density
118 m_beta_CC_grad, //!< Clausius-Clapeyron gradient
119 m_melting_point_temp; //!< for water, 273.15 K
120 EnthalpyConverter::Ptr m_EC;
121
122 double softness_paterson_budd(double T_pa) const;
123
124 //! regularizing length
125 double m_schoofLen;
126 //! regularizing velocity
127 double m_schoofVel;
128 //! regularization parameter for @f$ \gamma @f$
129 double m_schoofReg;
130
131 //! @f$ (1 - n) / (2n) @f$; used to compute viscosity
132 double m_viscosity_power;
133 //! @f$ - 1 / n @f$; used to compute hardness
134 double m_hardness_power;
135
136 //! Paterson-Budd softness, cold case
137 double m_A_cold;
138 //! Paterson-Budd softness, warm case
139 double m_A_warm;
140 //! Activation energy, cold case
141 double m_Q_cold;
142 //! Activation energy, warm case
143 double m_Q_warm;
144 //! critical temperature (cold -- warm transition)
145 double m_crit_temp;
146
147 //! acceleration due to gravity
148 double m_standard_gravity;
149 //! ideal gas constant
150 double m_ideal_gas_constant;
151 //! flow enhancement factor
152 double m_e;
153 //! flow enhancement factor for interglacial periods
154 double m_e_interglacial;
155 //! power law exponent
156 double m_n;
157 };
158
159 double averaged_hardness(const FlowLaw &ice,
160 double ice_thickness,
161 int kbelowH,
162 const double *zlevels,
163 const double *enthalpy);
164
165 void averaged_hardness_vec(const FlowLaw &ice,
166 const IceModelVec2S &ice_thickness,
167 const IceModelVec3 &enthalpy,
168 IceModelVec2S &result);
169
170 // Helper functions:
171 bool FlowLawUsesGrainSize(const FlowLaw &flow_law);
172
173 } // end of namespace rheology
174 } // end of namespace pism
175
176 #endif // __flowlaws_hh