trun.py - pism-exp-gsw - ice stream and sediment transport experiments
(HTM) git clone git://src.adamsgaard.dk/pism-exp-gsw
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
trun.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()