CISM exercise II: run diagnostic test cases
In this exercise, we will run a few of CISM's higher-order test cases, which span a wide range of flow regimes. From the top level directory where you build the code, change to into the higher-order tests directory:
You will see the following list of subdirectories
dome/ # parabolic shaped dome with simple boundary conditions ismip-hom/ # ISMIP-HOM test suite ross/ # Ross ice shelf test case shelf/ # ice shelf test cases on simplified domains
along with a few other files (note that the comments after the hash marks above have been added here). While you may want to look over the tests/higher-order/README file at some point, most of the necessary information from that file is contained on this page. While you are welcome to explore any of the test cases on your own (most of them can be run by simply following the instructions in the README files within each subdirectory), for this exercise we will pick a few representative examples that can be run relatively quickly.
All of the test cases use the simple_glide executable, which is built from the simple_glide.F90 driver in example-drivers/simple_glide/src/. To run any of the tests, you will need to change into the respective subdirectory and either copy the simple_glide executable to that directory
cp ../../../example-drivers/simple_glide/src/simple_gide ./
or make a virtual link to the actual executable file
ln -s ../../../example-drivers/simple_glide/src/simple_glide ./
The Dome test case
This is a very simple test case, simulating the three-dimensional flow field within an isothermal, parabolic shaped dome with no-slip basal boundary conditions and zero flux lateral boundary conditions. You can execute the test with
The call to the python script first builds an input netCDF file in the output/ subdirectory and executes simple_glide. In general, whenever simple_glide is executed, it expects to be followed by the name of a text file with the ".config" extension. If such a file is not specified, you will usually see something like
Enter name of GLIDE configuration file to be read
Here, the python script knows to pass the included "dome.config" script to simple_glide on its own. While there are numerous default settings in the code, in general it will need a configuration file of some sort to specify various things like grid size and spacing, various solver options, boundary conditions, etc. A good way to get a feel for what these options are and what parts of the code they trigger is to look in the ".config" file, find an options (e.g. "evolution = 3"), and then "grep" for that option in the file glide_types.F90 in the libglide/ subdirectory.
As the code runs, you will see some output to the screen like
(dH/dt using incremental remapping) time = 0.0000000 Running Payne/Price higher-order dynamics solver iter # resid (L2 norm) target resid 1 223.257 0.100000E-03 2 223.150 0.100000E-03 3 216.909 0.100000E-03 4 203.843 0.100000E-03 5 180.486 0.100000E-03 6 149.276 0.100000E-03 7 116.333 0.100000E-03 ... 39 0.341224E-03 0.100000E-03 40 0.226718E-03 0.100000E-03 41 0.150639E-03 0.100000E-03 42 0.100090E-03 0.100000E-03 43 0.665039E-04 0.100000E-03
The output you see here is fairly standard. It tells us the following information
- Which solver we are using to evolve the ice thickness (if at all). Here, we see that we are using incremental remapping (which we will discuss further later on). However, looking in the ".config" file we see that the start and end times are identical, so no geometric evolution will take place; we are simply after a diagnostic solution here.
- The current time step we are solving for.
- The dynamics scheme we are using (here, the Blatter-Pattyn equations as formulated and solved by Payne and Price).
- the non-linear iteration number, the residual (the L2 norm of the vector r = Ax - b), and the target residual
Note that when the residual is less than or equal to the target residual, the nonlinear iterations are halted and we have a converged solution. Here it took 43 iterations to arrive at a converged solution.
You can look at the model output using any convenient netCDF file viewer. A python-based netCDF file viewer, viewNetCDF.py is included in top level of the tests/higher-order subdirectory. Another common viewer installed on many machines is NCVIEW. Using the latter, your output for the variable velnorm at level "0" (that is, sigma coordinate level 0, which is the upper surface of the ice sheet), should look something like the following:
Admittedly, this is not a very exciting test case. However, as a sanity check it should confirm whether or not the model is working as expected and shows that, for a very "shallow-ice" like test case, the model indeed reproduces something that looks very much like shallow-ice flow.
The Confined Shelf test case
Now we will go all the way to the other end of the spectrum and demonstrate that the exact same model can also accurately reproduce ice shelf flow. In the follow idealized ice shelf test case, there has been no change at all in the governing equations of the model. The only thing that has changed is the geometry and the boundary conditions of the test problem. Instead of a parabolic dome with no basal slip we now have a flat, floating slab of ice with free-basal slip, zero-flux boundary conditions on three sides, and open-ocean (ice shelf) boundary conditions on the fourth side.
To run the test case, change into the tests/higher-order/shelf/ subdirectory. There are two idealized ice shelf tests cases here, confined-shelf.py and circular-shelf.py. To run the confined shelf test case, proceed with a similar set of steps as when running the dome test case. First, cp or link to the simple_glide executable, then run the test,
As in the previous test case, you should see gradually decreasing residuals as screen output. When the test completes, examine the output with
A color contour plot of CISM output (made in Matlab) is shown below. The black and white contour plot shows output for the same experiment using an SSA model. It is from page 7 of the EISMINT (European Ice Sheet Model InTercomparison) ice shelf intercomparison project documentation.
The ISMIP-HOM test cases
ISMIP-HOM test suite problems. The tests we'll run are for 3d models, so the domain and boundary conditions vary in the x and y directions (i.e. in map plane). For test A, the topography varies periodically in x and y, and for test C, the basal traction varies periodically in x and y. While the amplitude of the variations is the same for all tests, the wavelength is decreased by a factor of two for each successive test. For λ=160 km, the velocities solutions essentially look like that from a shallow ice model. Halving λ to 80 km, then to 40, 20, 10, and finally 5 km, the higher-order components of the stress balance become successively more important to the velocity solution. Figures 1 and 2 below shows relevant input data for each of the two experiments for λ = 80km. Here, in the interest of time, we will only run tests for the first three wavelengths in the series (160, 80, and 40 km).
Figure 1: ISMIP-HOM test A input (periodic basal roughness with no sliding); ice thickness, basal topography, and surface elevation. The basal boundary condition is no slip and the lateral boundary conditions are periodic velocities in x and y.
Figure 2: ISMIP-HOM test C input (sliding according to periodic basal traction); ice thickness,bBetasquared, and surface elevation. Sliding takes place along the basal boundary according to a "betasquared" (traction) type sliding law. The lateral boundary conditions are periodic velocities in x and y (NOTE: In this experiment we have a slab of constant thickness on an inclined plane with only the sliding properties changing along/across the domain. Also, the surface slope is ~10x smaller than for experiment A).
To set up the experiments, we will use some configuration files and python scripts developed by Tim Bocek and Jesse Johnson (also, see this link). These set the correct flags, so that Glimmer/CISM calls the necessary subroutines, and construct the necessary input netCDF files.
First, we need to change into the correct directory where the test scripts and configuration files live. Assuming that you are starting in the directory from the directory glimmer-cism-ho-wo, type
cd tests/ISMIP-HOM/; ls -l
to change into the appropriate directory and list its contents. The files with a .config extension are read by the appropriate python scripts to construct the appropriate fields for the input netCDF file. The .config files are also read by the model at run time, as they specify the values for various flags (including calls to the HO subroutines rather than the shallow ice dynamics routines). The files with a .py extension are the relevant python scripts.
Let's set up test cases A and C for a domain length of 160 km. First check the configuration file to make sure that the domain length, number of grid spaces, and the grid spacing give the correct input values.
emacs ishom.a.config &
[grid] upn = 11 ewn = 51 nsn = 51 dew = 3200 dns = 3200
for the grid variables. Note that 51 x 3200 = 163.2 km, so that our domain will actually have one 3.2 km grid space extra in each of the x and y directions. This is necessary in order to implement the periodic boundary conditions at the lateral domain boundaries. If we now type
python ismip_hom_a.py ishom.a.config
we will generate the necessary netCDF input file for the experiment. Using NCVIEW, you can look at the input data fields and make sure that they are the correct lateral dimensions. To run the model using these input data, we need to execute simple_glide. If that file is not in your path you can copy it into the test directory as follows
cp ../../src/fortran/simple_glide ./
and then execute it by typing
You will be prompted for the relevant configuration file, which is of course ishom.a.config. Another way to do this would have been to use the linux/unix pipe command
echo ishom.a.config | ./simple_glide - or - echo ishom.a.config | simple_glide
(the latter if "simple_glide" is already in your path)
Either way, after responding to the prompt, you should see some model output that looks something like this:
Running Payne/Price higher-order dynamics solver iter # uvel resid vvel resid target resid 2 1.00000 1.00000 0.100000E-04 3 0.143388 0.215689E-02 0.100000E-04 4 0.392197E-01 0.902012E-03 0.100000E-04 5 0.655156E-01 0.745786E-03 0.100000E-04 6 0.503367E-01 0.465037E-03 0.100000E-04 7 0.344782E-02 0.329053E-03 0.100000E-04 8 0.100065E-01 0.256138E-03 0.100000E-04 9 0.163779E-01 0.190435E-03 0.100000E-04 10 0.866196E-02 0.136968E-03 0.100000E-04 11 0.549490E-02 0.104118E-03 0.100000E-04 12 0.546429E-02 0.798739E-04 0.100000E-04
At this point, you know that the higher-order dynamics routine is working on a solution. The 1st column tells you which "outer loop" iteration you are on (that is, iteration on the effective viscosity - the "inner loop" iteration is the conjugate gradient iterative solution to the matrix inversion, and that output is normally suppressed). The 2nd and 3rd columns display the x (uvel) and y (vvel) residuals (the normalized, maximum change in the velocity field between the current and previous iterations) and the last column shows the target residual, at which point the solution is considered to be converged. When the model stops iterating it will create an output netCDF file that we can evaluate. In this case, the file name is ishom.a.out.nc.
- HAVING PROBLEMS?
If you find that something weird is happening like your residuals are all over the place, the model just crashes, etc., trying replacing your configuration files with these default configuration files for ISMIP-HOM test.
Plotting model output
NOTE: This part of the exercise is NOT working at the moment. You can, however, still use NCVIEW to create a profile of the velocity, which you can then compare with the plots shown here.
We can do a quick evaluation of ishom.a.out.nc using NCVIEW or some other netCDF file viewing utility. However, what we really want to do is compare our model solutions with the ISMIP-HOM benchmark solutions, so we can see how are model is doing relative to other models that took place in the benchmark exercise (see Pattyn et. al (2008) for a detailed discussion of the results). To do that, we'll use some of handy python test-suite scripts we mentioned above.
First, we need to generate a text file of output data, from our model result, which will be compared with other model results from the benchmarking study. To do that, we type
./formatData.py a ishom.a.out.nc glm1a160.txt
The script formatData.py reads the netCDF output file and generates the text file glm1a160.txt (the "...a160..." denotes test A for a domain length of 160 km). We then type
./createVisuals.py --exp=a --size=160 -or- ./createVisuals.py -ea -s160
which uses some python modules to make a nice Matlab style figure (see below). That figure will have a ".png" extension and can be found in
Now, do the same set of steps but decrease the domain wavelength to 80 and then 40 km. To make it easy for you, some grid parameters that work for this are
[grid] upn = 11 ewn = 51 nsn = 51 dew = 1600 dns = 1600
for the 80 km test and
[grid] upn = 11 ewn = 51 nsn = 51 dew = 800 dns = 800
for the 40 km test. Output files for these tests will be generated in the same way as above,
./formatData.py a ishom.a.out.nc glm1a080.txt - and - ./formatData.py a ishom.a.out.nc glm1a040.txt
To plot all of the three ouputs for test A on one plot, use
./createVisuals.py --exp=a --size=40,80,160 - or - ./createVisuals.py -ea -s40,80,160
which should give you something that looks similar to Figure 3.
Figure 3: Higher-order model output for ISMIP-HOM test A with domain lengths of 160, 80, and 40 km. Solid black line is output from the current model and the colored, shaded regions represent the standard deviation of other models participating in the benchmarking study (see HERE for a more detailed description of these plots).
Now go through the same set of steps for test case C (again, with wavelengths of 160, 80, and 40 km). You should get a figure that looks something like Figure 4.
Figure 4: Higher-order model output for ISMIP-HOM test C with domain lengths of 160, 80, and 40 km. Solid black line is output from the current model and the colored, shaded regions represent the standard deviation of other models participating in the benchmarking study (see HERE for a more detailed description of these plots).
- Try adjusting the horizontal and vertical grid spacing to see how it affects the results and/or model performance. For example, for the 80km tests, decrease the number of horizontal grid cells by a factor of two and increase the grid spacing by a factor of two,
[grid] upn = 11 ewn = 26 nsn = 26 dew = 3200 dns = 3200
How much faster does the model converge on a solution? Does the output still fall within the standard deviation given by the benchmarks? How What happens if the vertical resolution is doubled?
- Compare higher-order and 0-order solutions for test A with the 80 km domain length. To do this, in ishom.a.config, set the diagnostic_run flag to 0 instead of 1, rebuild the ishom.a.nc file using the python script (as done above), and re-run the model. When the model has finished running, examine ishom.a.out.nc using NCVIEW. Click on the variable uvelhom to make a colormap of the higher-order x component of velocity at time 1 (as shown in figure below).
Figure 4: Using NCVIEW to plot output of "ishom.a.out.nc" to compare higher-order and SIA solutions. Note that the value of current time is 2001, not 2000.
Click somewhere on the image to get a 2d velocity profile (choose x0 under Xaxis). Next, pick the variable uvel (the velocity from the SIA model) and do the same thing. When comparing the two profiles you should see something like in Figure 5.
Figure 5: Comparison of higher-order (top) and SIA (bottom) velocity profiles. For the same model domain, the HO velocities are ~25% slower due to the influence of horizontal-stress gradients, which the SIA model does not "feel" at all.
- Do the same for ISMIP-HOM test A for the 40 km domain. You should notice that, as the magnitude of the higher-order velocities continues to decrease with decreasing domain length, those for the SIA model do not. Why is this?
- Compare the values of the variable vvel (the across-flow velocity calculated from the SIA model) and vvelhom (the across-flow velocity calculated from the higher-order model) at time 200001. Can you explain the differences?