torographic_precipitation.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
---
torographic_precipitation.py (5851B)
---
1 #!/usr/bin/env python3
2 import numpy as np
3
4 import PISM
5
6 # silence initialization messages
7 PISM.Context().log.set_threshold(1)
8
9 def triangle_ridge_grid(dx=5e4, dy=5e4):
10 "Allocate the grid for the synthetic geometry test."
11 x_min, x_max = -100e3, 100e3
12 y_min, y_max = -100e3, 100e3
13
14 x0 = (x_max + x_min) / 2.0
15 y0 = (y_max + y_min) / 2.0
16
17 Lx = (x_max - x_min) / 2.0
18 Ly = (y_max - y_min) / 2.0
19 Mx = int((x_max - x_min) / dx)
20 My = int((y_max - y_min) / dy)
21
22 return PISM.IceGrid_Shallow(PISM.Context().ctx,
23 Lx, Ly,
24 x0, y0,
25 Mx, My,
26 PISM.CELL_CORNER, PISM.NOT_PERIODIC)
27
28 def triangle_ridge(x, A=500.0, d=50e3):
29 "Create the 'triangle ridge' topography"
30 return np.maximum(A * (1 - np.fabs(x) / d), 0)
31
32 def triangle_ridge_exact(x, u, Cw, tau, A=500.0, d=50e3):
33 """Exact precipitation for the triangle ridge setup.
34
35 See equations 44, 45, 46 in Smith and Barstad (2004).
36 """
37 assert d > 0
38
39 C = Cw * u * A / d
40 Ut = u * tau
41
42 xc = Ut * np.log(2 - np.exp(-d / Ut))
43
44 def P(x):
45 if x < 0 and x >= -d:
46 return C * (1.0 - np.exp(-(x + d) / Ut))
47 elif x >= 0 and x <= xc:
48 return C * (np.exp(-x / Ut) * (2 - np.exp(-d / Ut)) - 1)
49 else:
50 return 0
51
52 try:
53 return np.array([P(t) for t in x])
54 except TypeError:
55 return P(x)
56
57 def run_model(grid, orography):
58 "Run the PISM implementation of the model to compare to the Python version."
59
60 model = PISM.AtmosphereOrographicPrecipitation(grid, PISM.AtmosphereUniform(grid))
61 geometry = PISM.Geometry(grid)
62
63 with PISM.vec.Access(nocomm=geometry.ice_thickness):
64 for i,j in grid.points():
65 geometry.ice_thickness[i, j] = orography[j, i]
66
67 geometry.bed_elevation.set(0.0)
68 geometry.sea_level_elevation.set(0.0)
69 geometry.ice_area_specific_volume.set(0.0)
70
71 # compute surface elevation from ice thickness and bed elevation
72 geometry.ensure_consistency(0)
73
74 model.init(geometry)
75 model.update(geometry, 0, 1)
76
77 config = PISM.Context().config
78 water_density = config.get_number("constants.fresh_water.density")
79
80 # convert from kg / (m^2 s) to mm/s
81 return model.mean_precipitation().numpy() / (1e-3 * water_density)
82
83 def max_error(spacing, wind_direction):
84 # Set conversion time to zero (we could set fallout time to zero instead: it does not
85 # matter which one is zero)
86
87 wind_speed = 15
88
89 config = PISM.Context().config
90
91 # set wind speed and direction
92 config.set_number("atmosphere.orographic_precipitation.wind_speed", wind_speed)
93 config.set_number("atmosphere.orographic_precipitation.wind_direction", wind_direction)
94
95 # set conversion time to zero
96 config.set_number("atmosphere.orographic_precipitation.conversion_time", 0.0)
97 # eliminate the effect of airflow dynamics
98 config.set_number("atmosphere.orographic_precipitation.water_vapor_scale_height", 0.0)
99 # eliminate the effect of the Coriolis force
100 config.set_number("atmosphere.orographic_precipitation.coriolis_latitude", 0.0)
101
102 # get constants needed to compute the exact solution
103 tau = config.get_number("atmosphere.orographic_precipitation.fallout_time")
104 Theta_m = config.get_number("atmosphere.orographic_precipitation.moist_adiabatic_lapse_rate")
105 rho_Sref = config.get_number("atmosphere.orographic_precipitation.reference_density")
106 gamma = config.get_number("atmosphere.orographic_precipitation.lapse_rate")
107 Cw = rho_Sref * Theta_m / gamma
108
109 if wind_direction == 90 or wind_direction == 270:
110 # east or west
111 grid = triangle_ridge_grid(dx=spacing)
112 t = np.array(grid.x())
113
114 h = triangle_ridge(t)
115 orography = np.tile(h, (grid.My(), 1))
116
117 P = run_model(grid, orography)
118 P = P[grid.My() // 2, :]
119 else:
120 # north or south
121 grid = triangle_ridge_grid(dy=spacing)
122 t = np.array(grid.y())
123
124 h = triangle_ridge(t)
125 orography = np.tile(h, (grid.Mx(), 1)).T
126
127 P = run_model(grid, orography)
128 P = P[:, grid.Mx() // 2]
129
130 if wind_direction == 0 or wind_direction == 90:
131 P_exact = triangle_ridge_exact(-t, wind_speed, Cw, tau)
132 else:
133 P_exact = triangle_ridge_exact(t, wind_speed, Cw, tau)
134
135 return np.max(np.fabs(P - P_exact))
136
137 def convergence_rate(dxs, error, wind_direction, plot):
138 errors = [error(dx, wind_direction) for dx in dxs]
139
140 p = np.polyfit(np.log10(dxs), np.log10(errors), 1)
141
142 if plot:
143 import pylab as plt
144
145 direction = {0 : "north", 90 : "east", 180 : "south", 270 : "west"}
146
147 def log_plot(x, y, style, label):
148 plt.plot(np.log10(x), np.log10(y), style, label=label)
149 plt.xticks(np.log10(x), x)
150
151 def log_fit_plot(x, p, label):
152 plt.plot(np.log10(x), np.polyval(p, np.log10(x)), label=label)
153
154 plt.figure()
155 plt.title("Precipitation errors (wind direction: {})".format(direction[wind_direction]))
156 log_fit_plot(dxs, p, "polynomial fit (dx^{:1.4})".format(p[0]))
157 log_plot(dxs, errors, 'o', "errors")
158 plt.legend()
159 plt.grid()
160 plt.xlabel("grid spacing (meters)")
161 plt.ylabel("log10(error)")
162 plt.show()
163
164 return p[0]
165
166 def ltop_test(dxs=[2000, 1000, 500], plot=False):
167 "Orographic precipitation (triangle ridge test case)"
168
169 assert convergence_rate(dxs, max_error, 0, plot) > 1.99
170 assert convergence_rate(dxs, max_error, 90, plot) > 1.99
171 assert convergence_rate(dxs, max_error, 180, plot) > 1.99
172 assert convergence_rate(dxs, max_error, 270, plot) > 1.99
173
174 if __name__ == "__main__":
175 ltop_test(dxs=[2000, 1000, 500, 250, 125], plot=True)