vonMisesCalving.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