tvonMisesCalving.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tvonMisesCalving.cc (7520B)
---
1 /* Copyright (C) 2016, 2017, 2018, 2019, 2020 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/stressbalance/StressBalance.hh"
26 #include "pism/rheology/FlowLawFactory.hh"
27 #include "pism/rheology/FlowLaw.hh"
28 #include "pism/geometry/Geometry.hh"
29
30 namespace pism {
31 namespace calving {
32
33 vonMisesCalving::vonMisesCalving(IceGrid::ConstPtr grid,
34 std::shared_ptr<const rheology::FlowLaw> flow_law)
35 : StressCalving(grid, 2),
36 m_flow_law(flow_law) {
37
38 if (m_config->get_flag("calving.vonmises_calving.use_custom_flow_law")) {
39 EnthalpyConverter::Ptr EC = grid->ctx()->enthalpy_converter();
40 rheology::FlowLawFactory factory("calving.vonmises_calving.", m_config, EC);
41 m_flow_law = factory.create();
42 }
43
44 m_calving_rate.metadata().set_name("vonmises_calving_rate");
45 m_calving_rate.set_attrs("diagnostic",
46 "horizontal calving rate due to von Mises calving",
47 "m s-1", "m year-1", "", 0);
48
49 m_calving_threshold.create(m_grid, "vonmises_calving_threshold", WITHOUT_GHOSTS);
50
51 m_calving_threshold.set_attrs("diagnostic",
52 "threshold used by the 'von Mises' calving method",
53 "Pa", "Pa",
54 "", 0); // no standard name
55 m_calving_threshold.set_time_independent(true);
56
57 }
58
59 vonMisesCalving::~vonMisesCalving() {
60 // empty
61 }
62
63 void vonMisesCalving::init() {
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, CRITICAL);
79 } else {
80 m_log->message(2,
81 " von Mises calving threshold: %3.3f Pa.\n", sigma_max);
82 }
83
84 if (fabs(m_grid->dx() - m_grid->dy()) / std::min(m_grid->dx(), m_grid->dy()) > 1e-2) {
85 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
86 "-calving vonmises_calving using a non-square grid cell is not implemented (yet);\n"
87 "dx = %f, dy = %f, relative difference = %f",
88 m_grid->dx(), m_grid->dy(),
89 fabs(m_grid->dx() - m_grid->dy()) / std::max(m_grid->dx(), m_grid->dy()));
90 }
91
92 m_strain_rates.set(0.0);
93 }
94
95 //! \brief Uses principal strain rates to apply the "von Mises" calving method
96 /*!
97 See equation (4) in [@ref Morlighem2016].
98 */
99 void vonMisesCalving::update(const IceModelVec2CellType &cell_type,
100 const IceModelVec2S &ice_thickness,
101 const IceModelVec2V &ice_velocity,
102 const IceModelVec3 &ice_enthalpy) {
103
104 using std::max;
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
112 stressbalance::compute_2D_principal_strain_rates(ice_velocity, m_cell_type,
113 m_strain_rates);
114 m_strain_rates.update_ghosts();
115
116 IceModelVec::AccessList list{&ice_enthalpy, &ice_thickness, &m_cell_type, &ice_velocity,
117 &m_strain_rates, &m_calving_rate, &m_calving_threshold};
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 * (PetscSqr(max(0.0, eigen1)) +
179 PetscSqr(max(0.0, eigen2))));
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
193 const IceModelVec2S& vonMisesCalving::threshold() const {
194 return m_calving_threshold;
195 }
196
197 DiagnosticList vonMisesCalving::diagnostics_impl() const {
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