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()