PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
UNO.cc
Go to the documentation of this file.
1 /* Copyright (C) 2020, 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 <cmath>
21 
22 #include "pism/geometry/UNO.hh"
23 #include "pism/geometry/flux_limiter.hh"
24 
25 #include "pism/util/array/CellType.hh"
26 
27 /*! References:
28  *
29  * [Li2008] J.-G. Li, “Upstream Nonoscillatory Advection Schemes,” Monthly Weather Review,
30  * vol. 136, Art. no. 12, Dec. 2008.
31  */
32 
33 namespace pism {
34 
35 // 1D Upwind Non-Oscillatory (UNO) advection schemes and related approximations
36 namespace uno {
37 
38 // the difference quotient (equation 2 in [Li2008])
39 static inline double difference_quotient(const double *x, const double *f, size_t A, size_t B) {
40  return (f[A] - f[B]) / (x[A] - x[B]);
41 }
42 
43 // index of the "upstream" cell
44 static inline size_t upstream(double v, size_t j) {
45  return v > 0.0 ? j - 1 : j + 2;
46 }
47 
48 // index of the "central" cell
49 static inline size_t central(double v, size_t j) {
50  return v > 0.0 ? j : j + 1;
51 }
52 
53 // index of the "downstream" cell
54 static inline size_t downstream(double v, size_t j) {
55  return v > 0.0 ? j + 1 : j;
56 }
57 
58 // the sign
59 static inline double sign(double x) {
60  return x >= 0.0 ? 1.0 : -1.0;
61 }
62 
63 // first order upwinding
64 static inline double psi_upwind1(const double *x,
65  const double *y,
66  size_t j,
67  double velocity,
68  double dx,
69  double dt) {
70  (void) x;
71  (void) dx;
72  (void) dt;
73 
74  return y[central(velocity, j)];
75 }
76 
77 // Lax-Wendroff
78 static inline double psi_lax_wendroff(const double *x,
79  const double *y,
80  size_t j,
81  double velocity,
82  double dx,
83  double dt) {
84  size_t
85  c = central(velocity, j),
86  d = downstream(velocity, j);
87 
88  double gc = difference_quotient(x, y, d, c);
89 
90  // equation 3 in [Li2008]
91  return y[c] + 0.5 * sign(velocity) * (dx - std::abs(velocity) * dt) * gc;
92 }
93 
94 // Fromm
95 static inline double psi_fromm(const double *x,
96  const double *y,
97  size_t j,
98  double velocity,
99  double dx,
100  double dt) {
101  size_t
102  u = upstream(velocity, j),
103  c = central(velocity, j),
104  d = downstream(velocity, j);
105 
106  double gc = difference_quotient(x, y, d, u);
107 
108  // equation 3 in [Li2008]
109  return y[c] + 0.5 * sign(velocity) * (dx - std::abs(velocity) * dt) * gc;
110 }
111 
112 // UNO2
113 static inline double psi_uno2(const double *x,
114  const double *y,
115  size_t j,
116  double velocity,
117  double dx,
118  double dt) {
119  size_t
120  u = upstream(velocity, j),
121  c = central(velocity, j),
122  d = downstream(velocity, j);
123 
124  double
125  gdc = difference_quotient(x, y, d, c),
126  gcu = difference_quotient(x, y, c, u);
127 
128  // equation 5 in [Li2008]
129  double gc = sign(gdc) * std::min(std::abs(gdc), std::abs(gcu));
130 
131  // equation 3 in [Li2008]
132  return y[c] + 0.5 * sign(velocity) * (dx - std::abs(velocity) * dt) * gc;
133 }
134 
135 // UNO3
136 static inline double psi_uno3(const double *x,
137  const double *y,
138  size_t j,
139  double velocity,
140  double dx,
141  double dt) {
142  size_t
143  u = upstream(velocity, j),
144  c = central(velocity, j),
145  d = downstream(velocity, j);
146 
147  double
148  gdc = difference_quotient(x, y, d, c),
149  gcu = difference_quotient(x, y, c, u),
150  gdu = difference_quotient(x, y, d, u);
151 
152  double abs_u = std::abs(velocity);
153  double sgn_u = sign(velocity);
154 
155  // equation 10 in [Li2008]
156  double gc = 0.0;
157  // NOLINTBEGIN(readability-magic-numbers)
158  if (std::abs(gdc - gcu) < 1.2 * std::abs(gdu)) {
159  gc = gdc - ((dx + abs_u * dt) / (1.5 * sgn_u)) * ((gdc - gcu) / (x[d] - x[u]));
160  } else if (gdc * gcu > 0.0) {
161  gc = 2 * sign(gdc) * std::min(std::abs(gdc), std::abs(gcu));
162  } else {
163  gc = sign(gdc) * std::min(std::abs(gdc), std::abs(gcu)); // UNO2
164  }
165  // NOLINTEND(readability-magic-numbers)
166 
167  return y[c] + 0.5 * sgn_u * (dx - abs_u * dt) * gc;
168 }
169 
170 } // end of namespace uno
171 
172 const array::Scalar& UNO::x() const {
173  return m_x;
174 }
175 
176 UNO::UNO(std::shared_ptr<const Grid> grid, UNOType type)
177  : m_q(grid, "interface_fluxes"),
178  m_q_limited(grid, "limited_interface_fluxes"),
179  m_v_ghosted(grid, "velocity"),
180  m_x_ghosted(grid, "old_state"),
181  m_x(grid, "new_state")
182 {
183  switch (type) {
184  case PISM_UNO_UPWIND1:
186  break;
189  break;
190  case PISM_UNO_FROMM:
192  break;
193  case PISM_UNO_2:
195  break;
196  default:
197  case PISM_UNO_3:
199  break;
200  }
201 }
202 
203 /*!
204  * Compute staggered grid (cell interface) fluxes given the domain extent (cell type),
205  * regular grid velocities, and the current distribution of the advected quantity.
206  */
208  const array::Vector1 &velocity,
209  const array::Scalar2 &x_old,
210  double dt,
211  array::Staggered &result) const {
212 
213  auto grid = cell_type.grid();
214 
215  double
216  dx = grid->dx(),
217  dy = grid->dy();
218 
219  array::AccessScope scope{&cell_type, &velocity, &x_old, &result};
220 
221  // temporary storage for values needed by stencil computations
222  double coords[4];
223  double values[4];
224 
225  for (auto p = grid->points(); p; p.next()) {
226  const int i = p.i(), j = p.j();
227 
228  // velocities through east and north cell interfaces:
229  double vx, vy;
230  {
231  Vector2d
232  V = velocity(i, j),
233  V_e = velocity(i + 1, j),
234  V_n = velocity(i, j + 1);
235 
236  double
237  W = static_cast<double>(cell_type.icy(i, j)),
238  W_n = static_cast<double>(cell_type.icy(i, j + 1)),
239  W_e = static_cast<double>(cell_type.icy(i + 1, j));
240 
241  vx = (W * V.u + W_e * V_e.u) / std::max(W + W_e, 1.0);
242  vy = (W * V.v + W_n * V_n.v) / std::max(W + W_n, 1.0);
243  }
244 
245  // east
246  {
247  // gather inputs to the 1D UNO code: we need four numbers centered around i+1/2,
248  // i.e. [i-1, i, i+1, i+2].
249  double x0 = i > 0 ? grid->x(i - 1) : grid->x(0) - dx;
250  for (int k = 0; k < 4; ++k) {
251  coords[k] = x0 + static_cast<double>(k) * dx;
252  values[k] = x_old((i - 1) + k, j);
253  }
254 
255  // use the 1D "mid-flux" approximation to get the flux
256  result(i, j, 0) = vx * m_approx(coords, values, 1, vx, dx, dt);
257  }
258 
259  // north
260  {
261  // gather inputs to the 1D UNO code: we need four numbers centered around j+1/2,
262  // i.e. [j-1, j, j+1, j+2].
263  double y0 = j > 0 ? grid->y(j - 1) : grid->y(0) - dy;
264  for (int k = 0; k < 4; ++k) {
265  coords[k] = y0 + k * dy;
266  values[k] = x_old(i, (j - 1) + k);
267  }
268 
269  // use the 1D "mid-flux" approximation to get the flux
270  result(i, j, 1) = vy * m_approx(coords, values, 1, vy, dy, dt);
271  }
272  } // end of the loop over grid points
273 }
274 
275 /*!
276  * Perform an explicit step using the flux form of the advection equation.
277  *
278  * @param[in] dt time step length
279  * @param[in] flux cell interface fluxes
280  * @param[in] x_old current state
281  * @param[out] result new state
282  */
283 static void step(double dt,
284  const array::Staggered1 &flux,
285  const array::Scalar &x_old,
286  array::Scalar &result) {
287 
288  auto grid = result.grid();
289 
290  double
291  dx = grid->dx(),
292  dy = grid->dy();
293 
294  array::AccessScope scope{&flux, &x_old, &result};
295 
296  for (auto p = grid->points(); p; p.next()) {
297  const int i = p.i(), j = p.j();
298 
299  const auto Q = flux.star(i, j);
300 
301  result(i, j) = x_old(i, j) - dt * ((Q.e - Q.w) / dx + (Q.n - Q.s) / dy);
302  }
303 }
304 
305 void UNO::update(double dt,
306  const array::CellType1 &cell_type,
307  const array::Scalar &x,
308  const array::Vector &velocity,
309  bool nonnegative) {
310 
311  // make ghosted copies:
312  m_v_ghosted.copy_from(velocity);
314 
316 
317  m_q.update_ghosts();
318 
319  // limit fluxes to preserve non-negativity
320  if (nonnegative) {
323  }
324 
325  step(dt, m_q, x, m_x);
326 }
327 
328 } // end of namespace pism
void compute_interface_fluxes(const array::CellType1 &cell_type, const array::Vector1 &velocity, const array::Scalar2 &x_old, double dt, array::Staggered &result) const
Definition: UNO.cc:207
array::Staggered1 m_q
Definition: UNO.hh:71
const array::Scalar & x() const
Definition: UNO.cc:172
void update(double dt, const array::CellType1 &cell_type, const array::Scalar &x, const array::Vector &velocity, bool nonnegative=false)
Definition: UNO.cc:305
MidFluxApproximation m_approx
Definition: UNO.hh:68
array::Staggered1 m_q_limited
Definition: UNO.hh:71
array::Scalar m_x
Definition: UNO.hh:83
array::Scalar2 m_x_ghosted
Definition: UNO.hh:80
array::Vector1 m_v_ghosted
Definition: UNO.hh:74
UNO(std::shared_ptr< const Grid > grid, UNOType type)
Definition: UNO.cc:176
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
void copy_from(const Array2D< T > &source)
Definition: Array2D.hh:73
std::shared_ptr< const Grid > grid() const
Definition: Array.cc:132
void update_ghosts()
Updates ghost points.
Definition: Array.cc:693
bool icy(int i, int j) const
Definition: CellType.hh:42
stencils::Star< double > star(int i, int j) const
Returns the values at interfaces of the cell i,j using the staggered grid.
Definition: Staggered.hh:73
void copy_from(const array::Staggered &input)
Definition: Staggered.cc:42
A class for storing and accessing internal staggered-grid 2D fields. Uses dof=2 storage....
Definition: Staggered.hh:35
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:165
double min(const array::Scalar &input)
Finds minimum over all the values in an array::Scalar object. Ignores ghosts.
Definition: Scalar.cc:193
static size_t downstream(double v, size_t j)
Definition: UNO.cc:54
static double psi_lax_wendroff(const double *x, const double *y, size_t j, double velocity, double dx, double dt)
Definition: UNO.cc:78
static double difference_quotient(const double *x, const double *f, size_t A, size_t B)
Definition: UNO.cc:39
static double psi_fromm(const double *x, const double *y, size_t j, double velocity, double dx, double dt)
Definition: UNO.cc:95
static size_t upstream(double v, size_t j)
Definition: UNO.cc:44
static double psi_uno3(const double *x, const double *y, size_t j, double velocity, double dx, double dt)
Definition: UNO.cc:136
static size_t central(double v, size_t j)
Definition: UNO.cc:49
static double psi_uno2(const double *x, const double *y, size_t j, double velocity, double dx, double dt)
Definition: UNO.cc:113
static double psi_upwind1(const double *x, const double *y, size_t j, double velocity, double dx, double dt)
Definition: UNO.cc:64
static double sign(double x)
Definition: UNO.cc:59
static void step(double dt, const array::Staggered1 &velocity, const array::Scalar1 &x_old, array::Scalar &x)
Definition: MPDATA2.cc:281
static const double k
Definition: exactTestP.cc:42
void make_nonnegative_preserving(double dt, const array::Scalar1 &x, const array::Staggered1 &flux, array::Staggered &result)
Definition: flux_limiter.cc:75
UNOType
Definition: UNO.hh:35
@ PISM_UNO_2
Definition: UNO.hh:35
@ PISM_UNO_LAX_WENDROFF
Definition: UNO.hh:35
@ PISM_UNO_UPWIND1
Definition: UNO.hh:35
@ PISM_UNO_3
Definition: UNO.hh:35
@ PISM_UNO_FROMM
Definition: UNO.hh:35