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