Debugged and extended the Dukowicz slab test case

This commit modifies the run and plot scripts for the Dukowicz slab test case,
as described in Section 5 of this paper:

   J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a
   more efficient computational solution. The Cryosphere, 6, 21-34,
   https://doi.org/10.5194/tc-6-21-2012.

The test case consists of an ice slab of uniform thickness moving down an
inclined plane by a combination of sliding and shearing.
Analytic Stokes and first-order velocity solutions exist for all values of Glen's exponent n >= 1.
The solutions for n = 1 are derived in Dukowicz (2012), and solutions for n > 1
are derived in an unpublished manuscript by Dukowicz (2013).
These solutions can be compared to those simulated by CISM.

The original scripts, runSlab.py and plotSlab.py, were written by Matt Hoffman
with support for n = 1.  They came with warnings that the test is not supported.

The test is now supported, and the scripts include some new features:
* The user may specify any value of n >= 1 (not necessarily an integer).
  The tests assume which_ho_efvs = 2 (nonlinear viscosity) with flow_law = 0 (constant).
* Physics parameters are no longer hard-coded.  The user can enter the ice thickness,
  beta, viscosity coefficient (mu_n), and slope angle (theta) on the command line.
* The user can specify time parameters dt (the dynamic time step) and nt (number of steps).
  The previous version did not support transient runs.
* The user can specify a small thickness perturbation dh, which is added to the initial
  uniform thickness via random sampling from a Gaussian distribution.
  The perturbation will grow or decay, depending on the solver stability for given dx and dt.

For n = 1, the viscosity coefficient mu_1 has a default value of 1.e6 Pa yr in the relation
mu = mu_1 * eps((1-n)/n), where eps is the effective strain rate.

For n > 1, the user can specify a coefficient mu_n; otherwise the run script computes mu_n
such that the basal and surface speeds are nearly the same as for an n = 1 case with the
mu_1 = 1.e6 Pa yr and the same values of thickness, beta, and theta.

(Note: There is a subtle difference between the Dukowicz and CISM definitions of the
effective strain rate; the Dukowicz value is twice as large. Later, it might be helpful
to make the Dukowicz convention consistent with CISM.)

I modified the plotting script, plotSlab.py, to look in the config and output files
for physics parameters that are no longer hardwired.
I slightly modified the analytic formulas to give the correct solution for non-integer n.

This script creates two plots.  The first plot shows excellent agreement between higher-order
CISM solutions and the analytic solution for small values of the slope angle theta.
For steep slopes, the answers diverge as expected.

For the second plot, the extent of the y-axis is wrong. This remains to be fixed.

I also added a new script, stabilitySlab.py, to carry out stability tests as described in:

     Robinson, A., D. Goldberg, and W. H. Lipscomb, A comparison of the performance
     of depth-integrated ice-dynamics solvers, to be submitted to The Cryosphere.

The idea is that for a given set of physics parameters and stress-balance approximation
(DIVA, L1L2, etc.), the script launches multiple CISM runs at a range of grid resolutions.
At each grid resolution, the script determines the maximum stable time step.
A run is deemed stable when the standard deviation of an initial small thickness perturbation
is reduced over the course of 100 time steps.  A run is unstable if the standard deviation
increases or if the model aborts (usually with a CFL violation).

I have run the stability script for several solvers (DIVA, L1L2, SIA, SSA) for each of
two physical cases: one with thick shearing ice and one with thin sliding ice.
Each suite of experiments runs in a few minutes on 4 Macbook cores for solvers other than BP.
As expected, DIVA and SSA are much more stable than L1L2 and SIA.

This and the previous commit correspond to the CISM code and scripts used for
the initial submission by Robinson et al. (2021).

tests/slab/plotSlab.py

@@ -1,24 +1,22 @@
 #!/usr/bin/env python2
 
-
 """
 This script plots the results of an experiment with an ice "slab" on an
 inclined plane. Test case is described in sections 5.1-2 of:
-    J.K. Dukoqicz, 2012. Reformulating the full-Stokes ice sheet model for a
+    J.K. Dukowicz, 2012. Reformulating the full-Stokes ice sheet model for a
     more efficient computational solution. The Cryosphere, 6, 21-34.
     www.the-cryosphere.net/6/21/2012/
 
 Blatter-Pattyn First-order solution is described in J.K. Dukowicz, manuscript
 in preparation.
 """
 #FIXME: Manuscript likely not in prep anymore -- JHK, 08/07/2015
+#       Not published as of July 2021 -- WHL
 
 # Written by Matt Hoffman, Dec. 16, 2013
 # Reconfigured by Joseph H Kennedy at ORNL on August 7, 2015 to work with the regression testing
 #     NOTE: Did not adjust inner workings except where needed.
-
-
-#NOTE: this script is assuming n=3, but more general solutions are available.
+# Revised by William Lipscomb in 2021 to support more options, including general values of Glen's n.
 
 import os
 import sys
@@ -28,8 +26,12 @@
 import matplotlib.pyplot as plt
 
 from netCDF import *
-from math import tan, pi, sin, cos
-from runSlab import n, rhoi, grav, theta, beta, efvs, thickness  # Get the values used to run the experiment
+from math import tan, pi, sin, cos, atan
+
+# Get hard-coded parameters from the run script.
+from runSlab import rhoi, grav
+
+import ConfigParser
 
 import argparse
 parser = argparse.ArgumentParser(description=__doc__,
@@ -46,16 +48,15 @@
         help="The tests output file you would like to plot. If a path is" \
         +"passed via this option, the -o/--output-dir option will be ignored.")
 
+parser.add_argument('-c','--config-file',
+        help="The configure file used to set up the test case and run CISM")
+
 
 # ===========================================================
 # Define some variables and functions used in the main script
 # ===========================================================
 
-# Calculate scales from Ducowicz unpub. man.
-eta = beta * thickness * efvs**-n * (rhoi * grav * thickness)**(n-1) 
-velscale = (rhoi * grav * thickness / efvs)**n * thickness
-thetar = theta * pi/180.0  # theta in radians
-
+#WHL args.output-file with a hyphen?
 def get_in_file():
     if args.output_file:
         out_d, out_f = os.path.split(args.output_file)
@@ -76,7 +77,7 @@ def get_in_file():
             newest = max(matching, key=os.path.getmtime)
             print("\nWARNING: MULTIPLE *.out.nc FILES DETECTED!")
             print(  "==========================================")
-            print(  "Ploting the most recently modified file in the output directory:")
+            print(  "Plotting the most recently modified file in the output directory:")
             print(  "    "+newest)
             print(  "To plot another file, specify it with the -f/--outfile option.\n")
 
@@ -94,6 +95,25 @@ def get_in_file():
 
     return filein
 
+def split_file_name(file_name):
+    """
+    Get the root name, size, and number of processors from an out.nc filename.
+    #WHL - Adapted from plotISMIP_HOM.py
+    """
+    root = ''
+    size = ''
+    proc = ''
+
+    file_details = file_name.replace('.out.nc','') .split('.')
+#    print(file_details)
+#    print('len = ' + str(len(file_details)))
+
+    if len(file_details) > 2:
+        proc = '.'+file_details[2]
+    size = '.'+file_details[1]
+    root = file_details[0]
+
+    return (root, size, proc)
 
 # =========================
 # Actual script starts here
@@ -103,10 +123,7 @@ def main():
     Plot the slab test results.
     """
 
-    print("WARNING: THIS TEST CASE IS IN DEVELOPMENT. USE AT YOUR OWN RISK!")
-
-
-    filein = get_in_file()    
+    filein = get_in_file()
 
     # Get needed variables from the output file
     x1 = filein.variables['x1'][:]
@@ -120,28 +137,96 @@ def main():
     # use integer floor division operator to get an index close to the center 
     xp = len(x0)//2
     yp = len(y0)//2
-    #yp = 15
-    #xp = 15
     # =====================================================================
-    print 'Using x index of '+str(xp)+'='+str(x0[xp])
-    print 'Using y index of '+str(yp)+'='+str(y0[yp])
+
+    print('Using x index of '+str(xp)+'='+str(x0[xp]))
+    print('Using y index of '+str(yp)+'='+str(y0[yp]))
 
     thk = filein.variables['thk'][:]
     if netCDF_module == 'Scientific.IO.NetCDF':
-       thk = thk * filein.variables['thk'].scale_factor
+        thk = thk * filein.variables['thk'].scale_factor
     topg = filein.variables['topg'][:]
     if netCDF_module == 'Scientific.IO.NetCDF':
-       topg = topg * filein.variables['topg'].scale_factor
+        topg = topg * filein.variables['topg'].scale_factor
     uvel = filein.variables['uvel'][:]
     if netCDF_module == 'Scientific.IO.NetCDF':
-       uvel = uvel * filein.variables['uvel'].scale_factor
-
+        uvel = uvel * filein.variables['uvel'].scale_factor
+    beta_2d = filein.variables['beta'][:]
+    if netCDF_module == 'Scientific.IO.NetCDF':
+        beta_2d = beta_2d * filein.variables['beta'].scale_factor
+
+    # Get the name of the config file
+    # If not entered on the command line, then construct from the output filename
+
+    if not args.config_file:
+        root, size, proc = split_file_name(args.output_file)
+        args.config_file = root + size + proc + '.config'
+
+    configpath = os.path.join(args.output_dir, args.config_file)
+    print('configpath = ' + configpath)
+
+    # Get gn and default_flwa from the config file
+
+    try:
+        config_parser = ConfigParser.SafeConfigParser()
+        config_parser.read( configpath )
+
+        gn = float(config_parser.get('parameters','n_glen'))
+        flwa = float(config_parser.get('parameters', 'default_flwa'))
+
+    except ConfigParser.Error as error:
+        print("Error parsing " + args.config )
+        print("   "),
+        print(error)
+        sys.exit(1)
+
+    # Derive the viscosity constant mu_n from flwa
+    # This expression is derived in the comments on flwa in runSlab.py.
+    mu_n = 1.0 / (2.0**((1.0+gn)/(2.0*gn)) * flwa**(1.0/gn))
+
+    # Get the ice thickness from the output file.
+    # If thickness = constant (i.e., the optional perturbation dh = 0), it does not matter where we sample.
+    # Note: In general, this thickness will differ from the baseline 'thk' that is used in runSlab.py
+    #        to create the input file.
+    #       This is because the baseline value is measured perpendicular to the sloped bed,
+    #        whereas the CISM value is in the vertical direction, which is not perpendicular to the bed.
+    thickness = thk[0,yp,xp]
+
+    # Get beta from the output file.
+    # Since beta = constant, it does not matter where we sample.
+    beta = beta_2d[0,yp,xp]
+
+    # Derive theta from the output file as atan(slope(topg))
+    # Since the slope is constant, it does not matter where we sample.
+    slope = (topg[0,yp,xp] - topg[0,yp,xp+1]) / (x0[xp+1] - x0[xp])
+    thetar = atan(slope)
+    theta = thetar * 180.0/pi
+
+    # Compute the dimensionless parameter eta and the velocity scale,
+    # which appear in the scaled velocity solution.
+    eta = (beta * thickness / mu_n**gn) * (rhoi * grav * thickness)**(gn-1)
+    velscale = (rhoi * grav * thickness / mu_n)**gn * thickness
+
+    print('gn   = ' + str(gn))
+    print('rhoi = ' + str(rhoi))
+    print('grav = ' + str(grav))
+    print('thck = ' + str(thickness))
+    print('mu_n = ' + str(mu_n))
+    print('flwa = ' + str(flwa))
+    print('beta = ' + str(beta))
+    print('eta  = '  + str(eta))
+    print('theta= ' + str(theta))
+    print('velscale = ' + str(velscale))
 
     # === Plot the results at the given location ===
     # Note we are not plotting like in Fig 3 of paper.
     # That figure plotted a profile against zprime.
     # It seemed more accurate to plot a profile against z to avoid interpolating model results (analytic solution can be calculated anywhere).
-    # Also, the analytic solution calculates the bed-parallel u velocity, but CISM calculates u as parallel to the geoid, so we need to transform the analytic solution to the CISM coordinate system.
+    # Also, the analytic solution calculates the bed-parallel u velocity, but CISM calculates u as parallel to the geoid,
+    #  so we need to transform the analytic solution to the CISM coordinate system.
+
+    #WHL - I think the analytic solution is actually for u(z'), which is not bed-parallel.
+    #      The bed-parallel solution would be u'(z'), with w'(z') = 0.
 
     fig = plt.figure(1, facecolor='w', figsize=(12, 6))
 
@@ -151,24 +236,23 @@ def main():
     x = (x0-x0[xp]) / thickness
     # calculate rotated zprime coordinates for this column (we assume the solution truly is spatially uniform)
     zprime = x[xp] * sin(thetar) + z * cos(thetar)
-    #print 'zprime', zprime
 
     # Calculate analytic solution for x-component of velocity (eq. 39 in paper) for the CISM-column
-    #uvelStokesAnalyticScaled =  sin(theta * pi/180.0) * cos(theta * pi/180.0) * (0.5 * zprime**2 - zprime - 1.0/eta)
-    uvelStokesAnalyticScaled = (-1)**n * 2**((1.0-n)/2.0) * sin(thetar)**n * cos(thetar) / (n+1) \
-                    * ( (zprime - 1.0)**(n+1) - (-1.0)**(n+1) ) + sin(thetar) * cos(thetar) / eta
+    uvelStokesAnalyticScaled = sin(thetar) * cos(thetar) / eta   \
+        - 2**((1.0-gn)/2.0) * sin(thetar)**gn * cos(thetar) / (gn+1) * ( (1.0 - zprime)**(gn+1) - 1.0 )
 
-    # Calculate the BP FO solution for x-component of velocity (Ducowicz, in prep. paper, Eq.30, n=3)
-    #uvelFOAnalyticScaled = (tan(theta * pi/180.0))**3 / (8.0 * (1.0 + 3.0 * (sin(theta * pi/180.0)**2))**2) \
-    uvelFOAnalyticScaled = (-1)**n * 2**((1.0-n)/2.0) * tan(thetar)**n /  \
-                           ( (n + 1) * (1.0 + 3.0 * sin(thetar)**2)**((n+1.0)/2.0) )  \
-                           * ( (zprime - 1.0)**(n+1) - (-1.0)**(n+1) ) + tan(thetar) / eta
+    # Calculate the BP FO solution for x-component of velocity (Dukowicz, in prep. paper, Eq.30, n=3)
+    uvelFOAnalyticScaled = + tan(thetar) / eta  \
+                           -  2**((1.0-gn)/2.0) * tan(thetar)**gn /  \
+                           ( (gn + 1) * (1.0 + 3.0 * sin(thetar)**2)**((gn+1.0)/2.0) )  \
+                           * ( (1.0 - zprime)**(gn+1) - 1.0 )
 
     ### 1. Plot as nondimensional variables
     # Plot analytic solution 
     fig.add_subplot(1,2,1)
     plt.plot(uvelStokesAnalyticScaled, z, '-kx', label='Analytic Stokes')
     plt.plot(uvelFOAnalyticScaled, z, '-ko', label='Analytic FO')
+
     # Plot model results
     plt.plot(uvel[0,:,yp,xp] / velscale, z, '--ro', label='CISM') 
     plt.ylim((-0.05, 1.05))
@@ -191,7 +275,16 @@ def main():
     plt.title('Velocity profile at x=' + str(x0[xp]) + ' m, y=' + str(y0[yp]) + ' m\n(Unscaled coordinates)')
 
     #################
+#    print('y0_min:')
+#    print(y0.min())
+#    print('y0_max:')
+#    print(y0.max())
+
     # Now plot maps to show if the velocities vary over the domain (they should not)
+    # For some reason, the y-axis has a greater extent than the range (y0.min, y0.max).
+    #TODO - Fix the y-axis extent.  Currently, the extent is too large for small values of ny.
+    #TODO - Plot the thickness relative to the initial thickness.
+
     fig = plt.figure(2, facecolor='w', figsize=(12, 6))
     fig.add_subplot(1,2,1)
     uvelDiff = uvel[0,0,:,:] - uvel[0,0,yp,xp]
@@ -224,14 +317,11 @@ def main():
     #plt.plot(level, tan(thetar)**3 / (8.0 * (1.0 + 3.0 * sin(thetar)**2)**2) * (1.0 - (level-1.0)**4 ) + tan(thetar)/eta, 'b--' , label='nonlinear fo')
     #plt.ylim((0.0, 0.04)); plt.xlabel("z'"); plt.ylabel('u'); plt.legend()
 
-
     plt.draw()
     plt.show()
 
     filein.close()
 
-    print("WARNING: THIS TEST CASE IS IN DEVELOPMENT. USE AT YOUR OWN RISK!")
-
 # Run only if this is being run as a script.
 if __name__=='__main__':
 
@@ -240,4 +330,3 @@ def main():
 
     # run the script
     sys.exit(main())
-