Difference between revisions of "CISM exercise II: run diagnostic test cases"

From Interactive System for Ice sheet Simulation
Jump to: navigation, search
(The Confined Shelf test case)
(Introduction)
 
(58 intermediate revisions by 2 users not shown)
Line 1: Line 1:
 
== Introduction ==
 
== Introduction ==
  
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:
+
In this exercise, we will run a few of CISM's higher-order test cases, which span a wide range of flow regimes. Since the test cases are small and can be run on a single processor, we will run them "interactively" (that is, without submitting the jobs to the queue). First, make a directory on the scratch space to run the code and store the output, for example,
  
  cd tests/higher-order/
+
  mkdir /scratch2/username/cism
  
You will see the following list of subdirectories
+
(where "username" is your username). Now copy the relevant directories there. From within the "CISM-LANL-4-2011/tests/" subdirectory,
 +
 
 +
cp -r higher-order /scratch2/username/cism/
 +
 
 +
This subdirectory contains various test cases for the newer, "higher-order" dynamical core in CISM. Within the "higher-order" directory you will see the following subdirectories,
  
 
  dome/                # parabolic shaped dome with simple boundary conditions
 
  dome/                # parabolic shaped dome with simple boundary conditions
Line 12: Line 16:
 
  shelf/                # ice shelf test cases on simplified domains
 
  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.
+
along with a few other files. You may want to look over the tests/higher-order/README file at some point but 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.
 +
 
 +
Before we can run the code we first need to start an interactive session using one node,
 +
 
 +
msub -I -A s11_cesm
 +
 
 +
Once your prompt has been returned to you, make sure that you are still in the directory you want to be in on your scratch space. Since we've started up a new session you will need to re-source the environment script,
  
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
+
source /usr/projects/cesm/cism/cism-serial-env
  
cp ../../../example-drivers/simple_glide/src/simple_gide ./
+
CISM uses Python to read and write input/output netCDF files. Here, we need a fairly specific set of Python toolboxes. To make sure that you have access to the correct version of Python, type "python". You should see the following first few lines:
  
or make a virtual link to the actual executable file
+
Enthought Python Distribution -- www.enthought.com
 +
Version: 7.0-2 (64-bit)
 +
Python 2.7.1 |EPD 7.0-2 (64-bit)| (r271:86832, Nov 29 2010, 13:51:37)
  
  ln -s ../../../example-drivers/simple_glide/src/simple_glide ./
+
If instead you see,
 +
 
 +
  Python 2.4.3 (#1, Sep  8 2010, 11:37:47)
 +
[GCC 4.1.2 20080704 (Red Hat 4.1.2-48)] on linux2
 +
 
 +
let us know so that we can alter a few things and give you access to the correct version of Python. To escape out of Python, use <CTRL><D>.
  
 
== The Dome test case ==
 
== 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
+
This is a very simple test case, simulating the three-dimensional flow field within an isothermal, 3d, parabolic shaped dome with no-slip basal boundary conditions and zero flux lateral boundary conditions. To run the test case, first change into the "dome" subdirectory,
 +
 
 +
cd /scratch2/username/cism/higher-order/dome/
 +
 
 +
Now make a virtual link to the executable file you built in the first exercise. 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 link to it from your scratch space type,
 +
 
 +
ln -s /path/to/your/project/space/CISM-LANL-4-2011/example-drivers/simple_glide/src/simple_glide ./
 +
 
 +
You can then execute the test with
  
 
  python dome.py
 
  python dome.py
Line 32: Line 57:
 
  Enter name of GLIDE configuration file to be read
 
  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.
+
Here, the "dome.config" script from within this directory is passed to '''simple_glide''' by the "dome.py" python script.  
 +
 
 +
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 for an option you want to understand (e.g. "evolution = 3") and then "grep" for that option in the file ''glide_types.F90'' in the ''libglide/'' subdirectory from within the main directory where you built the code.
 +
 
 +
[[Image:dome-grab.png|thumb|left|700px|'''Dome test case:''' Screen grab of model output from dome test case using NCVIEW. Shown are the velocity magnitude in map view (color plot), the surface speed across the ice dome, and a vertical velocity profile from approximately half way between the ice divide (center) and the margin (click for higher-resolution image).]]
  
 
As the code runs, you will see some output to the screen like
 
As the code runs, you will see some output to the screen like
Line 62: Line 91:
 
# The current time step we are solving for.
 
# 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 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
+
# the non-linear iteration number, the current 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.
+
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 (i.e. we have the answer). Here it took 43 nonlinear 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:
+
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 (including the machine used here) is ''NCVIEW''. To examine the output file using NCVIEW, type
  
[[Image:dome-grab.png|thumb|left|700px|Screen grab of model output from dome test case using NCVIEW. Shown are the velocity magnitude in map view (color plot), the surface speed across the ice dome, and a vertical velocity profile from approximately half way between the ice divide (center) and the margin (click for higher-resolution image).]]
+
ncview output/dome.out.nc
 +
 
 +
Your output for the variable ''velnorm'' (the ice speed) at level "0" (sigma coordinate level 0, which is the upper surface of the ice sheet), should look something like what is shown in the figure labeled '''Dome test case'''. Take a minute to play around with the different buttons on NCVIEW to see what they do. You can step through the vertical levels of any of the 3d model output fields, click on the color contour plots to obtain 2d profiles, change the color scheme, and for variables that change in time, run simple "movies" showing a variables evolution over time (we will do this later).
  
 
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.
 
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.
Line 74: Line 105:
 
== The Confined Shelf test case ==
 
== 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.   
+
[[Image:confined-shelf.png|thumb|right|500px|'''Confined shelf test case:''' Along (right) and across (left) flow velocities for the confined shelf test case. The upper panel (color) shows results using CISM and the lower panel shows results from the [http://homepages.vub.ac.be/~phuybrec/eismint/iceshelf.html EISMINT ice shelf intercomparison project] for the same test case. Solid black contour lines are the same in both plots (click for higher-resolution image).]]
 +
 
 +
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 following 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,
+
To run the test case, change from the "higher-order/dome/" subdirectory into the "/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, link to the '''simple_glide''' executable as you did before, then run the test with,
  
 
  python confined-shelf.py
 
  python confined-shelf.py
Line 84: Line 117:
 
  ncview output/confined-shelf.out.nc
 
  ncview output/confined-shelf.out.nc
  
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 [[shallow-shelf approxmation|SSA]] model. It is from [http://homepages.vub.ac.be/~phuybrec/eismint/shelf-descr.pdf page 7] of the [http://homepages.vub.ac.be/~phuybrec/eismint.html EISMINT] (European Ice Sheet Model InTercomparison) [http://homepages.vub.ac.be/~phuybrec/eismint/iceshelf.html ice shelf intercomparison project] documentation.
+
A color contour plot of CISM output (made in Matlab) is shown in the figure labeled '''Confined shelf test case'''. The black and white contour plot shows output for the same experiment using an [[shallow-shelf approxmation|SSA]] (ice shelf) model. It is from experiment 3 (page 7) of the [http://homepages.vub.ac.be/~phuybrec/eismint.html EISMINT] (European Ice Sheet Model InTercomparison) [http://homepages.vub.ac.be/~phuybrec/eismint/iceshelf.html ice shelf intercomparison project] [http://homepages.vub.ac.be/~phuybrec/eismint/shelf-descr.pdf documentation]. We will return to this experiment and add some additional complexity to it in a later exercise.
 +
 
 +
If there is time, you may also want to try running the "circular-shelf" experiment, which demonstrates that CISM can also implement an accurate ice shelf boundary condition for an ice shelf front with a non-trivial shape in map view (i.e. one for which the shelf-front normal vectors are not parallel to coordinate directions).
  
 
== The ISMIP-HOM test cases ==
 
== The ISMIP-HOM test cases ==
  
 +
[[Image:ismiphom.a.jpg|thumb|left|500px|'''ISMIP-HOM A Setup: '''Doubly periodic basal roughness with no sliding. Ice thickness, basal topography, and surface elevation are shown. At the lateral boundaries, velocities are doubly periodic (click for higher-resolution image).]]
  
[http://homepages.ulb.ac.be/~fpattyn/ismip/ 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 &lambda;=160 km, the velocities solutions essentially look like that from a shallow ice model. Halving &lambda; 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 &lambda; = 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).
+
[[Image:ismiphom.c.jpg|thumb|left|500px|'''ISMIP-HOM C Setup: '''Doubly periodic basal traction coefficient with sliding. Ice thickness, basal topography, and surface elevation are shown. At the lateral boundaries, velocities are doubly periodic (click for higher-resolution image).]]
  
 +
The last set of diagnostic problems we will look at are from the [http://homepages.ulb.ac.be/~fpattyn/ismip/ ISMIP-HOM] test suite, which was specifically designed for "higher-order" models and nicely demonstrates the difference between higher-order and 0-order (or "shallow ice") models. While the test suite includes a total of 6 tests we will look only at the tests for diagnostic solutions on idealized, three-dimensional domains. Each of these tests (A and C) includes a subset of 6 tests for a range of domain lengths.
  
'''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''.
+
Both tests consist of a uniformly sloping slab of ice with periodic lateral velocities in the ''x'' and ''y'' directions (i.e. in map plane). For test A, the basal topography varies periodically in ''x'' and ''y'' directions and the there is a no-slip basal boundary condition. For test C, basal sliding is allowed. The basal traction coefficient varies periodically in ''x'' and ''y'' and the thickness is uniform throughout the domain. While the ''amplitude'' of the variations (topography in A and traction coefficient in C) is the same for all tests, the ''wavelength'', &lambda;, is decreased by a factor of two for each successive test. For &lambda;=160 km, the velocity solutions are essentially equal to those from a 0-order shallow ice model. When halving &lambda; 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 show relevant input data for each of the two experiments for &lambda; = 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).
  
[[Image:ismiphom.a.jpg]]
+
=== Running the model test cases ===
  
 +
To run the experiments, we will use some python scripts developed by colleagues at the University of Montana (also, see this [[Validation and Verification|link]]). As with the other test cases, several Python scripts set up the necessary netCDF input and output files. The python scripts simplify things here by allowing you to run and plot the results from multiple tests and multiple domain wavelengths sequentially. In addition, they plot CISM output relative to the model means and standard deviations from the actual benchmark study of [http://www.the-cryosphere.net/2/95/2008/tc-2-95-2008.html Pattyn et. al (2008)]. This is a great convenience (as anyone who has ever done this on their own will attest to!). First, move into the ''tests/higher-order/ismip-hom'' subdirectory. To execute test A, for &lambda;=160, 80, and 40 km, type
  
'''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).
+
python runISMIPHOM.py --exp=a --size=160,80,40
  
[[Image:ismiphom.c.jpg]]
+
As with the other test cases above, you should see some screen output showing model residuals decreasing as the nonlinear iterations proceed (***Did you remember to make a virtual link to the '''simple_glide''' executable?***?).
  
 +
=== Plotting model output ===
  
 +
To compare CISM output from your model runs with that from the ISMIP-HOM benchmark study of [http://www.the-cryosphere.net/2/95/2008/tc-2-95-2008.html Pattyn et. al (2008)], we will execute the python plotting scripts in a similar manner. For test A, type
  
To set up the experiments, we will use some configuration files and python scripts developed by Tim Bocek and Jesse Johnson (also, see this [[Validation and Verification|link]]). These set the correct flags, so that Glimmer/CISM calls the necessary subroutines, and construct the necessary input netCDF files.
+
python plotISMIPHOM.py --exp=a --size=160,80,40
  
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
+
Your output figure will have a ".png" extension and will be placed in the ''output/'' subdirectory. It should look something like the figure here labeled '''ISMIP-HOM A Output'''. There is a simple image viewer on the cluster where we have been running the code (the classic XV). To use it to view the output form you test case type,
  
  cd tests/ISMIP-HOM/; ls -l
+
  xv output/ISMIP-HOM-A-glm1.png
  
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.
+
[[Image:ISMIP-HOM-A-glm1.png|thumb|right|500px|'''ISMIP-HOM A Output: '''CISM output for experiment A run for wavelengths of 160, 80, and 40 km. Black solid lines are CISM output, blue (red) dashed and shaded areas show the mean and standard deviation from all first-order (Stokes) models participating in the benchmark exercise. Model velocity is shown on the vertical axis and normalized distance along the domain, at ''y''=0.25*&lambda;, is shown on the horizontal axis (click for higher-resolution image).]]
  
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.  
+
[[Image:ISMIP-HOM-C-glm1.png|thumb|right|500px|'''ISMIP-HOM C Output: '''CISM output for experiment C run for wavelengths of 160, 80, and 40 km. Black solid lines are CISM output, blue (red) dashed and shaded areas show the mean and standard deviation from all first-order (Stokes) models participating in the benchmark exercise. Model velocity is shown on the vertical axis and normalized distance along the domain, at ''y''=0.25*&lambda;, is shown on the horizontal axis  (click for higher-resolution image).]]
  
emacs ishom.a.config &
+
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 like the figure labeled '''ISMIP-HOM-C Output'''.
  
gives
+
For additional information on running and plotting results for the ISMIP-HOM test suite, see the README file in the ''tests/higher-order/ismip-hom'' subdirectory.
  
[grid]
+
== Additional Exercises ==
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
+
=== ISMIP-HOM A with shallow-ice dynamics ===
  
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
+
To clarify the importance of the higher-order stresses in the model velocity solutions, it is instructive to go back and re-run one of the above tests using the shallow-ice model. To do this, we first need to edit some of the configuration file options in the file ''ishom.a.config''. Copy the original file to a backup version first (e.g. cp ishom.a.config ishom.a.config.orig). Now open ''ishom.a.config'' with your favorite editor (e.g. VI or Emacs) and look for the following sections:
  
  cp ../../src/fortran/simple_glide ./
+
  [options]
 +
flow_law = 2            # constant and uniform rate factor
 +
periodic_ew = 1        # doubly periodic lateral boundary conditions
 +
periodic_ns = 1
 +
evolution = 3
  
and then execute it by typing
+
[ho_options]
 +
diagnostic_scheme = 1  # Payne/Price 1st-order dynamics
 +
which_ho_babc = 4      # no-slip basal boundary conditions       
 +
which_ho_efvs = 0      # nonlinear eff. visc. w/ n=3
 +
which_ho_sparse = 1    # use SLAP GMRES for linear solver
  
./simple_glide
+
To implement 0-order shallow ice dynamics rather than first-order dynamics, change the following flags in the ''options'' and ''ho_options'' sections,
  
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
+
[options]
 +
evolution = 0          # now SIA dynamics!
  
  echo ishom.a.config | ./simple_glide    - or -    echo ishom.a.config | simple_glide         
+
  [ho_options]
(the latter if "simple_glide" is already in your path)
+
diagnostic_scheme = 0  # now SIA dynamics!
  
Either way, after responding to the prompt, you ''should'' see some model output that looks something like this:
+
Now re-run the ISMIP-HOM test case
  
  Running Payne/Price higher-order dynamics solver
+
  python runISMIPHOM.py --exp=a, --size=160,80,40
 
+
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'''.
+
You won't see any output, but you will probably notice that the model gets through both of these tests much more quickly than when using the first-order stress balance (one good thing about shallow ice dynamics, they are computationally cheap!). When the model is done running, plot the results again,
  
 +
python plotISMIPHOM.py --exp=a, --size=160,80,40
  
*'''HAVING PROBLEMS?'''  
+
Your results should look something like what is shown in the figure labeled '''ISMIP-HOM-A SIA Output'''. Can you explain why the velocity field from the shallow ice model is identical despite the change in the wavelength of the basal topography for the three experiments? As far as the shallow-ice model is concerned, these three domains are all identical because the flow rate is controlled only by the local slope and ice thickness (which is the same despite the different wavelengths of the basal topography).
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 ==
+
[[Image:ishom-a-sia.png|thumb|left|500px|'''ISMIP-HOM A SIA Output: '''CISM output for experiment A, run for wavelengths of 160, 80, and 40 km '''using shallow-ice dyanmics'''. Black solid lines are CISM output, blue (red) dashed and shaded areas show the mean and standard deviation from all first-order (Stokes) models participating in the benchmark exercise. Model velocity is shown on the vertical axis and normalized distance along the domain, at ''y''=0.25*&lambda;, is shown on the horizontal axis  (click for higher-resolution image).]]
  
 +
=== ISMIP-HOM A: Newton versus Picard ===
  
 +
The increasing difficulty of the ISMIP-HOM experiments as the domain wavelength decreases provides a good opportunity to demonstrate the differences between handling the model nonlinear with a Picard versus a Newton iteration (as discussed in more detail on the [[CISM higher-order dynamics: model solution|model solution]] page). Because the current Newton solver in CISM is still under development, we have to make a few more simplifications here to compare the two. In particular, we have not yet implemented periodic boundary conditions in the Newton iteration in which case we will have to turn these "off" for the Picard iteration as well. As above, we need to edit a few sections in your your '''original''' '''ishom.a.config''' file (which was hopefully copied before you made the previous edits). First, to run the test cases using Picard, change the periodic flags in the ''options'' section of your .config file as follows:
  
'''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.'''
+
[options]
 +
periodic_ew = 0        # Now zero-velocity rather than doubly periodic!
 +
periodic_ns = 0
  
 +
No re-run the test with the Picard iteration and the new boundary conditions. To get an approximate total time for the run and to dump the output to a .txt file (in case we want to plot it later on), use
  
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 [http://www.the-cryosphere.net/2/95/2008/tc-2-95-2008.html 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.
+
nohup time python runISMIPHOM.py -e a -s 40,20,10 > picard-log.txt &
 
+
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
+
tests/ISMIP-HOM/
+
 
+
 
+
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 [[Validation and Verification#Higher Order Isothermal Flow|HERE]] for a more detailed description of these plots).
+
 
+
[[Image:ISMIP-HOM-A-glm1.png]]
+
 
+
 
+
 
+
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 [[Validation and Verification#Higher Order Isothermal Flow|HERE]] for a more detailed description of these plots).
+
 
+
[[Image:ISMIP-HOM-C-glm1.png]]
+
 
+
== Additional Exercises ==
+
 
+
* 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]
+
Also notice that we are now running a few of the shorter wavelength test cases in order to work the model a little bit harder.
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?
+
Do the same but using the Newton iteration instead. For this case we need to add an additional flag to the ''ho_options'' section of the .config file:
  
* 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).
+
[ho_options]
 +
which_ho_nonlinear = 1  # add this flag to call JFNK for nonlinear iteration rather than Picard!
  
 +
To re-run the test case, time the model run and save the output, use
  
'''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.
+
nohup time python runISMIPHOM.py -e a -s 40,20,10 > newton-log.txt &
  
[[Image: ncviewHO.png]]
+
Note that the "nohup" command sets your job to running in the background so that you can do other things while you wait for it to complete. To check the status of your job, type
  
 +
jobs -l
  
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.
+
If your job is still running you will see something like
  
 +
72195 Running    nohup time python runISMIPHOM.py -e a -s 40,20,10 > newton-log.txt &
  
'''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.
+
where the first number is the job ID.
  
[[Image: ishoma-80km-HOvsSIA.jpg]]
+
When your jobs have both completed, look in the ''scratch/'' subdirectory for your log.txt files. You should have a record of the iteration count for each job and also a record of the total time to run the job at the very bottom. You should notice that it takes ~4x fewer nonlinear iterations to reach converged solutions. The Newton iteration in this version of the code has not been optimized yet, so there is an additional "cost" associated with using it. Thus, you may notice that the overall savings in computational time is only about ~25% relative to Picard. In other developmental versions of the code for which the Newton iteration has been better optimized the computational time savings is usually a factor of ~2-5x (e.g. see the figure in [[CISM higher-order dynamics: model solution#Newton-based Methods for Solutions of the Non-linear System|'''this''']]) section.
  
  
* 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?
+
Go to the [[CISM exercise III: add a module to CISM and run prognostic test case|third set of exercises]].
  
* 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?
+
[[International Workshop on Ice Flow Modeling, Beijing Normal University (March 21-25, 2011)|Return to main course page]]

Latest revision as of 16:17, 9 May 2011

Contents

Introduction

In this exercise, we will run a few of CISM's higher-order test cases, which span a wide range of flow regimes. Since the test cases are small and can be run on a single processor, we will run them "interactively" (that is, without submitting the jobs to the queue). First, make a directory on the scratch space to run the code and store the output, for example,

mkdir /scratch2/username/cism

(where "username" is your username). Now copy the relevant directories there. From within the "CISM-LANL-4-2011/tests/" subdirectory,

cp -r higher-order /scratch2/username/cism/

This subdirectory contains various test cases for the newer, "higher-order" dynamical core in CISM. Within the "higher-order" directory you will see the following 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. You may want to look over the tests/higher-order/README file at some point but 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.

Before we can run the code we first need to start an interactive session using one node,

msub -I -A s11_cesm

Once your prompt has been returned to you, make sure that you are still in the directory you want to be in on your scratch space. Since we've started up a new session you will need to re-source the environment script,

source /usr/projects/cesm/cism/cism-serial-env 

CISM uses Python to read and write input/output netCDF files. Here, we need a fairly specific set of Python toolboxes. To make sure that you have access to the correct version of Python, type "python". You should see the following first few lines:

Enthought Python Distribution -- www.enthought.com
Version: 7.0-2 (64-bit)
Python 2.7.1 |EPD 7.0-2 (64-bit)| (r271:86832, Nov 29 2010, 13:51:37) 

If instead you see,

Python 2.4.3 (#1, Sep  8 2010, 11:37:47)
[GCC 4.1.2 20080704 (Red Hat 4.1.2-48)] on linux2

let us know so that we can alter a few things and give you access to the correct version of Python. To escape out of Python, use <CTRL><D>.

The Dome test case

This is a very simple test case, simulating the three-dimensional flow field within an isothermal, 3d, parabolic shaped dome with no-slip basal boundary conditions and zero flux lateral boundary conditions. To run the test case, first change into the "dome" subdirectory,

cd /scratch2/username/cism/higher-order/dome/

Now make a virtual link to the executable file you built in the first exercise. 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 link to it from your scratch space type,

ln -s /path/to/your/project/space/CISM-LANL-4-2011/example-drivers/simple_glide/src/simple_glide ./

You can then execute the test with

python dome.py

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 "dome.config" script from within this directory is passed to simple_glide by the "dome.py" python script.

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 for an option you want to understand (e.g. "evolution = 3") and then "grep" for that option in the file glide_types.F90 in the libglide/ subdirectory from within the main directory where you built the code.

Dome test case: Screen grab of model output from dome test case using NCVIEW. Shown are the velocity magnitude in map view (color plot), the surface speed across the ice dome, and a vertical velocity profile from approximately half way between the ice divide (center) and the margin (click for higher-resolution image).

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

  1. 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.
  2. The current time step we are solving for.
  3. The dynamics scheme we are using (here, the Blatter-Pattyn equations as formulated and solved by Payne and Price).
  4. the non-linear iteration number, the current 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 (i.e. we have the answer). Here it took 43 nonlinear 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 (including the machine used here) is NCVIEW. To examine the output file using NCVIEW, type

ncview output/dome.out.nc

Your output for the variable velnorm (the ice speed) at level "0" (sigma coordinate level 0, which is the upper surface of the ice sheet), should look something like what is shown in the figure labeled Dome test case. Take a minute to play around with the different buttons on NCVIEW to see what they do. You can step through the vertical levels of any of the 3d model output fields, click on the color contour plots to obtain 2d profiles, change the color scheme, and for variables that change in time, run simple "movies" showing a variables evolution over time (we will do this later).

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

Confined shelf test case: Along (right) and across (left) flow velocities for the confined shelf test case. The upper panel (color) shows results using CISM and the lower panel shows results from the EISMINT ice shelf intercomparison project for the same test case. Solid black contour lines are the same in both plots (click for higher-resolution image).

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 following 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 from the "higher-order/dome/" subdirectory into the "/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, link to the simple_glide executable as you did before, then run the test with,

python confined-shelf.py

As in the previous test case, you should see gradually decreasing residuals as screen output. When the test completes, examine the output with

ncview output/confined-shelf.out.nc

A color contour plot of CISM output (made in Matlab) is shown in the figure labeled Confined shelf test case. The black and white contour plot shows output for the same experiment using an SSA (ice shelf) model. It is from experiment 3 (page 7) of the EISMINT (European Ice Sheet Model InTercomparison) ice shelf intercomparison project documentation. We will return to this experiment and add some additional complexity to it in a later exercise.

If there is time, you may also want to try running the "circular-shelf" experiment, which demonstrates that CISM can also implement an accurate ice shelf boundary condition for an ice shelf front with a non-trivial shape in map view (i.e. one for which the shelf-front normal vectors are not parallel to coordinate directions).

The ISMIP-HOM test cases

ISMIP-HOM A Setup: Doubly periodic basal roughness with no sliding. Ice thickness, basal topography, and surface elevation are shown. At the lateral boundaries, velocities are doubly periodic (click for higher-resolution image).
ISMIP-HOM C Setup: Doubly periodic basal traction coefficient with sliding. Ice thickness, basal topography, and surface elevation are shown. At the lateral boundaries, velocities are doubly periodic (click for higher-resolution image).

The last set of diagnostic problems we will look at are from the ISMIP-HOM test suite, which was specifically designed for "higher-order" models and nicely demonstrates the difference between higher-order and 0-order (or "shallow ice") models. While the test suite includes a total of 6 tests we will look only at the tests for diagnostic solutions on idealized, three-dimensional domains. Each of these tests (A and C) includes a subset of 6 tests for a range of domain lengths.

Both tests consist of a uniformly sloping slab of ice with periodic lateral velocities in the x and y directions (i.e. in map plane). For test A, the basal topography varies periodically in x and y directions and the there is a no-slip basal boundary condition. For test C, basal sliding is allowed. The basal traction coefficient varies periodically in x and y and the thickness is uniform throughout the domain. While the amplitude of the variations (topography in A and traction coefficient in C) is the same for all tests, the wavelength, λ, is decreased by a factor of two for each successive test. For λ=160 km, the velocity solutions are essentially equal to those from a 0-order shallow ice model. When 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 show 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).

Running the model test cases

To run the experiments, we will use some python scripts developed by colleagues at the University of Montana (also, see this link). As with the other test cases, several Python scripts set up the necessary netCDF input and output files. The python scripts simplify things here by allowing you to run and plot the results from multiple tests and multiple domain wavelengths sequentially. In addition, they plot CISM output relative to the model means and standard deviations from the actual benchmark study of Pattyn et. al (2008). This is a great convenience (as anyone who has ever done this on their own will attest to!). First, move into the tests/higher-order/ismip-hom subdirectory. To execute test A, for λ=160, 80, and 40 km, type

python runISMIPHOM.py --exp=a --size=160,80,40

As with the other test cases above, you should see some screen output showing model residuals decreasing as the nonlinear iterations proceed (***Did you remember to make a virtual link to the simple_glide executable?***?).

Plotting model output

To compare CISM output from your model runs with that from the ISMIP-HOM benchmark study of Pattyn et. al (2008), we will execute the python plotting scripts in a similar manner. For test A, type

python plotISMIPHOM.py --exp=a --size=160,80,40

Your output figure will have a ".png" extension and will be placed in the output/ subdirectory. It should look something like the figure here labeled ISMIP-HOM A Output. There is a simple image viewer on the cluster where we have been running the code (the classic XV). To use it to view the output form you test case type,

xv output/ISMIP-HOM-A-glm1.png
ISMIP-HOM A Output: CISM output for experiment A run for wavelengths of 160, 80, and 40 km. Black solid lines are CISM output, blue (red) dashed and shaded areas show the mean and standard deviation from all first-order (Stokes) models participating in the benchmark exercise. Model velocity is shown on the vertical axis and normalized distance along the domain, at y=0.25*λ, is shown on the horizontal axis (click for higher-resolution image).
ISMIP-HOM C Output: CISM output for experiment C run for wavelengths of 160, 80, and 40 km. Black solid lines are CISM output, blue (red) dashed and shaded areas show the mean and standard deviation from all first-order (Stokes) models participating in the benchmark exercise. Model velocity is shown on the vertical axis and normalized distance along the domain, at y=0.25*λ, is shown on the horizontal axis (click for higher-resolution image).

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 like the figure labeled ISMIP-HOM-C Output.

For additional information on running and plotting results for the ISMIP-HOM test suite, see the README file in the tests/higher-order/ismip-hom subdirectory.

Additional Exercises

ISMIP-HOM A with shallow-ice dynamics

To clarify the importance of the higher-order stresses in the model velocity solutions, it is instructive to go back and re-run one of the above tests using the shallow-ice model. To do this, we first need to edit some of the configuration file options in the file ishom.a.config. Copy the original file to a backup version first (e.g. cp ishom.a.config ishom.a.config.orig). Now open ishom.a.config with your favorite editor (e.g. VI or Emacs) and look for the following sections:

[options]
flow_law = 2            # constant and uniform rate factor
periodic_ew = 1         # doubly periodic lateral boundary conditions
periodic_ns = 1
evolution = 3
[ho_options]
diagnostic_scheme = 1   # Payne/Price 1st-order dynamics
which_ho_babc = 4       # no-slip basal boundary conditions        
which_ho_efvs = 0       # nonlinear eff. visc. w/ n=3
which_ho_sparse = 1     # use SLAP GMRES for linear solver

To implement 0-order shallow ice dynamics rather than first-order dynamics, change the following flags in the options and ho_options sections,

[options]
evolution = 0           # now SIA dynamics!
[ho_options]
diagnostic_scheme = 0   # now SIA dynamics!

Now re-run the ISMIP-HOM test case

python runISMIPHOM.py --exp=a, --size=160,80,40

You won't see any output, but you will probably notice that the model gets through both of these tests much more quickly than when using the first-order stress balance (one good thing about shallow ice dynamics, they are computationally cheap!). When the model is done running, plot the results again,

python plotISMIPHOM.py --exp=a, --size=160,80,40

Your results should look something like what is shown in the figure labeled ISMIP-HOM-A SIA Output. Can you explain why the velocity field from the shallow ice model is identical despite the change in the wavelength of the basal topography for the three experiments? As far as the shallow-ice model is concerned, these three domains are all identical because the flow rate is controlled only by the local slope and ice thickness (which is the same despite the different wavelengths of the basal topography).

ISMIP-HOM A SIA Output: CISM output for experiment A, run for wavelengths of 160, 80, and 40 km using shallow-ice dyanmics. Black solid lines are CISM output, blue (red) dashed and shaded areas show the mean and standard deviation from all first-order (Stokes) models participating in the benchmark exercise. Model velocity is shown on the vertical axis and normalized distance along the domain, at y=0.25*λ, is shown on the horizontal axis (click for higher-resolution image).

ISMIP-HOM A: Newton versus Picard

The increasing difficulty of the ISMIP-HOM experiments as the domain wavelength decreases provides a good opportunity to demonstrate the differences between handling the model nonlinear with a Picard versus a Newton iteration (as discussed in more detail on the model solution page). Because the current Newton solver in CISM is still under development, we have to make a few more simplifications here to compare the two. In particular, we have not yet implemented periodic boundary conditions in the Newton iteration in which case we will have to turn these "off" for the Picard iteration as well. As above, we need to edit a few sections in your your original ishom.a.config file (which was hopefully copied before you made the previous edits). First, to run the test cases using Picard, change the periodic flags in the options section of your .config file as follows:

[options]
periodic_ew = 0         # Now zero-velocity rather than doubly periodic!
periodic_ns = 0 

No re-run the test with the Picard iteration and the new boundary conditions. To get an approximate total time for the run and to dump the output to a .txt file (in case we want to plot it later on), use

nohup time python runISMIPHOM.py -e a -s 40,20,10 > picard-log.txt &

Also notice that we are now running a few of the shorter wavelength test cases in order to work the model a little bit harder.

Do the same but using the Newton iteration instead. For this case we need to add an additional flag to the ho_options section of the .config file:

[ho_options]
which_ho_nonlinear = 1  # add this flag to call JFNK for nonlinear iteration rather than Picard!

To re-run the test case, time the model run and save the output, use

nohup time python runISMIPHOM.py -e a -s 40,20,10 > newton-log.txt &

Note that the "nohup" command sets your job to running in the background so that you can do other things while you wait for it to complete. To check the status of your job, type

jobs -l

If your job is still running you will see something like

72195 Running     nohup time python runISMIPHOM.py -e a -s 40,20,10 > newton-log.txt &

where the first number is the job ID.

When your jobs have both completed, look in the scratch/ subdirectory for your log.txt files. You should have a record of the iteration count for each job and also a record of the total time to run the job at the very bottom. You should notice that it takes ~4x fewer nonlinear iterations to reach converged solutions. The Newton iteration in this version of the code has not been optimized yet, so there is an additional "cost" associated with using it. Thus, you may notice that the overall savings in computational time is only about ~25% relative to Picard. In other developmental versions of the code for which the Newton iteration has been better optimized the computational time savings is usually a factor of ~2-5x (e.g. see the figure in this) section.


Go to the third set of exercises.

Return to main course page