PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
flux_limiter.cc
Go to the documentation of this file.
1 /* Copyright (C) 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/geometry/flux_limiter.hh"
21 
22 #include <cassert>
23 #include <algorithm>
24 #include <limits>
25 
26 #include "pism/util/array/Scalar.hh"
27 #include "pism/util/array/Staggered.hh"
28 #include "pism/util/Context.hh"
29 #include "pism/util/Logger.hh"
30 #include "pism/util/pism_utilities.hh" // GlobalSum()
31 
32 namespace pism {
33 
34 /*! References:
35  *
36  * [Smolarkiewicz1989] P. K. Smolarkiewicz, “Comment on “A Positive Definite Advection
37  * Scheme Obtained by Nonlinear Renormalization of the Advective Fluxes”,” Monthly Weather
38  * Review, vol. 117, Art. no. 11, 1989.
39  *
40  * [Zalesak1979] S. T. Zalesak, “Fully multidimensional flux-corrected transport
41  * algorithms for fluids,” Journal of Computational Physics, vol. 31, Art. no. 3, Jun.
42  * 1979.
43  */
44 
45 namespace details {
46 
47 // positive part
48 static inline double pp(double x) {
49  return std::max(x, 0.0);
50 }
51 
52 // negative part
53 //
54 // Note the negative sign!
55 static inline double np(double x) {
56  return -std::min(x, 0.0);
57 }
58 
59 // Total flux out of a cell over a time step dt
60 static inline double flux_out(const stencils::Star<double> &u, double dx, double dy, double dt) {
61  return dt * ((pp(u.e) + np(u.w)) / dx + (pp(u.n) + np(u.s)) / dy);
62 }
63 
64 } // end of namespace details
65 
66 /*! Limit fluxes to preserve non-negativity of a transported quantity.
67  *
68  * This method is described in [Smolarkiewicz1989].
69  *
70  * It is based on the [Zalesak1979] flux corrected transport limiter, but for the
71  * "regular" flux instead of the "anti-diffusive" flux and with a different limiting
72  * criterion (non-negativity instead of monotonicity).
73  *
74  */
76  const array::Scalar1 &x,
77  const array::Staggered1 &flux,
78  array::Staggered &result) {
79 
80  using details::pp;
81  using details::np;
82  using details::flux_out;
83 
84  auto grid = result.grid();
85 
86  double
87  eps = std::numeric_limits<double>::epsilon(),
88  dx = grid->dx(),
89  dy = grid->dy();
90 
91  array::AccessScope list{&flux, &x, &result};
92 
93  // flux divergence
94  auto div = [dx, dy](const stencils::Star<double> &Q) {
95  return (Q.e - Q.w) / dx + (Q.n - Q.s) / dy;
96  };
97 
98  int limiter_count = 0;
99 
100  for (auto p = grid->points(); p; p.next()) {
101  const int i = p.i(), j = p.j();
102 
103  auto Q = flux.star(i, j);
104  auto Q_n = flux.star(i, j + 1);
105  auto Q_e = flux.star(i + 1, j);
106 
107  const double
108  div_Q = div(Q),
109  div_Q_e = div(Q_e),
110  div_Q_n = div(Q_n);
111 
112  if ((div_Q <= 0.0 or x(i, j) - dt * div_Q >= eps) and
113  (div_Q_e <= 0.0 or x(i + 1, j) - dt * div_Q_e >= eps) and
114  (div_Q_n <= 0.0 or x(i, j + 1) - dt * div_Q_n >= eps)) {
115  // No need to limit fluxes: total fluxes out of cells (i, j), (i + 1, j), (i, j + 1)
116  // may be able to create a negative thickness, but fluxes *into* these cells make up for it
117  //
118  // Without this check the limiter is a little over-eager and may affect results in
119  // areas where mass conservation is not an issue.
120  result(i, j, 0) = Q.e;
121  result(i, j, 1) = Q.n;
122  continue;
123  }
124 
125  limiter_count += 1;
126 
127  // compute total amounts moved *out* of the current cell and its north and east
128  // neighbors over the course of the time step dt
129  //
130  // see equation (A4) in [Smolarkiewicz1989]
131  //
132  // note that we can compute all these using the width=1 stencil because of the way
133  // PISM's staggered grid is set up
134  double F_out = flux_out(Q, dx, dy, dt);
135  double F_out_n = flux_out(Q_n, dx, dy, dt);
136  double F_out_e = flux_out(Q_e, dx, dy, dt);
137 
138  // amounts moved through the eastern and northern cell faces
139  double F_e = Q.e * dt / dx;
140  double F_n = Q.n * dt / dy;
141 
142  // Maximum amounts the current cell and its neighbors can lose while maintaining
143  // non-negativity
144  //
145  // Note: we limit total amounts so that
146  //
147  // - if a cell value X is below eps, the flux is zero
148  //
149  // - otherwise the total flux out of a cell can remove at most (X - eps) over the
150  // course of a time step
151  //
152  // This is needed to avoid small negative values resulting from rounding errors.
153  double X_ij = pp(x(i, j) - eps);
154  double X_e = pp(x(i + 1, j) - eps);
155  double X_n = pp(x(i, j + 1) - eps);
156 
157  // limit total amounts (see equation (10) in [Smolarkiewicz1989])
158  double F_e_limited = std::max(std::min(F_e, (pp(F_e) / F_out) * X_ij),
159  (-np(F_e) / F_out_e) * X_e);
160 
161  assert(x(i, j) - F_e_limited >= 0);
162  assert(x(i + 1, j) + F_e_limited >= 0);
163 
164  double F_n_limited = std::max(std::min(F_n, (pp(F_n) / F_out) * X_ij),
165  (-np(F_n) / F_out_n) * X_n);
166 
167  assert(x(i, j) - F_n_limited >= 0);
168  assert(x(i, j + 1) + F_n_limited >= 0);
169 
170  // convert back to fluxes:
171  result(i, j, 0) = F_e_limited * dx / dt;
172  result(i, j, 1) = F_n_limited * dy / dt;
173  }
174 
175  limiter_count = GlobalSum(grid->com, limiter_count);
176  if (limiter_count > 0) {
177  auto log = grid->ctx()->log();
178  log->message(2, "limited ice flux at %d locations\n", limiter_count);
179  }
180 }
181 
182 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
Definition: Array.hh:65
std::shared_ptr< const Grid > grid() const
Definition: Array.cc:132
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
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 double pp(double x)
Definition: flux_limiter.cc:48
static double np(double x)
Definition: flux_limiter.cc:55
static double flux_out(const stencils::Star< double > &u, double dx, double dy, double dt)
Definition: flux_limiter.cc:60
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
void make_nonnegative_preserving(double dt, const array::Scalar1 &x, const array::Staggered1 &flux, array::Staggered &result)
Definition: flux_limiter.cc:75
Star stencil points (in the map-plane).
Definition: stencils.hh:30