tFractureDensity.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
---
tFractureDensity.cc (17423B)
---
1 /* Copyright (C) 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 "FractureDensity.hh"
21 #include "pism/geometry/Geometry.hh"
22 #include "pism/stressbalance/StressBalance.hh"
23 #include "pism/util/pism_options.hh"
24 #include "pism/util/pism_utilities.hh"
25
26 namespace pism {
27
28 FractureDensity::FractureDensity(IceGrid::ConstPtr grid,
29 std::shared_ptr<const rheology::FlowLaw> flow_law)
30 : Component(grid),
31 m_density(grid, "fracture_density", WITH_GHOSTS, 1),
32 m_density_new(grid, "new_fracture_density", WITHOUT_GHOSTS),
33 m_growth_rate(grid, "fracture_growth_rate", WITHOUT_GHOSTS),
34 m_healing_rate(grid, "fracture_healing_rate", WITHOUT_GHOSTS),
35 m_flow_enhancement(grid, "fracture_flow_enhancement", WITHOUT_GHOSTS),
36 m_age(grid, "fracture_age", WITH_GHOSTS, 1),
37 m_age_new(grid, "new_fracture_age", WITHOUT_GHOSTS),
38 m_toughness(grid, "fracture_toughness", WITHOUT_GHOSTS),
39 m_strain_rates(grid, "strain_rates", WITHOUT_GHOSTS,
40 2, // stencil width
41 2), // dof
42 m_deviatoric_stresses(grid, "sigma", WITHOUT_GHOSTS,
43 0, // stencil width
44 3), // dof
45 m_velocity(grid, "ghosted_velocity", WITH_GHOSTS, 1),
46 m_flow_law(flow_law) {
47
48 m_density.set_attrs("model_state", "fracture density in ice shelf", "1", "1", "", 0);
49 m_density.metadata().set_number("valid_max", 1.0);
50 m_density.metadata().set_number("valid_min", 0.0);
51
52 m_growth_rate.set_attrs("model_state", "fracture growth rate", "second-1", "second-1", "", 0);
53 m_growth_rate.metadata().set_number("valid_min", 0.0);
54
55 m_healing_rate.set_attrs("model_state", "fracture healing rate", "second-1", "second-1", "", 0);
56
57 m_flow_enhancement.set_attrs("model_state", "fracture-induced flow enhancement", "", "", "", 0);
58
59 m_age.set_attrs("model_state", "age since fracturing", "seconds", "seconds", "", 0);
60 m_age.metadata().set_string("glaciological_units", "years");
61
62 m_toughness.set_attrs("model_state", "fracture toughness", "Pa", "Pa", "", 0);
63
64 m_strain_rates.metadata(0).set_name("eigen1");
65 m_strain_rates.set_attrs("internal",
66 "major principal component of horizontal strain-rate",
67 "second-1", "second-1", "", 0);
68
69 m_strain_rates.metadata(1).set_name("eigen2");
70 m_strain_rates.set_attrs("internal",
71 "minor principal component of horizontal strain-rate",
72 "second-1", "second-1", "", 1);
73
74 m_deviatoric_stresses.metadata(0).set_name("sigma_xx");
75 m_deviatoric_stresses.set_attrs("internal",
76 "deviatoric stress in x direction",
77 "Pa", "Pa", "", 0);
78
79 m_deviatoric_stresses.metadata(1).set_name("sigma_yy");
80 m_deviatoric_stresses.set_attrs("internal",
81 "deviatoric stress in y direction",
82 "Pa", "Pa", "", 1);
83
84 m_deviatoric_stresses.metadata(2).set_name("sigma_xy");
85 m_deviatoric_stresses.set_attrs("internal",
86 "deviatoric shear stress",
87 "Pa", "Pa", "", 2);
88 }
89
90 FractureDensity::~FractureDensity() {
91 // empty
92 }
93
94 void FractureDensity::restart(const File &input_file, int record) {
95 m_log->message(2, "* Restarting the fracture density model from %s...\n",
96 input_file.filename().c_str());
97
98 m_density.read(input_file, record);
99 m_age.read(input_file, record);
100
101 regrid("Fracture density model", m_density, REGRID_WITHOUT_REGRID_VARS);
102 regrid("Fracture density model", m_age, REGRID_WITHOUT_REGRID_VARS);
103 }
104
105 void FractureDensity::bootstrap(const File &input_file) {
106 m_log->message(2, "* Bootstrapping the fracture density model from %s...\n",
107 input_file.filename().c_str());
108
109 m_density.regrid(input_file, OPTIONAL, 0.0);
110 m_age.regrid(input_file, OPTIONAL, 0.0);
111 }
112
113 void FractureDensity::initialize(const IceModelVec2S &density,
114 const IceModelVec2S &age) {
115 m_density.copy_from(density);
116 m_age.copy_from(age);
117 }
118
119 void FractureDensity::initialize() {
120 m_density.set(0.0);
121 m_age.set(0.0);
122 }
123
124 void FractureDensity::define_model_state_impl(const File &output) const {
125 m_density.define(output);
126 m_age.define(output);
127 }
128
129 void FractureDensity::write_model_state_impl(const File &output) const {
130 m_density.write(output);
131 m_age.write(output);
132 }
133
134 void FractureDensity::update(double dt,
135 const Geometry &geometry,
136 const IceModelVec2V &velocity,
137 const IceModelVec2S &hardness,
138 const IceModelVec2S &bc_mask) {
139 const double
140 dx = m_grid->dx(),
141 dy = m_grid->dy();
142 const int
143 Mx = m_grid->Mx(),
144 My = m_grid->My();
145
146 IceModelVec2S
147 &D = m_density,
148 &D_new = m_density_new,
149 &A = m_age,
150 &A_new = m_age_new;
151
152 m_velocity.copy_from(velocity);
153
154 stressbalance::compute_2D_principal_strain_rates(m_velocity,
155 geometry.cell_type,
156 m_strain_rates);
157
158 stressbalance::compute_2D_stresses(*m_flow_law,
159 m_velocity,
160 hardness,
161 geometry.cell_type,
162 m_deviatoric_stresses);
163
164 IceModelVec::AccessList list{&m_velocity, &m_strain_rates, &m_deviatoric_stresses,
165 &D, &D_new, &geometry.cell_type, &bc_mask, &A, &A_new,
166 &m_growth_rate, &m_healing_rate, &m_flow_enhancement,
167 &m_toughness};
168
169 D_new.copy_from(D);
170
171 //options
172 /////////////////////////////////////////////////////////
173 double soft_residual = options::Real("-fracture_softening", "soft_residual", 1.0);
174 // assume linear response function: E_fr = (1-(1-soft_residual)*phi) -> 1-phi
175 //
176 // more: T. Albrecht, A. Levermann; Fracture-induced softening for
177 // large-scale ice dynamics; (2013), The Cryosphere Discussions 7;
178 // 4501-4544; DOI:10.5194/tcd-7-4501-2013
179
180 // get four options for calculation of fracture density.
181 // 1st: fracture growth constant gamma
182 // 2nd: fracture initiation stress threshold sigma_cr
183 // 3rd: healing rate constant gamma_h
184 // 4th: healing strain rate threshold
185 // more: T. Albrecht, A. Levermann; Fracture field for large-scale
186 // ice dynamics; (2012), Journal of Glaciology, Vol. 58, No. 207,
187 // 165-176, DOI: 10.3189/2012JoG11J191.
188
189 double gamma = 1.0, initThreshold = 7.0e4, gammaheal = 0.0, healThreshold = 2.0e-10;
190 {
191 options::RealList fractures("-fracture_parameters",
192 "gamma, initThreshold, gammaheal, healThreshold",
193 {gamma, initThreshold, gammaheal, healThreshold});
194 if (fractures->size() != 4) {
195 throw RuntimeError(PISM_ERROR_LOCATION, "option -fracture_parameters requires exactly 4 arguments");
196 }
197 gamma = fractures[0];
198 initThreshold = fractures[1];
199 gammaheal = fractures[2];
200 healThreshold = fractures[3];
201 }
202
203 m_log->message(3, "PISM-PIK INFO: fracture density is found with parameters:\n"
204 " gamma=%.2f, sigma_cr=%.2f, gammah=%.2f, healing_cr=%.1e and soft_res=%f \n",
205 gamma, initThreshold, gammaheal, healThreshold, soft_residual);
206
207 bool do_fracground = m_config->get_flag("fracture_density.include_grounded_ice");
208
209 double fdBoundaryValue = m_config->get_number("fracture_density.phi0");
210
211 bool constant_healing = m_config->get_flag("fracture_density.constant_healing");
212
213 bool fracture_weighted_healing = m_config->get_flag("fracture_density.fracture_weighted_healing");
214
215 bool max_shear_stress = m_config->get_flag("fracture_density.max_shear_stress");
216
217 bool lefm = m_config->get_flag("fracture_density.lefm");
218
219 bool constant_fd = m_config->get_flag("fracture_density.constant_fd");
220
221 bool fd2d_scheme = m_config->get_flag("fracture_density.fd2d_scheme");
222
223 for (Points p(*m_grid); p; p.next()) {
224 const int i = p.i(), j = p.j();
225
226 double tempFD = 0.0;
227
228 double u = m_velocity(i, j).u;
229 double v = m_velocity(i, j).v;
230
231 if (fd2d_scheme) {
232 if (u >= dx * v / dy and v >= 0.0) { //1
233 tempFD = u * (D(i, j) - D(i - 1, j)) / dx + v * (D(i - 1, j) - D(i - 1, j - 1)) / dy;
234 } else if (u <= dx * v / dy and u >= 0.0) { //2
235 tempFD = u * (D(i, j - 1) - D(i - 1, j - 1)) / dx + v * (D(i, j) - D(i, j - 1)) / dy;
236 } else if (u >= -dx * v / dy and u <= 0.0) { //3
237 tempFD = -u * (D(i, j - 1) - D(i + 1, j - 1)) / dx + v * (D(i, j) - D(i, j - 1)) / dy;
238 } else if (u <= -dx * v / dy and v >= 0.0) { //4
239 tempFD = -u * (D(i, j) - D(i + 1, j)) / dx + v * (D(i + 1, j) - D(i + 1, j - 1)) / dy;
240 } else if (u <= dx * v / dy and v <= 0.0) { //5
241 tempFD = -u * (D(i, j) - D(i + 1, j)) / dx - v * (D(i + 1, j) - D(i + 1, j + 1)) / dy;
242 } else if (u >= dx * v / dy and u <= 0.0) { //6
243 tempFD = -u * (D(i, j + 1) - D(i + 1, j + 1)) / dx - v * (D(i, j) - D(i, j + 1)) / dy;
244 } else if (u <= -dx * v / dy and u >= 0.0) { //7
245 tempFD = u * (D(i, j + 1) - D(i - 1, j + 1)) / dx - v * (D(i, j) - D(i, j + 1)) / dy;
246 } else if (u >= -dx * v / dy and v <= 0.0) { //8
247 tempFD = u * (D(i, j) - D(i - 1, j)) / dx - v * (D(i - 1, j) - D(i - 1, j + 1)) / dy;
248 } else {
249 m_log->message(3, "######### missing case of angle %f of %f and %f at %d, %d \n",
250 atan(v / u) / M_PI * 180., u * 3e7, v * 3e7, i, j);
251 }
252 } else {
253 tempFD += u * (u < 0 ? D(i + 1, j) - D(i, j) : D(i, j) - D(i - 1, j)) / dx;
254 tempFD += v * (v < 0 ? D(i, j + 1) - D(i, j) : D(i, j) - D(i, j - 1)) / dy;
255 }
256
257 D_new(i, j) -= tempFD * dt;
258
259 //sources /////////////////////////////////////////////////////////////////
260 ///von mises criterion
261
262 double
263 txx = m_deviatoric_stresses(i, j, 0),
264 tyy = m_deviatoric_stresses(i, j, 1),
265 txy = m_deviatoric_stresses(i, j, 2),
266 T1 = 0.5 * (txx + tyy) + sqrt(0.25 * PetscSqr(txx - tyy) + PetscSqr(txy)), //Pa
267 T2 = 0.5 * (txx + tyy) - sqrt(0.25 * PetscSqr(txx - tyy) + PetscSqr(txy)), //Pa
268 sigmat = sqrt(PetscSqr(T1) + PetscSqr(T2) - T1 * T2);
269
270
271 ///max shear stress criterion (more stringent than von mises)
272 if (max_shear_stress) {
273 double maxshear = fabs(T1);
274 maxshear = std::max(maxshear, fabs(T2));
275 maxshear = std::max(maxshear, fabs(T1 - T2));
276
277 sigmat = maxshear;
278 }
279
280 ///lefm mixed-mode criterion
281 if (lefm) {
282 double sigmamu = 0.1; //friction coefficient between crack faces
283
284 double sigmac = 0.64 / M_PI; //initial crack depth 20cm
285
286 double sigmabetatest, sigmanor, sigmatau, Kone, Ktwo, KSI, KSImax = 0.0, sigmatetanull;
287
288 for (int l = 46; l <= 90; ++l) { //optimize for various precursor angles beta
289 sigmabetatest = l * M_PI / 180.0;
290
291 //rist_sammonds99
292 sigmanor = 0.5 * (T1 + T2) - (T1 - T2) * cos(2 * sigmabetatest);
293 sigmatau = 0.5 * (T1 - T2) * sin(2 * sigmabetatest);
294 //shayam_wu90
295 if (sigmamu * sigmanor < 0.0) { //compressive case
296 if (fabs(sigmatau) <= fabs(sigmamu * sigmanor)) {
297 sigmatau = 0.0;
298 } else {
299 if (sigmatau > 0) { //coulomb friction opposing sliding
300 sigmatau += (sigmamu * sigmanor);
301 } else {
302 sigmatau -= (sigmamu * sigmanor);
303 }
304 }
305 }
306
307 //stress intensity factors
308 Kone = sigmanor * sqrt(M_PI * sigmac); //normal
309 Ktwo = sigmatau * sqrt(M_PI * sigmac); //shear
310
311 if (Ktwo == 0.0) {
312 sigmatetanull = 0.0;
313 } else { //eq15 in hulbe_ledoux10 or eq15 shayam_wu90
314 sigmatetanull = -2.0 * atan((sqrt(PetscSqr(Kone) + 8.0 * PetscSqr(Ktwo)) - Kone) / (4.0 * Ktwo));
315 }
316
317 KSI = cos(0.5 * sigmatetanull) *
318 (Kone * cos(0.5 * sigmatetanull) * cos(0.5 * sigmatetanull) - 0.5 * 3.0 * Ktwo * sin(sigmatetanull));
319 // mode I stress intensity
320
321 KSImax = std::max(KSI, KSImax);
322 }
323 sigmat = KSImax;
324 }
325
326 //////////////////////////////////////////////////////////////////////////////
327
328 //fracture density
329 double fdnew = gamma * (m_strain_rates(i, j, 0) - 0.0) * (1 - D_new(i, j));
330 if (sigmat > initThreshold) {
331 D_new(i, j) += fdnew * dt;
332 }
333
334 //healing
335 double fdheal = gammaheal * (m_strain_rates(i, j, 0) - healThreshold);
336 if (geometry.cell_type.icy(i, j)) {
337 if (constant_healing) {
338 fdheal = gammaheal * (-healThreshold);
339 if (fracture_weighted_healing) {
340 D_new(i, j) += fdheal * dt * (1 - D(i, j));
341 } else {
342 D_new(i, j) += fdheal * dt;
343 }
344 } else if (m_strain_rates(i, j, 0) < healThreshold) {
345 if (fracture_weighted_healing) {
346 D_new(i, j) += fdheal * dt * (1 - D(i, j));
347 } else {
348 D_new(i, j) += fdheal * dt;
349 }
350 }
351 }
352
353 // bounding
354 D_new(i, j) = pism::clip(D_new(i, j), 0.0, 1.0);
355
356 if (geometry.cell_type.icy(i, j)) {
357 //fracture toughness
358 m_toughness(i, j) = sigmat;
359
360 // fracture growth rate
361 if (sigmat > initThreshold) {
362 m_growth_rate(i, j) = fdnew;
363 //m_growth_rate(i,j)=gamma*(vPrinStrain1(i,j)-0.0)*(1-D_new(i,j));
364 } else {
365 m_growth_rate(i, j) = 0.0;
366 }
367
368 // fracture healing rate
369 if (geometry.cell_type.icy(i, j)) {
370 if (constant_healing or (m_strain_rates(i, j, 0) < healThreshold)) {
371 if (fracture_weighted_healing) {
372 m_healing_rate(i, j) = fdheal * (1 - D(i, j));
373 } else {
374 m_healing_rate(i, j) = fdheal;
375 }
376 } else {
377 m_healing_rate(i, j) = 0.0;
378 }
379 }
380
381 // fracture age since fracturing occurred
382 {
383 auto a = A.star(i, j);
384 A_new(i, j) -= dt * u * (u < 0 ? a.e - a.ij : a.ij - a.w) / dx;
385 A_new(i, j) -= dt * v * (v < 0 ? a.n - a.ij : a.ij - a.s) / dy;
386 A_new(i, j) += dt;
387 if (sigmat > initThreshold) {
388 A_new(i, j) = 0.0;
389 }
390 }
391
392 // additional flow enhancement due to fracture softening
393 double phi_exp = 3.0; //flow_law->exponent();
394 double softening = pow((1.0 - (1.0 - soft_residual) * D_new(i, j)), -phi_exp);
395 if (geometry.cell_type.icy(i, j)) {
396 m_flow_enhancement(i, j) = 1.0 / pow(softening, 1 / 3.0);
397 } else {
398 m_flow_enhancement(i, j) = 1.0;
399 }
400 }
401
402 // boundary condition
403 if (geometry.cell_type.grounded(i, j) and not do_fracground) {
404
405 if (bc_mask(i, j) > 0.5) {
406 D_new(i, j) = fdBoundaryValue;
407
408 {
409 A_new(i, j) = 0.0;
410 m_growth_rate(i, j) = 0.0;
411 m_healing_rate(i, j) = 0.0;
412 m_flow_enhancement(i, j) = 1.0;
413 m_toughness(i, j) = 0.0;
414 }
415 }
416 }
417
418 // ice free regions and boundary of computational domain
419 if (geometry.cell_type.ice_free(i, j) or i == 0 or j == 0 or i == Mx - 1 or j == My - 1) {
420
421 D_new(i, j) = 0.0;
422
423 {
424 A_new(i, j) = 0.0;
425 m_growth_rate(i, j) = 0.0;
426 m_healing_rate(i, j) = 0.0;
427 m_flow_enhancement(i, j) = 1.0;
428 m_toughness(i, j) = 0.0;
429 }
430 }
431
432 if (constant_fd) { // no fd evolution
433 D_new(i, j) = D(i, j);
434 }
435 }
436
437 A_new.update_ghosts(A);
438 D_new.update_ghosts(D);
439 }
440
441 DiagnosticList FractureDensity::diagnostics_impl() const {
442 return {{"fracture_density", Diagnostic::wrap(m_density)},
443 {"fracture_growth_rate", Diagnostic::wrap(m_growth_rate)},
444 {"fracture_healing_rate", Diagnostic::wrap(m_healing_rate)},
445 {"fracture_flow_enhancement", Diagnostic::wrap(m_flow_enhancement)},
446 {"fracture_age", Diagnostic::wrap(m_age)},
447 {"fracture_toughness", Diagnostic::wrap(m_toughness)}
448 };
449 }
450
451 const IceModelVec2S& FractureDensity::density() const {
452 return m_density;
453 }
454
455 const IceModelVec2S& FractureDensity::growth_rate() const {
456 return m_growth_rate;
457 }
458
459 const IceModelVec2S& FractureDensity::healing_rate() const {
460 return m_healing_rate;
461 }
462
463 const IceModelVec2S& FractureDensity::flow_enhancement() const {
464 return m_flow_enhancement;
465 }
466
467 const IceModelVec2S& FractureDensity::age() const {
468 return m_age;
469 }
470
471 const IceModelVec2S& FractureDensity::toughness() const {
472 return m_toughness;
473 }
474
475 } // end of namespace pism