tbeddef_lc_elastic.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tbeddef_lc_elastic.py (4543B)
---
1 #!/usr/bin/env python3
2
3 from unittest import TestCase
4
5 import numpy as np
6 import scipy.integrate
7 import scipy.signal
8
9 import PISM
10
11 # silence models' initialization messages
12 PISM.Context().log.set_threshold(1)
13
14 class LingleClarkElastic(TestCase):
15 @staticmethod
16 def lrm(Mx, My, dx, dy):
17 "Compute the load response matrix, taking advantage of its symmetry."
18
19 ge = PISM.greens_elastic()
20
21 def ge_integrand(eta, xi, dx, dy, p, q):
22 xi_shift = p * dx - xi
23 eta_shift = q * dy - eta
24 r = np.sqrt(xi_shift * xi_shift + eta_shift * eta_shift)
25
26 return ge(r)
27
28 def f(dx, dy, p, q):
29 return scipy.integrate.dblquad(ge_integrand,
30 -dx/2.0, dx/2.0,
31 lambda x: -dy / 2.0, lambda x: dy / 2.0,
32 (dx, dy, p, q), epsrel=1e-8)[0]
33
34 Mx2 = int(Mx) // 2
35 My2 = int(My) // 2
36
37 a = np.zeros((My, Mx))
38
39 # top half
40 for j in range(My2 + 1):
41 # top left quarter
42 for i in range(Mx2 + 1):
43 p = Mx2 - i
44 q = My2 - j
45
46 a[j, i] = f(dx, dy, p, q)
47
48 # top right quarter
49 for i in range(Mx2 + 1, Mx):
50 a[j, i] = a[j, 2 * Mx2 - i]
51
52 # bottom half
53 for j in range(My2 + 1, My):
54 for i in range(Mx):
55 a[j, i] = a[2 * My2 - j, i]
56
57 return a
58
59 @staticmethod
60 def run_model(grid):
61 geometry = PISM.Geometry(grid)
62
63 bed_model = PISM.LingleClark(grid)
64
65 bed_uplift = PISM.IceModelVec2S(grid, "uplift", PISM.WITHOUT_GHOSTS)
66
67 # start with a flat bed, no ice, and no uplift
68 geometry.bed_elevation.set(0.0)
69 geometry.ice_thickness.set(0.0)
70 geometry.sea_level_elevation.set(-1000.0) # everything is grounded
71 geometry.ensure_consistency(0.0)
72
73 bed_uplift.set(0.0)
74
75 bed_model.bootstrap(geometry.bed_elevation, bed_uplift, geometry.ice_thickness,
76 geometry.sea_level_elevation)
77
78 Mx2 = int(grid.Mx()) // 2
79 My2 = int(grid.My()) // 2
80
81 # add the load
82 with PISM.vec.Access(nocomm=geometry.ice_thickness):
83 for (i, j) in grid.points():
84 # if i == Mx2 and j == My2:
85 if abs(i - Mx2) < 2 and abs(j - My2) < 2:
86 geometry.ice_thickness[i, j] = 1000.0
87
88 # dt of zero disables the viscous part of the model, so all we get is the elastic
89 # response
90 bed_model.step(geometry.ice_thickness, geometry.sea_level_elevation, 0)
91
92 return (geometry.ice_thickness.numpy(),
93 bed_model.total_displacement().numpy(),
94 bed_model.elastic_load_response_matrix().numpy())
95
96 def setUp(self):
97 self.ctx = PISM.Context()
98
99 self.elastic = self.ctx.config.get_flag("bed_deformation.lc.elastic_model")
100 self.ctx.config.set_flag("bed_deformation.lc.elastic_model", True)
101
102 self.size_factor = self.ctx.config.get_number("bed_deformation.lc.grid_size_factor")
103 self.ctx.config.set_number("bed_deformation.lc.grid_size_factor", 2)
104
105 Lx = 2000e3
106 Mx = 11
107 My = 2 * Mx # non-square grid
108
109 self.grid = PISM.IceGrid.Shallow(self.ctx.ctx, Lx, Lx,
110 0, 0, Mx, My, PISM.CELL_CORNER, PISM.NOT_PERIODIC)
111
112 self.H, self.db_pism, self.lrm_pism = self.run_model(self.grid)
113
114 def convolution_test(self):
115 "Compare PISM's FFTW-based convolution to scipy.signal.fftconvolve()"
116 rho = self.ctx.config.get_number("constants.ice.density")
117
118 db_scipy = scipy.signal.fftconvolve(rho * self.H, self.lrm_pism, mode="same")
119
120 np.testing.assert_allclose(self.db_pism, db_scipy, rtol=1e-12)
121
122 def lrm_test(self):
123 "Compare PISM's load response matrix to the one computed using scipy.integrate.dblquad()"
124
125 dx = self.grid.dx()
126 dy = self.grid.dy()
127
128 Ny, Nx = self.lrm_pism.shape
129 lrm_python = self.lrm(Nx, Ny, dx, dy)
130
131 # This is a crappy relative tolerance. Oh well...
132 np.testing.assert_allclose(self.lrm_pism, lrm_python, rtol=1e-2)
133
134 def tearDown(self):
135 # reset configuration parameters
136 self.ctx.config.set_flag("bed_deformation.lc.elastic_model", self.elastic)
137 self.ctx.config.set_number("bed_deformation.lc.grid_size_factor", self.size_factor)