tplace all initialization syntax in Makefile - pism-exp-gsw - ice stream and sediment transport experiments
git clone git://src.adamsgaard.dk/pism-exp-gsw
Log
Files
Refs
README
LICENSE
---
commit 8b985b054a316b0cc9f72041ba2f3fc9ae284f33
parent 83a341c191e31e46e6410801802061491a48f445
Author: Anders Damsgaard 
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
t@@ -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
t@@ -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
t@@ -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()