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)