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)