Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.1-cd005eec8 committed by Constantine Khrulev on 2025-03-07
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BlatterMod.cc
Go to the documentation of this file.
1/* Copyright (C) 2020, 2021, 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 <algorithm> // std::min, std::max
21
22#include "pism/stressbalance/blatter/BlatterMod.hh"
23
24#include "pism/rheology/FlowLawFactory.hh"
25
26#include "pism/geometry/Geometry.hh"
27#include "pism/stressbalance/StressBalance.hh" // Inputs
28#include "pism/util/pism_utilities.hh" // GlobalMax
29
30namespace pism {
31namespace stressbalance {
32
33BlatterMod::BlatterMod(std::shared_ptr<Blatter> solver)
34 : SSB_Modifier(solver->grid()),
35 m_solver(solver) {
36
37 rheology::FlowLawFactory ice_factory("stress_balance.blatter.", m_config, m_EC);
38 ice_factory.remove(ICE_GOLDSBY_KOHLSTEDT);
39 m_flow_law = ice_factory.create();
40}
41
43 // empty
44}
45
46/*!
47 * Post-process ice velocity computed by the Blatter solver.
48 *
49 * - transfers velocity from the sigma grid onto PISM's vertical grid
50 *
51 * - estimates the maximum diffusivity used to compute the time step restriction
52 */
53void BlatterMod::update(const array::Vector &sliding_velocity,
54 const Inputs &inputs,
55 bool full_update) {
56 (void) sliding_velocity;
57 (void) full_update;
58
59 // transfer velocity onto the PISM vertical grid, setting m_u and m_v
61
62 // estimate max diffusivity to use in adaptive time stepping
64 inputs.geometry->ice_thickness,
66
68}
69
70/*!
71 * Copy ice velocity from the sigma vertical grid onto PISM's vertical grid.
72 *
73 * Uses constant extrapolation above the ice surface.
74 */
75void BlatterMod::transfer(const array::Scalar &ice_thickness) {
76
77 auto u_sigma = m_solver->velocity_u_sigma();
78 auto v_sigma = m_solver->velocity_v_sigma();
79
80 const auto &zlevels = m_u.levels();
81 int Mz = zlevels.size();
82
83 array::AccessScope list{&m_u, &m_v, u_sigma.get(), v_sigma.get(), &ice_thickness};
84
85 for (auto p = m_grid->points(); p; p.next()) {
86 const int i = p.i(), j = p.j();
87
88 auto *u = m_u.get_column(i, j);
89 auto *v = m_v.get_column(i, j);
90
91 double H = ice_thickness(i, j);
92
93 if (H > 0.0) {
94 for (int k = 0; k < Mz; ++k) {
95 double sigma = std::min(zlevels[k] / H, 1.0);
96
97 u[k] = u_sigma->interpolate(i, j, sigma);
98 v[k] = v_sigma->interpolate(i, j, sigma);
99 }
100 } else {
101 m_u.set_column(i, j, 0.0);
102 m_v.set_column(i, j, 0.0);
103 }
104 }
105
108}
109
110/*!
111 * Estimate max SIA-type diffusivity assuming that `Q = -D \nabla h`.
112 */
114 const array::Scalar &ice_thickness,
115 const array::Scalar1 &surface) {
116 const double eps = 1e-3;
117 double
118 dx = m_grid->dx(),
119 dy = m_grid->dy();
120
121 array::AccessScope list{&velocity, &ice_thickness, &surface};
122
123 m_D_max = 0.0;
124 for (auto p = m_grid->points(); p; p.next()) {
125 const int i = p.i(), j = p.j();
126
127 auto h = surface.star(i, j);
128 auto H = ice_thickness(i, j);
129
130 if (H > 0.0) {
131 Vector2d grad_h = {(h.e - h.w) / (2.0 * dx),
132 (h.n - h.s) / (2.0 * dy)};
133
134 double D = H * velocity(i, j).magnitude() / (grad_h.magnitude() + eps);
135
136 m_D_max = std::max(D, m_D_max);
137 }
138 }
139
141}
142
143} // end of namespace stressbalance
144} // end of namespace pism
#define ICE_GOLDSBY_KOHLSTEDT
const Config::ConstPtr m_config
configuration database used by this component
Definition Component.hh:158
const std::shared_ptr< const Grid > m_grid
grid used by this component
Definition Component.hh:156
array::Scalar2 ice_surface_elevation
Definition Geometry.hh:57
array::Scalar2 ice_thickness
Definition Geometry.hh:51
double magnitude() const
Magnitude.
Definition Vector2d.hh:40
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:64
stencils::Star< T > star(int i, int j) const
Definition Array2D.hh:79
void set_column(int i, int j, double c)
Set all values of scalar quantity to given a single value in a particular column.
Definition Array3D.cc:50
double * get_column(int i, int j)
Definition Array3D.cc:121
const std::vector< double > & levels() const
Definition Array.cc:139
void set(double c)
Result: v[j] <- c for all j.
Definition Array.cc:629
void update_ghosts()
Updates ghost points.
Definition Array.cc:615
std::shared_ptr< FlowLaw > create()
void remove(const std::string &name)
std::shared_ptr< Blatter > m_solver
Definition BlatterMod.hh:45
void update(const array::Vector &sliding_velocity, const Inputs &inputs, bool full_update)
Definition BlatterMod.cc:53
void transfer(const array::Scalar &ice_thickness)
Definition BlatterMod.cc:75
BlatterMod(std::shared_ptr< Blatter > solver)
Definition BlatterMod.cc:33
void compute_max_diffusivity(const array::Vector &velocity, const array::Scalar &ice_thickness, const array::Scalar1 &surface)
EnthalpyConverter::Ptr m_EC
std::shared_ptr< rheology::FlowLaw > m_flow_law
Shallow stress balance modifier (such as the non-sliding SIA).
void GlobalMax(MPI_Comm comm, double *local, double *result, int count)
static const double k
Definition exactTestP.cc:42