PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
WeertmanSliding.cc
Go to the documentation of this file.
1 /* Copyright (C) 2018, 2021, 2022, 2023 PISM Authors
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 
20 #include "pism/stressbalance/WeertmanSliding.hh"
21 
22 #include "pism/rheology/FlowLawFactory.hh"
23 #include "pism/geometry/Geometry.hh"
24 #include "pism/stressbalance/StressBalance.hh"
25 
26 namespace pism {
27 namespace stressbalance {
28 
29 WeertmanSliding::WeertmanSliding(std::shared_ptr<const Grid> grid)
30  : ShallowStressBalance(grid) {
31  // Use the SIA flow law.
32  rheology::FlowLawFactory ice_factory("stress_balance.sia.", m_config, m_EC);
33  m_flow_law = ice_factory.create();
34 }
35 
37  m_log->message(2, "* Initializing Weertman-style basal sliding...\n");
38 }
39 
40 
41 /*!
42  * Compute sliding velocity using a Weertman-style parameterization from [@ref Tomkin2007],
43  * equation 5.
44  *
45  * @f[ u_s = \frac{2 A_s \beta_c (\rho g H)^{n}}{N - P} \cdot |\nabla h|^{n-1} \cdot \nabla h,
46  * @f]
47  *
48  * where
49  *
50  * - @f$ A_s @f$ is the sliding parameter,
51  * - @f$ \beta_c @f$ is a parameter capturing the effect of the presence of valley walls
52  * (set to 1 in this implementation),
53  * - @f$ \rho @f$ is the ice density,
54  * - @f$ H @f$ is the ice thickness,
55  * - @f$ h @f$ is the ice surface elevation ,
56  * - @f$ n @f$ is the flow law exponent (usually 3),
57  * - @f$ g @f$ is the acceleration due to gravity,
58  * - @f$ N @f$ is the ice overburden pressure,
59  * - @f$ P @f$ is the basal water pressure.
60  *
61  * With these modifications and noting that @f$ N = \rho g H @f$, the formula above
62  * becomes
63  *
64  * @f[ u_s = \frac{2 A_s}{1 - k} \cdot (N |\nabla h|)^{n-1} \cdot \nabla h,
65  * @f]
66  *
67  * where
68  * - @f$ N = \rho g H @f$,
69  * - @f$ P = k N @f$ (we assume that basal water pressure is a given fraction of overburden)
70  *
71  * This parameterization is used for areas of grounded ice where the base of the ice is
72  * temperate.
73  */
74 void WeertmanSliding::update(const Inputs &inputs, bool full_update) {
75 
76  (void) full_update;
77 
78  const array::Scalar &H = inputs.geometry->ice_thickness;
79  const array::Scalar &h = inputs.geometry->ice_surface_elevation;
80  const auto &cell_type = inputs.geometry->cell_type;
81  const array::Array3D &enthalpy = *inputs.enthalpy;
82 
83  double n = m_flow_law->exponent();
84  double A_s = m_config->get_number("stress_balance.weertman_sliding.A", "Pa-3 s-1 m-2");
85  double k = m_config->get_number("stress_balance.weertman_sliding.k");
86 
87  array::AccessScope list{&m_velocity, &H, &h, &enthalpy, &cell_type};
88 
89  ParallelSection loop(m_grid->com);
90  try {
91  for (auto p = m_grid->points(); p; p.next()) {
92  const int i = p.i(), j = p.j();
93 
94  double
95  P_o = m_EC->pressure(H(i, j)),
96  E_base = enthalpy.get_column(i, j)[0];
97 
98  if (not m_EC->is_temperate(E_base, P_o) or cell_type.ocean(i, j)) {
99  m_velocity(i, j) = 0.0;
100  continue;
101  }
102 
103  // Note: we may need to decide if we should use one-sided FD at ice margins.
104  Vector2d grad_h = {diff_x_p(h, i, j), diff_y_p(h, i, j)};
105 
106  // Note: this could be optimized by computing this instead
107  // 2 * A_s / (1 - k) * pow(P * P * (h_x * h_x + h_y * h_y), (n - 1) / 2) * grad_h;
108  // ... but I'm not sure we need to and the current code is cleaner.
109  m_velocity(i, j) = - 2.0 * A_s / (1.0 - k) * pow(P_o * grad_h.magnitude(), n - 1) * grad_h;
110  }
111  } catch (...) {
112  loop.failed();
113  }
114  loop.check();
115 
116 }
117 
118 } // end of namespace stressbalance
119 } // end of namespace pism
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:158
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:162
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition: Component.hh:156
array::Scalar2 ice_surface_elevation
Definition: Geometry.hh:57
array::CellType2 cell_type
Definition: Geometry.hh:55
array::Scalar2 ice_thickness
Definition: Geometry.hh:51
void failed()
Indicates a failure of a parallel section.
double magnitude() const
Magnitude.
Definition: Vector2d.hh:40
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2d.hh:29
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
double * get_column(int i, int j)
Definition: Array3D.cc:120
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: Array3D.hh:33
std::shared_ptr< FlowLaw > create()
const array::Array3D * enthalpy
std::shared_ptr< rheology::FlowLaw > m_flow_law
Shallow stress balance (such as the SSA).
WeertmanSliding(std::shared_ptr< const Grid > g)
virtual void update(const Inputs &inputs, bool full_update)
#define n
Definition: exactTestM.c:37
double diff_x_p(const array::Scalar &array, int i, int j)
Returns the x-derivative at i,j approximated using centered finite differences. Respects grid periodi...
Definition: Scalar.cc:108
double diff_y_p(const array::Scalar &array, int i, int j)
Returns the y-derivative at i,j approximated using centered finite differences. Respects grid periodi...
Definition: Scalar.cc:129
static const double k
Definition: exactTestP.cc:42