tbeddef_lc_viscous.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_viscous.py (8490B)
---
1 #!/usr/bin/env python3
2
3 """Solves equation (16) of BuelerLingleBrown to obtain the steady
4 state viscous plate deflection corresponding to a given disc load.
5
6 Used as a verification (and regression) test for LingleClarkSerial::bootstrap().
7 """
8
9 import PISM
10 import numpy as np
11 from PISM.util import convert
12
13 config = PISM.Context().config
14 log = PISM.Context().log
15
16 config.set_number("bed_deformation.lc.grid_size_factor", 2)
17 config.set_flag("bed_deformation.lc.elastic_model", False)
18
19 # constants
20 standard_gravity = config.get_number("constants.standard_gravity")
21 ice_density = config.get_number("constants.ice.density")
22 mantle_density = config.get_number("bed_deformation.mantle_density")
23 mantle_viscosity = config.get_number("bed_deformation.mantle_viscosity")
24 lithosphere_flexural_rigidity = config.get_number("bed_deformation.lithosphere_flexural_rigidity")
25
26 # disc load parameters
27 disc_radius = convert(1000, "km", "m")
28 disc_thickness = 1000.0 # meters
29 # domain size
30 Lx = 2 * disc_radius
31 # time to use for the comparison
32 T = convert(1e6, "years", "second")
33 t_final = convert(20000, "years", "second")
34 dt = convert(500, "years", "second")
35
36
37 def deflection(time, radius, disc_thickness, disc_radius):
38 """Compute the viscous plate deflection. See formula (17) in 'Fast
39 computation of a viscoelastic deformable Earth model for ice-sheet
40 simulations' by Bueler, Lingle, and Brown, 2007.
41 """
42 return PISM.viscDisc(time, disc_thickness, disc_radius, radius,
43 mantle_density, ice_density, standard_gravity,
44 lithosphere_flexural_rigidity, mantle_viscosity)
45
46
47 def exact(dics_radius, disc_thickness, t, L, N):
48 "Evaluate the exact solution at N points."
49 r = np.linspace(0, L, N)
50 z = [deflection(t, rr, disc_thickness, disc_radius) for rr in r]
51 return (r, z)
52
53
54 def modeled_time_dependent(dics_radius, disc_thickness, t_end, L, Nx, dt):
55 "Use the LingleClark class to compute plate deflection."
56
57 Ny = Nx
58 Mx = int(2 * Nx - 1)
59 My = int(2 * Ny - 1)
60
61 ctx = PISM.Context().ctx
62 grid = PISM.IceGrid.Shallow(ctx, L, L, 0, 0, Mx, My, PISM.CELL_CORNER, PISM.NOT_PERIODIC)
63
64 bed_model = PISM.LingleClark(grid)
65
66 ice_thickness = PISM.IceModelVec2S(grid, "thk", PISM.WITHOUT_GHOSTS)
67
68 bed = PISM.IceModelVec2S(grid, "topg", PISM.WITHOUT_GHOSTS)
69
70 bed_uplift = PISM.IceModelVec2S(grid, "uplift", PISM.WITHOUT_GHOSTS)
71
72 sea_level = PISM.IceModelVec2S(grid, "sea_level", PISM.WITHOUT_GHOSTS)
73
74 # start with a flat bed, no ice, and no uplift
75 bed.set(0.0)
76 bed_uplift.set(0.0)
77 ice_thickness.set(0.0)
78 sea_level.set(-1000.0)
79
80 bed_model.bootstrap(bed, bed_uplift, ice_thickness, sea_level)
81
82 # add the disc load
83 with PISM.vec.Access(nocomm=ice_thickness):
84 for (i, j) in grid.points():
85 r = PISM.radius(grid, i, j)
86 if r <= disc_radius:
87 ice_thickness[i, j] = disc_thickness
88
89 t = 0.0
90 while t < t_end:
91 # make sure we don't step past t_end
92 if t + dt > t_end:
93 dt = t_end - t
94
95 bed_model.step(ice_thickness, sea_level, dt)
96
97 t += dt
98 log.message(2, ".")
99 log.message(2, "\n")
100
101 # extract half of the x grid
102 r = grid.x()[Nx-1:]
103
104 # extract values along the x direction (along the radius of the disc)
105 z = bed_model.bed_elevation().numpy()[Ny-1, Nx-1:]
106
107 return r, z
108
109
110 def modeled_steady_state(dics_radius, disc_thickness, time, L, Nx):
111 "Use the LingleClark class to compute plate deflection."
112
113 Ny = Nx
114 Mx = int(2 * Nx - 1)
115 My = int(2 * Ny - 1)
116
117 ctx = PISM.Context().ctx
118 grid = PISM.IceGrid.Shallow(ctx, L, L, 0, 0, Mx, My, PISM.CELL_CORNER, PISM.NOT_PERIODIC)
119
120 bed_model = PISM.LingleClark(grid)
121
122 ice_thickness = PISM.IceModelVec2S(grid, "thk", PISM.WITHOUT_GHOSTS)
123
124 bed = PISM.IceModelVec2S(grid, "topg", PISM.WITHOUT_GHOSTS)
125 bed.set(0.0)
126
127 bed_uplift = PISM.IceModelVec2S(grid, "uplift", PISM.WITHOUT_GHOSTS)
128 bed_uplift.set(0.0)
129
130 sea_level = PISM.IceModelVec2S(grid, "sea_level", PISM.WITHOUT_GHOSTS)
131 sea_level.set(-1000.0)
132
133 with PISM.vec.Access(nocomm=ice_thickness):
134 for (i, j) in grid.points():
135 r = PISM.radius(grid, i, j)
136 if r <= disc_radius:
137 ice_thickness[i, j] = disc_thickness
138 else:
139 ice_thickness[i, j] = 0.0
140
141 bed_model.bootstrap(bed, bed_uplift, ice_thickness, sea_level)
142
143 # extract half of the x grid
144 r = grid.x()[Nx-1:]
145
146 # extract values along the x direction (along the radius of the disc)
147 z = bed_model.total_displacement().numpy()[Ny-1, Nx-1:]
148
149 return r, z
150
151
152 def compare_steady_state(N):
153 "Compare eact and modeled results."
154 r_exact, z_exact = exact(disc_radius, disc_thickness, T, Lx, N)
155 r, z = modeled_steady_state(disc_radius, disc_thickness, T, Lx, N)
156
157 diff_origin = np.fabs(z_exact[0] - z[0])
158 diff_average = np.average(np.fabs(z_exact - z))
159 diff_max = np.max(np.fabs(z_exact - z))
160
161 return diff_origin, diff_max, diff_average
162
163
164 def compare_time_dependent(N):
165 "Compare exact and modeled results."
166 r_exact, z_exact = exact(disc_radius, disc_thickness, t_final, Lx, N)
167
168 dx = r_exact[1] - r_exact[0]
169 log.message(2, "N = {}, dx = {} km\n".format(N, dx/1000.0))
170
171 r, z = modeled_time_dependent(disc_radius, disc_thickness, t_final, Lx, N, dt)
172
173 diff_origin = np.fabs(z_exact[0] - z[0])
174 diff_average = np.average(np.fabs(z_exact - z))
175 diff_max = np.max(np.fabs(z_exact - z))
176
177 return diff_origin, diff_max, diff_average, dx
178
179
180 def time_dependent_test():
181 "Time dependent bed deformation (disc load)"
182 diff = np.array([compare_time_dependent(n)[:3] for n in [34, 67]])
183
184 stored = [[0.04099917, 5.05854, 0.93909436],
185 [0.05710513, 4.14329508, 0.71246272]]
186
187 return np.testing.assert_almost_equal(diff, stored)
188
189
190 def steady_state_test():
191 "Steady state bed deformation (disc load)"
192 Ns = 10 * np.arange(1, 5) + 1
193 diff = np.array([compare_steady_state(n) for n in Ns])
194
195 stored = [[ 0.0399697, 15.71882867, 3.80458833],
196 [ 0.04592036, 11.43876195, 1.94967725],
197 [ 0.04357962, 9.7207298, 1.76262896],
198 [ 0.04019595, 7.71929661, 1.38746767]]
199
200 return np.testing.assert_almost_equal(diff, stored)
201
202
203 def verify_steady_state():
204 "Set up a grid refinement study and produce convergence plots."
205
206 Ns = 101 + 10 * np.arange(0, 10)
207
208 diff = np.array([compare_steady_state(n) for n in Ns])
209
210 plt.figure()
211 plt.title("Steady state")
212
213 d = np.log10(diff)
214 log_n = np.log10(1.0 / Ns)
215 for j, label in enumerate(["origin", "max", "average"]):
216 p = np.polyfit(log_n, d[:, j], 1)
217 plt.plot(log_n, d[:, j], marker='o', label="{}, O(dx^{:.2})".format(label, p[0]))
218 plt.plot(log_n, np.polyval(p, log_n), ls="--")
219
220 plt.legend()
221 plt.grid(True)
222 plt.xlabel("log10(1/N)")
223 plt.ylabel("log10(error)")
224 plt.title("Convergence rates for the steady-state problem")
225 plt.show()
226
227
228 def verify_time_dependent():
229 "Set up a spatial grid refinement study and produce convergence plots."
230
231 dxs = [15, 30, 60, 125, 250, 500]
232 # in km, same as in BuelerLingleBrown, figure 4
233 #
234 # Note that we compute max and average errors differently here, so
235 # convergence rates produced by this scrips are not directly
236 # comparable to the paper.
237
238 Ns = [int(Lx / (1000 * dx)) + 1 for dx in dxs]
239
240 diff = np.array([compare_time_dependent(n) for n in Ns])
241
242 plt.figure()
243 plt.title("Time-dependent")
244
245 d = np.log10(diff)
246 dx = diff[:, 3] / 1000.0 # convert to km
247 log_dx = np.log10(dx)
248 for j, label in enumerate(["origin", "max", "average"]):
249 p = np.polyfit(log_dx, d[:, j], 1)
250 plt.plot(log_dx, d[:, j], marker='o', label="{}, O(dx^{:.2})".format(label, p[0]))
251 plt.plot(log_dx, np.polyval(p, log_dx), ls="--")
252 plt.xticks(log_dx, ["{:.0f}".format(x) for x in dx])
253
254 plt.legend()
255 plt.grid(True)
256 plt.xlabel("dx, km")
257 plt.ylabel("log10(error)")
258 plt.title("Convergence rates for the time-dependent problem")
259 plt.show()
260
261 if __name__ == "__main__":
262 import pylab as plt
263 log.message(2, " Creating convergence plots (spatial refinement)...\n")
264 log.message(2, " 1. Steady state problem...\n")
265 verify_steady_state()
266 log.message(2, " 2. Time-dependent problem...\n")
267 verify_time_dependent()