22 #include "pism/geometry/UNO.hh"
23 #include "pism/geometry/flux_limiter.hh"
25 #include "pism/util/array/CellType.hh"
40 return (f[A] - f[B]) / (x[A] - x[B]);
44 static inline size_t upstream(
double v,
size_t j) {
45 return v > 0.0 ? j - 1 : j + 2;
49 static inline size_t central(
double v,
size_t j) {
50 return v > 0.0 ? j : j + 1;
55 return v > 0.0 ? j + 1 : j;
59 static inline double sign(
double x) {
60 return x >= 0.0 ? 1.0 : -1.0;
91 return y[c] + 0.5 *
sign(velocity) * (dx - std::abs(velocity) * dt) * gc;
109 return y[c] + 0.5 *
sign(velocity) * (dx - std::abs(velocity) * dt) * gc;
129 double gc =
sign(gdc) *
std::min(std::abs(gdc), std::abs(gcu));
132 return y[c] + 0.5 *
sign(velocity) * (dx - std::abs(velocity) * dt) * gc;
152 double abs_u = std::abs(velocity);
153 double sgn_u =
sign(velocity);
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));
163 gc =
sign(gdc) *
std::min(std::abs(gdc), std::abs(gcu));
167 return y[c] + 0.5 * sgn_u * (dx - abs_u * dt) * gc;
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")
213 auto grid = cell_type.
grid();
225 for (
auto p = grid->points(); p; p.next()) {
226 const int i = p.i(), j = p.j();
233 V_e = velocity(i + 1, j),
234 V_n = velocity(i, j + 1);
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));
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);
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);
256 result(i, j, 0) = vx *
m_approx(coords, values, 1, vx, dx, dt);
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);
270 result(i, j, 1) = vy *
m_approx(coords, values, 1, vy, dy, dt);
288 auto grid = result.
grid();
296 for (
auto p = grid->points(); p; p.next()) {
297 const int i = p.i(), j = p.j();
299 const auto Q = flux.
star(i, j);
301 result(i, j) = x_old(i, j) - dt * ((Q.e - Q.w) / dx + (Q.n - Q.s) / dy);
void compute_interface_fluxes(const array::CellType1 &cell_type, const array::Vector1 &velocity, const array::Scalar2 &x_old, double dt, array::Staggered &result) const
const array::Scalar & x() const
void update(double dt, const array::CellType1 &cell_type, const array::Scalar &x, const array::Vector &velocity, bool nonnegative=false)
MidFluxApproximation m_approx
array::Staggered1 m_q_limited
array::Scalar2 m_x_ghosted
array::Vector1 m_v_ghosted
UNO(std::shared_ptr< const Grid > grid, UNOType type)
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Makes sure that we call begin_access() and end_access() for all accessed array::Arrays.
void copy_from(const Array2D< T > &source)
std::shared_ptr< const Grid > grid() const
void update_ghosts()
Updates ghost points.
bool icy(int i, int j) const
stencils::Star< double > star(int i, int j) const
Returns the values at interfaces of the cell i,j using the staggered grid.
void copy_from(const array::Staggered &input)
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 size_t downstream(double v, size_t j)
static double psi_lax_wendroff(const double *x, const double *y, size_t j, double velocity, double dx, double dt)
static double difference_quotient(const double *x, const double *f, size_t A, size_t B)
static double psi_fromm(const double *x, const double *y, size_t j, double velocity, double dx, double dt)
static size_t upstream(double v, size_t j)
static double psi_uno3(const double *x, const double *y, size_t j, double velocity, double dx, double dt)
static size_t central(double v, size_t j)
static double psi_uno2(const double *x, const double *y, size_t j, double velocity, double dx, double dt)
static double psi_upwind1(const double *x, const double *y, size_t j, double velocity, double dx, double dt)
static double sign(double x)
static void step(double dt, const array::Staggered1 &velocity, const array::Scalar1 &x_old, array::Scalar &x)
void make_nonnegative_preserving(double dt, const array::Scalar1 &x, const array::Staggered1 &flux, array::Staggered &result)