tSSA.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
---
tSSA.cc (13982B)
---
1 // Copyright (C) 2004--2019 Constantine Khroulev, Ed Bueler, Jed Brown, Torsten Albrecht
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 #include "SSA.hh"
20 #include "pism/basalstrength/basal_resistance.hh"
21 #include "pism/util/EnthalpyConverter.hh"
22 #include "pism/rheology/FlowLawFactory.hh"
23 #include "pism/util/Mask.hh"
24 #include "pism/util/Vars.hh"
25 #include "pism/util/error_handling.hh"
26 #include "pism/util/io/File.hh"
27 #include "pism/util/pism_options.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/util/IceModelVec2CellType.hh"
30 #include "pism/stressbalance/StressBalance.hh"
31 #include "pism/geometry/Geometry.hh"
32
33 #include "SSA_diagnostics.hh"
34
35 namespace pism {
36 namespace stressbalance {
37
38 SSAStrengthExtension::SSAStrengthExtension(const Config &config) {
39 m_min_thickness = config.get_number("stress_balance.ssa.strength_extension.min_thickness");
40 m_constant_nu = config.get_number("stress_balance.ssa.strength_extension.constant_nu");
41 }
42
43 //! Set strength = (viscosity times thickness).
44 /*! Determines nu by input strength and current min_thickness. */
45 void SSAStrengthExtension::set_notional_strength(double my_nuH) {
46 if (my_nuH <= 0.0) {
47 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "nuH must be positive, got %f", my_nuH);
48 }
49 m_constant_nu = my_nuH / m_min_thickness;
50 }
51
52 //! Set minimum thickness to trigger use of extension.
53 /*! Preserves strength (nuH) by also updating using current nu. */
54 void SSAStrengthExtension::set_min_thickness(double my_min_thickness) {
55 if (my_min_thickness <= 0.0) {
56 throw RuntimeError::formatted(PISM_ERROR_LOCATION, "min_thickness must be positive, got %f",
57 my_min_thickness);
58 }
59 double nuH = m_constant_nu * m_min_thickness;
60 m_min_thickness = my_min_thickness;
61 m_constant_nu = nuH / m_min_thickness;
62 }
63
64 //! Returns strength = (viscosity times thickness).
65 double SSAStrengthExtension::get_notional_strength() const {
66 return m_constant_nu * m_min_thickness;
67 }
68
69 //! Returns minimum thickness to trigger use of extension.
70 double SSAStrengthExtension::get_min_thickness() const {
71 return m_min_thickness;
72 }
73
74
75 SSA::SSA(IceGrid::ConstPtr g)
76 : ShallowStressBalance(g)
77 {
78 strength_extension = new SSAStrengthExtension(*m_config);
79
80 const unsigned int WIDE_STENCIL = m_config->get_number("grid.max_stencil_width");
81
82 // grounded_dragging_floating integer mask
83 m_mask.create(m_grid, "ssa_mask", WITH_GHOSTS, WIDE_STENCIL);
84 m_mask.set_attrs("diagnostic", "ice-type (ice-free/grounded/floating/ocean) integer mask",
85 "", "", "", 0);
86 std::vector<double> mask_values(4);
87 mask_values[0] = MASK_ICE_FREE_BEDROCK;
88 mask_values[1] = MASK_GROUNDED;
89 mask_values[2] = MASK_FLOATING;
90 mask_values[3] = MASK_ICE_FREE_OCEAN;
91 m_mask.metadata().set_numbers("flag_values", mask_values);
92 m_mask.metadata().set_string("flag_meanings",
93 "ice_free_bedrock grounded_ice floating_ice ice_free_ocean");
94
95 m_taud.create(m_grid, "taud", WITHOUT_GHOSTS);
96 m_taud.set_attrs("diagnostic",
97 "X-component of the driving shear stress at the base of ice",
98 "Pa", "Pa", "", 0);
99 m_taud.set_attrs("diagnostic",
100 "Y-component of the driving shear stress at the base of ice",
101 "Pa", "Pa", "", 1);
102
103
104 // override velocity metadata
105 m_velocity.metadata(0).set_name("u_ssa");
106 m_velocity.metadata(0).set_string("long_name", "SSA model ice velocity in the X direction");
107
108 m_velocity.metadata(1).set_name("v_ssa");
109 m_velocity.metadata(1).set_string("long_name", "SSA model ice velocity in the Y direction");
110
111 m_velocity_global.create(m_grid, "bar", WITHOUT_GHOSTS);
112
113 m_da = m_velocity_global.dm();
114
115 {
116 rheology::FlowLawFactory ice_factory("stress_balance.ssa.", m_config, m_EC);
117 ice_factory.remove(ICE_GOLDSBY_KOHLSTEDT);
118 m_flow_law = ice_factory.create();
119 }
120 }
121
122 SSA::~SSA() {
123 if (strength_extension != NULL) {
124 delete strength_extension;
125 strength_extension = NULL;
126 }
127 }
128
129
130 //! \brief Initialize a generic regular-grid SSA solver.
131 void SSA::init_impl() {
132
133 ShallowStressBalance::init_impl();
134
135 m_log->message(2, "* Initializing the SSA stress balance...\n");
136 m_log->message(2,
137 " [using the %s flow law]\n", m_flow_law->name().c_str());
138
139 InputOptions opts = process_input_options(m_grid->com, m_config);
140
141 // Check if PISM is being initialized from an output file from a previous run
142 // and read the initial guess (unless asked not to).
143 if (opts.type == INIT_RESTART) {
144 if (m_config->get_flag("stress_balance.ssa.read_initial_guess")) {
145 File input_file(m_grid->com, opts.filename, PISM_GUESS, PISM_READONLY);
146 bool u_ssa_found = input_file.find_variable("u_ssa");
147 bool v_ssa_found = input_file.find_variable("v_ssa");
148 unsigned int start = input_file.nrecords() - 1;
149
150 if (u_ssa_found and v_ssa_found) {
151 m_log->message(3, "Reading u_ssa and v_ssa...\n");
152
153 m_velocity.read(input_file, start);
154 }
155 }
156 } else {
157 m_velocity.set(0.0); // default initial guess
158 }
159 }
160
161 //! \brief Update the SSA solution.
162 void SSA::update(const Inputs &inputs, bool full_update) {
163
164 // update the cell type mask using the ice-free thickness threshold for stress balance
165 // computations
166 {
167 const double H_threshold = m_config->get_number("stress_balance.ice_free_thickness_standard");
168 GeometryCalculator gc(*m_config);
169 gc.set_icefree_thickness(H_threshold);
170
171 gc.compute_mask(inputs.geometry->sea_level_elevation,
172 inputs.geometry->bed_elevation,
173 inputs.geometry->ice_thickness,
174 m_mask);
175 }
176
177 if (full_update) {
178 solve(inputs);
179 compute_basal_frictional_heating(m_velocity,
180 *inputs.basal_yield_stress,
181 m_mask,
182 m_basal_frictional_heating);
183 }
184 }
185
186 /*!
187 * Compute the weight used to determine if the difference between locations `i,j` and `n`
188 * (neighbor) should be used in the computation of the surface gradient in
189 * SSA::compute_driving_stress().
190 *
191 * We avoid differencing across
192 *
193 * - ice margins if stress boundary condition at ice margins (CFBC) is active
194 * - grounding lines
195 * - ice margins next to ice free locations above the surface elevation of the ice (fjord
196 * walls, nunataks, headwalls)
197 */
198 static int weight(bool margin_bc,
199 int M_ij, int M_n,
200 double h_ij, double h_n,
201 int N_ij, int N_n) {
202 using mask::grounded;
203 using mask::icy;
204 using mask::floating_ice;
205 using mask::ice_free;
206 using mask::ice_free_ocean;
207
208 // grounding lines and calving fronts
209 if ((grounded(M_ij) and floating_ice(M_n)) or
210 (floating_ice(M_ij) and grounded(M_n)) or
211 (floating_ice(M_ij) and ice_free_ocean(M_n))) {
212 return 0;
213 }
214
215 // fjord walls, nunataks, headwalls
216 if ((icy(M_ij) and ice_free(M_n) and h_n > h_ij) or
217 (ice_free(M_ij) and icy(M_n) and h_ij > h_n)) {
218 return 0;
219 }
220
221 // This condition has to match the one used to implement the calving front stress
222 // boundary condition in SSAFD::assemble_rhs().
223 if (margin_bc and
224 ((icy(M_ij) and ice_free(M_n)) or
225 (ice_free(M_ij) and icy(M_n)))) {
226 return 0;
227 }
228
229 // boundaries of the "no model" area
230 if ((N_ij == 0 and N_n == 1) or (N_ij == 1 and N_n == 0)) {
231 return 0;
232 }
233
234 return 1;
235 }
236
237 //! \brief Compute the gravitational driving stress.
238 /*!
239 Computes the gravitational driving stress at the base of the ice:
240 \f[ \tau_d = - \rho g H \nabla h \f]
241
242 If configuration parameter `sia.surface_gradient_method` = `eta` then the surface
243 gradient \f$\nabla h\f$ is computed by the gradient of the transformed variable
244 \f$\eta= H^{(2n+2)/n}\f$ (frequently, \f$\eta= H^{8/3}\f$). The idea is that
245 this quantity is more regular at ice sheet margins, and so we get a better
246 surface gradient. When the thickness at a grid point is very small (below \c
247 minThickEtaTransform in the procedure), the formula is slightly modified to
248 give a lower driving stress. The transformation is not used in floating ice.
249 */
250 void SSA::compute_driving_stress(const IceModelVec2S &ice_thickness,
251 const IceModelVec2S &surface_elevation,
252 const IceModelVec2CellType &cell_type,
253 const IceModelVec2Int *no_model_mask,
254 IceModelVec2V &result) const {
255
256 bool cfbc = m_config->get_flag("stress_balance.calving_front_stress_bc");
257 bool surface_gradient_inward = m_config->get_flag("stress_balance.ssa.compute_surface_gradient_inward");
258
259 double
260 dx = m_grid->dx(),
261 dy = m_grid->dy();
262
263 IceModelVec::AccessList list{&surface_elevation, &cell_type, &ice_thickness, &result};
264
265 if (no_model_mask) {
266 list.add(*no_model_mask);
267 }
268
269 for (Points p(*m_grid); p; p.next()) {
270 const int i = p.i(), j = p.j();
271
272 const double pressure = m_EC->pressure(ice_thickness(i, j)); // FIXME issue #15
273 if (pressure <= 0.0) {
274 result(i, j) = 0.0;
275 continue;
276 }
277
278 // Special case for verification tests.
279 if (surface_gradient_inward) {
280 double
281 h_x = surface_elevation.diff_x_p(i, j),
282 h_y = surface_elevation.diff_y_p(i, j);
283 result(i, j) = - pressure * Vector2(h_x, h_y);
284 continue;
285 }
286
287 // To compute the x-derivative we use
288 //
289 // * away from the grounding line, ice margins, and no_model mask transitions -- 2nd
290 // order centered difference
291 //
292 // * at the grounded cell near the grounding line -- 1st order
293 // one-sided difference using the grounded neighbor
294 //
295 // * at the floating cell near the grounding line -- 1st order
296 // one-sided difference using the floating neighbor
297 //
298 // All these cases can be combined by writing h_x as the weighted
299 // average of one-sided differences, with weights of 0 if a finite
300 // difference is not used and 1 if it is.
301 //
302 // The y derivative is handled the same way.
303
304 auto M = cell_type.int_star(i, j);
305 auto h = surface_elevation.star(i, j);
306 StarStencil<int> N(0);
307
308 if (no_model_mask) {
309 N = no_model_mask->int_star(i, j);
310 }
311
312 // x-derivative
313 double h_x = 0.0;
314 {
315 double
316 west = weight(cfbc, M.ij, M.w, h.ij, h.w, N.ij, N.w),
317 east = weight(cfbc, M.ij, M.e, h.ij, h.e, N.ij, N.e);
318
319 if (east + west > 0) {
320 h_x = 1.0 / ((west + east) * dx) * (west * (h.ij - h.w) + east * (h.e - h.ij));
321 } else {
322 h_x = 0.0;
323 }
324 }
325
326 // y-derivative
327 double h_y = 0.0;
328 {
329 double
330 south = weight(cfbc, M.ij, M.s, h.ij, h.s, N.ij, N.s),
331 north = weight(cfbc, M.ij, M.n, h.ij, h.n, N.ij, N.n);
332
333 if (north + south > 0) {
334 h_y = 1.0 / ((south + north) * dy) * (south * (h.ij - h.s) + north * (h.n - h.ij));
335 } else {
336 h_y = 0.0;
337 }
338 }
339
340 result(i, j) = - pressure * Vector2(h_x, h_y);
341 }
342 }
343
344 std::string SSA::stdout_report() const {
345 return m_stdout_ssa;
346 }
347
348
349 //! \brief Set the initial guess of the SSA velocity.
350 void SSA::set_initial_guess(const IceModelVec2V &guess) {
351 m_velocity.copy_from(guess);
352 }
353
354 const IceModelVec2V& SSA::driving_stress() const {
355 return m_taud;
356 }
357
358
359 void SSA::define_model_state_impl(const File &output) const {
360 m_velocity.define(output);
361 }
362
363 void SSA::write_model_state_impl(const File &output) const {
364 m_velocity.write(output);
365 }
366
367 DiagnosticList SSA::diagnostics_impl() const {
368 DiagnosticList result = ShallowStressBalance::diagnostics_impl();
369
370 // replace these diagnostics
371 result["taud"] = Diagnostic::Ptr(new SSA_taud(this));
372 result["taud_mag"] = Diagnostic::Ptr(new SSA_taud_mag(this));
373
374 return result;
375 }
376
377 SSA_taud::SSA_taud(const SSA *m)
378 : Diag<SSA>(m) {
379
380 // set metadata:
381 m_vars = {SpatialVariableMetadata(m_sys, "taud_x"),
382 SpatialVariableMetadata(m_sys, "taud_y")};
383
384 set_attrs("X-component of the driving shear stress at the base of ice", "",
385 "Pa", "Pa", 0);
386 set_attrs("Y-component of the driving shear stress at the base of ice", "",
387 "Pa", "Pa", 1);
388
389 for (auto &v : m_vars) {
390 v.set_string("comment",
391 "this is the driving stress used by the SSA solver");
392 }
393 }
394
395 IceModelVec::Ptr SSA_taud::compute_impl() const {
396
397 IceModelVec2V::Ptr result(new IceModelVec2V(m_grid, "result", WITHOUT_GHOSTS));
398 result->metadata(0) = m_vars[0];
399 result->metadata(1) = m_vars[1];
400
401 result->copy_from(model->driving_stress());
402
403 return result;
404 }
405
406 SSA_taud_mag::SSA_taud_mag(const SSA *m)
407 : Diag<SSA>(m) {
408
409 // set metadata:
410 m_vars = {SpatialVariableMetadata(m_sys, "taud_mag")};
411
412 set_attrs("magnitude of the driving shear stress at the base of ice", "",
413 "Pa", "Pa", 0);
414 m_vars[0].set_string("comment",
415 "this is the magnitude of the driving stress used by the SSA solver");
416 }
417
418 IceModelVec::Ptr SSA_taud_mag::compute_impl() const {
419
420 // Allocate memory:
421 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "taud_mag", WITHOUT_GHOSTS));
422 result->metadata() = m_vars[0];
423
424 result->set_to_magnitude(model->driving_stress());
425
426 return result;
427 }
428
429 } // end of namespace stressbalance
430 } // end of namespace pism