commit 8b985b054a316b0cc9f72041ba2f3fc9ae284f33
parent 83a341c191e31e46e6410801802061491a48f445
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 29 Nov 2021 13:46:07 +0100
place all initialization syntax in Makefile
Diffstat:
| M | Makefile | | | 105 | ++++++++++++++++++++++++++++++++++++++++++------------------------------------- |
| M | plot-time-series.py | | | 16 | ++++++++-------- |
| D | run.py | | | 383 | ------------------------------------------------------------------------------- |
3 files changed, 64 insertions(+), 440 deletions(-)
diff --git a/Makefile b/Makefile
@@ -1,59 +1,66 @@
NPROC != awk '/^cpu cores/ {print $$4; exit}' /proc/cpuinfo
SLSERIES=sealvl.nc
+T_END=1e3
all: \
- ABC1_1a_M1_A1-flux.pdf\
- ex_1a_M1_A1_1d-profile.pdf\
- #ABC1_1a_M3_A1-flux.pdf\
- #ABC1_3a_M1_A1-flux.pdf\
- #ABC1_3a_M3_A1-flux.pdf\
-
-ABC1_1a_M1_A1-flux.pdf: ABC1_1a_M1_A1.nc plot.py
- ./plot.py ABC1_1a_M1_A1.nc
-
-ABC1_1a_M3_A1-flux.pdf: ABC1_1a_M3_A1.nc plot.py
- ./plot.py ABC1_1a_M3_A1.nc
-
-ABC1_3a_M1_A1-flux.pdf: ABC1_3a_M1_A1.nc plot.py
- ./plot.py ABC1_3a_M1_A1.nc
-
-ABC1_3a_M3_A1-flux.pdf: ABC1_3a_M3_A1.nc plot.py
- ./plot.py ABC1_3a_M3_A1.nc
-
-ex_1a_M1_A1_1d-profile.pdf: ex_1a_M1_A1_1d.nc
- ./plot-time-series.py ex_1a_M1_A1_1d.nc
-
-ex_1a_M1_A1_1d.nc: ex_ABC1_1a_M1_A1.nc
- flowline.py -o $@ --collapse -d y ex_ABC1_1a_M1_A1.nc
-
-ABC1_1a_M1_A1.nc: experiment-1a-mode-1.sh
- time sh experiment-1a-mode-1.sh ${NPROC} | tee out.1a-mode-1
-
-ABC1_1a_M3_A1.nc: experiment-1a-mode-3.sh
- time sh experiment-1a-mode-3.sh ${NPROC} | tee out.1a-mode-3
-
-ABC1_3a_M1_A1.nc: experiment-3a-mode-1.sh
- time sh experiment-3a-mode-1.sh ${NPROC} | tee out.3a-mode-1
-
-ABC1_3a_M3_A1.nc: experiment-3a-mode-3.sh
- time sh experiment-3a-mode-3.sh ${NPROC} | tee out.3a-mode-3
-
-experiment-3a-mode-1.sh: run.py ${SLSERIES}
- ./run.py -e 3a --mode=1 > $@
-
-experiment-3a-mode-3.sh: run.py ${SLSERIES}
- ./run.py -e 3a --mode=3 --Mx=1201 > $@
-
-experiment-1a-mode-1.sh: run.py ${SLSERIES}
- ./run.py -e 1a --mode=1 > $@
-
-experiment-1a-mode-3.sh: run.py ${SLSERIES}
- ./run.py -e 1a --mode=3 --Mx=1201 > $@
+ ex_deltaSL-linear_1d-profile.pdf\
+
+ex_deltaSL-linear_1d-profile.pdf: ex_deltaSL-linear_1d.nc
+ ./plot-time-series.py ex_deltaSL-linear_1d.nc
+
+ex_deltaSL-linear_1d.nc: deltaSL-linear.nc
+ flowline.py -o $@ --collapse -d y ex_deltaSL-linear.nc
+
+deltaSL-linear.nc: init-linear.nc ${SLSERIES}
+ mpiexec -n ${NPROC} \
+ pismr -i init-linear.nc\
+ -options_left\
+ -sea_level constant,delta_sl -ocean_delta_sl_file ${SLSERIES}\
+ -extra_file ex_$@\
+ -extra_times 0:50:${T_END}\
+ -extra_vars thk,topg,effbwp,velbar_mag,flux_mag,mask,dHdt,usurf,hardav,velbase_mag,nuH,tauc,taud,taub,flux_divergence,cell_grounded_fraction\
+ -ts_file ts_$@\
+ -ts_times 0:50:${T_END}\
+ -ys 0 -ye ${T_END}\
+ -backup_size big\
+ -o $@\
+ -o_order zyx
${SLSERIES}: sealvl.py Makefile
- printf '0\t0\n30000\t0\n' | ./sealvl.py
+ printf '0\t0\n${T_END}\t200\n' | ./sealvl.py
+
+init-linear.nc: MISMIP_boot_1a_M1_A1.nc
+ mpiexec -n ${NPROC} \
+ pismr -i MISMIP_boot_1a_M1_A1.nc\
+ -bootstrap\
+ -Mx 301 -My 3 -Mz 15\
+ -Lz 5000\
+ -stress_balance ssa+sia\
+ -gradient eta\
+ -front_retreat_file MISMIP_boot_1a_M1_A1.nc\
+ -cfbc -part_grid -kill_icebergs -subgl\
+ -ys 0 -ye 10e3\
+ -options_left\
+ -plastic_phi 24\
+ -till_cohesion 0\
+ -hydrology routing\
+ -geothermal_flux 70e-3\
+ -sia_e 4.5 -ssa_e 0.512\
+ -bed_def lc\
+ -stress_balance.sia.max_diffusivity 1e4\
+ -extra_file ex_$@\
+ -extra_times 0:50:3e4\
+ -extra_vars thk,topg,effbwp,velbar_mag,flux_mag,mask,dHdt,usurf,hardav,velbase_mag,nuH,tauc,taud,taub,flux_divergence,cell_grounded_fraction\
+ -ts_file ts_$@\
+ -ts_times 0:50:3e4\
+ -backup_size big\
+ -o $@\
+ -o_order zyx
+
+MISMIP_boot_1a_M1_A1.nc: prepare.py MISMIP.py
+ ./prepare.py -o $@ -e 1a -m 1
clean:
- rm -f out.* experiment-*.sh *.nc *.pdf ${SLSERIES}
+ rm -f *.nc *.nc~ *.pdf ${SLSERIES} SSAFD_kspdivergederror.petsc*
.PHONY: all clean
\ No newline at end of file
diff --git a/plot-time-series.py b/plot-time-series.py
@@ -86,24 +86,24 @@ def plot_profile(in_file, out_file):
# mask out ice-free areas
usurf = np.ma.array(usurf, mask=mask == 4)
- # compute the lower surface elevation
- lsrf = topg.copy()
- lsrf[mask == 3] = -MISMIP.rho_i() / MISMIP.rho_w() * thk[mask == 3]
- lsrf = np.ma.array(lsrf, mask=mask == 4)
-
# convert x to kilometers
x /= 1e3
figure(1)
ax = subplot(111)
- for i in range(0, thk.shape[0]):
+ for i in range(1, thk.shape[0]):
+ # compute the lower surface elevation
+ lsrf = topg[i].copy()
+ thk_ = thk[i].copy()
+ lsrf[mask[i] == 3] = -MISMIP.rho_i() / MISMIP.rho_w() * thk_[mask[i] == 3]
+ lsrf = np.ma.array(lsrf, mask=mask[i] == 4)
# modeled grounding line position
- #xg_PISM = find_grounding_line(x, topg[i], thk[i], mask[i])
+ #xg_PISM = find_grounding_line(x, lsrf, thk_, mask[i])
#plot(x, np.zeros_like(x), ls='dotted', color='red')
icecolor = cm.cool(i / thk.shape[0])
plot(x, topg[i], color='black')
plot(x, usurf[i], color=icecolor, label='{}'.format(i))
- plot(x, lsrf[i], color=icecolor)
+ plot(x, lsrf, color=icecolor)
xlabel('distance from the divide, km')
ylabel('elevation, m')
diff --git a/run.py b/run.py
@@ -1,383 +0,0 @@
-#!/usr/bin/env python3
-import MISMIP
-
-# This scripts generates bash scripts that run MISMIP experiments and generates
-# all the necessary input files.
-#
-# Run run.py > my_new_mismip.sh and use that.
-
-try:
- from netCDF4 import Dataset as NC
-except:
- print("netCDF4 is not installed!")
- sys.exit(1)
-
-import sys
-
-# The "standard" preamble used in many PISM scripts:
-preamble = '''
-if [ -n "${SCRIPTNAME:+1}" ] ; then
- echo "[SCRIPTNAME=$SCRIPTNAME (already set)]"
- echo ""
-else
- SCRIPTNAME="#(mismip.sh)"
-fi
-
-echo
-echo "# =================================================================================="
-echo "# MISMIP experiments"
-echo "# =================================================================================="
-echo
-
-set -e # exit on error
-
-NN=2 # default number of processors
-if [ $# -gt 0 ] ; then # if user says "mismip.sh 8" then NN = 8
- NN="$1"
-fi
-
-echo "$SCRIPTNAME NN = $NN"
-
-# set MPIDO if using different MPI execution command, for example:
-# $ export PISM_MPIDO="aprun -n "
-if [ -n "${PISM_MPIDO:+1}" ] ; then # check if env var is already set
- echo "$SCRIPTNAME PISM_MPIDO = $PISM_MPIDO (already set)"
-else
- PISM_MPIDO="mpiexec -n "
- echo "$SCRIPTNAME PISM_MPIDO = $PISM_MPIDO"
-fi
-
-# check if env var PISM_DO was set (i.e. PISM_DO=echo for a 'dry' run)
-if [ -n "${PISM_DO:+1}" ] ; then # check if env var DO is already set
- echo "$SCRIPTNAME PISM_DO = $PISM_DO (already set)"
-else
- PISM_DO=""
-fi
-
-# prefix to pism (not to executables)
-if [ -n "${PISM_BIN:+1}" ] ; then # check if env var is already set
- echo "$SCRIPTNAME PISM_BIN = $PISM_BIN (already set)"
-else
- PISM_BIN="" # just a guess
- echo "$SCRIPTNAME PISM_BIN = $PISM_BIN"
-fi
-'''
-
-
-class Experiment:
-
- "A MISMIP experiment."
- experiment = ""
- mode = 1
- model = 1
- semianalytic = True
- Mx = 151
- My = 3
- Mz = 15
- initials = "ABC"
- executable = "$PISM_DO $PISM_MPIDO $NN ${PISM_BIN}pismr"
-
- def __init__(self, experiment, model=1, mode=1, Mx=None, Mz=15, semianalytic=True,
- initials="ABC", executable=None):
- self.model = model
- self.mode = mode
- self.experiment = experiment
- self.initials = initials
- self.semianalytic = semianalytic
-
- if executable:
- self.executable = executable
-
- if mode == 3:
- self.Mx = Mx
- else:
- self.Mx = 2 * MISMIP.N(self.mode) + 1
-
- self.My = 3
-
- if self.experiment == "2b":
- self.Lz = 7000
- else:
- self.Lz = 6000
- self.Lz *= 2.0
-
- def physics_options(self, input_file, step):
- "Options corresponding to modeling choices."
- config_filename = self.config(step)
-
- #options = ["-energy none", # isothermal setup; allows selecting cold-mode flow laws
- #"-ssa_flow_law isothermal_glen", # isothermal setup
- #"-yield_stress constant",
- #"-tauc %e" % MISMIP.C(self.experiment),
- #"-pseudo_plastic",
- #"-gradient eta",
- #"-pseudo_plastic_q %e" % MISMIP.m(self.experiment),
- #"-pseudo_plastic_uthreshold %e" % MISMIP.secpera(),
- #"-front_retreat_file %s" % input_file, # prescribe the maximum ice extent
- #"-config_override %s" % config_filename,
- #"-ssa_method fd",
- #"-cfbc", # calving front boundary conditions
- #"-part_grid", # sub-grid front motion parameterization
- #"-ssafd_ksp_rtol 1e-7",
- #"-ys 0",
- #"-ye %d" % MISMIP.run_length(self.experiment, step),
- #"-options_left",
- #]
-
-
- options = ["-stress_balance ssa+sia",
- #"-ssa_flow_law isothermal_glen", # isothermal setup
- #"-yield_stress constant",
- #"-tauc %e" % MISMIP.C(self.experiment),
- #"-pseudo_plastic",
- "-gradient eta",
- #"-pseudo_plastic_q %e" % MISMIP.m(self.experiment),
- #"-pseudo_plastic_uthreshold %e" % MISMIP.secpera(),
- #"-yield_stress mohr_coulomb",
- #"-till_flux",
- "-front_retreat_file %s" % input_file, # prescribe the maximum ice extent
- "-config_override %s" % config_filename,
- #"-ssa_method fd",
- "-cfbc", # PIK: calving front boundary conditions
- "-part_grid", # PIK: sub-grid front motion parameterization
- "-kill_icebergs", # PIK: https://pism-docs.org/sphinx/manual/modeling-choices/marine/pik.html#iceberg-removal
- "-subgl", # PIK: https://pism-docs.org/sphinx/manual/modeling-choices/marine/pik.html#sub-grid-treatment-of-the-grounding-line-position
- #"-ssafd_ksp_rtol 1e-7",
- "-ys 0",
- "-ye 1e3",
- #"-ye %d" % MISMIP.run_length(self.experiment, step),
- "-options_left",
- "-hydrology routing",
- "-plastic_phi 24", # Tulaczyk et al 2000
- "-till_cohesion 0", # Tulaczyk et al 2000
- "-geothermal_flux 70e-3", # Feldmann and Levermann 2017
- "-sia_e 4.5", # Martin et al 2010
- "-ssa_e 0.512", # Martin et al 2010
- "-bed_def lc", # Lingle and Clark 1985, Bueler et al 2007
- "-stress_balance.sia.max_diffusivity 1e4",
- "-sea_level constant,delta_sl -ocean_delta_sl_file sealvl.nc",
- ]
-
- #options = [,
- #"-yield_stress mohr_coulomb",
- #"-Lz 1e4",
- ##"-till_flux",
- #"-cfbc", # calving front boundary conditions
- #"-part_grid", # sub-grid front motion parameterization
- #"-front_retreat_file %s" % input_file, # prescribe the maximum ice extent
- #"-ys 0",
- #"-ye %d" % MISMIP.run_length(self.experiment, step),
- #"-options_left",
- #]
-
- return options
-
- def config(self, step):
- '''Generates a config file containing flags and parameters
- for a particular experiment and step.
-
- This takes care of flags and parameters that *cannot* be set using
- command-line options. (We try to use command-line options whenever we can.)
- '''
-
- filename = "MISMIP_conf_%s_A%s.nc" % (self.experiment, step)
-
- nc = NC(filename, 'w', format="NETCDF3_CLASSIC")
-
- var = nc.createVariable("pism_overrides", 'i')
-
- attrs = {"geometry.update.use_basal_melt_rate": "no",
- "stress_balance.ssa.compute_surface_gradient_inward": "no",
- "flow_law.isothermal_Glen.ice_softness": MISMIP.A(self.experiment, step),
- "constants.ice.density": MISMIP.rho_i(),
- "constants.sea_water.density": MISMIP.rho_w(),
- "bootstrapping.defaults.geothermal_flux": 0.0,
- "stress_balance.ssa.Glen_exponent": MISMIP.n(),
- "constants.standard_gravity": MISMIP.g(),
- "ocean.sub_shelf_heat_flux_into_ice": 0.0,
- }
-
- if self.model != 1:
- attrs["stress_balance.sia.bed_smoother.range"] = 0.0
-
- for name, value in attrs.items():
- var.setncattr(name, value)
-
- nc.close()
-
- return filename
-
- def bootstrap_options(self, step):
- boot_filename = "MISMIP_boot_%s_M%s_A%s.nc" % (self.experiment, self.mode, step)
-
- import prepare
- prepare.pism_bootstrap_file(boot_filename, self.experiment, step, self.mode, N=self.Mx,
- semianalytical_profile=self.semianalytic)
-
- options = ["-i %s -bootstrap" % boot_filename,
- "-Mx %d" % self.Mx,
- "-My %d" % self.My,
- "-Mz %d" % self.Mz,
- "-Lz %d" % self.Lz]
-
- return options, boot_filename
-
- def output_options(self, step):
- output_file = self.output_filename(self.experiment, step)
- extra_file = "ex_" + output_file
- ts_file = "ts_" + output_file
-
- options = ["-extra_file %s" % extra_file,
- "-extra_times 0:50:3e4",
- "-extra_vars thk,topg,velbar_mag,flux_mag,mask,dHdt,usurf,hardav,velbase_mag,nuH,tauc,taud,taub,flux_divergence,cell_grounded_fraction",
- "-ts_file %s" % ts_file,
- "-ts_times 0:50:3e4",
- "-backup_size big",
- "-o %s" % output_file,
- "-o_order zyx",
- ]
-
- return output_file, options
-
- def output_filename(self, experiment, step):
- return "%s%s_%s_M%s_A%s.nc" % (self.initials, self.model, experiment, self.mode, step)
-
- def options(self, step, input_file=None):
- '''Generates a string of PISM options corresponding to a MISMIP experiment.'''
-
- if input_file is None:
- input_options, input_file = self.bootstrap_options(step)
- else:
- input_options = ["-i %s" % input_file]
-
- physics = self.physics_options(input_file, step)
-
- output_file, output_options = self.output_options(step)
-
- return output_file, (input_options + physics + output_options)
-
- def run_step(self, step, input_file=None):
- out, opts = self.options(step, input_file)
- print('echo "# Step %s-%s"' % (self.experiment, step))
- print("%s %s" % (self.executable, ' '.join(opts)))
- print('echo "Done."\n')
-
- return out
-
- def run(self, step=None):
- print('echo "# Experiment %s"' % self.experiment)
-
- if self.experiment in ('1a', '1b'):
- # bootstrap
- input_file = None
- # steps 1 to 9
- steps = list(range(1, 10))
-
- if self.experiment in ('2a', '2b'):
- # start from step 9 of the corresponding experiment 1
- input_file = self.output_filename(self.experiment.replace("2", "1"), 9)
- # steps 8 to 1
- steps = list(range(8, 0, -1))
-
- if self.experiment == '3a':
- # bootstrap
- input_file = None
- # steps 1 to 13
- steps = list(range(1, 14))
-
- if self.experiment == '3b':
- # bootstrap
- input_file = None
- # steps 1 to 15
- steps = list(range(1, 16))
-
- if step is not None:
- input_file = None
- steps = [step]
-
- for step in steps:
- input_file = self.run_step(step, input_file)
-
-
-def run_mismip(initials, executable, semianalytic):
- Mx = 601
- models = (1, 2)
- modes = (1, 2, 3)
- experiments = ('1a', '1b', '2a', '2b', '3a', '3b')
-
- print(preamble)
-
- for model in models:
- for mode in modes:
- for experiment in experiments:
- e = Experiment(experiment,
- initials=initials,
- executable=executable,
- model=model, mode=mode, Mx=Mx,
- semianalytic=semianalytic)
- e.run()
-
-
-if __name__ == "__main__":
- from optparse import OptionParser
-
- parser = OptionParser()
-
- parser.usage = "%prog [options]"
- parser.description = "Creates a script running MISMIP experiments."
- parser.add_option("--initials", dest="initials", type="string",
- help="Initials (3 letters)", default="ABC")
- parser.add_option("-e", "--experiment", dest="experiment", type="string",
- default='1a',
- help="MISMIP experiments (one of '1a', '1b', '2a', '2b', '3a', '3b')")
- parser.add_option("-s", "--step", dest="step", type="int",
- help="MISMIP step number")
- parser.add_option("-u", "--uniform_thickness",
- action="store_false", dest="semianalytic", default=True,
- help="Use uniform 10 m ice thickness")
- parser.add_option("-a", "--all",
- action="store_true", dest="all", default=False,
- help="Run all experiments")
- parser.add_option("-m", "--mode", dest="mode", type="int", default=1,
- help="MISMIP grid mode")
- parser.add_option("--Mx", dest="Mx", type="int", default=601,
- help="Custom grid size; use with --mode=3")
- parser.add_option("--model", dest="model", type="int", default=1,
- help="Models: 1 - SSA only; 2 - SIA+SSA")
- parser.add_option("--executable", dest="executable", type="string",
- help="Executable to run, e.g. 'mpiexec -n 4 pismr'")
-
- (opts, args) = parser.parse_args()
-
- if opts.all:
- run_mismip(opts.initials, opts.executable, opts.semianalytic)
- exit(0)
-
- def escape(arg):
- if arg.find(" ") >= 0:
- parts = arg.split("=")
- return "%s=\"%s\"" % (parts[0], ' '.join(parts[1:]))
- else:
- return arg
-
- arg_list = [escape(a) for a in sys.argv]
-
- print("#!/bin/bash")
- print("# This script was created by examples/mismip/run.py. The command was:")
- print("# %s" % (' '.join(arg_list)))
-
- if opts.executable is None:
- print(preamble)
-
- e = Experiment(opts.experiment,
- initials=opts.initials,
- executable=opts.executable,
- model=opts.model,
- mode=opts.mode,
- Mx=opts.Mx,
- semianalytic=opts.semianalytic)
-
- if opts.step is not None:
- e.run(opts.step)
- else:
- e.run()
