Updated the Dukowicz slab test case #37
Open
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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.
whlipscomb
force-pushed
the
lipscomb/update_slab
branch
from
f91e55c
to
bd9b1f4
Compare
|
@Katetc, this PR can be removed, since the slab test changes were brought into main as part of the basal_physics4 branch. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This commit modifies the run and plot scripts for the Dukowicz slab test case, as described in Section 5 of Dukowicz (2012, The Cryosphere). 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 n >= 1, and can be compared to those simulated by CISM. The original scripts contained warnings that the test was not fully supported. With these changes, the test is supported for CISM users.
This commit also adds a new script to run the solver stability tests described in Section 3.4 of Robinson et al. (2021, submitted to The Cryosphere), and updates the README file.
In most cases, answers are unchanged. To support the slab test, Glen's n is now a user-configurable variable. This might lead to roundoff-level changes on some machines.