tLingleClarkSerial.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
---
tLingleClarkSerial.cc (18276B)
---
1 // Copyright (C) 2004-2009, 2011, 2013, 2014, 2015, 2016, 2017, 2018, 2019 Ed Bueler and Constantine Khroulev
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 <cassert>
20 #include <cmath> // sqrt
21 #include <fftw3.h>
22 #include <gsl/gsl_math.h> // M_PI
23
24 #include "matlablike.hh"
25 #include "greens.hh"
26 #include "LingleClarkSerial.hh"
27
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/util/ConfigInterface.hh"
30 #include "pism/util/error_handling.hh"
31 #include "pism/util/petscwrappers/Vec.hh"
32 #include "pism/util/fftw_utilities.hh"
33
34 namespace pism {
35 namespace bed {
36
37 /*!
38 * @param[in] config configuration database
39 * @param[in] include_elastic include elastic deformation component
40 * @param[in] Mx grid size in the X direction
41 * @param[in] My grid size in the Y direction
42 * @param[in] dx grid spacing in the X direction
43 * @param[in] dy grid spacing in the Y direction
44 * @param[in] Nx extended grid size in the X direction
45 * @param[in] Ny extended grid size in the Y direction
46 */
47 LingleClarkSerial::LingleClarkSerial(Logger::ConstPtr log,
48 const Config &config,
49 bool include_elastic,
50 int Mx, int My,
51 double dx, double dy,
52 int Nx, int Ny)
53 : m_log(log) {
54
55 // set parameters
56 m_include_elastic = include_elastic;
57
58 if (include_elastic) {
59 // check if the extended grid is large enough (it has to be at least twice the size of
60 // the physical grid so that the load in one corner of the domain affects the grid
61 // point in the opposite corner).
62
63 if (config.get_number("bed_deformation.lc.grid_size_factor") < 2) {
64 throw RuntimeError::formatted(PISM_ERROR_LOCATION,
65 "bed_deformation.lc.elastic_model"
66 " requires bed_deformation.lc.grid_size_factor > 1");
67 }
68 }
69
70 // grid parameters
71 m_Mx = Mx;
72 m_My = My;
73 m_dx = dx;
74 m_dy = dy;
75 m_Nx = Nx;
76 m_Ny = Ny;
77
78 m_load_density = config.get_number("constants.ice.density");
79 m_mantle_density = config.get_number("bed_deformation.mantle_density");
80 m_eta = config.get_number("bed_deformation.mantle_viscosity");
81 m_D = config.get_number("bed_deformation.lithosphere_flexural_rigidity");
82
83 m_standard_gravity = config.get_number("constants.standard_gravity");
84
85 // derive more parameters
86 m_Lx = 0.5 * (m_Nx - 1.0) * m_dx;
87 m_Ly = 0.5 * (m_Ny - 1.0) * m_dy;
88 m_i0_offset = (Nx - Mx) / 2;
89 m_j0_offset = (Ny - My) / 2;
90
91 // memory allocation
92 PetscErrorCode ierr = 0;
93
94 // total displacement
95 ierr = VecCreateSeq(PETSC_COMM_SELF, m_Mx * m_My, m_U.rawptr());;
96 PISM_CHK(ierr, "VecCreateSeq");
97
98 // elastic displacement
99 ierr = VecCreateSeq(PETSC_COMM_SELF, m_Mx * m_My, m_Ue.rawptr());;
100 PISM_CHK(ierr, "VecCreateSeq");
101
102 // viscous displacement
103 ierr = VecCreateSeq(PETSC_COMM_SELF, m_Nx * m_Ny, m_Uv.rawptr());
104 PISM_CHK(ierr, "VecCreateSeq");
105
106 // setup fftw stuff: FFTW builds "plans" based on observed performance
107 m_fftw_input = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * m_Nx * m_Ny);
108 m_fftw_output = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * m_Nx * m_Ny);
109 m_loadhat = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * m_Nx * m_Ny);
110 m_lrm_hat = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * m_Nx * m_Ny);
111
112 clear_fftw_array(m_fftw_input, m_Nx, m_Ny);
113 m_dft_forward = fftw_plan_dft_2d(m_Nx, m_Ny, m_fftw_input, m_fftw_output,
114 FFTW_FORWARD, FFTW_ESTIMATE);
115 m_dft_inverse = fftw_plan_dft_2d(m_Nx, m_Ny, m_fftw_input, m_fftw_output,
116 FFTW_BACKWARD, FFTW_ESTIMATE);
117
118 // Note: FFTW is weird. If a malloc() call fails it will just call
119 // abort() on you without giving you a chance to recover or tell the
120 // user what happened. This is why we don't check return values of
121 // fftw_malloc() and fftw_plan_dft_2d() calls here...
122 //
123 // (Constantine Khroulev, February 1, 2015)
124
125 precompute_coefficients();
126 }
127
128 LingleClarkSerial::~LingleClarkSerial() {
129 fftw_destroy_plan(m_dft_forward);
130 fftw_destroy_plan(m_dft_inverse);
131 fftw_free(m_fftw_input);
132 fftw_free(m_fftw_output);
133 fftw_free(m_loadhat);
134 fftw_free(m_lrm_hat);
135 }
136
137 /*!
138 * Return total displacement.
139 */
140 Vec LingleClarkSerial::total_displacement() const {
141 return m_U;
142 }
143
144 /*!
145 * Return viscous plate displacement.
146 */
147 Vec LingleClarkSerial::viscous_displacement() const {
148 return m_Uv;
149 }
150
151 /*!
152 * Return elastic plate displacement.
153 */
154 Vec LingleClarkSerial::elastic_displacement() const {
155 return m_Ue;
156 }
157
158 void LingleClarkSerial::compute_load_response_matrix(fftw_complex *output) {
159
160 FFTWArray LRM(output, m_Nx, m_Ny);
161
162 greens_elastic G;
163 ge_data ge_data {m_dx, m_dy, 0, 0, &G};
164
165 int Nx2 = m_Nx / 2;
166 int Ny2 = m_Ny / 2;
167
168 // Top half
169 for (int j = 0; j <= Ny2; ++j) {
170 // Top left quarter
171 for (int i = 0; i <= Nx2; ++i) {
172 ge_data.p = Nx2 - i;
173 ge_data.q = Ny2 - j;
174
175 LRM(i, j) = dblquad_cubature(ge_integrand,
176 -m_dx / 2, m_dx / 2,
177 -m_dy / 2, m_dy / 2,
178 1.0e-8, &ge_data);
179 }
180
181 // Top right quarter
182 //
183 // Note: Nx2 = m_Nx / 2 (using integer division!), so
184 //
185 // - If m_Nx is even then 2 * Nx2 == m_Nx. So i < m_Nx implies i < 2 * Nx2 and
186 // 2 * Nx2 - i > 0.
187 //
188 // - If m_Nx is odd then 2 * Nx2 == m_Nx - 1 or m_Nx == 2 * Nx2 + 1. So i < m_Nx
189 // implies i < 2 * Nx2 + 1, which is the same as 2 * Nx2 - i > -1 or
190 // 2 * Nx2 - i >= 0.
191 //
192 // Also, i == Nx2 + 1 gives 2 * Nx2 - i == Nx2 - 1
193 //
194 // So, in both cases (even and odd) 0 <= 2 * Nx2 - i <= Nx2 - 1.
195 //
196 // This means that LRM(2 * Nx2 - i, j) will not use indexes that are out of bounds
197 // *and* will only use values computed in the for loop above.
198 for (int i = Nx2 + 1; i < m_Nx; ++i) {
199 assert(2 * Nx2 - i >= 0);
200 LRM(i, j) = LRM(2 * Nx2 - i, j);
201 }
202 } // End of the loop over the top half
203
204 // Bottom half
205 //
206 // See the comment above the "top right quarter" loop.
207 for (int j = Ny2 + 1; j < m_Ny; ++j) {
208 for (int i = 0; i < m_Nx; ++i) {
209 assert(2 * Ny2 - j >= 0);
210 LRM(i, j) = LRM(i, 2 * Ny2 - j);
211 }
212 }
213 }
214
215 /**
216 * Pre-compute coefficients used by the model.
217 */
218 void LingleClarkSerial::precompute_coefficients() {
219
220 // Coefficients for Fourier spectral method Laplacian
221 // MATLAB version: cx=(pi/Lx)*[0:Nx/2 Nx/2-1:-1:1]
222 m_cx = fftfreq(m_Nx, m_Lx / (m_Nx * M_PI));
223 m_cy = fftfreq(m_Ny, m_Ly / (m_Ny * M_PI));
224
225 // compare geforconv.m
226 if (m_include_elastic) {
227 m_log->message(2, " computing spherical elastic load response matrix ...");
228 {
229 compute_load_response_matrix(m_fftw_input);
230 // Compute fft2(LRM) and save it in m_lrm_hat
231 fftw_execute(m_dft_forward);
232 copy_fftw_array(m_fftw_output, m_lrm_hat, m_Nx, m_Ny);
233 }
234 m_log->message(2, " done\n");
235 }
236 }
237
238 /*!
239 * Solve
240 *
241 * @f$ 2 \nu |\nabla| \diff{u}{t} + \rho_r g U + D\nabla^4 U = \sigma_{zz}@f$
242 *
243 * for @f$ U @f$, treating @f$ \diff{u}{t} @f$ and @f$ \sigma_{zz} @f$ as known.
244 *
245 * @param[in] load_thickness load thickness, meters
246 * @param[in] bed_uplift bed uplift, m/second
247 *
248 * Here `load_thickness` is used to compute the load @f$ \sigma_{zz} @f$ and `bed_uplift` is
249 * @f$ \diff{u}{t} @f$ itself.
250 *
251 */
252 void LingleClarkSerial::uplift_problem(Vec load_thickness, Vec bed_uplift,
253 Vec output) {
254
255 // Compute fft2(-load_density * g * load_thickness)
256 {
257 clear_fftw_array(m_fftw_input, m_Nx, m_Ny);
258 set_real_part(load_thickness, - m_load_density * m_standard_gravity,
259 m_Mx, m_My, m_Nx, m_Ny, m_i0_offset, m_j0_offset,
260 m_fftw_input);
261 fftw_execute(m_dft_forward);
262 // Save fft2(-load_density * g * load_thickness) in loadhat.
263 copy_fftw_array(m_fftw_output, m_loadhat, m_Nx, m_Ny);
264 }
265
266 // fft2(uplift)
267 {
268 clear_fftw_array(m_fftw_input, m_Nx, m_Ny);
269 set_real_part(bed_uplift, 1.0, m_Mx, m_My, m_Nx, m_Ny, m_i0_offset, m_j0_offset,
270 m_fftw_input);
271 fftw_execute(m_dft_forward);
272 }
273
274 {
275 FFTWArray
276 u0_hat(m_fftw_input, m_Nx, m_Ny),
277 load_hat(m_loadhat, m_Nx, m_Ny),
278 uplift_hat(m_fftw_output, m_Nx, m_Ny);
279
280 for (int i = 0; i < m_Nx; i++) {
281 for (int j = 0; j < m_Ny; j++) {
282 const double
283 C = m_cx[i]*m_cx[i] + m_cy[j]*m_cy[j],
284 A = - 2.0 * m_eta * sqrt(C),
285 B = m_mantle_density * m_standard_gravity + m_D * C * C;
286
287 u0_hat(i, j) = (load_hat(i, j) + A * uplift_hat(i, j)) / B;
288 }
289 }
290 }
291
292 fftw_execute(m_dft_inverse);
293 get_real_part(m_fftw_output, 1.0 / (m_Nx * m_Ny), m_Nx, m_Ny, m_Nx, m_Ny, 0, 0, output);
294
295 tweak(load_thickness, output, m_Nx, m_Ny, 0.0);
296 }
297
298 /*! Initialize using provided load thickness and the bed uplift rate.
299 *
300 * Here we solve:
301 *
302 * rho_r g U + D grad^4 U = -rho g H - 2 eta |grad| uplift
303 *
304 * Compare equation (16) in Bueler, Lingle, Brown (2007) "Fast computation of a viscoelastic
305 * deformable Earth model for ice sheet simulations", Ann. Glaciol. 46, 97--105
306 *
307 * @note Probable sign error in eqn (16)?: load "rho g H" should be "- rho g H"]
308 *
309 * This initialization method is used to "bootstrap" the model. It should not be used to re-start a
310 * stopped modeling run.
311 *
312 * @param[in] thickness load thickness, meters
313 * @param[in] uplift initial bed uplift on the PISM grid
314 *
315 * Sets m_Uv, m_Ue, m_U.
316 */
317 void LingleClarkSerial::bootstrap(Vec thickness, Vec uplift) {
318
319 // compute viscous displacement
320 uplift_problem(thickness, uplift, m_Uv);
321
322 if (m_include_elastic) {
323 compute_elastic_response(thickness, m_Ue);
324 } else {
325 PetscErrorCode ierr = VecSet(m_Ue, 0.0); PISM_CHK(ierr, "VecSet");
326 }
327
328 update_displacement(m_Uv, m_Ue, m_U);
329 }
330
331 /*!
332 * Initialize using provided plate displacement.
333 *
334 * @param[in] viscous_displacement initial viscous plate displacement (meters) on the extended grid
335 * @param[in] elastic_displacement initial viscous plate displacement (meters) on the regular grid
336 *
337 * Sets m_Uv, m_Ue, m_U.
338 */
339 void LingleClarkSerial::init(Vec viscous_displacement,
340 Vec elastic_displacement) {
341 PetscErrorCode ierr = 0;
342
343 ierr = VecCopy(viscous_displacement, m_Uv); PISM_CHK(ierr, "VecCopy");
344
345 if (m_include_elastic) {
346 ierr = VecCopy(elastic_displacement, m_Ue); PISM_CHK(ierr, "VecCopy");
347 } else {
348 ierr = VecSet(m_Ue, 0.0); PISM_CHK(ierr, "VecSet");
349 }
350
351 update_displacement(m_Uv, m_Ue, m_U);
352 }
353
354 /*!
355 * Perform a time step.
356 *
357 * @param[in] dt time step length
358 * @param[in] H load thickness on the physical (Mx*My) grid
359 */
360 void LingleClarkSerial::step(double dt, Vec H) {
361 // solves:
362 // (2 eta |grad| U^{n+1}) + (dt/2) * (rho_r g U^{n+1} + D grad^4 U^{n+1})
363 // = (2 eta |grad| U^n) - (dt/2) * (rho_r g U^n + D grad^4 U^n) - dt * rho g H_start
364 // where U=plate displacement; see equation (7) in
365 // Bueler, Lingle, Brown (2007) "Fast computation of a viscoelastic
366 // deformable Earth model for ice sheet simulations", Ann. Glaciol. 46, 97--105
367
368 // Compute viscous displacement if dt > 0 and bypass this computation if dt == 0.
369 //
370 // This makes it easier to test the elastic part of the model.
371 if (dt > 0.0) {
372 // Non-zero time step: include the viscous part of the model.
373
374 // Compute fft2(-load_density * g * dt * H)
375 {
376 clear_fftw_array(m_fftw_input, m_Nx, m_Ny);
377 set_real_part(H,
378 - m_load_density * m_standard_gravity * dt,
379 m_Mx, m_My, m_Nx, m_Ny, m_i0_offset, m_j0_offset,
380 m_fftw_input);
381 fftw_execute(m_dft_forward);
382
383 // Save fft2(-load_density * g * H * dt) in loadhat.
384 copy_fftw_array(m_fftw_output, m_loadhat, m_Nx, m_Ny);
385 }
386
387 // Compute fft2(u).
388 // no need to clear fftw_input: all values are overwritten
389 {
390 set_real_part(m_Uv, 1.0, m_Nx, m_Ny, m_Nx, m_Ny, 0, 0, m_fftw_input);
391 fftw_execute(m_dft_forward);
392 }
393
394 // frhs = right.*fft2(uun) + fft2(dt*sszz);
395 // uun1 = real(ifft2(frhs./left));
396 {
397 FFTWArray input(m_fftw_input, m_Nx, m_Ny),
398 u_hat(m_fftw_output, m_Nx, m_Ny), load_hat(m_loadhat, m_Nx, m_Ny);
399 for (int i = 0; i < m_Nx; i++) {
400 for (int j = 0; j < m_Ny; j++) {
401 const double
402 C = m_cx[i]*m_cx[i] + m_cy[j]*m_cy[j],
403 part1 = 2.0 * m_eta * sqrt(C),
404 part2 = (dt / 2.0) * (m_mantle_density * m_standard_gravity + m_D * C * C),
405 A = part1 - part2,
406 B = part1 + part2;
407
408 input(i, j) = (load_hat(i, j) + A * u_hat(i, j)) / B;
409 }
410 }
411 }
412
413 fftw_execute(m_dft_inverse);
414 get_real_part(m_fftw_output, 1.0 / (m_Nx * m_Ny), m_Nx, m_Ny, m_Nx, m_Ny, 0, 0, m_Uv);
415
416 // Now tweak. (See the "correction" in section 5 of BuelerLingleBrown.)
417 //
418 // Here 1e16 approximates t = \infty.
419 tweak(H, m_Uv, m_Nx, m_Ny, 1e16);
420 } else {
421 // zero time step: viscous displacement is zero
422 PetscErrorCode ierr = VecSet(m_Uv, 0.0); PISM_CHK(ierr, "VecSet");
423 }
424
425 // now compute elastic response if desired
426 if (m_include_elastic) {
427 compute_elastic_response(H, m_Ue);
428 }
429
430 update_displacement(m_Uv, m_Ue, m_U);
431 }
432
433 /*!
434 * Compute elastic response to the load H
435 *
436 * @param[in] H load thickness (ice equivalent meters)
437 * @param[out] dE elastic plate displacement
438 */
439 void LingleClarkSerial::compute_elastic_response(Vec H, Vec dE) {
440
441 // Compute fft2(load_density * H)
442 //
443 // Note that here the load is placed in the corner of the array on the extended grid
444 // (offsets i0 and j0 are zero).
445 {
446 clear_fftw_array(m_fftw_input, m_Nx, m_Ny);
447 set_real_part(H, m_load_density, m_Mx, m_My, m_Nx, m_Ny, 0, 0, m_fftw_input);
448 fftw_execute(m_dft_forward);
449 }
450
451 // fft2(m_response_matrix) * fft2(load_density*H)
452 //
453 // Compute the product of Fourier transforms of the LRM and the load. This uses C++'s
454 // native support for complex arithmetic.
455 {
456 FFTWArray
457 input(m_fftw_input, m_Nx, m_Ny),
458 LRM_hat(m_lrm_hat, m_Nx, m_Ny),
459 load_hat(m_fftw_output, m_Nx, m_Ny);
460 for (int i = 0; i < m_Nx; i++) {
461 for (int j = 0; j < m_Ny; j++) {
462 input(i, j) = LRM_hat(i, j) * load_hat(i, j);
463 }
464 }
465 }
466
467 // Compute the inverse transform and extract the elastic response.
468 //
469 // Here the offsets are:
470 // i0 = m_Nx / 2,
471 // j0 = m_Ny / 2.
472 fftw_execute(m_dft_inverse);
473 get_real_part(m_fftw_output, 1.0 / (m_Nx * m_Ny), m_Mx, m_My, m_Nx, m_Ny,
474 m_Nx/2, m_Ny/2, dE);
475 }
476
477 /*! Compute total displacement by combining viscous and elastic contributions.
478 *
479 * @param[in] V viscous displacement
480 * @param[in] dE elastic displacement
481 * @param[out] dU total displacement
482 */
483 void LingleClarkSerial::update_displacement(Vec Uv, Vec Ue, Vec U) {
484 // PISM grid
485 petsc::VecArray2D
486 u(U, m_Mx, m_My),
487 u_elastic(Ue, m_Mx, m_My);
488 // extended grid
489 petsc::VecArray2D
490 u_viscous(Uv, m_Nx, m_Ny, m_i0_offset, m_j0_offset);
491
492 for (int i = 0; i < m_Mx; i++) {
493 for (int j = 0; j < m_My; j++) {
494 u(i, j) = u_viscous(i, j) + u_elastic(i, j);
495 }
496 }
497 }
498
499
500 /*!
501 * Modify the plate displacement to correct for the effect of imposing periodic boundary conditions
502 * at a finite distance.
503 *
504 * See Section 5 in [@ref BuelerLingleBrown].
505 *
506 * @param[in] load_thickness thickness of the load (used to compute the corresponding disc volume)
507 * @param[in,out] U viscous plate displacement
508 * @param[in] Nx grid size
509 * @param[in] Ny grid size
510 * @param[in] time time, seconds (usually 0 or a large number approximating \infty)
511 */
512 void LingleClarkSerial::tweak(Vec load_thickness, Vec U, int Nx, int Ny, double time) {
513 PetscErrorCode ierr = 0;
514 petsc::VecArray2D u(U, Nx, Ny);
515
516 // find average value along "distant" boundary of [-Lx, Lx]X[-Ly, Ly]
517 // note domain is periodic, so think of cut locus of torus (!)
518 // (will remove it: uun1=uun1-(sum(uun1(1, :))+sum(uun1(:, 1)))/(2*N);)
519 double average = 0.0;
520 for (int i = 0; i < Nx; i++) {
521 average += u(i, 0);
522 }
523
524 for (int j = 0; j < Ny; j++) {
525 average += u(0, j);
526 }
527
528 average /= (double) (Nx + Ny);
529
530 double shift = 0.0;
531
532 if (time > 0.0) {
533 // tweak continued: replace far field with value for an equivalent disc load which has
534 // R0=Lx*(2/3)=L/3 (instead of 1000km in MATLAB code: H0 = dx*dx*sum(sum(H))/(pi*1e6^2); %
535 // trapezoid rule)
536 const double L_average = (m_Lx + m_Ly) / 2.0;
537 const double R = L_average * (2.0 / 3.0);
538
539 double H_sum = 0.0;
540 ierr = VecSum(load_thickness, &H_sum); PISM_CHK(ierr, "VecSum");
541
542 // compute disc thickness by dividing its volume by the area
543 const double H = (H_sum * m_dx * m_dy) / (M_PI * R * R);
544
545 shift = viscDisc(time, // time in seconds
546 H, // disc thickness
547 R, // disc radius
548 L_average, // compute deflection at this radius
549 m_mantle_density, m_load_density, // mantle and load densities
550 m_standard_gravity, //
551 m_D, // flexural rigidity
552 m_eta); // mantle viscosity
553 }
554
555 ierr = VecShift(U, shift - average); PISM_CHK(ierr, "VecShift");
556 }
557
558 } // end of namespace bed
559 } // end of namespace pism