ttest_29.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 --- ttest_29.py (4147B) --- 1 #!/usr/bin/env python3 2 3 import subprocess 4 import shutil 5 import shlex 6 import os 7 from sys import exit, argv 8 from netCDF4 import Dataset as NC 9 import numpy as np 10 11 12 def process_arguments(): 13 from argparse import ArgumentParser 14 parser = ArgumentParser() 15 parser.add_argument("PISM_PATH") 16 parser.add_argument("MPIEXEC") 17 parser.add_argument("PISM_SOURCE_DIR") 18 19 return parser.parse_args() 20 21 22 def copy_input(opts): 23 shutil.copy(os.path.join(opts.PISM_SOURCE_DIR, "test/test_hydrology/inputforP_regression.nc"), ".") 24 25 26 def generate_config(): 27 """Generates the config file with custom ice softness and hydraulic conductivity.""" 28 29 print("generating testPconfig.nc ...") 30 31 nc = NC("testPconfig.nc", 'w') 32 pism_overrides = nc.createVariable("pism_overrides", 'b') 33 34 attrs = { 35 "constants.standard_gravity": 9.81, 36 "constants.standard_gravity_doc": "m s-2; = g; acceleration due to gravity on Earth geoid", 37 38 "constants.fresh_water.density": 1000.0, 39 "constants.fresh_water.density_doc": "kg m-3; = rhow", 40 41 "flow_law.isothermal_Glen.ice_softness": 3.1689e-24, 42 "flow_law.isothermal_Glen.ice_softness_doc": "Pa-3 s-1; ice softness; NOT DEFAULT", 43 44 "hydrology.hydraulic_conductivity": 1.0e-2 / (1000.0 * 9.81), 45 "hydrology.hydraulic_conductivity_doc": "= k; NOT DEFAULT", 46 47 "hydrology.tillwat_max": 0.0, 48 "hydrology.tillwat_max_doc": "m; turn off till water mechanism", 49 50 "hydrology.thickness_power_in_flux": 1.0, 51 "hydrology.thickness_power_in_flux_doc": "; = alpha in notes", 52 53 "hydrology.gradient_power_in_flux": 2.0, 54 "hydrology.gradient_power_in_flux_doc": "; = beta in notes", 55 56 "hydrology.roughness_scale": 1.0, 57 "hydrology.roughness_scale_doc": "m; W_r in notes; roughness scale", 58 59 "hydrology.regularizing_porosity": 0.01, 60 "hydrology.regularizing_porosity_doc": "[pure]; phi_0 in notes", 61 62 "basal_yield_stress.model": "constant", 63 "basal_yield_stress.model_doc": "only the constant yield stress model works without till", 64 65 "basal_yield_stress.constant.value": 1e6, 66 "basal_yield_stress.constant.value_doc": "set default to 'high tauc'" 67 } 68 69 for k, v in list(attrs.items()): 70 pism_overrides.setncattr(k, v) 71 72 nc.close() 73 74 75 def run_pism(opts): 76 cmd = "%s %s/pismr -config_override testPconfig.nc -i inputforP_regression.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" % ( 77 opts.MPIEXEC, opts.PISM_PATH, 21, 21) 78 79 print(cmd) 80 subprocess.call(shlex.split(cmd)) 81 82 83 def check_drift(file1, file2): 84 nc1 = NC(file1) 85 nc2 = NC(file2) 86 87 stored_drift = {'bwat_max': 0.023524626576411189, 88 'bwp_max': 79552.478734239354, 89 'bwp_avg': 6261.1642337484445, 90 'bwat_avg': 0.0034449380393343091} 91 92 drift = {} 93 for name in ("bwat", "bwp"): 94 var1 = nc1.variables[name] 95 var2 = nc2.variables[name] 96 diff = np.abs(np.squeeze(var1[:]) - np.squeeze(var2[:])) 97 98 drift["%s_max" % name] = np.max(diff) 99 drift["%s_avg" % name] = np.average(diff) 100 101 print("drift = ", drift) 102 print("stored_drift = ", stored_drift) 103 104 for name in list(drift.keys()): 105 rel_diff = np.abs(stored_drift[name] - drift[name]) / stored_drift[name] 106 107 if rel_diff > 1e-3: 108 print("Stored and computed drifts in %s differ: %f != %f" % (name, stored_drift[name], drift[name])) 109 exit(1) 110 111 112 def cleanup(): 113 for fname in ("inputforP_regression.nc", "testPconfig.nc", "end.nc"): 114 os.remove(fname) 115 116 117 if __name__ == "__main__": 118 opts = process_arguments() 119 120 print("Copying input files...") 121 copy_input(opts) 122 123 print("Generating the -config_override file...") 124 generate_config() 125 126 print("Running PISM...") 127 run_pism(opts) 128 129 print("Checking the drift...") 130 check_drift("inputforP_regression.nc", "end.nc") 131 132 print("Cleaning up...") 133 cleanup()