22 #include "pism/geometry/Geometry.hh"
23 #include "pism/util/pism_utilities.hh"
24 #include "pism/geometry/Geometry.hh"
29 namespace diagnostics {
46 const int i = p.i(), j = p.j();
48 auto P = psi.
star(i, j);
56 if (domain_mask(i, j) > 0.5 and v_e <= 0.0 and v_w >= 0.0 and v_n <= 0.0 and v_s >= 0.0) {
76 {&ice_thickness, &bed_elevation, &cell_type, &sea_level_elevation,
77 &water_flux, &result};
80 grid_spacing = 0.5 * (grid->dx() + grid->dy()),
84 const int i = p.i(), j = p.j();
86 if (cell_type.
icy(i, j)) {
96 double water_depth =
std::max(sea_level_elevation(i, j) - bed_elevation(i, j), 0.0),
97 submerged_front_area = water_depth * grid_spacing;
99 if (submerged_front_area > 0.0) {
100 auto Q_sg = water_flux(i, j) * grid_spacing;
101 auto q_sg = Q_sg /
std::max(submerged_front_area, eps);
117 m_potential(grid,
"hydraulic_potential",
WITH_GHOSTS, 1),
120 m_W(grid,
"remaining_water_thickness",
WITH_GHOSTS, 1),
123 m_domain_mask(grid,
"domain_mask",
WITH_GHOSTS, 1),
126 m_adjustment(grid,
"hydraulic_potential_adjustment",
WITHOUT_GHOSTS),
129 m_potential.
set_attrs(
"diagnostic",
"estimate of the steady state hydraulic potential in the steady hydrology model",
136 "scaled water thickness in the steady state hydrology model"
137 " (has no physical meaning)",
141 "m s-1",
"m s-1",
"", 0);
145 m_Q.
set_attrs(
"diagnostic",
"steady state water flux",
"m2 s-1",
"m2 s-1",
"", 0);
147 m_q_sg.
set_attrs(
"diagnostic",
"x-component of the effective water velocity in the steady-state hydrology model",
148 "m s-1",
"m day-1",
"", 0);
149 m_q_sg.
set_attrs(
"diagnostic",
"y-component of the effective water velocity in the steady-state hydrology model",
150 "m s-1",
"m day-1",
"", 1);
152 m_sinks.
set_attrs(
"diagnostic",
"map of sinks in the domain (for debugging)",
"",
"",
"", 0);
155 "potential adjustment needed to fill sinks"
156 " when computing an estimate of the steady-state hydraulic potential",
164 m_tau =
m_config->get_number(
"hydrology.steady.input_rate_scaling");
177 bool recompute_potential) {
181 cell_area =
m_grid->cell_area(),
184 dt = 0.5 / (u_max /
m_dx + v_max /
m_dy),
185 volume_ratio =
m_config->get_number(
"hydrology.steady.volume_ratio");
187 const int n_iterations =
m_config->get_number(
"hydrology.steady.n_iterations");
189 if (recompute_potential) {
211 double volume_0 = 0.0;
216 const int i = p.i(), j = p.j();
219 m_W(i, j) =
m_tau * water_input_rate(i, j);
224 volume_0 +=
m_W(i, j);
236 if (volume_0 == 0.0) {
243 int step_counter = 0;
247 for (step_counter = 0; step_counter < n_iterations; ++step_counter) {
251 const int i = p.i(), j = p.j();
257 q_n = v.n * (v.n >= 0.0 ? w.ij : w.n),
258 q_e = v.e * (v.e >= 0.0 ? w.ij : w.e),
259 q_s = v.s * (v.s >= 0.0 ? w.s : w.ij),
260 q_w = v.w * (v.w >= 0.0 ? w.w : w.ij),
261 divQ = (q_e - q_w) /
m_dx + (q_n - q_s) /
m_dy;
265 m_tmp(i, j) = w.ij + dt * (- divQ);
270 if (
m_tmp(i, j) < -eps) {
276 m_Qsum(i, j, 0) += dt * q_e;
277 m_Qsum(i, j, 1) += dt * q_n;
280 volume +=
m_tmp(i, j);
286 if (volume / volume_0 <= volume_ratio) {
289 m_log->message(3,
"%04d V = %f\n", step_counter, volume / volume_0);
292 m_log->message(3,
"Emptying problem: stopped after %d iterations. V = %f\n",
293 step_counter, volume / volume_0);
295 double epsilon = volume / volume_0;
316 g =
m_config->get_number(
"constants.standard_gravity"),
317 rho_i =
m_config->get_number(
"constants.ice.density"),
318 rho_w =
m_config->get_number(
"constants.fresh_water.density");
323 const int i = p.i(), j = p.j();
325 result(i, j) = rho_i *
g * H(i, j) + rho_w *
g * b(i, j);
337 double delta =
m_config->get_number(
"hydrology.steady.potential_delta");
339 int n_iterations =
m_config->get_number(
"hydrology.steady.potential_n_iterations");
340 int step_counter = 0;
342 int n_sinks_remaining = 0;
347 for (step_counter = 0; step_counter < n_iterations; ++step_counter) {
349 n_sinks_remaining = 0;
351 const int i = p.i(), j = p.j();
353 if (domain_mask(i, j) > 0.5) {
354 auto P = result.
star(i, j);
357 v_n = - (P.n - P.ij),
358 v_e = - (P.e - P.ij),
359 v_s = - (P.ij - P.s),
360 v_w = - (P.ij - P.w);
362 if (v_e <= 0.0 and v_w >= 0.0 and v_n <= 0.0 and v_s >= 0.0) {
365 psi_new(i, j) = 0.25 * (P.n + P.e + P.s + P.w) + delta;
367 psi_new(i, j) = result(i, j);
370 psi_new(i, j) = result(i, j);
379 if (step_counter == 0) {
380 n_sinks = n_sinks_remaining;
383 if (n_sinks_remaining == 0) {
384 m_log->message(3,
"Emptying problem: filled %d sinks after %d iterations.\n",
385 n_sinks - n_sinks_remaining, step_counter);
390 if (n_sinks_remaining > 0) {
391 m_log->message(2,
"WARNING: %d sinks left.\n", n_sinks_remaining);
396 static double K(
double psi_x,
double psi_y,
double speed,
double epsilon) {
410 const int i = p.i(), j = p.j();
412 for (
int o = 0; o < 2; ++o) {
415 psi_x = (psi(i + 1, j) - psi(i, j)) /
m_dx;
416 psi_y = (psi(i + 1, j + 1) + psi(i, j + 1) - psi(i + 1, j - 1) - psi(i, j - 1)) / (4.0 *
m_dy);
419 psi_x = (psi(i + 1, j + 1) + psi(i + 1, j) - psi(i - 1, j + 1) - psi(i - 1, j)) / (4.0 *
m_dx);
420 psi_y = (psi(i, j + 1) - psi(i, j)) /
m_dy;
424 auto M = domain_mask.
star(i, j);
426 if (M.ij == 0 and M.e == 0) {
427 result(i, j, 0) = 0.0;
430 if (M.ij == 0 and M.n == 0) {
431 result(i, j, 1) = 0.0;
448 list.
add(*no_model_mask);
452 const int i = p.i(), j = p.j();
460 if (no_model_mask and no_model_mask->
as_int(i, j) == 1) {
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
const Config::ConstPtr m_config
configuration database used by this component
const Logger::ConstPtr m_log
logger (for easy access)
const IceGrid::ConstPtr m_grid
grid used by this component
A class defining a common interface for most PISM sub-models.
static Ptr wrap(const T &input)
IceModelVec2CellType cell_type
IceModelVec2S bed_elevation
IceModelVec2S sea_level_elevation
IceModelVec2S ice_thickness
std::shared_ptr< const IceGrid > ConstPtr
bool icy(int i, int j) const
bool ice_free_ocean(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
int as_int(int i, int j) const
stencils::Star< int > star(int i, int j) const
A simple class "hiding" the fact that the mask is stored as floating-point scalars (instead of intege...
void add(double alpha, const IceModelVec2S &x)
void copy_from(const IceModelVec2S &source)
stencils::Star< double > star(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.
A class for storing and accessing internal staggered-grid 2D fields. Uses dof=2 storage....
void update_ghosts()
Updates ghost points.
void set_attrs(const std::string &pism_intent, const std::string &long_name, const std::string &units, const std::string &glaciological_units, const std::string &standard_name, unsigned int component)
Sets NetCDF attributes of an IceModelVec object.
void set(double c)
Result: v[j] <- c for all j.
void scale(double alpha)
Result: v <- v * alpha. Calls VecScale.
IceGrid::ConstPtr grid() const
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
const IceModelVec2V & flux() const
virtual void compute_raw_potential(const IceModelVec2S &ice_thickness, const IceModelVec2S &ice_bottom_surface, IceModelVec2S &result) const
void update(const Geometry &geometry, const IceModelVec2Int *no_model_mask, const IceModelVec2S &water_input_rate, bool recompute_potential=true)
IceModelVec2S m_adjustment
const IceModelVec2Int & sinks() const
IceModelVec2S m_bottom_surface
DiagnosticList diagnostics() const
IceModelVec2Int m_domain_mask
EmptyingProblem(IceGrid::ConstPtr g)
const IceModelVec2V & effective_water_velocity() const
const IceModelVec2S & potential() const
IceModelVec2S m_potential
void compute_velocity(const IceModelVec2S &hydraulic_potential, const IceModelVec2Int &mask, IceModelVec2Stag &result) const
void compute_potential(const IceModelVec2S &ice_thickness, const IceModelVec2S &ice_bottom_surface, const IceModelVec2Int &domain_mask, IceModelVec2S &result)
const IceModelVec2S & remaining_water_thickness() const
void compute_mask(const IceModelVec2CellType &cell_type, const IceModelVec2Int *no_model_mask, IceModelVec2Int &result) const
const IceModelVec2S & adjustment() const
#define PISM_ERROR_LOCATION
static void effective_water_velocity(const Geometry &geometry, const IceModelVec2V &water_flux, IceModelVec2V &result)
static void compute_sinks(const IceModelVec2Int &domain_mask, const IceModelVec2S &psi, IceModelVec2S &result)
static double K(double psi_x, double psi_y, double speed, double epsilon)
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
std::map< std::string, Diagnostic::Ptr > DiagnosticList
void ice_bottom_surface(const Geometry &geometry, IceModelVec2S &result)
void GlobalSum(MPI_Comm comm, double *local, double *result, int count)
void staggered_to_regular(const IceModelVec2CellType &cell_type, const IceModelVec2Stag &input, bool include_floating_ice, IceModelVec2S &result)