Loading [MathJax]/extensions/tex2jax.js
PISM, A Parallel Ice Sheet Model 2.2.2-d6b3a29ca committed by Constantine Khrulev on 2025-03-28
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LingleClarkSerial.cc
Go to the documentation of this file.
1// Copyright (C) 2004-2009, 2011, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2023 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 "pism/earth/matlablike.hh"
25#include "pism/earth/greens.hh"
26#include "pism/earth/LingleClarkSerial.hh"
27
28#include "pism/util/ConfigInterface.hh"
29#include "pism/util/error_handling.hh"
30#include "pism/util/petscwrappers/Vec.hh"
31#include "pism/util/fftw_utilities.hh"
32
33namespace pism {
34namespace bed {
35
36/*!
37 * @param[in] config configuration database
38 * @param[in] include_elastic include elastic deformation component
39 * @param[in] Mx grid size in the X direction
40 * @param[in] My grid size in the Y direction
41 * @param[in] dx grid spacing in the X direction
42 * @param[in] dy grid spacing in the Y direction
43 * @param[in] Nx extended grid size in the X direction
44 * @param[in] Ny extended grid size in the Y direction
45 */
47 const Config &config,
48 bool include_elastic,
49 int Mx, int My,
50 double dx, double dy,
51 int Nx, int Ny)
52 : m_t_infty(1e16), // around 317 million years
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) {
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
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
126}
127
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 */
141 return m_U;
142}
143
144/*!
145 * Return viscous plate displacement.
146 */
150
151/*!
152 * Return elastic plate displacement.
153 */
157
159
160 FFTWArray LRM(output, m_Nx, m_Ny);
161
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 */
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 {
230 // Compute fft2(LRM) and save it in m_lrm_hat
231 fftw_execute(m_dft_forward);
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 */
253 petsc::Vec &bed_uplift,
254 petsc::Vec &output) {
255
256 // Compute fft2(-load_density * g * load_thickness)
257 {
262 fftw_execute(m_dft_forward);
263 // Save fft2(-load_density * g * load_thickness) in loadhat.
265 }
266
267 // fft2(uplift)
268 {
270 set_real_part(bed_uplift, 1.0, m_Mx, m_My, m_Nx, m_Ny, m_i0_offset, m_j0_offset,
272 fftw_execute(m_dft_forward);
273 }
274
275 {
277 u0_hat(m_fftw_input, m_Nx, m_Ny),
278 load_hat(m_loadhat, m_Nx, m_Ny),
279 uplift_hat(m_fftw_output, m_Nx, m_Ny);
280
281 for (int i = 0; i < m_Nx; i++) {
282 for (int j = 0; j < m_Ny; j++) {
283 const double
284 C = m_cx[i]*m_cx[i] + m_cy[j]*m_cy[j],
285 A = - 2.0 * m_eta * sqrt(C),
287
288 u0_hat(i, j) = (load_hat(i, j) + A * uplift_hat(i, j)) / B;
289 }
290 }
291 }
292
293 fftw_execute(m_dft_inverse);
294 get_real_part(m_fftw_output, 1.0 / (m_Nx * m_Ny), m_Nx, m_Ny, m_Nx, m_Ny, 0, 0, output);
295
296 tweak(load_thickness, output, m_Nx, m_Ny);
297}
298
299/*! Initialize using provided load thickness and the bed uplift rate.
300 *
301 * Here we solve:
302 *
303 * rho_r g U + D grad^4 U = -rho g H - 2 eta |grad| uplift
304 *
305 * Compare equation (16) in Bueler, Lingle, Brown (2007) "Fast computation of a viscoelastic
306 * deformable Earth model for ice sheet simulations", Ann. Glaciol. 46, 97--105
307 *
308 * @note Probable sign error in eqn (16)?: load "rho g H" should be "- rho g H"]
309 *
310 * This initialization method is used to "bootstrap" the model. It should not be used to re-start a
311 * stopped modeling run.
312 *
313 * @param[in] thickness load thickness, meters
314 * @param[in] uplift initial bed uplift on the PISM grid
315 *
316 * Sets m_Uv, m_Ue, m_U.
317 */
319
320 // compute viscous displacement
321 uplift_problem(thickness, uplift, m_Uv);
322
323 if (m_include_elastic) {
324 compute_elastic_response(thickness, m_Ue);
325 } else {
326 PetscErrorCode ierr = VecSet(m_Ue, 0.0); PISM_CHK(ierr, "VecSet");
327 }
328
330}
331
332/*!
333 * Initialize using provided plate displacement.
334 *
335 * @param[in] viscous_displacement initial viscous plate displacement (meters) on the extended grid
336 * @param[in] elastic_displacement initial viscous plate displacement (meters) on the regular grid
337 *
338 * Sets m_Uv, m_Ue, m_U.
339 */
340void LingleClarkSerial::init(petsc::Vec &viscous_displacement,
341 petsc::Vec &elastic_displacement) {
342 PetscErrorCode ierr = 0;
343
344 ierr = VecCopy(viscous_displacement, m_Uv); PISM_CHK(ierr, "VecCopy");
345
346 if (m_include_elastic) {
347 ierr = VecCopy(elastic_displacement, m_Ue); PISM_CHK(ierr, "VecCopy");
348 } else {
349 ierr = VecSet(m_Ue, 0.0); PISM_CHK(ierr, "VecSet");
350 }
351
353}
354
355/*!
356 * Perform a time step.
357 *
358 * @param[in] dt time step length
359 * @param[in] H load thickness on the physical (Mx*My) grid
360 */
362 // solves:
363 // (2 eta |grad| U^{n+1}) + (dt/2) * (rho_r g U^{n+1} + D grad^4 U^{n+1})
364 // = (2 eta |grad| U^n) - (dt/2) * (rho_r g U^n + D grad^4 U^n) - dt * rho g H_start
365 // where U=plate displacement; see equation (7) in
366 // Bueler, Lingle, Brown (2007) "Fast computation of a viscoelastic
367 // deformable Earth model for ice sheet simulations", Ann. Glaciol. 46, 97--105
368
369 // Compute viscous displacement if dt > 0 and bypass this computation if dt == 0.
370 //
371 // This makes it easier to test the elastic part of the model.
372 if (dt > 0.0) {
373 // Non-zero time step: include the viscous part of the model.
374
375 // Compute fft2(-load_density * g * dt * H)
376 {
382 fftw_execute(m_dft_forward);
383
384 // Save fft2(-load_density * g * H * dt) in loadhat.
386 }
387
388 // Compute fft2(u).
389 // no need to clear fftw_input: all values are overwritten
390 {
392 fftw_execute(m_dft_forward);
393 }
394
395 // frhs = right.*fft2(uun) + fft2(dt*sszz);
396 // uun1 = real(ifft2(frhs./left));
397 {
399 u_hat(m_fftw_output, m_Nx, m_Ny), load_hat(m_loadhat, m_Nx, m_Ny);
400 for (int i = 0; i < m_Nx; i++) {
401 for (int j = 0; j < m_Ny; j++) {
402 const double
403 C = m_cx[i]*m_cx[i] + m_cy[j]*m_cy[j],
404 part1 = 2.0 * m_eta * sqrt(C),
405 part2 = (dt / 2.0) * (m_mantle_density * m_standard_gravity + m_D * C * C),
406 A = part1 - part2,
407 B = part1 + part2;
408
409 input(i, j) = (load_hat(i, j) + A * u_hat(i, j)) / B;
410 }
411 }
412 }
413
414 fftw_execute(m_dft_inverse);
415 get_real_part(m_fftw_output, 1.0 / (m_Nx * m_Ny), m_Nx, m_Ny, m_Nx, m_Ny, 0, 0, m_Uv);
416
417 // Now tweak. (See the "correction" in section 5 of BuelerLingleBrown.)
418 tweak(H, m_Uv, m_Nx, m_Ny);
419 } else {
420 // zero time step: viscous displacement is zero
421 PetscErrorCode ierr = VecSet(m_Uv, 0.0); PISM_CHK(ierr, "VecSet");
422 }
423
424 // now compute elastic response if desired
425 if (m_include_elastic) {
427 }
428
430}
431
432/*!
433 * Compute elastic response to the load H
434 *
435 * @param[in] H load thickness (ice equivalent meters)
436 * @param[out] dE elastic plate displacement
437 */
439
440 // Compute fft2(load_density * H)
441 //
442 // Note that here the load is placed in the corner of the array on the extended grid
443 // (offsets i0 and j0 are zero).
444 {
447 fftw_execute(m_dft_forward);
448 }
449
450 // fft2(m_response_matrix) * fft2(load_density*H)
451 //
452 // Compute the product of Fourier transforms of the LRM and the load. This uses C++'s
453 // native support for complex arithmetic.
454 {
456 input(m_fftw_input, m_Nx, m_Ny),
457 LRM_hat(m_lrm_hat, m_Nx, m_Ny),
458 load_hat(m_fftw_output, m_Nx, m_Ny);
459 for (int i = 0; i < m_Nx; i++) {
460 for (int j = 0; j < m_Ny; j++) {
461 input(i, j) = LRM_hat(i, j) * load_hat(i, j);
462 }
463 }
464 }
465
466 // Compute the inverse transform and extract the elastic response.
467 //
468 // Here the offsets are:
469 // i0 = m_Nx / 2,
470 // j0 = m_Ny / 2.
471 fftw_execute(m_dft_inverse);
473 m_Nx/2, m_Ny/2, dE);
474}
475
476/*! Compute total displacement by combining viscous and elastic contributions.
477 *
478 * @param[in] V viscous displacement
479 * @param[in] dE elastic displacement
480 * @param[out] dU total displacement
481 */
483 // PISM grid
485 u(U, m_Mx, m_My),
486 u_elastic(Ue, m_Mx, m_My);
487 // extended grid
489 u_viscous(Uv, m_Nx, m_Ny, m_i0_offset, m_j0_offset);
490
491 for (int i = 0; i < m_Mx; i++) {
492 for (int j = 0; j < m_My; j++) {
493 u(i, j) = u_viscous(i, j) + u_elastic(i, j);
494 }
495 }
496}
497
498
499/*!
500 * Modify the plate displacement to correct for the effect of imposing periodic boundary conditions
501 * at a finite distance.
502 *
503 * See Section 5 in [@ref BuelerLingleBrown].
504 *
505 * @param[in] load_thickness thickness of the load (used to compute the corresponding disc volume)
506 * @param[in,out] U viscous plate displacement
507 * @param[in] Nx grid size
508 * @param[in] Ny grid size
509 */
510void LingleClarkSerial::tweak(petsc::Vec &load_thickness, petsc::Vec &U, int Nx, int Ny) {
511 PetscErrorCode ierr = 0;
512 petsc::VecArray2D u(U, Nx, Ny);
513
514 // find average value along "distant" boundary of [-Lx, Lx]X[-Ly, Ly]
515 // note domain is periodic, so think of cut locus of torus (!)
516 // (will remove it: uun1=uun1-(sum(uun1(1, :))+sum(uun1(:, 1)))/(2*N);)
517 double average = 0.0;
518 for (int i = 0; i < Nx; i++) {
519 average += u(i, 0);
520 }
521
522 for (int j = 0; j < Ny; j++) {
523 average += u(0, j);
524 }
525
526 average /= (double) (Nx + Ny);
527
528 double shift = 0.0;
529
530 {
531 // tweak continued: replace far field with value for an equivalent disc load which has
532 // R0=Lx*(2/3)=L/3 (instead of 1000km in MATLAB code: H0 = dx*dx*sum(sum(H))/(pi*1e6^2); %
533 // trapezoid rule)
534 const double L_average = (m_Lx + m_Ly) / 2.0;
535 const double R = L_average * (2.0 / 3.0);
536
537 double H_sum = 0.0;
538 ierr = VecSum(load_thickness, &H_sum); PISM_CHK(ierr, "VecSum");
539
540 // compute disc thickness by dividing its volume by the area
541 const double H = (H_sum * m_dx * m_dy) / (M_PI * R * R);
542
543 shift = viscDisc(m_t_infty, // time in seconds
544 H, // disc thickness
545 R, // disc radius
546 L_average, // compute deflection at this radius
547 m_mantle_density, m_load_density, // mantle and load densities
549 m_D, // flexural rigidity
550 m_eta); // mantle viscosity
551 }
552
553 ierr = VecShift(U, shift - average); PISM_CHK(ierr, "VecShift");
554}
555
556} // end of namespace bed
557} // end of namespace pism
double get_number(const std::string &name, UseFlag flag=REMEMBER_THIS_USE) const
A class for storing and accessing PISM configuration flags and parameters.
std::shared_ptr< const Logger > ConstPtr
Definition Logger.hh:46
static RuntimeError formatted(const ErrorLocation &location, const char format[],...) __attribute__((format(printf
build a RuntimeError with a formatted message
T * rawptr()
Definition Wrapper.hh:39
void uplift_problem(petsc::Vec &load_thickness, petsc::Vec &bed_uplift, petsc::Vec &output)
void update_displacement(petsc::Vec &V, petsc::Vec &dE, petsc::Vec &dU)
LingleClarkSerial(Logger::ConstPtr log, const Config &config, bool include_elastic, int Mx, int My, double dx, double dy, int Nx, int Ny)
void step(double dt_seconds, petsc::Vec &H)
const petsc::Vec & viscous_displacement() const
double m_D
lithosphere flexural rigidity
void bootstrap(petsc::Vec &thickness, petsc::Vec &uplift)
void compute_elastic_response(petsc::Vec &H, petsc::Vec &dE)
double m_mantle_density
mantle density
void compute_load_response_matrix(fftw_complex *output)
void init(petsc::Vec &viscous_displacement, petsc::Vec &elastic_displacement)
const petsc::Vec & elastic_displacement() const
void tweak(petsc::Vec &load_thickness, petsc::Vec &U, int Nx, int Ny)
const petsc::Vec & total_displacement() const
double m_eta
mantle viscosity
double m_load_density
load density (for computing load from its thickness)
Wrapper around VecGetArray2d and VecRestoreArray2d.
Definition Vec.hh:55
#define PISM_CHK(errcode, name)
#define PISM_ERROR_LOCATION
const double G
Definition exactTestK.c:40
double dblquad_cubature(integrand f, double ax, double bx, double ay, double by, double reqRelError, void *fdata)
Definition matlablike.cc:22
double ge_integrand(unsigned ndim, const double *args, void *data)
The integrand in defining the elastic Green's function from the [Farrell] earth model.
Definition greens.cc:30
double viscDisc(double t, double H0, double R0, double r, double rho, double rho_ice, double grav, double D, double eta)
Actually compute the response of the viscous half-space model in LingleClark, to a disc load.
Definition greens.cc:135
void clear_fftw_array(fftw_complex *input, int Nx, int Ny)
Fill input with zeros.
void copy_fftw_array(fftw_complex *source, fftw_complex *destination, int Nx, int Ny)
Copy source to destination.
void get_real_part(fftw_complex *input, double normalization, int Mx, int My, int Nx, int Ny, int i0, int j0, petsc::Vec &output)
Get the real part of input and put it in output.
void set_real_part(petsc::Vec &input, double normalization, int Mx, int My, int Nx, int Ny, int i0, int j0, fftw_complex *output)
std::vector< double > fftfreq(int M, double normalization)
Parameters used to access elastic Green's function from the [Farrell] earth model.
Definition greens.hh:56