Updated the Dukowicz slab test case #37

Until now, the exponent n in the Glen flow law has been set in glimmer_physcon.F90 as an integer parameter, gn = 3. With this commit, n is renamed n_glen (a more greppable name) for use in Glissade. It is declared in glimmer_physcon.F90 as real(dp) with default value 3.0d0. Since n_glen is no longer a parameter, it can now be set in the config file like other physical 'constants' (e.g., rhoi and rhoo) that are not truly constant, but can take different values in different models or benchmarking experiments. To avoid changing answers in old Glide code, I retained the integer parameter gn = 3 in glimmer_physcon.F90. This parameter is now used only in the Glide dycore (e.g., glide_velo.F90). In Glissade, I replaced gn with n_glen in several places: (1) In subroutine glissade_interior_dissipation_sia, the dissipation factor includes n_glen. Note: n_glen does not appear explicitly in the 1st-order dissipation, which is proportional to tau_eff^2/efvs. (2) In glissade_velo_sia, n_glen appears in the vertical integral for the velocity. (3) In glissade_velo_higher, flwafact = flwa**(-1./n_glen) where flwa = A. (4) In glissade_velo_higher, the exponent p_effstr = (1.d0 - n_glen)/n_glen is used to compute effective_viscosity for BP, DIVA, or L1L2. (5) In glissade_velo_higher, subroutine compute_3d_velocity_l1l2 has a factor proportional to tau**((n_glen-1.)/2.) in the vertical integral. For (1) and (2), n_glen was previously assumed to be an integer. Switching it to real(dp) is answer-changing at the machine roundoff level for runs using the local SIA solver in glissade_velo_sia.F90. For (3), (4) and (5), n_glen replaces the equivalent expression real(gn,dp). For this reason, answers are BFB when using the BP, DIVA or L1L2 solver. In glissade_basal_traction, n_glen appeared in the expression for beta in the Coulomb sliding option, HO_BABC_COULOMB_FRICTION. Here, I replaced n_glen with powerlaw_m (also = 3.0d0 by default) to be consistent with the expressions for beta in the Schoof and Tsai laws. In glimmer_scales, Glen's n appears in the expressions for the scaling parameters vis0 and vis_scale, Here, I defined a local integer parameter gn = 3 so that the scales are unchanged. It is now possible for the user to specify arbitrary values of n_glen in tests such as the slab test. Another minor change: I added some logic to the subroutine that computes L1L2 velocities. For which_ho_efvs = 2 = HO_EFVS_NONLINEAR, the effective viscosity (efvs) is computed from the effective strain rate using an equation from Perego et al. (2012). But for option 0 (efvs = constant) and option 1 (efvs = multiple of flow factor), this strain rate equation in the code does not apply. For these options, I added an alternative that computes velocity in terms of the strain-rate-independent efvs. This allows us to use L1L2 for problems with constant efvs (e.g., the slab test).

The code was calling subroutine init_isostasy with isostasy = 0 = ISOSTASY_NONE. This subroutine is now called only if isostasy = 1 = ISOSTASY_COMPUTE.

I modified glissade.F90 to abort cleanly with a call to glide_finalise after an advective CFL error. This is done only when the user does *not* specify adaptive subcycling. The clean abort allows the new slabStability script to keep going, launching a new run with a shorter timestep. In subroutine glissade_flow_factor of glissade_therm.F90, I removed the FLWA_INPUT option (option 3 of whichflwa). This option is redundant with option 0, FLWA_CONST_FLWA, in which the user can set default_flwa in the parameters section of the config file, and otherwise CISM uses default_flwa = 1.0e-16 Pa^-n yr^-1. We probably should rename default_flwa to constant_flwa, but leaving it for now to avoid breaking config files in test cases. Note: This option was added by Matt Hoffman in 2014. I am unaware of test cases that use this option (most have flow_law = 0 or 2), but if there are any, we will need to fix them by switching to whichflwa = 0. In subroutine glissade_therm_driver of glissade_therm.F90, I increased the threshold for column energy conservation errors from 1.0d-8 to 1.0d-7 W/m^2. For very small timesteps of ~1.0e-6 yr, the smaller threshold can be exceeded because of roundoff errors. In subroutine glissade_check_cfl of glissade_transport.F90, I modified the fatal abort for large CFL violations (advective CFL > 10). Now, CISM aborts for CFL > 10 only when adaptive_cfl_threshold > 0, i.e. transport subcycling is enabled. This prevents excessive subcycling for egregious CFL violations. If adaptive_cfl_threshold = 0. (the default), i.e. transport subcycling is not enabled, then the code exits cleanly (with a call to glide_finalise) in glissade.F90 when advective CFL > 1. I added a TODO note to set the default value of geot (the geothermal heat flux) to 0 instead of 0.05 W/m^2. With the default value, simple prognostic tests like the dome are not mass-conserving. Not changing just yet because answers will change for the dome test.

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).

Rewrote the slab README file to describe the new command line options for runSlab.py, and the new script stabilitySlab.py.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Updated the Dukowicz slab test case #37

Updated the Dukowicz slab test case #37

Commits on Sep 18, 2021

Updated the Dukowicz slab test case #37

Are you sure you want to change the base?

Updated the Dukowicz slab test case #37

Commits on Sep 18, 2021