trunTestP.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
---
trunTestP.py (9630B)
---
1 #!/usr/bin/env python3
2 from sys import argv, stderr, exit
3 import subprocess
4 import numpy as np
5 from argparse import ArgumentParser
6
7 try:
8 from PISM import exactP
9 except:
10 stderr.write("Please build PISM's Python bindings'.\n")
11 exit(1)
12
13 try:
14 from PISMNC import PISMDataset
15 except:
16 subprocess.call("ln -sf ../../util/PISMNC.py", shell=True)
17 from PISMNC import PISMDataset
18
19
20 def parse_options():
21 stderr.write("reading options ...\n")
22
23 parser = ArgumentParser()
24 parser.description = "Test P (verification of '-hydrology distributed')."
25 parser.add_argument("--pism_path", dest="PISM_PATH", default=".")
26 parser.add_argument("--mpiexec", dest="MPIEXEC", default="")
27 parser.add_argument("--Mx", dest="Mx",
28 help="Horizontal grid size. Default corresponds to a 1km grid.", type=int, default=51)
29 parser.add_argument("--keep", dest="keep", action="store_true", help="Keep the generated PISM input file.")
30
31 return parser.parse_args()
32
33
34 def generate_config():
35 """Generates the config file with custom ice softness and hydraulic conductivity."""
36
37 stderr.write("generating testPconfig.nc ...\n")
38
39 nc = PISMDataset("testPconfig.nc", 'w')
40 pism_overrides = nc.createVariable("pism_overrides", 'b')
41
42 attrs = {
43 "flow_law.isothermal_Glen.ice_softness": 3.1689e-24,
44 "flow_law.isothermal_Glen.ice_softness_doc": "Pa-3 s-1; ice softness; NOT DEFAULT",
45
46 "hydrology.hydraulic_conductivity": 1.0e-2 / (1000.0 * 9.81),
47 "hydrology.hydraulic_conductivity_doc": "= k; NOT DEFAULT",
48
49 "hydrology.regularizing_porosity": 0.01,
50 "hydrology.regularizing_porosity_doc": "[pure]; phi_0 in notes",
51
52 "hydrology.tillwat_max": 0.0,
53 "hydrology.tillwat_max_doc": "m; turn off till water mechanism",
54
55 "hydrology.thickness_power_in_flux": 1.0,
56 "hydrology.thickness_power_in_flux_doc": "; = alpha in notes",
57
58 "hydrology.gradient_power_in_flux": 2.0,
59 "hydrology.gradient_power_in_flux_doc": "; = beta in notes",
60
61 "hydrology.roughness_scale": 1.0,
62 "hydrology.roughness_scale_doc": "m; W_r in notes; roughness scale",
63
64 "basal_yield_stress.model": "constant",
65 "basal_yield_stress.model_doc": "only the constant yield stress model works without till",
66
67 "basal_yield_stress.constant.value": 1e6,
68 "basal_yield_stress.constant.value_doc": "set default to 'high tauc'",
69 }
70 keys = list(attrs.keys())
71 keys.sort()
72 for k in keys:
73 pism_overrides.setncattr(k, attrs[k])
74
75 nc.close()
76
77
78 def report_drift(name, file1, file2, xx, yy, doshow=False):
79 "Report on the difference between two files."
80 nc1 = PISMDataset(file1)
81 nc2 = PISMDataset(file2)
82
83 var1 = nc1.variables[name]
84 var2 = nc2.variables[name]
85 diff = np.abs(np.squeeze(var1[:]) - np.squeeze(var2[:]))
86
87 rr = np.sqrt(xx ** 2 + yy ** 2)
88 diff[rr >= 0.89 * 25000.0] = 0.0
89
90 if (doshow):
91 import matplotlib.pyplot as plt
92 plt.pcolormesh(xx, yy, diff)
93 plt.axis('equal')
94 plt.axis('tight')
95 plt.colorbar()
96 plt.show()
97
98 #stderr.write("Drift in %s: average = %f, max = %f [%s]" % (name, np.average(diff), np.max(diff), var1.units) + "\n")
99 return np.average(diff), np.max(diff)
100
101
102 def create_grid(Mx):
103 Lx = 25.0e3 # outside L = 22.5 km
104
105 x = np.linspace(-Lx, Lx, Mx)
106 xx, yy = np.meshgrid(x, x)
107 return x, x, xx, yy
108
109
110 def radially_outward(mag, x, y):
111 """return components of a vector field V(x,y) which is radially-outward from
112 the origin and has magnitude mag"""
113 r = np.sqrt(x * x + y * y)
114 if r == 0.0:
115 return (0.0, 0.0)
116 vx = mag * x / r
117 vy = mag * y / r
118 return (vx, vy)
119
120
121 def compute_sorted_radii(xx, yy):
122 stderr.write("sorting radial variable ...\n")
123
124 Mx = xx.shape[0]
125 # create 1D array of tuples (r,j,k), sorted by r-value
126 dtype = [('r', float), ('j', int), ('k', int)]
127 rr = np.empty((Mx, Mx), dtype=dtype)
128
129 for j in range(Mx):
130 for k in range(Mx):
131 rr[j, k] = (np.sqrt(xx[j, k] ** 2 + yy[j, k] ** 2), j, k)
132
133 r = np.sort(rr.flatten(), order='r')
134
135 return np.flipud(r)
136
137
138 def generate_pism_input(x, y, xx, yy):
139 stderr.write("calling exactP() ...\n")
140
141 EPS_ABS = 1.0e-12
142 EPS_REL = 1.0e-15
143
144 p = exactP(r[:]['r'], EPS_ABS, EPS_REL, 1)
145 h_r, magvb_r, W_r, P_r = p.h, p.magvb, p.W, p.P
146
147 stderr.write("creating gridded variables ...\n")
148 # put on grid
149 h = np.zeros_like(xx)
150 W = np.zeros_like(xx)
151 P = np.zeros_like(xx)
152
153 magvb = np.zeros_like(xx)
154 ussa = np.zeros_like(xx)
155 vssa = np.zeros_like(xx)
156
157 for n, pt in enumerate(r):
158 j = pt['j']
159 k = pt['k']
160 h[j, k] = h_r[n] # ice thickness in m
161 magvb[j, k] = magvb_r[n] # sliding speed in m s-1
162 ussa[j, k], vssa[j, k] = radially_outward(magvb[j, k], xx[j, k], yy[j, k])
163 W[j, k] = W_r[n] # water thickness in m
164 P[j, k] = P_r[n] # water pressure in Pa
165
166 stderr.write("creating inputforP.nc ...\n")
167
168 nc = PISMDataset("inputforP.nc", 'w')
169 nc.create_dimensions(x, y, time_dependent=True, use_time_bounds=True)
170
171 nc.define_2d_field("thk", time_dependent=False,
172 attrs={"long_name": "ice thickness",
173 "units": "m",
174 "valid_min": 0.0,
175 "standard_name": "land_ice_thickness"})
176 nc.define_2d_field("topg", time_dependent=False,
177 attrs={"long_name": "bedrock topography",
178 "units": "m",
179 "standard_name": "bedrock_altitude"})
180 nc.define_2d_field("climatic_mass_balance", time_dependent=False,
181 attrs={"long_name": "climatic mass balance for -surface given",
182 "units": "kg m-2 year-1",
183 "standard_name": "land_ice_surface_specific_mass_balance"})
184 nc.define_2d_field("ice_surface_temp", time_dependent=False,
185 attrs={"long_name": "ice surface temp (K) for -surface given",
186 "units": "Kelvin",
187 "valid_min": 0.0})
188 nc.define_2d_field("bmelt", time_dependent=False,
189 attrs={"long_name": "basal melt rate",
190 "units": "m year-1",
191 "standard_name": "land_ice_basal_melt_rate"})
192
193 nc.define_2d_field("bwat", time_dependent=False,
194 attrs={"long_name": "thickness of basal water layer",
195 "units": "m",
196 "valid_min": 0.0})
197 nc.define_2d_field("bwp", time_dependent=False,
198 attrs={"long_name": "water pressure in basal water layer",
199 "units": "Pa",
200 "valid_min": 0.0})
201
202 nc.define_2d_field("bc_mask", time_dependent=False,
203 attrs={"long_name": "if =1, apply u_ssa_bc and v_ssa_bc as sliding velocity"})
204 nc.define_2d_field("u_ssa_bc", time_dependent=False,
205 attrs={"long_name": "x-component of prescribed sliding velocity",
206 "units": "m s-1"})
207 nc.define_2d_field("v_ssa_bc", time_dependent=False,
208 attrs={"long_name": "y-component of prescribed sliding velocity",
209 "units": "m s-1"})
210
211 Phi0 = 0.20 # 20 cm/year basal melt rate
212 T_surface = 260 # ice surface temperature, K
213
214 variables = {"topg": np.zeros_like(xx),
215 "climatic_mass_balance": np.zeros_like(xx),
216 "ice_surface_temp": np.ones_like(xx) + T_surface,
217 "bmelt": np.zeros_like(xx) + Phi0,
218 "thk": h,
219 "bwat": W,
220 "bwp": P,
221 "bc_mask": np.ones_like(xx),
222 "u_ssa_bc": ussa,
223 "v_ssa_bc": vssa}
224
225 for name in list(variables.keys()):
226 nc.write(name, variables[name])
227
228 nc.history = subprocess.list2cmdline(argv)
229 nc.close()
230 stderr.write("NetCDF file %s written\n" % "inputforP.nc")
231
232
233 def run_pism(opts):
234 stderr.write("Testing: Test P verification of '-hydrology distributed'.\n")
235
236 cmd = "%s %s/pismr -config_override testPconfig.nc -i inputforP.nc -bootstrap -Mx %d -My %d -Mz 11 -Lz 4000 -hydrology distributed -report_mass_accounting -y 0.08333333333333 -max_dt 0.01 -no_mass -energy none -stress_balance ssa+sia -ssa_dirichlet_bc -o end.nc" % (
237 opts.MPIEXEC, opts.PISM_PATH, opts.Mx, opts.Mx)
238
239 stderr.write(cmd + "\n")
240 subprocess.call(cmd, shell=True)
241
242 # high-res and parallel example:
243 # ./runTestP.py --pism_path=../../build --mpiexec="mpiexec -n 4" --Mx=201
244 # example which should suffice for regression:
245 # ./runTestP.py --pism_path=../../build --Mx=21
246
247
248 if __name__ == "__main__":
249
250 opts = parse_options()
251
252 x, y, xx, yy = create_grid(opts.Mx)
253
254 r = compute_sorted_radii(xx, yy)
255
256 generate_config()
257
258 generate_pism_input(x, y, xx, yy)
259
260 run_pism(opts)
261
262 (bwatav, bwatmax) = report_drift("bwat", "inputforP.nc", "end.nc", xx, yy, doshow=False)
263 (bwpav, bwpmax) = report_drift("bwp", "inputforP.nc", "end.nc", xx, yy, doshow=False)
264
265 print("NUMERICAL ERRORS:")
266 print("%d %f %f %f %f\n" % (opts.Mx, bwatav, bwatmax, bwpav, bwpmax))
267
268 # cleanup:
269 if opts.keep == False:
270 subprocess.call("rm testPconfig.nc inputforP.nc end.nc", shell=True)