tFlowLaw.cc - 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.cc (9375B)
---
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 #include "FlowLaw.hh"
20 #include "pism/util/pism_utilities.hh"
21 #include "pism/util/EnthalpyConverter.hh"
22 #include "pism/util/pism_options.hh"
23 #include "pism/util/iceModelVec.hh"
24
25 #include "pism/util/ConfigInterface.hh"
26 #include "pism/util/IceGrid.hh"
27
28 #include "pism/util/error_handling.hh"
29
30 namespace pism {
31 namespace rheology {
32
33 FlowLaw::FlowLaw(const std::string &prefix, const Config &config,
34 EnthalpyConverter::Ptr ec)
35 : m_EC(ec), m_e(1) {
36
37 if (not m_EC) {
38 throw RuntimeError(PISM_ERROR_LOCATION, "EC is NULL in FlowLaw::FlowLaw()");
39 }
40
41 m_standard_gravity = config.get_number("constants.standard_gravity");
42 m_ideal_gas_constant = config.get_number("constants.ideal_gas_constant");
43
44 m_rho = config.get_number("constants.ice.density");
45 m_beta_CC_grad = config.get_number("constants.ice.beta_Clausius_Clapeyron") * m_rho * m_standard_gravity;
46 m_melting_point_temp = config.get_number("constants.fresh_water.melting_point_temperature");
47 m_e = config.get_number(prefix + "enhancement_factor");
48 m_e_interglacial = config.get_number(prefix + "enhancement_factor_interglacial");
49 m_n = config.get_number(prefix + "Glen_exponent");
50 m_viscosity_power = (1.0 - m_n) / (2.0 * m_n);
51 m_hardness_power = -1.0 / m_n;
52
53 m_A_cold = config.get_number("flow_law.Paterson_Budd.A_cold");
54 m_A_warm = config.get_number("flow_law.Paterson_Budd.A_warm");
55 m_Q_cold = config.get_number("flow_law.Paterson_Budd.Q_cold");
56 m_Q_warm = config.get_number("flow_law.Paterson_Budd.Q_warm");
57 m_crit_temp = config.get_number("flow_law.Paterson_Budd.T_critical");
58 m_schoofLen = config.get_number("flow_law.Schoof_regularizing_length", "m"); // convert to meters
59 m_schoofVel = config.get_number("flow_law.Schoof_regularizing_velocity", "m second-1"); // convert to m second-1
60 m_schoofReg = PetscSqr(m_schoofVel/m_schoofLen);
61 }
62
63 FlowLaw::~FlowLaw() {
64 // empty
65 }
66
67 std::string FlowLaw::name() const {
68 return m_name;
69 }
70
71 EnthalpyConverter::Ptr FlowLaw::EC() const {
72 return m_EC;
73 }
74
75 double FlowLaw::exponent() const {
76 return m_n;
77 }
78
79 double FlowLaw::enhancement_factor() const {
80 return m_e;
81 }
82
83 double FlowLaw::enhancement_factor_interglacial() const {
84 return m_e_interglacial;
85 }
86
87 //! Return the softness parameter A(T) for a given temperature T.
88 /*! This is not a natural part of all FlowLaw instances. */
89 double FlowLaw::softness_paterson_budd(double T_pa) const {
90 const double A = T_pa < m_crit_temp ? m_A_cold : m_A_warm;
91 const double Q = T_pa < m_crit_temp ? m_Q_cold : m_Q_warm;
92
93 return A * exp(-Q / (m_ideal_gas_constant * T_pa));
94 }
95
96 //! The flow law itself.
97 double FlowLaw::flow(double stress, double enthalpy,
98 double pressure, double gs) const {
99 return this->flow_impl(stress, enthalpy, pressure, gs);
100 }
101
102 double FlowLaw::flow_impl(double stress, double enthalpy,
103 double pressure, double /* gs */) const {
104 return softness(enthalpy, pressure) * pow(stress, m_n-1);
105 }
106
107 void FlowLaw::flow_n(const double *stress, const double *enthalpy,
108 const double *pressure, const double *grainsize,
109 unsigned int n, double *result) const {
110 this->flow_n_impl(stress, enthalpy, pressure, grainsize, n, result);
111 }
112
113 void FlowLaw::flow_n_impl(const double *stress, const double *enthalpy,
114 const double *pressure, const double *grainsize,
115 unsigned int n, double *result) const {
116 for (unsigned int k = 0; k < n; ++k) {
117 result[k] = this->flow(stress[k], enthalpy[k], pressure[k], grainsize[k]);
118 }
119 }
120
121
122 double FlowLaw::softness(double E, double p) const {
123 return this->softness_impl(E, p);
124 }
125
126 double FlowLaw::hardness(double E, double p) const {
127 return this->hardness_impl(E, p);
128 }
129
130 void FlowLaw::hardness_n(const double *enthalpy, const double *pressure,
131 unsigned int n, double *result) const {
132 this->hardness_n_impl(enthalpy, pressure, n, result);
133 }
134
135 void FlowLaw::hardness_n_impl(const double *enthalpy, const double *pressure,
136 unsigned int n, double *result) const {
137 for (unsigned int k = 0; k < n; ++k) {
138 result[k] = this->hardness(enthalpy[k], pressure[k]);
139 }
140 }
141
142 double FlowLaw::hardness_impl(double E, double p) const {
143 return pow(softness(E, p), m_hardness_power);
144 }
145
146 //! \brief Computes the regularized effective viscosity and its derivative with respect to the
147 //! second invariant \f$ \gamma \f$.
148 /*!
149 *
150 * @f{align*}{
151 * \nu &= \frac{1}{2} B \left( \epsilon + \gamma \right)^{(1-n)/(2n)},\\
152 * \diff{\nu}{\gamma} &= \frac{1}{2} B \cdot \frac{1-n}{2n} \cdot \left(\epsilon + \gamma \right)^{(1-n)/(2n) - 1}, \\
153 * &= \frac{1-n}{2n} \cdot \frac{1}{2} B \left( \epsilon + \gamma \right)^{(1-n)/(2n)} \cdot \frac{1}{\epsilon + \gamma}, \\
154 * &= \frac{1-n}{2n} \cdot \frac{\nu}{\epsilon + \gamma}.
155 * @f}
156 * Here @f$ \gamma @f$ is the second invariant
157 * @f{align*}{
158 * \gamma &= \frac{1}{2} D_{ij} D_{ij}\\
159 * &= \frac{1}{2}\, ((u_x)^2 + (v_y)^2 + (u_x + v_y)^2 + \frac{1}{2}\, (u_y + v_x)^2) \\
160 * @f}
161 * and
162 * @f[ D_{ij}(\mathbf{u}) = \frac{1}{2}\left(\diff{u_{i}}{x_{j}} + \diff{u_{j}}{x_{i}}\right). @f]
163 *
164 * Either one of \c nu and \c dnu can be NULL if the corresponding output is not needed.
165 *
166 * \param[in] B ice hardness
167 * \param[in] gamma the second invariant
168 * \param[out] nu effective viscosity
169 * \param[out] dnu derivative of \f$ \nu \f$ with respect to \f$ \gamma \f$
170 */
171 void FlowLaw::effective_viscosity(double B, double gamma,
172 double *nu, double *dnu) const {
173 const double
174 my_nu = 0.5 * B * pow(m_schoofReg + gamma, m_viscosity_power);
175
176 if (PetscLikely(nu != NULL)) {
177 *nu = my_nu;
178 }
179
180 if (PetscLikely(dnu != NULL)) {
181 *dnu = m_viscosity_power * my_nu / (m_schoofReg + gamma);
182 }
183 }
184
185 void averaged_hardness_vec(const FlowLaw &ice,
186 const IceModelVec2S &thickness,
187 const IceModelVec3 &enthalpy,
188 IceModelVec2S &result) {
189
190 const IceGrid &grid = *thickness.grid();
191
192 IceModelVec::AccessList list{&thickness, &result, &enthalpy};
193
194 ParallelSection loop(grid.com);
195 try {
196 for (Points p(grid); p; p.next()) {
197 const int i = p.i(), j = p.j();
198
199 // Evaluate column integrals in flow law at every quadrature point's column
200 double H = thickness(i,j);
201 const double *enthColumn = enthalpy.get_column(i, j);
202 result(i,j) = averaged_hardness(ice, H, grid.kBelowHeight(H),
203 &(grid.z()[0]), enthColumn);
204 }
205 } catch (...) {
206 loop.failed();
207 }
208 loop.check();
209
210 result.update_ghosts();
211 }
212
213 //! Computes vertical average of `B(E, p)` ice hardness, namely @f$\bar B(E, p)@f$.
214 /*!
215 * See comment for hardness(). Note `E[0], ..., E[kbelowH]` must be valid.
216 */
217 double averaged_hardness(const FlowLaw &ice,
218 double thickness, int kbelowH,
219 const double *zlevels,
220 const double *enthalpy) {
221 double B = 0;
222
223 EnthalpyConverter &EC = *ice.EC();
224
225 // Use trapezoidal rule to integrate from 0 to zlevels[kbelowH]:
226 if (kbelowH > 0) {
227 double
228 p0 = EC.pressure(thickness),
229 E0 = enthalpy[0],
230 h0 = ice.hardness(E0, p0); // ice hardness at the left endpoint
231
232 for (int i = 1; i <= kbelowH; ++i) { // note the "1" and the "<="
233 const double
234 p1 = EC.pressure(thickness - zlevels[i]), // pressure at the right endpoint
235 E1 = enthalpy[i], // enthalpy at the right endpoint
236 h1 = ice.hardness(E1, p1); // ice hardness at the right endpoint
237
238 // The trapezoid rule sans the "1/2":
239 B += (zlevels[i] - zlevels[i-1]) * (h0 + h1);
240
241 h0 = h1;
242 }
243 }
244
245 // Add the "1/2":
246 B *= 0.5;
247
248 // use the "rectangle method" to integrate from
249 // zlevels[kbelowH] to thickness:
250 double
251 depth = thickness - zlevels[kbelowH],
252 p = EC.pressure(depth);
253
254 B += depth * ice.hardness(enthalpy[kbelowH], p);
255
256 // Now B is an integral of ice hardness; next, compute the average:
257 if (thickness > 0) {
258 B = B / thickness;
259 } else {
260 B = 0;
261 }
262
263 return B;
264 }
265
266 bool FlowLawUsesGrainSize(const FlowLaw &flow_law) {
267 static const double gs[] = {1e-4, 1e-3, 1e-2, 1}, s=1e4, E=400000, p=1e6;
268 double ref = flow_law.flow(s, E, p, gs[0]);
269 for (int i=1; i<4; i++) {
270 if (flow_law.flow(s, E, p, gs[i]) != ref) {
271 return true;
272 }
273 }
274 return false;
275 }
276
277 } // end of namespace rheology
278 } // end of namespace pism