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