tShallowStressBalance.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
---
tShallowStressBalance.cc (12804B)
---
1 // Copyright (C) 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 Constantine Khroulev and Ed Bueler
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 "ShallowStressBalance.hh"
20 #include "pism/basalstrength/basal_resistance.hh"
21 #include "pism/rheology/FlowLawFactory.hh"
22
23 #include "pism/util/Mask.hh"
24 #include "pism/util/Vars.hh"
25 #include "pism/util/error_handling.hh"
26 #include "pism/util/pism_options.hh"
27 #include "pism/util/IceModelVec2CellType.hh"
28
29 #include "SSB_diagnostics.hh"
30
31 namespace pism {
32 namespace stressbalance {
33
34 //! Evaluate the margin pressure difference term in the calving-front BC.
35 //
36 // Units: (kg / m3) * (m / s2) * m2 = Pa m
37 double margin_pressure_difference(bool shelf, bool dry_mode, double H, double bed,
38 double sea_level, double rho_ice, double rho_ocean,
39 double g) {
40 if (shelf) {
41 // floating shelf
42 return 0.5 * rho_ice * g * (1.0 - (rho_ice / rho_ocean)) * H * H;
43 } else {
44 // grounded terminus
45 if (bed >= sea_level or dry_mode) {
46 return 0.5 * rho_ice * g * H * H;
47 } else {
48 return 0.5 * rho_ice * g * (H * H - (rho_ocean / rho_ice) * pow(sea_level - bed, 2.0));
49 }
50 }
51 }
52
53 using pism::mask::ice_free;
54
55 ShallowStressBalance::ShallowStressBalance(IceGrid::ConstPtr g)
56 : Component(g), m_basal_sliding_law(NULL), m_flow_law(NULL), m_EC(g->ctx()->enthalpy_converter()) {
57
58 const unsigned int WIDE_STENCIL = m_config->get_number("grid.max_stencil_width");
59
60 if (m_config->get_flag("basal_resistance.pseudo_plastic.enabled") == true) {
61 m_basal_sliding_law = new IceBasalResistancePseudoPlasticLaw(*m_config);
62 } else {
63 m_basal_sliding_law = new IceBasalResistancePlasticLaw(*m_config);
64 }
65
66 m_velocity.create(m_grid, "bar", WITH_GHOSTS, WIDE_STENCIL); // components ubar, vbar
67 m_velocity.set_attrs("model_state",
68 "thickness-advective ice velocity (x-component)",
69 "m s-1", "m s-1", "", 0);
70 m_velocity.set_attrs("model_state",
71 "thickness-advective ice velocity (y-component)",
72 "m s-1", "m s-1", "", 1);
73
74 m_basal_frictional_heating.create(m_grid, "bfrict", WITHOUT_GHOSTS);
75 m_basal_frictional_heating.set_attrs("diagnostic",
76 "basal frictional heating",
77 "W m-2", "mW m-2", "", 0);
78 }
79
80 ShallowStressBalance::~ShallowStressBalance() {
81 delete m_basal_sliding_law;
82 }
83
84 void ShallowStressBalance::init() {
85 this->init_impl();
86 }
87
88 void ShallowStressBalance::init_impl() {
89 // empty
90 }
91
92 std::string ShallowStressBalance::stdout_report() const {
93 return "";
94 }
95
96 std::shared_ptr<const rheology::FlowLaw> ShallowStressBalance::flow_law() const {
97 return m_flow_law;
98 }
99
100 EnthalpyConverter::Ptr ShallowStressBalance::enthalpy_converter() const {
101 return m_EC;
102 }
103
104 const IceBasalResistancePlasticLaw* ShallowStressBalance::sliding_law() const {
105 return m_basal_sliding_law;
106 }
107
108 //! \brief Get the thickness-advective 2D velocity.
109 const IceModelVec2V& ShallowStressBalance::velocity() const {
110 return m_velocity;
111 }
112
113 //! \brief Get the basal frictional heating (for the adaptive energy time-stepping).
114 const IceModelVec2S& ShallowStressBalance::basal_frictional_heating() {
115 return m_basal_frictional_heating;
116 }
117
118
119 DiagnosticList ShallowStressBalance::diagnostics_impl() const {
120 DiagnosticList result = {
121 {"beta", Diagnostic::Ptr(new SSB_beta(this))},
122 {"taub", Diagnostic::Ptr(new SSB_taub(this))},
123 {"taub_mag", Diagnostic::Ptr(new SSB_taub_mag(this))},
124 {"taud", Diagnostic::Ptr(new SSB_taud(this))},
125 {"taud_mag", Diagnostic::Ptr(new SSB_taud_mag(this))}
126 };
127
128 if(m_config->get_flag("output.ISMIP6")) {
129 result["strbasemag"] = Diagnostic::Ptr(new SSB_taub_mag(this));
130 }
131
132 return result;
133 }
134
135
136 ZeroSliding::ZeroSliding(IceGrid::ConstPtr g)
137 : ShallowStressBalance(g) {
138
139 // Use the SIA flow law.
140 rheology::FlowLawFactory ice_factory("stress_balance.sia.", m_config, m_EC);
141 m_flow_law = ice_factory.create();
142 }
143
144 ZeroSliding::~ZeroSliding() {
145 // empty
146 }
147
148 //! \brief Update the trivial shallow stress balance object.
149 void ZeroSliding::update(const Inputs &inputs, bool full_update) {
150 (void) inputs;
151
152 if (full_update) {
153 m_velocity.set(0.0);
154 m_basal_frictional_heating.set(0.0);
155 }
156 }
157
158 //! \brief Compute the basal frictional heating.
159 /*!
160 Ice shelves have zero basal friction heating.
161
162 \param[in] V *basal* sliding velocity
163 \param[in] tauc basal yield stress
164 \param[in] mask (used to determine if floating or grounded)
165 \param[out] result
166 */
167 void ShallowStressBalance::compute_basal_frictional_heating(const IceModelVec2V &V,
168 const IceModelVec2S &tauc,
169 const IceModelVec2CellType &mask,
170 IceModelVec2S &result) const {
171
172 IceModelVec::AccessList list{&V, &result, &tauc, &mask};
173
174 for (Points p(*m_grid); p; p.next()) {
175 const int i = p.i(), j = p.j();
176
177 if (mask.ocean(i,j)) {
178 result(i,j) = 0.0;
179 } else {
180 const double
181 C = m_basal_sliding_law->drag(tauc(i,j), V(i,j).u, V(i,j).v),
182 basal_stress_x = - C * V(i,j).u,
183 basal_stress_y = - C * V(i,j).v;
184 result(i,j) = - basal_stress_x * V(i,j).u - basal_stress_y * V(i,j).v;
185 }
186 }
187 }
188
189
190 SSB_taud::SSB_taud(const ShallowStressBalance *m)
191 : Diag<ShallowStressBalance>(m) {
192
193 // set metadata:
194 m_vars = {SpatialVariableMetadata(m_sys, "taud_x"),
195 SpatialVariableMetadata(m_sys, "taud_y")};
196
197 set_attrs("X-component of the driving shear stress at the base of ice", "",
198 "Pa", "Pa", 0);
199 set_attrs("Y-component of the driving shear stress at the base of ice", "",
200 "Pa", "Pa", 1);
201
202 for (auto &v : m_vars) {
203 v.set_string("comment",
204 "this field is purely diagnostic (not used by the model)");
205 }
206 }
207
208 /*!
209 * The driving stress computed here is not used by the model, so this
210 * implementation intentionally does not use the eta-transformation or special
211 * cases at ice margins.
212 */
213 IceModelVec::Ptr SSB_taud::compute_impl() const {
214
215 IceModelVec2V::Ptr result(new IceModelVec2V);
216 result->create(m_grid, "result", WITHOUT_GHOSTS);
217 result->metadata(0) = m_vars[0];
218 result->metadata(1) = m_vars[1];
219
220 const IceModelVec2S *thickness = m_grid->variables().get_2d_scalar("land_ice_thickness");
221 const IceModelVec2S *surface = m_grid->variables().get_2d_scalar("surface_altitude");
222
223 double standard_gravity = m_config->get_number("constants.standard_gravity"),
224 ice_density = m_config->get_number("constants.ice.density");
225
226 IceModelVec::AccessList list{surface, thickness, result.get()};
227
228 for (Points p(*m_grid); p; p.next()) {
229 const int i = p.i(), j = p.j();
230
231 double pressure = ice_density * standard_gravity * (*thickness)(i,j);
232 if (pressure <= 0.0) {
233 (*result)(i,j).u = 0.0;
234 (*result)(i,j).v = 0.0;
235 } else {
236 (*result)(i,j).u = - pressure * surface->diff_x_p(i,j);
237 (*result)(i,j).v = - pressure * surface->diff_y_p(i,j);
238 }
239 }
240
241 return result;
242 }
243
244 SSB_taud_mag::SSB_taud_mag(const ShallowStressBalance *m)
245 : Diag<ShallowStressBalance>(m) {
246
247 // set metadata:
248 m_vars = {SpatialVariableMetadata(m_sys, "taud_mag")};
249
250 set_attrs("magnitude of the gravitational driving stress at the base of ice", "",
251 "Pa", "Pa", 0);
252 m_vars[0].set_string("comment",
253 "this field is purely diagnostic (not used by the model)");
254 }
255
256 IceModelVec::Ptr SSB_taud_mag::compute_impl() const {
257 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "taud_mag", WITHOUT_GHOSTS));
258 result->metadata(0) = m_vars[0];
259
260 IceModelVec2V::Ptr taud = IceModelVec2V::ToVector(SSB_taud(model).compute());
261
262 result->set_to_magnitude(*taud);
263
264 return result;
265 }
266
267 SSB_taub::SSB_taub(const ShallowStressBalance *m)
268 : Diag<ShallowStressBalance>(m) {
269 // set metadata:
270 m_vars = {SpatialVariableMetadata(m_sys, "taub_x"),
271 SpatialVariableMetadata(m_sys, "taub_y")};
272
273 set_attrs("X-component of the shear stress at the base of ice", "",
274 "Pa", "Pa", 0);
275 set_attrs("Y-component of the shear stress at the base of ice", "",
276 "Pa", "Pa", 1);
277
278 for (auto &v : m_vars) {
279 v.set_string("comment",
280 "this field is purely diagnostic (not used by the model)");
281 }
282 }
283
284
285 IceModelVec::Ptr SSB_taub::compute_impl() const {
286
287 IceModelVec2V::Ptr result(new IceModelVec2V);
288 result->create(m_grid, "result", WITHOUT_GHOSTS);
289 result->metadata() = m_vars[0];
290 result->metadata(1) = m_vars[1];
291
292 const IceModelVec2V &velocity = model->velocity();
293 const IceModelVec2S *tauc = m_grid->variables().get_2d_scalar("tauc");
294 const IceModelVec2CellType &mask = *m_grid->variables().get_2d_cell_type("mask");
295
296 const IceBasalResistancePlasticLaw *basal_sliding_law = model->sliding_law();
297
298 IceModelVec::AccessList list{tauc, &velocity, &mask, result.get()};
299 for (Points p(*m_grid); p; p.next()) {
300 const int i = p.i(), j = p.j();
301
302 if (mask.grounded_ice(i,j)) {
303 double beta = basal_sliding_law->drag((*tauc)(i,j), velocity(i,j).u, velocity(i,j).v);
304 (*result)(i,j).u = - beta * velocity(i,j).u;
305 (*result)(i,j).v = - beta * velocity(i,j).v;
306 } else {
307 (*result)(i,j).u = 0.0;
308 (*result)(i,j).v = 0.0;
309 }
310 }
311
312 return result;
313 }
314
315 SSB_taub_mag::SSB_taub_mag(const ShallowStressBalance *m)
316 : Diag<ShallowStressBalance>(m) {
317
318 auto ismip6 = m_config->get_flag("output.ISMIP6");
319
320 // set metadata:
321 m_vars = {SpatialVariableMetadata(m_sys, ismip6 ? "strbasemag" : "taub_mag")};
322
323 set_attrs("magnitude of the basal shear stress at the base of ice",
324 "land_ice_basal_drag", // ISMIP6 "standard" name
325 "Pa", "Pa", 0);
326 m_vars[0].set_string("comment",
327 "this field is purely diagnostic (not used by the model)");
328 }
329
330 IceModelVec::Ptr SSB_taub_mag::compute_impl() const {
331 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "taub_mag", WITHOUT_GHOSTS));
332 result->metadata(0) = m_vars[0];
333
334 IceModelVec2V::Ptr taub = IceModelVec2V::ToVector(SSB_taub(model).compute());
335
336 result->set_to_magnitude(*taub);
337
338 return result;
339 }
340
341 /**
342 * Shallow stress balance class that reads `u` and `v` fields from a
343 * file and holds them constant.
344 *
345 * The only use I can think of right now is testing.
346 */
347 PrescribedSliding::PrescribedSliding(IceGrid::ConstPtr g)
348 : ZeroSliding(g) {
349 // empty
350 }
351
352 PrescribedSliding::~PrescribedSliding() {
353 // empty
354 }
355
356 void PrescribedSliding::update(const Inputs &inputs, bool full_update) {
357 (void) inputs;
358 if (full_update) {
359 m_basal_frictional_heating.set(0.0);
360 }
361 }
362
363 void PrescribedSliding::init_impl() {
364 ShallowStressBalance::init_impl();
365
366 auto input_filename = m_config->get_string("stress_balance.prescribed_sliding.file");
367
368 if (input_filename.empty()) {
369 throw RuntimeError(PISM_ERROR_LOCATION,
370 "stress_balance.prescribed_sliding.file is required.");
371 }
372
373 m_velocity.regrid(input_filename, CRITICAL);
374 }
375
376 SSB_beta::SSB_beta(const ShallowStressBalance *m)
377 : Diag<ShallowStressBalance>(m) {
378 // set metadata:
379 m_vars = {SpatialVariableMetadata(m_sys, "beta")};
380
381 set_attrs("basal drag coefficient", "", "Pa s / m", "Pa s / m", 0);
382 }
383
384 IceModelVec::Ptr SSB_beta::compute_impl() const {
385 IceModelVec2S::Ptr result(new IceModelVec2S(m_grid, "beta", WITHOUT_GHOSTS));
386 result->metadata(0) = m_vars[0];
387
388 const IceModelVec2S *tauc = m_grid->variables().get_2d_scalar("tauc");
389
390 const IceBasalResistancePlasticLaw *basal_sliding_law = model->sliding_law();
391
392 const IceModelVec2V &velocity = model->velocity();
393
394 IceModelVec::AccessList list{tauc, &velocity, result.get()};
395 for (Points p(*m_grid); p; p.next()) {
396 const int i = p.i(), j = p.j();
397
398 (*result)(i,j) = basal_sliding_law->drag((*tauc)(i,j), velocity(i,j).u, velocity(i,j).v);
399 }
400
401 return result;
402 }
403
404 } // end of namespace stressbalance
405 } // end of namespace pism