trigid_rotation.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
---
trigid_rotation.py (5114B)
---
1 import numpy as np
2 import PISM
3
4 """This is a half-baked test for the mass transport code that sets up
5 a rotational velocity field and then steps forward in time for one revolution.
6
7 It is not ready to be used as a regression or a verification test. All
8 I can see right now is that our mass transport scheme is very
9 diffusive, which is not a surprise."""
10
11 log = PISM.Context().log
12
13
14 def disc(thickness, x0, y0, H, R):
15 """Set ice thickness to H within the disc centered at (x0,y0) of
16 radius R and 0 elsewhere.
17 """
18
19 grid = thickness.grid()
20
21 with PISM.vec.Access(nocomm=thickness):
22 for (i, j) in grid.points():
23 x = grid.x(i)
24 y = grid.y(j)
25 d2 = (x - x0)**2 + (y - y0)**2
26 r2 = R**2
27 if d2 <= r2:
28 thickness[i, j] = H
29 else:
30 thickness[i, j] = 0.0
31
32 thickness.update_ghosts()
33
34
35 def set_ice_thickness(output, time):
36 """Exact solution at time time. Corresponds to a disc that is rotated
37 around the origin (one revolution per time unit)."""
38
39 grid = output.grid()
40
41 L = min(grid.Lx(), grid.Ly())
42
43 # 1 revolution per 1 time unit
44 phi = 2 * np.pi * time
45
46 M = np.matrix([[np.cos(phi), -np.sin(phi)],
47 [np.sin(phi), np.cos(phi)]])
48
49 # center coordinates at time 0
50 x0, y0 = 0.5 * L, 0.0
51
52 R = 0.25 * L
53 x, y = M * np.matrix([x0, y0]).T
54
55 disc(output, x, y, 1, R)
56
57
58 def set_velocity(v):
59 """Initialize the velocity field to a rigid rotation around the
60 origin.
61
62 """
63 grid = v.grid()
64
65 radial_velocity = 2 * np.pi
66
67 with PISM.vec.Access(nocomm=v):
68 for (i, j) in grid.points():
69 x = grid.x(i)
70 y = grid.y(j)
71
72 v[i, j].u = y * radial_velocity
73 v[i, j].v = -x * radial_velocity
74
75 v.update_ghosts()
76
77
78 def quiver(v, **kwargs):
79 a = v.numpy()
80 plt.quiver(a[:, :, 0], a[:, :, 1], **kwargs)
81
82
83 class MassTransport(object):
84
85 def __init__(self, grid, part_grid=False):
86
87 self.grid = grid
88
89 if part_grid:
90 grid.ctx().config().set_flag("geometry.part_grid.enabled", True)
91
92 self.v = PISM.IceModelVec2V(grid, "velocity", PISM.WITHOUT_GHOSTS)
93 self.Q = PISM.IceModelVec2Stag(grid, "Q", PISM.WITHOUT_GHOSTS)
94 self.v_bc_mask = PISM.IceModelVec2Int(grid, "v_bc_mask", PISM.WITHOUT_GHOSTS)
95 self.H_bc_mask = PISM.IceModelVec2Int(grid, "H_bc_mask", PISM.WITHOUT_GHOSTS)
96
97 self.ge = PISM.GeometryEvolution(grid)
98
99 self.geometry = PISM.Geometry(grid)
100
101 self.reset()
102
103 def reset(self):
104 geometry = self.geometry
105 # grid info
106 geometry.latitude.set(0.0)
107 geometry.longitude.set(0.0)
108 # environment
109 geometry.bed_elevation.set(-10.0)
110 geometry.sea_level_elevation.set(0.0)
111 # ice
112 set_ice_thickness(geometry.ice_thickness, time=1)
113 geometry.ice_area_specific_volume.set(0.0)
114
115 geometry.ensure_consistency(0.0)
116
117 set_velocity(self.v)
118
119 def plot_thickness(self, levels, title):
120 import pylab as plt
121 cm = plt.contour(self.grid.x(), self.grid.y(),
122 self.geometry.ice_thickness.numpy(), levels=levels)
123 plt.clabel(cm)
124 plt.grid()
125 plt.title(title)
126
127 def step(self, t_final, C=1):
128 geometry = self.geometry
129 t = 0.0
130 j = 0
131 while t < t_final:
132
133 dt = PISM.max_timestep_cfl_2d(geometry.ice_thickness,
134 geometry.cell_type,
135 self.v).dt_max.value() * C
136
137 if t + dt > t_final:
138 dt = t_final - t
139
140 log.message(2, "{}, {}\n".format(t, dt))
141
142 self.ge.flow_step(geometry, dt,
143 self.v,
144 self.Q,
145 self.v_bc_mask,
146 self.H_bc_mask)
147
148 geometry.ice_thickness.add(1.0, self.ge.thickness_change_due_to_flow())
149 geometry.ice_area_specific_volume.add(1.0, self.ge.area_specific_volume_change_due_to_flow())
150 geometry.ensure_consistency(0.0)
151
152 t += dt
153 j += 1
154
155 def volume(geometry):
156 cell_area = geometry.ice_thickness.grid().cell_area()
157 volume = ((geometry.ice_thickness.numpy() + geometry.ice_area_specific_volume.numpy()) *
158 cell_area)
159 return volume.sum()
160
161 def test():
162 ctx = PISM.Context()
163 Mx = int(ctx.config.get_number("grid.Mx"))
164 My = int(ctx.config.get_number("grid.My"))
165 grid = PISM.IceGrid_Shallow(ctx.ctx, 1, 1, 0, 0, Mx, My, PISM.CELL_CORNER, PISM.NOT_PERIODIC)
166
167 import pylab as plt
168 plt.rcParams['figure.figsize'] = (8.0, 8.0)
169
170 mt = MassTransport(grid)
171
172 mt.reset()
173 levels = np.linspace(0, 1, 11)
174
175 t = 0
176 dt = 0.333
177 C = 1.0
178
179 for j in [1, 2, 4, 3]:
180 plt.subplot(2, 2, j)
181 mt.plot_thickness(levels, "time={}, V={}".format(t, volume(mt.geometry)))
182
183 mt.step(dt, C)
184 t += dt
185
186 plt.show()
187
188
189 if __name__ == "__main__":
190 test()