tEigenCalving.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
---
tEigenCalving.cc (4923B)
---
1 /* Copyright (C) 2013, 2014, 2015, 2016, 2017, 2018, 2019 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 "EigenCalving.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/geometry/Geometry.hh"
27
28 namespace pism {
29 namespace calving {
30
31 EigenCalving::EigenCalving(IceGrid::ConstPtr grid)
32 : StressCalving(grid, 2) {
33
34 m_K = m_config->get_number("calving.eigen_calving.K");
35
36 m_calving_rate.metadata().set_name("eigen_calving_rate");
37 m_calving_rate.set_attrs("diagnostic",
38 "horizontal calving rate due to eigen-calving",
39 "m s-1", "m year-1", "", 0);
40 }
41
42 EigenCalving::~EigenCalving() {
43 // empty
44 }
45
46 void EigenCalving::init() {
47
48 m_log->message(2, "* Initializing the 'eigen-calving' mechanism...\n");
49
50 if (fabs(m_grid->dx() - m_grid->dy()) / std::min(m_grid->dx(), m_grid->dy()) > 1e-2) {
51 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "-calving eigen_calving using a non-square grid cell is not implemented (yet);\n"
52 "dx = %f, dy = %f, relative difference = %f",
53 m_grid->dx(), m_grid->dy(),
54 fabs(m_grid->dx() - m_grid->dy()) / std::max(m_grid->dx(), m_grid->dy()));
55 }
56
57 m_strain_rates.set(0.0);
58 }
59
60 //! \brief Uses principal strain rates to apply "eigencalving" with constant K.
61 /*!
62 See equation (26) in [\ref Winkelmannetal2011].
63 */
64 void EigenCalving::update(const IceModelVec2CellType &cell_type,
65 const IceModelVec2V &ice_velocity) {
66
67 // make a copy with a wider stencil
68 m_cell_type.copy_from(cell_type);
69
70 // Distance (grid cells) from calving front where strain rate is evaluated
71 int offset = m_stencil_width;
72
73 // eigenCalvOffset allows adjusting the transition from compressive to extensive flow
74 // regime
75 const double eigenCalvOffset = 0.0;
76
77 stressbalance::compute_2D_principal_strain_rates(ice_velocity, m_cell_type,
78 m_strain_rates);
79 m_strain_rates.update_ghosts();
80
81 IceModelVec::AccessList list{&m_cell_type, &m_calving_rate, &m_strain_rates};
82
83 // Compute the horizontal calving rate
84 for (Points pt(*m_grid); pt; pt.next()) {
85 const int i = pt.i(), j = pt.j();
86
87 // Find partially filled or empty grid boxes on the icefree ocean, which
88 // have floating ice neighbors after the mass continuity step
89 if (m_cell_type.ice_free_ocean(i, j) and m_cell_type.next_to_floating_ice(i, j)) {
90
91 // Average of strain-rate eigenvalues in adjacent floating grid cells to be used for
92 // eigen-calving:
93 double
94 eigen1 = 0.0,
95 eigen2 = 0.0;
96 {
97 int N = 0;
98 for (int p = -1; p < 2; p += 2) {
99 const int I = i + p * offset;
100 if (m_cell_type.floating_ice(I, j) and not m_cell_type.ice_margin(I, j)) {
101 eigen1 += m_strain_rates(I, j, 0);
102 eigen2 += m_strain_rates(I, j, 1);
103 N += 1;
104 }
105 }
106
107 for (int q = -1; q < 2; q += 2) {
108 const int J = j + q * offset;
109 if (m_cell_type.floating_ice(i, J) and not m_cell_type.ice_margin(i, J)) {
110 eigen1 += m_strain_rates(i, J, 0);
111 eigen2 += m_strain_rates(i, J, 1);
112 N += 1;
113 }
114 }
115
116 if (N > 0) {
117 eigen1 /= N;
118 eigen2 /= N;
119 }
120 }
121
122 // Calving law
123 //
124 // eigen1 * eigen2 has units [s^-2] and calving_rate_horizontal
125 // [m*s^1] hence, eigen_calving_K has units [m*s]
126 if (eigen2 > eigenCalvOffset and eigen1 > 0.0) {
127 // spreading in all directions
128 m_calving_rate(i, j) = m_K * eigen1 * (eigen2 - eigenCalvOffset);
129 } else {
130 m_calving_rate(i, j) = 0.0;
131 }
132
133 } else { // end of "if (ice_free_ocean and next_to_floating)"
134 m_calving_rate(i, j) = 0.0;
135 }
136 } // end of the loop over grid points
137 }
138
139 DiagnosticList EigenCalving::diagnostics_impl() const {
140 return {{"eigen_calving_rate", Diagnostic::wrap(m_calving_rate)}};
141 }
142
143 } // end of namespace calving
144 } // end of namespace pism