PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
vonMisesCalving.cc
Go to the documentation of this file.
1 /* Copyright (C) 2016, 2017, 2018, 2019, 2020, 2021 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 "vonMisesCalving.hh"
21 
22 #include "pism/util/IceGrid.hh"
23 #include "pism/util/error_handling.hh"
24 #include "pism/util/IceModelVec2CellType.hh"
25 #include "pism/util/IceModelVec2V.hh"
26 #include "pism/stressbalance/StressBalance.hh"
27 #include "pism/rheology/FlowLawFactory.hh"
28 #include "pism/rheology/FlowLaw.hh"
29 #include "pism/geometry/Geometry.hh"
30 #include "pism/util/Context.hh"
31 
32 namespace pism {
33 namespace calving {
34 
36  std::shared_ptr<const rheology::FlowLaw> flow_law)
37  : StressCalving(grid, 2),
38  m_calving_threshold(m_grid, "vonmises_calving_threshold", WITHOUT_GHOSTS),
39  m_flow_law(flow_law)
40 {
41 
42  if (m_config->get_flag("calving.vonmises_calving.use_custom_flow_law")) {
43  EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
44  rheology::FlowLawFactory factory("calving.vonmises_calving.", m_config, EC);
45  m_flow_law = factory.create();
46  }
47 
48  m_calving_rate.metadata().set_name("vonmises_calving_rate");
49  m_calving_rate.set_attrs("diagnostic",
50  "horizontal calving rate due to von Mises calving",
51  "m s-1", "m year-1", "", 0);
52 
53  m_calving_threshold.set_attrs("diagnostic",
54  "threshold used by the 'von Mises' calving method",
55  "Pa", "Pa",
56  "", 0); // no standard name
58 
59 }
60 
62 
63  m_log->message(2,
64  "* Initializing the 'von Mises calving' mechanism...\n");
65 
66  std::string threshold_file = m_config->get_string("calving.vonmises_calving.threshold_file");
67  double sigma_max = m_config->get_number("calving.vonmises_calving.sigma_max");
68 
69  m_calving_threshold.set(sigma_max);
70 
71  if (not threshold_file.empty()) {
72  m_log->message(2,
73  " Reading von Mises calving threshold from file '%s'...\n",
74  threshold_file.c_str());
75 
76  m_calving_threshold.regrid(threshold_file, CRITICAL);
77  } else {
78  m_log->message(2,
79  " von Mises calving threshold: %3.3f Pa.\n", sigma_max);
80  }
81 
82  if (fabs(m_grid->dx() - m_grid->dy()) / std::min(m_grid->dx(), m_grid->dy()) > 1e-2) {
84  "-calving vonmises_calving using a non-square grid cell is not implemented (yet);\n"
85  "dx = %f, dy = %f, relative difference = %f",
86  m_grid->dx(), m_grid->dy(),
87  fabs(m_grid->dx() - m_grid->dy()) / std::max(m_grid->dx(), m_grid->dy()));
88  }
89 
90  m_strain_rates.set(0.0);
91 }
92 
93 //! \brief Uses principal strain rates to apply the "von Mises" calving method
94 /*!
95  See equation (4) in [@ref Morlighem2016].
96 */
98  const IceModelVec2S &ice_thickness,
99  const IceModelVec2V &ice_velocity,
100  const IceModelVec3 &ice_enthalpy) {
101 
102  using std::max;
103  using std::sqrt;
104  using std::pow;
105 
106  // Distance (grid cells) from calving front where strain rate is evaluated
107  int offset = m_stencil_width;
108 
109  // make a copy with a wider stencil
110  m_cell_type.copy_from(cell_type);
111 
115 
116  IceModelVec::AccessList list{&ice_enthalpy, &ice_thickness, &m_cell_type, &ice_velocity,
118 
119  const double *z = &m_grid->z()[0];
120 
121  double glen_exponent = m_flow_law->exponent();
122 
123  for (Points pt(*m_grid); pt; pt.next()) {
124  const int i = pt.i(), j = pt.j();
125 
126  // Find partially filled or empty grid boxes on the icefree ocean, which
127  // have floating ice neighbors after the mass continuity step
128  if (m_cell_type.ice_free_ocean(i, j) and m_cell_type.next_to_ice(i, j)) {
129 
130  double
131  velocity_magnitude = 0.0,
132  hardness = 0.0;
133  // Average of strain-rate eigenvalues in adjacent floating grid cells.
134  double
135  eigen1 = 0.0,
136  eigen2 = 0.0;
137  {
138  int N = 0;
139  for (int p = -1; p < 2; p += 2) {
140  const int I = i + p * offset;
141  if (m_cell_type.icy(I, j)) {
142  velocity_magnitude += ice_velocity(I, j).magnitude();
143  {
144  double H = ice_thickness(I, j);
145  unsigned int k = m_grid->kBelowHeight(H);
146  hardness += averaged_hardness(*m_flow_law, H, k, &z[0], ice_enthalpy.get_column(I, j));
147  }
148  eigen1 += m_strain_rates(I, j, 0);
149  eigen2 += m_strain_rates(I, j, 1);
150  N += 1;
151  }
152  }
153 
154  for (int q = -1; q < 2; q += 2) {
155  const int J = j + q * offset;
156  if (m_cell_type.icy(i, J)) {
157  velocity_magnitude += ice_velocity(i, J).magnitude();
158  {
159  double H = ice_thickness(i, J);
160  unsigned int k = m_grid->kBelowHeight(H);
161  hardness += averaged_hardness(*m_flow_law, H, k, &z[0], ice_enthalpy.get_column(i, J));
162  }
163  eigen1 += m_strain_rates(i, J, 0);
164  eigen2 += m_strain_rates(i, J, 1);
165  N += 1;
166  }
167  }
168 
169  if (N > 0) {
170  eigen1 /= N;
171  eigen2 /= N;
172  hardness /= N;
173  velocity_magnitude /= N;
174  }
175  }
176 
177  // [\ref Morlighem2016] equation 6
178  const double effective_tensile_strain_rate = sqrt(0.5 * (pow(max(0.0, eigen1), 2) +
179  pow(max(0.0, eigen2), 2)));
180  // [\ref Morlighem2016] equation 7
181  const double sigma_tilde = sqrt(3.0) * hardness * pow(effective_tensile_strain_rate,
182  1.0 / glen_exponent);
183 
184  // Calving law [\ref Morlighem2016] equation 4
185  m_calving_rate(i, j) = velocity_magnitude * sigma_tilde / m_calving_threshold(i, j);
186 
187  } else { // end of "if (ice_free_ocean and next_to_floating)"
188  m_calving_rate(i, j) = 0.0;
189  }
190  } // end of loop over grid points
191 }
192 
194  return m_calving_threshold;
195 }
196 
198  return {{"vonmises_calving_rate", Diagnostic::wrap(m_calving_rate)},
199  {"vonmises_calving_threshold", Diagnostic::wrap(m_calving_threshold)}};
200 }
201 
202 } // end of namespace calving
203 } // end of namespace pism
Makes sure that we call begin_access() and end_access() for all accessed IceModelVecs.
Definition: iceModelVec.hh:59
IceGrid::ConstPtr grid() const
Definition: Component.cc:105
const Config::ConstPtr m_config
configuration database used by this component
Definition: Component.hh:138
const Logger::ConstPtr m_log
logger (for easy access)
Definition: Component.hh:142
const IceGrid::ConstPtr m_grid
grid used by this component
Definition: Component.hh:136
static Ptr wrap(const T &input)
Definition: Diagnostic.hh:155
std::shared_ptr< EnthalpyConverter > Ptr
std::shared_ptr< const IceGrid > ConstPtr
Definition: IceGrid.hh:233
bool next_to_ice(int i, int j) const
Ice-free margin (at least one of four neighbors has ice).
bool icy(int i, int j) const
bool ice_free_ocean(int i, int j) const
"Cell type" mask. Adds convenience methods to IceModelVec2Int.
void copy_from(const IceModelVec2S &source)
double * get_column(int i, int j)
A virtual class collecting methods common to ice and bedrock 3D fields.
Definition: iceModelVec.hh:404
void update_ghosts()
Updates ghost points.
Definition: iceModelVec.cc:669
SpatialVariableMetadata & metadata(unsigned int N=0)
Returns a reference to the SpatialVariableMetadata object containing metadata for the compoment N.
Definition: iceModelVec.cc:533
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.
Definition: iceModelVec.cc:399
void regrid(const std::string &filename, RegriddingFlag flag, double default_value=0.0)
Definition: iceModelVec.cc:838
void set(double c)
Result: v[j] <- c for all j.
Definition: iceModelVec.cc:683
void set_time_independent(bool flag)
Set the time independent flag for all variables corresponding to this IceModelVec instance.
Definition: iceModelVec.cc:175
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
void set_name(const std::string &name)
IceModelVec2CellType m_cell_type
An abstract class containing fields used by all stress-based calving methods.
std::shared_ptr< const rheology::FlowLaw > m_flow_law
void update(const IceModelVec2CellType &cell_type, const IceModelVec2S &ice_thickness, const IceModelVec2V &ice_velocity, const IceModelVec3 &ice_enthalpy)
Uses principal strain rates to apply the "von Mises" calving method.
DiagnosticList diagnostics_impl() const
vonMisesCalving(IceGrid::ConstPtr grid, std::shared_ptr< const rheology::FlowLaw > flow_law)
const IceModelVec2S & threshold() const
std::shared_ptr< FlowLaw > create()
#define PISM_ERROR_LOCATION
double averaged_hardness(const FlowLaw &ice, double thickness, int kbelowH, const double *zlevels, const double *enthalpy)
Computes vertical average of B(E, p) ice hardness, namely .
Definition: FlowLaw.cc:214
void compute_2D_principal_strain_rates(const IceModelVec2V &V, const IceModelVec2CellType &mask, IceModelVec3 &result)
Compute eigenvalues of the horizontal, vertically-integrated strain rate tensor.
double max(const IceModelVec2S &input)
Finds maximum over all the values in an IceModelVec2S object. Ignores ghosts.
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:117
static const double k
Definition: exactTestP.cc:45
@ CRITICAL
Definition: IO_Flags.hh:70
double min(const IceModelVec2S &input)
Finds minimum over all the values in an IceModelVec2S object. Ignores ghosts.
@ WITHOUT_GHOSTS
Definition: iceModelVec.hh:49
const int J[]
Definition: ssafd_code.cc:34
const int I[]
Definition: ssafd_code.cc:24