tWeertmanSliding.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
---
tWeertmanSliding.cc (4100B)
---
1 /* Copyright (C) 2018 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 "WeertmanSliding.hh"
21
22 #include "pism/rheology/FlowLawFactory.hh"
23 #include "pism/geometry/Geometry.hh"
24 #include "StressBalance.hh"
25
26 namespace pism {
27 namespace stressbalance {
28
29 WeertmanSliding::WeertmanSliding(IceGrid::ConstPtr 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
36 WeertmanSliding::~WeertmanSliding() {
37 // empty
38 }
39
40 void WeertmanSliding::init_impl() {
41 m_log->message(2, "* Initializing Weertman-style basal sliding...\n");
42 }
43
44
45 /*!
46 * Compute sliding velocity using a Weertman-style parameterization from [@ref Tomkin2007],
47 * equation 5.
48 *
49 * @f[ u_s = \frac{2 A_s \beta_c (\rho g H)^{n}}{N - P} \cdot |\nabla h|^{n-1} \cdot \nabla h,
50 * @f]
51 *
52 * where
53 *
54 * - @f$ A_s @f$ is the sliding parameter,
55 * - @f$ \beta_c @f$ is a parameter capturing the effect of the presence of valley walls
56 * (set to 1 in this implementation),
57 * - @f$ \rho @f$ is the ice density,
58 * - @f$ H @f$ is the ice thickness,
59 * - @f$ h @f$ is the ice surface elevation ,
60 * - @f$ n @f$ is the flow law exponent (usually 3),
61 * - @f$ g @f$ is the acceleration due to gravity,
62 * - @f$ N @f$ is the ice overburden pressure,
63 * - @f$ P @f$ is the basal water pressure.
64 *
65 * With these modifications and noting that @f$ N = \rho g H @f$, the formula above
66 * becomes
67 *
68 * @f[ u_s = \frac{2 A_s}{1 - k} \cdot (N |\nabla h|)^{n-1} \cdot \nabla h,
69 * @f]
70 *
71 * where
72 * - @f$ N = \rho g H @f$,
73 * - @f$ P = k N @f$ (we assume that basal water pressure is a given fraction of overburden)
74 *
75 * This parameterization is used for areas of grounded ice where the base of the ice is
76 * temperate.
77 */
78 void WeertmanSliding::update(const Inputs &inputs, bool full_update) {
79
80 (void) full_update;
81
82 const IceModelVec2S &H = inputs.geometry->ice_thickness;
83 const IceModelVec2S &h = inputs.geometry->ice_surface_elevation;
84 const IceModelVec2CellType &cell_type = inputs.geometry->cell_type;
85 const IceModelVec3 &enthalpy = *inputs.enthalpy;
86
87 double n = m_flow_law->exponent();
88 double A_s = m_config->get_number("stress_balance.weertman_sliding.A");
89 double k = m_config->get_number("stress_balance.weertman_sliding.k");
90
91 IceModelVec::AccessList list{&m_velocity, &H, &h, &enthalpy, &cell_type};
92
93 ParallelSection loop(m_grid->com);
94 try {
95 for (Points p(*m_grid); p; p.next()) {
96 const int i = p.i(), j = p.j();
97
98 double
99 P_o = m_EC->pressure(H(i, j)),
100 E_base = enthalpy(i, j, 0);
101
102 if (not m_EC->is_temperate(E_base, P_o) or cell_type.ocean(i, j)) {
103 m_velocity(i, j) = {0.0, 0.0};
104 continue;
105 }
106
107 // Note: we may need to decide if we should use one-sided FD at ice margins.
108 Vector2 grad_h = {h.diff_x_p(i, j), h.diff_y_p(i, j)};
109
110 // Note: this could be optimized by computing this instead
111 // 2 * A_s / (1 - k) * pow(P * P * (h_x * h_x + h_y * h_y), (n - 1) / 2) * grad_h;
112 // ... but I'm not sure we need to and the current code is cleaner.
113 m_velocity(i, j) = 2.0 * A_s / (1.0 - k) * pow(P_o * grad_h.magnitude(), n - 1) * grad_h;
114 }
115 } catch (...) {
116 loop.failed();
117 }
118 loop.check();
119
120 }
121
122 } // end of namespace stressbalance
123 } // end of namespace pism