PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
basal_resistance.cc
Go to the documentation of this file.
1 // Copyright (C) 2004-2017, 2019, 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 <cmath> // pow, sqrt
20 
21 #include "basal_resistance.hh"
22 
23 #include "pism/util/ConfigInterface.hh"
24 #include "pism/util/Logger.hh"
25 
26 namespace pism {
27 
28 static double square(double x) {
29  return x * x;
30 }
31 
32 /* Purely plastic */
33 
35  m_plastic_regularize = config.get_number("basal_resistance.plastic.regularization", "m second-1");
36 }
37 
39  int threshold,
40  units::System::Ptr system) const {
41  log.message(threshold, "Using purely plastic till with eps = %10.5e m year-1.\n",
42  units::convert(system, m_plastic_regularize, "m second-1", "m year-1"));
43 }
44 
45 
46 //! Compute the drag coefficient for the basal shear stress.
47 double IceBasalResistancePlasticLaw::drag(double tauc, double vx, double vy) const {
48  const double magreg2 = square(m_plastic_regularize) + square(vx) + square(vy);
49 
50  return tauc / sqrt(magreg2);
51 }
52 
53 //! Compute the drag coefficient and its derivative with respect to \f$ \alpha = \frac 1 2 (u_x^2 + u_y^2) \f$
54 /**
55  * @f{align*}{
56  * \beta &= \frac{\tau_{c}}{|\mathbf{u}|} = \tau_{c}\cdot (|\mathbf{u}|^{2})^{-1/2}\\
57  * \diff{\beta}{\frac12 |\mathbf{u}|^{2}} &= -\frac{1}{|\mathbf{u}|^{2}}\cdot \beta
58  * @f}
59  */
60 void IceBasalResistancePlasticLaw::drag_with_derivative(double tauc, double vx, double vy,
61  double *beta, double *dbeta) const {
62  const double magreg2 = square(m_plastic_regularize) + square(vx) + square(vy);
63 
64  *beta = tauc / sqrt(magreg2);
65 
66  if (dbeta) {
67  *dbeta = -1 * (*beta) / magreg2;
68  }
69 }
70 
71 /* Pseudo-plastic */
72 
75  m_pseudo_q = config.get_number("basal_resistance.pseudo_plastic.q");
76  m_pseudo_u_threshold = config.get_number("basal_resistance.pseudo_plastic.u_threshold", "m second-1");
77  m_sliding_scale_factor_reduces_tauc = config.get_number("basal_resistance.pseudo_plastic.sliding_scale_factor");
78 }
79 
81  int threshold,
82  units::System::Ptr system) const {
83 
84  if (m_pseudo_q == 1.0) {
85  log.message(threshold,
86  "Using linearly viscous till with u_threshold = %.2f m year-1.\n",
87  units::convert(system, m_pseudo_u_threshold, "m second-1", "m year-1"));
88  } else {
89  log.message(threshold,
90  "Using pseudo-plastic till with eps = %10.5e m year-1, q = %.4f,"
91  " and u_threshold = %.2f m year-1.\n",
92  units::convert(system, m_plastic_regularize, "m second-1", "m year-1"),
93  m_pseudo_q,
94  units::convert(system, m_pseudo_u_threshold, "m second-1", "m year-1"));
95  }
96 }
97 
98 //! Compute the drag coefficient for the basal shear stress.
99 /*!
100 
101  The basal shear stress term @f$ \tau_b @f$ in the SSA stress balance
102  for ice is minus the return value here times (vx,vy). Thus this
103  method computes the basal shear stress as
104 
105  @f[ \tau_b = - \frac{\tau_c}{|\mathbf{U}|^{1-q} U_{\mathtt{th}}^q} \mathbf{U} @f]
106 
107  where @f$ \tau_b=(\tau_{(b)x},\tau_{(b)y}) @f$ , @f$ U=(u,v) @f$ ,
108  @f$ q= @f$ `pseudo_q`, and @f$ U_{\mathtt{th}}= @f$
109  `pseudo_u_threshold`. Typical values for the constants are
110  @f$ q=0.25 @f$ and @f$ U_{\mathtt{th}} = 100 @f$ m year-1.
111 
112  The linearly-viscous till case pseudo_q = 1.0 is allowed, in which
113  case @f$ \beta = \tau_c/U_{\mathtt{th}} @f$ . The purely-plastic till
114  case pseudo_q = 0.0 is also allowed; note that there is still a
115  regularization with data member plastic_regularize.
116 
117  One can scale tauc if desired:
118 
119  A scale factor of @f$ A @f$ is intended to increase basal sliding rate
120  by @f$ A @f$ . It would have exactly this effect \e if the driving
121  stress were \e hypothetically completely held by the basal
122  resistance. Thus this scale factor is used to reduce (if
123  `-sliding_scale_factor_reduces_tauc` @f$ A @f$ with @f$ A > 1 @f$) or increase (if @f$ A <
124  1 @f$) the value of the (pseudo-) yield stress `tauc`. The concept
125  behind this is described at
126  [the SeaRISE wiki](http://websrv.cs.umt.edu/isis/index.php/Category_1:_Whole_Ice_Sheet#Initial_Experiment_-_E1_-_Increased_Basal_Lubrication).
127 
128  Specifically, the concept behind this mechanism is to suppose
129  equality of driving and basal shear stresses,
130 
131  @f[ \rho g H \nabla h = \frac{\tau_c}{|\mathbf{U}|^{1-q} U_{\mathtt{th}}^q} \mathbf{U}. @f]
132 
133  (*For emphasis:* The membrane stress held by the ice itself is
134  missing from this incomplete stress balance.) Thus the pseudo yield
135  stress @f$ \tau_c @f$ would be related to the sliding speed
136  @f$ |\mathbf{U}| @f$ by
137 
138  @f[ |\mathbf{U}| = \frac{C}{\tau_c^{1/q}} \f]
139 
140  for some (geometry-dependent) constant @f$ C @f$ . Multiplying
141  @f$ |\mathbf{U}| @f$ by @f$ A @f$ in this equation corresponds to
142  dividing @f$ \tau_c @f$ by @f$ A^q @f$ .
143 
144  Note that the mechanism has no effect whatsoever if @f$ q=0 @f$ , which
145  is the purely plastic case. In that case there is \e no direct
146  relation between the yield stress and the sliding velocity, and the
147  difference between the driving stress and the yield stress is
148  entirely held by the membrane stresses. (There is also no singular
149  mathematical operation as @f$ A^q = A^0 = 1 @f$ .)
150 */
151 double IceBasalResistancePseudoPlasticLaw::drag(double tauc, double vx, double vy) const {
152  const double magreg2 = square(m_plastic_regularize) + square(vx) + square(vy);
153 
156  return (tauc / Aq) * pow(magreg2, 0.5*(m_pseudo_q - 1)) * pow(m_pseudo_u_threshold, -m_pseudo_q);
157  } else {
158  return tauc * pow(magreg2, 0.5*(m_pseudo_q - 1)) * pow(m_pseudo_u_threshold, -m_pseudo_q);
159  }
160 }
161 
162 
163 //! Compute the drag coefficient and its derivative with respect to @f$ \alpha = \frac 1 2 (u_x^2 + u_y^2) @f$
164 /**
165  * @f{align*}{
166  * \beta &= \frac{\tau_{c}}{u_{\text{threshold}}^q}\cdot (|u|^{2})^{\frac{q-1}{2}} \\
167  * \diff{\beta}{\frac12 |\mathbf{u}|^{2}} &= \frac{\tau_{c}}{u_{\text{threshold}}^q}\cdot \frac{q-1}{2}\cdot (|\mathbf{u}|^{2})^{\frac{q-1}{2} - 1}\cdot 2 \\
168  * &= \frac{q-1}{|\mathbf{u}|^{2}}\cdot \beta(\mathbf{u}) \\
169  * @f}
170  */
171 void IceBasalResistancePseudoPlasticLaw::drag_with_derivative(double tauc, double vx, double vy,
172  double *beta, double *dbeta) const
173 {
174  const double magreg2 = square(m_plastic_regularize) + square(vx) + square(vy);
175 
178  *beta = (tauc / Aq) * pow(magreg2, 0.5*(m_pseudo_q - 1)) * pow(m_pseudo_u_threshold, -m_pseudo_q);
179  } else {
180  *beta = tauc * pow(magreg2, 0.5*(m_pseudo_q - 1)) * pow(m_pseudo_u_threshold, -m_pseudo_q);
181  }
182 
183  if (dbeta) {
184  *dbeta = (m_pseudo_q - 1) * (*beta) / magreg2;
185  }
186 
187 }
188 
189 /* Regularized Coulomb */
190 
192  : IceBasalResistancePlasticLaw(config) {
193  m_pseudo_q = config.get_number("basal_resistance.pseudo_plastic.q");
194  m_pseudo_u_threshold = config.get_number("basal_resistance.pseudo_plastic.u_threshold", "m second-1");
195  m_sliding_scale_factor_reduces_tauc = config.get_number("basal_resistance.pseudo_plastic.sliding_scale_factor");
196 }
197 
199  int threshold,
200  units::System::Ptr system) const {
201  log.message(threshold,
202  "Using regularized Coulomb till with eps = %10.5e m year-1, q = %.4f,"
203  " and u_threshold = %.2f m year-1.\n",
204  units::convert(system, m_plastic_regularize, "m second-1", "m year-1"),
205  m_pseudo_q,
206  units::convert(system, m_pseudo_u_threshold, "m second-1", "m year-1"));
207 }
208 
209 double IceBasalResistanceRegularizedLaw::drag(double tauc, double vx, double vy) const {
210  const double magreg2sqr = sqrt(square(m_plastic_regularize) + square(vx) + square(vy));
211 
214  return (tauc / Aq) * pow(magreg2sqr, (m_pseudo_q - 1)) * pow((magreg2sqr + m_pseudo_u_threshold), -m_pseudo_q);
215  } else {
216  return tauc * pow(magreg2sqr, (m_pseudo_q - 1)) * pow((magreg2sqr + m_pseudo_u_threshold), -m_pseudo_q);
217  }
218 }
219 
220 //! Compute the drag coefficient and its derivative with respect to @f$ \alpha = \frac 1 2 (u_x^2 + u_y^2) @f$
221 /**
222  * @f{align*}{
223  * \beta &= \frac{\tau_{c}}{ \left( (|u|^{2})^{\frac12} + u_{\text{threshold}} \right)^{q} }\cdot (|\mathbf{u}|^{2})^{\frac{q-1}{2}} \\
224  * \diff{\beta}{\frac12 |\mathbf{u}|^{2}} &= \frac{\tau_{c}}{ \left( (|\mathbf{u}|^{2})^{\frac12} + u_{\text{threshold}}\right)^{q} }\cdot \frac{q-1}{2}\cdot (|\mathbf{u}|^{2})^{\frac{q-1}{2} - 1}\cdot 2 - \frac{\tau_{c}}{ \left( (|\mathbf{u}|^{2})^{\frac12} + u_{\text{threshold}}\right)^{q+1} }\cdot \frac{2 \cdot q }{(|\mathbf{u}|^{2})^{\frac12}} \cdot (|\mathbf{u}|^{2})^{\frac{q-1}{2} } \\
225  * &= \left( \frac{q-1}{|\mathbf{u}|^{2}} - \frac{q}{|\mathbf{u}| \cdot ( |\mathbf{u}| + u_{\text{threshold}}) } \right) \cdot \beta(\mathbf{u})\\
226  * @f}
227  * Same parameters are used as in the pseudo-plastic case, namely @f$ q, u_{\text{threshold}}, A_q @f$.
228  * It should be noted, that in the original equation (3) in Zoet et al, 2020 the exponent @f$ q=1/p @f$ is used.
229  */
230 void IceBasalResistanceRegularizedLaw::drag_with_derivative(double tauc, double vx, double vy,
231  double *beta, double *dbeta) const {
232  const double magreg2 = square(m_plastic_regularize) + square(vx) + square(vy),
233  magreg2sqr = sqrt(magreg2);
234 
237  *beta = (tauc / Aq) * pow(magreg2sqr, (m_pseudo_q - 1)) * pow((magreg2sqr + m_pseudo_u_threshold), -m_pseudo_q);
238  } else {
239  *beta = tauc * pow(magreg2sqr, (m_pseudo_q - 1)) * pow((magreg2sqr + m_pseudo_u_threshold), -m_pseudo_q);
240  }
241 
242  if (dbeta) {
243  *dbeta = (((m_pseudo_q - 1) / magreg2) - (m_pseudo_q / (magreg2sqr * (magreg2sqr + m_pseudo_u_threshold)))) * (*beta);
244  }
245 }
246 
247 
248 } // end of namespace pism
double get_number(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
A class for storing and accessing PISM configuration flags and parameters.
virtual double drag(double tauc, double vx, double vy) const
Compute the drag coefficient for the basal shear stress.
virtual void print_info(const Logger &log, int threshold, units::System::Ptr system) const
IceBasalResistancePlasticLaw(const Config &config)
virtual void drag_with_derivative(double tauc, double vx, double vy, double *drag, double *ddrag) const
Compute the drag coefficient and its derivative with respect to .
Class containing physical constants and the constitutive relation describing till for SSA.
IceBasalResistancePseudoPlasticLaw(const Config &config)
virtual double drag(double tauc, double vx, double vy) const
Compute the drag coefficient for the basal shear stress.
virtual void print_info(const Logger &log, int threshold, units::System::Ptr system) const
virtual void drag_with_derivative(double tauc, double vx, double vy, double *drag, double *ddrag) const
Compute the drag coefficient and its derivative with respect to .
virtual double drag(double tauc, double vx, double vy) const
Compute the drag coefficient for the basal shear stress.
virtual void print_info(const Logger &log, int threshold, units::System::Ptr system) const
IceBasalResistanceRegularizedLaw(const Config &config)
virtual void drag_with_derivative(double tauc, double vx, double vy, double *drag, double *ddrag) const
Compute the drag coefficient and its derivative with respect to .
void message(int threshold, const char format[],...) const __attribute__((format(printf
Print a message to the log.
Definition: Logger.cc:50
A basic logging class.
Definition: Logger.hh:40
std::shared_ptr< System > Ptr
Definition: Units.hh:47
double convert(System::Ptr system, double input, const std::string &spec1, const std::string &spec2)
Convert a quantity from unit1 to unit2.
Definition: Units.cc:72
static double square(double x)