tmass_transport.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
---
tmass_transport.py (5661B)
---
1 import PISM
2 import numpy as np
3
4 """ Test mass transport code using a radially-symmetric setup in which
5 a disc expands uniformly in all directions. Given mass conservation,
6 the thickness of this disc outside of the fixed circular area (the
7 thickness Dirichlet B.C. area) is inversely proportional to the
8 distance from the center. Here we use this 'exact solution' to test
9 the symmetry of the produced ice thickness and linear convergence
10 towards the exact solution."""
11
12 log = PISM.Context().log
13
14
15 def disc(thickness, x0, y0, H0, R_inner, R_outer):
16 """Set ice thickness to H0 within the disc centered at (x0,y0) of
17 radius R_inner, C/R in an annulus R_inner < r <= R_outer and 0
18 elsewhere.
19
20 """
21
22 grid = thickness.grid()
23
24 R_inner_2 = R_inner**2
25 R_outer_2 = R_outer**2
26
27 C = H0 * R_inner
28
29 with PISM.vec.Access(nocomm=thickness):
30 for (i, j) in grid.points():
31 x = grid.x(i)
32 y = grid.y(j)
33 d2 = (x - x0)**2 + (y - y0)**2
34 if d2 <= R_inner_2:
35 thickness[i, j] = H0
36 elif d2 <= R_outer_2:
37 thickness[i, j] = C / np.sqrt(d2)
38 else:
39 thickness[i, j] = 0.0
40
41 thickness.update_ghosts()
42
43
44 def set_velocity(scalar_velocity, v):
45 """Initialize the velocity field to a rigid rotation around the
46 origin. This is slow, but it works.
47
48 """
49 grid = v.grid()
50
51 with PISM.vec.Access(nocomm=v):
52 for (i, j) in grid.points():
53 x = grid.x(i)
54 y = grid.y(j)
55
56 r = max(PISM.radius(grid, i, j), 0.001)
57
58 v[i, j].u = scalar_velocity * x / r
59 v[i, j].v = scalar_velocity * y / r
60
61 v.update_ghosts()
62
63
64 def run(Mx, My, t_final, part_grid, C=1.0):
65 "Test GeometryEvolution::step()"
66
67 ctx = PISM.Context().ctx
68
69 config = PISM.Context().config
70
71 config.set_flag("geometry.part_grid.enabled", part_grid)
72
73 grid = PISM.IceGrid_Shallow(ctx, 1, 1, 0, 0, Mx, My, PISM.CELL_CORNER, PISM.NOT_PERIODIC)
74
75 assert t_final <= 1.0
76
77 L = min(grid.Lx(), grid.Ly())
78 R_inner = 0.25 * L
79 spreading_velocity = 0.7
80 R_outer = R_inner + spreading_velocity * t_final
81
82 geometry = PISM.Geometry(grid)
83
84 v = PISM.IceModelVec2V(grid, "velocity", PISM.WITHOUT_GHOSTS)
85 Q = PISM.IceModelVec2Stag(grid, "Q", PISM.WITHOUT_GHOSTS)
86 v_bc_mask = PISM.IceModelVec2Int(grid, "v_bc_mask", PISM.WITHOUT_GHOSTS)
87 H_bc_mask = PISM.IceModelVec2Int(grid, "H_bc_mask", PISM.WITHOUT_GHOSTS)
88
89 ge = PISM.GeometryEvolution(grid)
90
91 # grid info
92 geometry.latitude.set(0.0)
93 geometry.longitude.set(0.0)
94 # environment
95 geometry.bed_elevation.set(-10.0)
96 geometry.sea_level_elevation.set(0.0)
97 # set initial ice thickness
98 disc(geometry.ice_thickness, 0, 0, 1, R_inner, R_inner)
99 geometry.ice_area_specific_volume.set(0.0)
100
101 geometry.ensure_consistency(0.0)
102
103 set_velocity(spreading_velocity, v)
104 v_bc_mask.set(0.0)
105 disc(H_bc_mask, 0, 0, 1, R_inner, R_inner)
106
107 profiling = ctx.profiling()
108 profiling.start()
109
110 t = 0.0
111 j = 0
112 profiling.stage_begin("ge")
113 while t < t_final:
114 cfl_data = PISM.max_timestep_cfl_2d(geometry.ice_thickness,
115 geometry.cell_type,
116 v)
117
118 dt = cfl_data.dt_max.value() * C
119
120 if t + dt > t_final:
121 dt = t_final - t
122
123 log.message(2, "{}, {}\n".format(t, dt))
124
125 profiling.begin("step")
126 ge.flow_step(geometry, dt,
127 v,
128 Q,
129 v_bc_mask,
130 H_bc_mask)
131 profiling.end("step")
132
133 profiling.begin("modify")
134 ge.apply_flux_divergence(geometry)
135 geometry.ensure_consistency(0.0)
136 profiling.end("modify")
137
138 t += dt
139 j += 1
140 profiling.stage_end("ge")
141
142 profiling.report("profiling_%d_%d.py" % (Mx, My))
143
144 return geometry
145
146
147 def average_error(N):
148 t_final = 1.0
149 C = 1.0
150
151 log.disable()
152 geometry = run(N, N, t_final, True, C)
153 log.enable()
154 # combine stuff stored as thickness and as area specific volume
155 geometry.ice_thickness.add(1.0, geometry.ice_area_specific_volume)
156
157 grid = geometry.ice_thickness.grid()
158
159 diff = PISM.IceModelVec2S(grid, "difference", PISM.WITHOUT_GHOSTS)
160 exact = PISM.IceModelVec2S(grid, "thk", PISM.WITHOUT_GHOSTS)
161
162 L = min(grid.Lx(), grid.Ly())
163 R_inner = 0.25 * L
164 spreading_velocity = 0.7
165 R_outer = R_inner + spreading_velocity * t_final
166
167 disc(exact, 0, 0, 1, R_inner, R_outer)
168
169 exact.add(-1.0, geometry.ice_thickness, diff)
170
171 # return the average error
172 return diff.norm(PISM.PETSc.NormType.N1) / (N*N)
173
174
175 def part_grid_convergence_test():
176 "Test that the error does go down as O(1/N)"
177
178 np.testing.assert_almost_equal([average_error(N) for N in [51, 101]],
179 [0.0338388, 0.0158498])
180
181
182 def part_grid_symmetry_test():
183 """The initial condition and the velocity fields are radially
184 symmetric, so the result should be too."""
185
186 N = 51
187
188 log.disable()
189 geometry = run(N, N, 1, True, 1.0)
190 log.enable()
191
192 # combine stuff stored as thickness and as area specific volume
193 geometry.ice_thickness.add(1.0, geometry.ice_area_specific_volume)
194
195 # convert ice thickness to a NumPy array on rank 0 -- that way we can use flipud() and
196 # fliplr().
197 H = geometry.ice_thickness.numpy()
198
199 np.testing.assert_almost_equal(H, np.flipud(H))
200 np.testing.assert_almost_equal(H, np.fliplr(H))
201 np.testing.assert_almost_equal(H, np.flipud(np.fliplr(H)))