run.py (13935B)
1 #!/usr/bin/env python3 2 import MISMIP 3 4 # This scripts generates bash scripts that run MISMIP experiments and generates 5 # all the necessary input files. 6 # 7 # Run run.py > my_new_mismip.sh and use that. 8 9 try: 10 from netCDF4 import Dataset as NC 11 except: 12 print("netCDF4 is not installed!") 13 sys.exit(1) 14 15 import sys 16 17 # The "standard" preamble used in many PISM scripts: 18 preamble = ''' 19 if [ -n "${SCRIPTNAME:+1}" ] ; then 20 echo "[SCRIPTNAME=$SCRIPTNAME (already set)]" 21 echo "" 22 else 23 SCRIPTNAME="#(mismip.sh)" 24 fi 25 26 echo 27 echo "# ==================================================================================" 28 echo "# MISMIP experiments" 29 echo "# ==================================================================================" 30 echo 31 32 set -e # exit on error 33 34 NN=2 # default number of processors 35 if [ $# -gt 0 ] ; then # if user says "mismip.sh 8" then NN = 8 36 NN="$1" 37 fi 38 39 echo "$SCRIPTNAME NN = $NN" 40 41 # set MPIDO if using different MPI execution command, for example: 42 # $ export PISM_MPIDO="aprun -n " 43 if [ -n "${PISM_MPIDO:+1}" ] ; then # check if env var is already set 44 echo "$SCRIPTNAME PISM_MPIDO = $PISM_MPIDO (already set)" 45 else 46 PISM_MPIDO="mpiexec -n " 47 echo "$SCRIPTNAME PISM_MPIDO = $PISM_MPIDO" 48 fi 49 50 # check if env var PISM_DO was set (i.e. PISM_DO=echo for a 'dry' run) 51 if [ -n "${PISM_DO:+1}" ] ; then # check if env var DO is already set 52 echo "$SCRIPTNAME PISM_DO = $PISM_DO (already set)" 53 else 54 PISM_DO="" 55 fi 56 57 # prefix to pism (not to executables) 58 if [ -n "${PISM_BIN:+1}" ] ; then # check if env var is already set 59 echo "$SCRIPTNAME PISM_BIN = $PISM_BIN (already set)" 60 else 61 PISM_BIN="" # just a guess 62 echo "$SCRIPTNAME PISM_BIN = $PISM_BIN" 63 fi 64 ''' 65 66 67 class Experiment: 68 69 "A MISMIP experiment." 70 experiment = "" 71 mode = 1 72 model = 1 73 semianalytic = True 74 Mx = 151 75 My = 3 76 Mz = 15 77 initials = "ABC" 78 executable = "$PISM_DO $PISM_MPIDO $NN ${PISM_BIN}pismr" 79 80 def __init__(self, experiment, model=1, mode=1, Mx=None, Mz=15, semianalytic=True, 81 initials="ABC", executable=None): 82 self.model = model 83 self.mode = mode 84 self.experiment = experiment 85 self.initials = initials 86 self.semianalytic = semianalytic 87 88 if executable: 89 self.executable = executable 90 91 if mode == 3: 92 self.Mx = Mx 93 else: 94 self.Mx = 2 * MISMIP.N(self.mode) + 1 95 96 self.My = 3 97 98 if self.experiment == "2b": 99 self.Lz = 7000 100 else: 101 self.Lz = 6000 102 self.Lz *= 2.0 103 104 def physics_options(self, input_file, step): 105 "Options corresponding to modeling choices." 106 config_filename = self.config(step) 107 108 #options = ["-energy none", # isothermal setup; allows selecting cold-mode flow laws 109 #"-ssa_flow_law isothermal_glen", # isothermal setup 110 #"-yield_stress constant", 111 #"-tauc %e" % MISMIP.C(self.experiment), 112 #"-pseudo_plastic", 113 #"-gradient eta", 114 #"-pseudo_plastic_q %e" % MISMIP.m(self.experiment), 115 #"-pseudo_plastic_uthreshold %e" % MISMIP.secpera(), 116 #"-front_retreat_file %s" % input_file, # prescribe the maximum ice extent 117 #"-config_override %s" % config_filename, 118 #"-ssa_method fd", 119 #"-cfbc", # calving front boundary conditions 120 #"-part_grid", # sub-grid front motion parameterization 121 #"-ssafd_ksp_rtol 1e-7", 122 #"-ys 0", 123 #"-ye %d" % MISMIP.run_length(self.experiment, step), 124 #"-options_left", 125 #] 126 127 128 options = ["-stress_balance ssa+sia", 129 #"-ssa_flow_law isothermal_glen", # isothermal setup 130 #"-yield_stress constant", 131 #"-tauc %e" % MISMIP.C(self.experiment), 132 #"-pseudo_plastic", 133 "-gradient eta", 134 #"-pseudo_plastic_q %e" % MISMIP.m(self.experiment), 135 #"-pseudo_plastic_uthreshold %e" % MISMIP.secpera(), 136 #"-yield_stress mohr_coulomb", 137 #"-till_flux", 138 "-front_retreat_file %s" % input_file, # prescribe the maximum ice extent 139 "-config_override %s" % config_filename, 140 #"-ssa_method fd", 141 "-cfbc", # PIK: calving front boundary conditions 142 "-part_grid", # PIK: sub-grid front motion parameterization 143 "-kill_icebergs", # PIK: https://pism-docs.org/sphinx/manual/modeling-choices/marine/pik.html#iceberg-removal 144 "-subgl", # PIK: https://pism-docs.org/sphinx/manual/modeling-choices/marine/pik.html#sub-grid-treatment-of-the-grounding-line-position 145 #"-ssafd_ksp_rtol 1e-7", 146 "-ys 0", 147 "-ye 1e3", 148 #"-ye %d" % MISMIP.run_length(self.experiment, step), 149 "-options_left", 150 "-hydrology routing", 151 "-plastic_phi 24", # Tulaczyk et al 2000 152 "-till_cohesion 0", # Tulaczyk et al 2000 153 "-geothermal_flux 70e-3", # Feldmann and Levermann 2017 154 "-sia_e 4.5", # Martin et al 2010 155 "-ssa_e 0.512", # Martin et al 2010 156 "-bed_def lc", # Lingle and Clark 1985, Bueler et al 2007 157 "-stress_balance.sia.max_diffusivity 1e4", 158 "-sea_level constant,delta_sl -ocean_delta_sl_file sealvl.nc", 159 ] 160 161 #options = [, 162 #"-yield_stress mohr_coulomb", 163 #"-Lz 1e4", 164 ##"-till_flux", 165 #"-cfbc", # calving front boundary conditions 166 #"-part_grid", # sub-grid front motion parameterization 167 #"-front_retreat_file %s" % input_file, # prescribe the maximum ice extent 168 #"-ys 0", 169 #"-ye %d" % MISMIP.run_length(self.experiment, step), 170 #"-options_left", 171 #] 172 173 return options 174 175 def config(self, step): 176 '''Generates a config file containing flags and parameters 177 for a particular experiment and step. 178 179 This takes care of flags and parameters that *cannot* be set using 180 command-line options. (We try to use command-line options whenever we can.) 181 ''' 182 183 filename = "MISMIP_conf_%s_A%s.nc" % (self.experiment, step) 184 185 nc = NC(filename, 'w', format="NETCDF3_CLASSIC") 186 187 var = nc.createVariable("pism_overrides", 'i') 188 189 attrs = {"geometry.update.use_basal_melt_rate": "no", 190 "stress_balance.ssa.compute_surface_gradient_inward": "no", 191 "flow_law.isothermal_Glen.ice_softness": MISMIP.A(self.experiment, step), 192 "constants.ice.density": MISMIP.rho_i(), 193 "constants.sea_water.density": MISMIP.rho_w(), 194 "bootstrapping.defaults.geothermal_flux": 0.0, 195 "stress_balance.ssa.Glen_exponent": MISMIP.n(), 196 "constants.standard_gravity": MISMIP.g(), 197 "ocean.sub_shelf_heat_flux_into_ice": 0.0, 198 } 199 200 if self.model != 1: 201 attrs["stress_balance.sia.bed_smoother.range"] = 0.0 202 203 for name, value in attrs.items(): 204 var.setncattr(name, value) 205 206 nc.close() 207 208 return filename 209 210 def bootstrap_options(self, step): 211 boot_filename = "MISMIP_boot_%s_M%s_A%s.nc" % (self.experiment, self.mode, step) 212 213 import prepare 214 prepare.pism_bootstrap_file(boot_filename, self.experiment, step, self.mode, N=self.Mx, 215 semianalytical_profile=self.semianalytic) 216 217 options = ["-i %s -bootstrap" % boot_filename, 218 "-Mx %d" % self.Mx, 219 "-My %d" % self.My, 220 "-Mz %d" % self.Mz, 221 "-Lz %d" % self.Lz] 222 223 return options, boot_filename 224 225 def output_options(self, step): 226 output_file = self.output_filename(self.experiment, step) 227 extra_file = "ex_" + output_file 228 ts_file = "ts_" + output_file 229 230 options = ["-extra_file %s" % extra_file, 231 "-extra_times 0:50:3e4", 232 "-extra_vars thk,topg,velbar_mag,flux_mag,mask,dHdt,usurf,hardav,velbase_mag,nuH,tauc,taud,taub,flux_divergence,cell_grounded_fraction", 233 "-ts_file %s" % ts_file, 234 "-ts_times 0:50:3e4", 235 "-backup_size big", 236 "-o %s" % output_file, 237 "-o_order zyx", 238 ] 239 240 return output_file, options 241 242 def output_filename(self, experiment, step): 243 return "%s%s_%s_M%s_A%s.nc" % (self.initials, self.model, experiment, self.mode, step) 244 245 def options(self, step, input_file=None): 246 '''Generates a string of PISM options corresponding to a MISMIP experiment.''' 247 248 if input_file is None: 249 input_options, input_file = self.bootstrap_options(step) 250 else: 251 input_options = ["-i %s" % input_file] 252 253 physics = self.physics_options(input_file, step) 254 255 output_file, output_options = self.output_options(step) 256 257 return output_file, (input_options + physics + output_options) 258 259 def run_step(self, step, input_file=None): 260 out, opts = self.options(step, input_file) 261 print('echo "# Step %s-%s"' % (self.experiment, step)) 262 print("%s %s" % (self.executable, ' '.join(opts))) 263 print('echo "Done."\n') 264 265 return out 266 267 def run(self, step=None): 268 print('echo "# Experiment %s"' % self.experiment) 269 270 if self.experiment in ('1a', '1b'): 271 # bootstrap 272 input_file = None 273 # steps 1 to 9 274 steps = list(range(1, 10)) 275 276 if self.experiment in ('2a', '2b'): 277 # start from step 9 of the corresponding experiment 1 278 input_file = self.output_filename(self.experiment.replace("2", "1"), 9) 279 # steps 8 to 1 280 steps = list(range(8, 0, -1)) 281 282 if self.experiment == '3a': 283 # bootstrap 284 input_file = None 285 # steps 1 to 13 286 steps = list(range(1, 14)) 287 288 if self.experiment == '3b': 289 # bootstrap 290 input_file = None 291 # steps 1 to 15 292 steps = list(range(1, 16)) 293 294 if step is not None: 295 input_file = None 296 steps = [step] 297 298 for step in steps: 299 input_file = self.run_step(step, input_file) 300 301 302 def run_mismip(initials, executable, semianalytic): 303 Mx = 601 304 models = (1, 2) 305 modes = (1, 2, 3) 306 experiments = ('1a', '1b', '2a', '2b', '3a', '3b') 307 308 print(preamble) 309 310 for model in models: 311 for mode in modes: 312 for experiment in experiments: 313 e = Experiment(experiment, 314 initials=initials, 315 executable=executable, 316 model=model, mode=mode, Mx=Mx, 317 semianalytic=semianalytic) 318 e.run() 319 320 321 if __name__ == "__main__": 322 from optparse import OptionParser 323 324 parser = OptionParser() 325 326 parser.usage = "%prog [options]" 327 parser.description = "Creates a script running MISMIP experiments." 328 parser.add_option("--initials", dest="initials", type="string", 329 help="Initials (3 letters)", default="ABC") 330 parser.add_option("-e", "--experiment", dest="experiment", type="string", 331 default='1a', 332 help="MISMIP experiments (one of '1a', '1b', '2a', '2b', '3a', '3b')") 333 parser.add_option("-s", "--step", dest="step", type="int", 334 help="MISMIP step number") 335 parser.add_option("-u", "--uniform_thickness", 336 action="store_false", dest="semianalytic", default=True, 337 help="Use uniform 10 m ice thickness") 338 parser.add_option("-a", "--all", 339 action="store_true", dest="all", default=False, 340 help="Run all experiments") 341 parser.add_option("-m", "--mode", dest="mode", type="int", default=1, 342 help="MISMIP grid mode") 343 parser.add_option("--Mx", dest="Mx", type="int", default=601, 344 help="Custom grid size; use with --mode=3") 345 parser.add_option("--model", dest="model", type="int", default=1, 346 help="Models: 1 - SSA only; 2 - SIA+SSA") 347 parser.add_option("--executable", dest="executable", type="string", 348 help="Executable to run, e.g. 'mpiexec -n 4 pismr'") 349 350 (opts, args) = parser.parse_args() 351 352 if opts.all: 353 run_mismip(opts.initials, opts.executable, opts.semianalytic) 354 exit(0) 355 356 def escape(arg): 357 if arg.find(" ") >= 0: 358 parts = arg.split("=") 359 return "%s=\"%s\"" % (parts[0], ' '.join(parts[1:])) 360 else: 361 return arg 362 363 arg_list = [escape(a) for a in sys.argv] 364 365 print("#!/bin/bash") 366 print("# This script was created by examples/mismip/run.py. The command was:") 367 print("# %s" % (' '.join(arg_list))) 368 369 if opts.executable is None: 370 print(preamble) 371 372 e = Experiment(opts.experiment, 373 initials=opts.initials, 374 executable=opts.executable, 375 model=opts.model, 376 mode=opts.mode, 377 Mx=opts.Mx, 378 semianalytic=opts.semianalytic) 379 380 if opts.step is not None: 381 e.run(opts.step) 382 else: 383 e.run()