20 #include "pism/geometry/flux_limiter.hh"
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"
48 static inline double pp(
double x) {
55 static inline double np(
double x) {
61 return dt * ((
pp(u.
e) +
np(u.
w)) / dx + (
pp(u.
n) +
np(u.
s)) / dy);
84 auto grid = result.
grid();
87 eps = std::numeric_limits<double>::epsilon(),
95 return (Q.e - Q.w) / dx + (Q.n - Q.s) / dy;
98 int limiter_count = 0;
100 for (
auto p = grid->points(); p; p.next()) {
101 const int i = p.i(), j = p.j();
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);
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)) {
120 result(i, j, 0) = Q.e;
121 result(i, j, 1) = Q.n;
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);
139 double F_e = Q.e * dt / dx;
140 double F_n = Q.n * dt / dy;
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);
159 (-
np(F_e) / F_out_e) * X_e);
161 assert(x(i, j) - F_e_limited >= 0);
162 assert(x(i + 1, j) + F_e_limited >= 0);
165 (-
np(F_n) / F_out_n) * X_n);
167 assert(x(i, j) - F_n_limited >= 0);
168 assert(x(i, j + 1) + F_n_limited >= 0);
171 result(i, j, 0) = F_e_limited * dx / dt;
172 result(i, j, 1) = F_n_limited * dy / dt;
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);
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
std::shared_ptr< const Grid > grid() const
stencils::Star< double > star(int i, int j) const
Returns the values at interfaces of the cell i,j using the staggered grid.
A class for storing and accessing internal staggered-grid 2D fields. Uses dof=2 storage....
double max(const array::Scalar &input)
Finds maximum over all the values in an array::Scalar object. Ignores ghosts.
double min(const array::Scalar &input)
Finds minimum over all the values in an array::Scalar object. Ignores ghosts.
static double pp(double x)
static double np(double x)
static double flux_out(const stencils::Star< double > &u, double dx, double dy, double dt)
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)
Star stencil points (in the map-plane).