# Difference between revisions of "ISMIP-HOM test suite exercise"

## Introduction

In this exercise, we will test out Glimmer/CISM's higher-order stress balance subroutines by running the model through a few of the 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.

## Running the test cases

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/src/fortran/, 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.

```vi ishom.a.config
```

gives

```[grid]
upn = 21
ewn = 51
nsn = 51
dew = 3200
dns = 3200

```

for the grid variables. Note that 51 x 3200 = 16.32 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 in your path, you can simple type the name to execute it. If not, you can copy it into the test directory as follows

```cp ../../srf/fortran/simple_glide ./
```

and then execute it by typing

```./simple_glide
```

Either way, you will be prompted for the relevant configuration file, which is of course "ishom.a.config". After responding to the prompt, you will see 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.

## Plotting model output

how to plot the model output and see how you did relative to others

## Extra Credit

Some other things to try if you have time ...

Have them play w/ the grid spacing to see how that affects results?

Can we get a 0-order solution for these as well?