PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
SSAFD_Regional.cc
Go to the documentation of this file.
1 /* Copyright (C) 2015, 2016, 2017, 2020, 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/regional/SSAFD_Regional.hh"
21 #include "pism/util/Vars.hh"
22 #include "pism/stressbalance/StressBalance.hh"
23 
24 namespace pism {
25 
26 namespace stressbalance {
27 
28 SSAFD_Regional::SSAFD_Regional(std::shared_ptr<const Grid> g)
29  : SSAFD(g) {
30 
31  m_h_stored = NULL;
32  m_H_stored = NULL;
33  m_no_model_mask = NULL;
34 }
35 
37 
38  SSAFD::init();
39 
40  m_log->message(2, " using the regional version of the SSA solver...\n");
41 
42  if (m_config->get_flag("stress_balance.ssa.dirichlet_bc")) {
43  m_log->message(2, " using stored SSA velocities as Dirichlet B.C. in the no_model_strip...\n");
44  }
45 }
46 
47 void SSAFD_Regional::update(const Inputs &inputs, bool full_update) {
51 
52  SSA::update(inputs, full_update);
53 
54  m_h_stored = NULL;
55  m_H_stored = NULL;
56  m_no_model_mask = NULL;
57 }
58 
59 static int weight(int M_ij, int M_n, double h_ij, double h_n) {
60  // fjord walls, nunataks, headwalls
61  if ((mask::icy(M_ij) and mask::ice_free(M_n) and h_n > h_ij) or
62  (mask::ice_free(M_ij) and mask::icy(M_n) and h_ij > h_n)) {
63  return 0;
64  }
65 
66  return 1;
67 }
68 
70  const array::Scalar1 &surface_elevation,
71  const array::CellType1 &cell_type,
72  const array::Scalar1 *no_model_mask,
73  array::Vector &result) const {
74 
75  SSAFD::compute_driving_stress(ice_thickness, surface_elevation, cell_type, no_model_mask, result);
76 
77  double
78  dx = m_grid->dx(),
79  dy = m_grid->dy();
80 
81  int
82  Mx = m_grid->Mx(),
83  My = m_grid->My();
84 
85  array::AccessScope list{&result, &cell_type, no_model_mask, m_h_stored, m_H_stored};
86 
87  for (auto p = m_grid->points(); p; p.next()) {
88  const int i = p.i(), j = p.j();
89 
90  auto M = no_model_mask->star(i, j);
91 
92  if (M.c == 0) {
93  // this grid point is in the modeled area so we don't need to modify the driving
94  // stress
95  continue;
96  }
97 
98  double pressure = m_EC->pressure((*m_H_stored)(i, j));
99  if (pressure <= 0) {
100  result(i, j) = 0.0;
101  continue;
102  }
103 
104  auto h = m_h_stored->star(i, j);
105  auto CT = cell_type.star(i, j);
106 
107  // x-derivative
108  double h_x = 0.0;
109  {
110  double
111  west = M.w == 1 and i > 0,
112  east = M.e == 1 and i < Mx - 1;
113 
114  // don't use differences spanning "cliffs"
115  west *= weight(CT.c, CT.w, h.c, h.w);
116  east *= weight(CT.c, CT.e, h.c, h.e);
117 
118  if (east + west > 0) {
119  h_x = 1.0 / ((west + east) * dx) * (west * (h.c - h.w) + east * (h.e - h.c));
120  } else {
121  h_x = 0.0;
122  }
123  }
124 
125  // y-derivative
126  double h_y = 0.0;
127  {
128  double
129  south = M.s == 1 and j > 0,
130  north = M.n == 1 and j < My - 1;
131 
132  // don't use differences spanning "cliffs"
133  south *= weight(CT.c, CT.s, h.c, h.s);
134  north *= weight(CT.c, CT.n, h.c, h.n);
135 
136  if (north + south > 0) {
137  h_y = 1.0 / ((south + north) * dy) * (south * (h.c - h.s) + north * (h.n - h.c));
138  } else {
139  h_y = 0.0;
140  }
141  }
142 
143  result(i, j) = - pressure * Vector2d(h_x, h_y);
144  } // end of the loop over grid points
145 }
146 
147 SSA * SSAFD_RegionalFactory(std::shared_ptr<const Grid> grid) {
148  return new SSAFD_Regional(grid);
149 }
150 
151 } // end of namespace stressbalance
152 
153 } // 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
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
stencils::Star< T > star(int i, int j) const
Definition: Array2D.hh:79
const array::Scalar * no_model_ice_thickness
const array::Scalar2 * no_model_surface_elevation
const array::Scalar2 * no_model_mask
SSAFD_Regional(std::shared_ptr< const Grid > g)
const array::Scalar * m_no_model_mask
void update(const Inputs &inputs, bool full_update)
Update the SSA solution.
const array::Scalar1 * m_h_stored
virtual void compute_driving_stress(const array::Scalar &ice_thickness, const array::Scalar1 &surface_elevation, const array::CellType1 &cell_type, const array::Scalar1 *no_model_mask, array::Vector &result) const
Compute the gravitational driving stress.
A version of the SSA stress balance with tweaks for outlet glacier simulations.
PISM's SSA solver: the finite difference implementation.
Definition: SSAFD.hh:35
virtual void update(const Inputs &inputs, bool full_update)
Update the SSA solution.
Definition: SSA.cc:152
virtual void compute_driving_stress(const array::Scalar &ice_thickness, const array::Scalar1 &surface_elevation, const array::CellType1 &cell_type, const array::Scalar1 *no_model_mask, array::Vector &result) const
Compute the gravitational driving stress.
Definition: SSA.cc:232
PISM's SSA solver.
Definition: SSA.hh:110
bool icy(int M)
Ice-filled cell (grounded or floating).
Definition: Mask.hh:48
bool ice_free(int M)
Ice-free cell (grounded or ocean).
Definition: Mask.hh:58
static int weight(int M_ij, int M_n, double h_ij, double h_n)
SSA * SSAFD_RegionalFactory(std::shared_ptr< const Grid > grid)
static const double g
Definition: exactTestP.cc:36