# Verification¶

Two types of errors may be distinguished: modeling errors and numerical errors. Modeling errors arise from not solving the right equations. Numerical errors result from not solving the equations right. The assessment of modeling errors is

validation, whereas the assessment of numerical errors is calledverification… Validation makes sense only after verification, otherwise agreement between measured and computed results may well be fortuitous.P. Wesseling, (2001)

Principles of Computational Fluid Dynamics, pp. 560–561 [128]

Verification is the essentially mathematical task of checking that the predictions of the numerical code are close to the predictions of a continuum model, the one which the numerical code claims to approximate. It is a crucial task for a code as complicated as PISM. In verification there is no comparison between model output and observations of nature. Instead, one compares exact solutions of the continuum model, in circumstances in which they are available, to their numerical approximations.

Reference [127] gives a broad discussion of verification and validation in computational fluid dynamics. See [26] and [23] for discussion of verification issues for the isothermal and thermomechanically coupled shallow ice approximation (SIA), respectively, and for exact solutions to these models, and [22], [38] for verification using an exact solution to the SSA equations for ice streams.

In PISM there is a separate executable `pismv`

which is used for SIA-related
verification, and there are additional scripts for SSA-related verification. The source
codes which are verified by `pismv`

are, however, exactly the same source files as are
run by the normal PISM executable `pismr`

. In technical terms, `pismv`

runs a derived
class of the PISM base class.

Test | Continuum model tested | Reference | Comments |
---|---|---|---|

A | isothermal SIA, steady, flat bed, constant accumulation | [26] | |

B | isothermal SIA, flat bed, zero accumulation | [26], [104] | similarity solution |

C | isothermal SIA, flat bed, growing accumulation | [26] | similarity solution |

D | isothermal SIA, flat bed, oscillating accumulation | [26] | uses compensatory accumulation |

E | isothermal SIA; as A, but with sliding in a sector | [26] | uses compensatory accumulation |

F | thermomechanically coupled SIA (mass and energy conservation), steady, flat bed | [69], [23] | uses compensatory accumulation and heating |

G | thermomechanically coupled SIA; as F but with oscillating accumulation | [69], [23] | ditto |

H | bed deformation coupled with isothermal SIA | [89] | joined similarity solution |

I | stream velocity computation using SSA (plastic till) | [38], [22] | |

J | shelf velocity computation using SSA | (source code) | |

K | pure conduction in ice and bedrock | [137] | |

L | isothermal SIA, steady, non-flat bed | (source code) | numerical ODE solution |

Test | Example invocation |
---|---|

A | `pismv -test A -Mx 61 -My 61 -Mz 11 -y 25000` |

B | `pismv -test B -Mx 61 -My 61 -Mz 11 -ys 422.45 -y 25000` |

C | `pismv -test C -Mx 61 -My 61 -Mz 11 -y 15208.0` |

D | `pismv -test D -Mx 61 -My 61 -Mz 11 -y 25000` |

E | `pismv -test E -Mx 61 -My 61 -Mz 11 -y 25000` |

F | `pismv -test F -Mx 61 -My 61 -Mz 61 -y 25000` |

G | `pismv -test G -Mx 61 -My 61 -Mz 61 -y 25000` |

H | `pismv -test H -Mx 61 -My 61 -Mz 11 -y 40034 -bed_def iso` |

I | `ssa_testi -ssa_method fd -Mx 5 -My 500 -ssa_rtol 1e-6 -ssafd_ksp_rtol 1e-11` |

J | `ssa_testj -ssa_method fd -Mx 60 -My 60 -ssafd_ksp_rtol 1e-12` |

K | `pismv -test K -Mx 6 -My 6 -Mz 401 -Mbz 101 -y 130000` |

L | `pismv -test L -Mx 61 -My 61 -Mz 31 -y 25000` |

Option | Description |
---|---|

`-test` |
Choose verification test by single character name; see Table 40. |

`-no_report` |
Do not report errors at the end of a verification run. |

`-eo` |
Only evaluate the exact solution; no numerical approximation at all. |

`-report_file` |
Save error report to a netCDF file. |

`-append` |
Append to a report file. |

Table 40 summarizes the many exact solutions currently available in PISM. Most
of these exact solutions are solutions of *free boundary problems* for partial
differential equations; only Tests A, E, J, K are fixed boundary value problems.

Table 41 shows how to run each of them on a coarse grids. Note that tests
I and J require special executables `ssa_testi,ssa_testj`

which are built with
configuration flag `Pism_BUILD_EXTRA_EXECS`

equal to `ON`

. Table 42
gives the special verification-related options of the `pismv`

executable.

Numerical errors are not, however, the dominant reasons why ice sheet models give
imperfect results. The largest sources of errors include those from using the wrong (e.g.
over-simplified or incorrectly-parameterized) continuum model, and from observational or
pre-processing errors present in input data. Our focus here on numerical errors has a
model-maintenance goal. It is *easier* to maintain code by quantitatively confirming that
it produces small errors in cases where those can be measured, rather than “eyeballing”
results to see that they are “right” according to human judgment.

The goal of verification is not generally to see that the error is zero at any particular resolution, or even to show that the error is small in a predetermined absolute sense. Rather the goals are

- to see that the error
*is*decreasing, - to measure the rate at which it decreases, and
- to develop a sense of the magnitude of numerical error before doing realistic ice sheet model runs.

Knowing the error decay rate may give a prediction of how fine a grid is necessary to achieve a desired smallness for the numerical error.

Therefore one must “go down” a grid refinement “path” and measure numerical error for each grid [127]. The refinement path is defined by a sequence of spatial grid cell sizes which decrease toward the refinement limit of zero size [138]. In PISM the timestep \(\Delta t\) is determined adaptively by a stability criterion (see section Understanding adaptive time-stepping). In PISM one specifies the number of grid points, thus the grid cell sizes because the overall dimensions of the computational box are normally fixed; see section Computational box. By “measuring the error for each grid” we mean computing a norm (or norms) of the difference between the numerical solution and the exact solution.

For a grid refinement path example, in tests of the thermomechanically-coupled SIA model one refines in three dimensions, and these runs produced Figures 13, 14, and 15 of [23]:

```
pismv -test G -max_dt 10.0 -y 25000 -Mx 61 -My 61 -Mz 61 -z_spacing equal
pismv -test G -max_dt 10.0 -y 25000 -Mx 91 -My 91 -Mz 91 -z_spacing equal
pismv -test G -max_dt 10.0 -y 25000 -Mx 121 -My 121 -Mz 121 -z_spacing equal
pismv -test G -max_dt 10.0 -y 25000 -Mx 181 -My 181 -Mz 181 -z_spacing equal
pismv -test G -max_dt 10.0 -y 25000 -Mx 241 -My 241 -Mz 241 -z_spacing equal
pismv -test G -max_dt 10.0 -y 25000 -Mx 361 -My 361 -Mz 361 -z_spacing equal
```

The last two runs require a supercomputer! In fact the \(361\times 361\times 361\) run
involves more than \(100\) million unknowns, updated at each of millions of time
steps. Appropriate use of parallelism (`mpiexec -n NN pismv`

) and of the `-skip`

modification to adaptive timestepping accelerates such fine-grid runs; see section
Understanding adaptive time-stepping.

Figures Fig. 25 through Fig. 29 in
Sample convergence plots show a sampling of the results of verifying PISM using the
tests described above. These figures were produced automatically using Python scripts
`test/vfnow.py`

and `test/vnreport.py`

. See section Utility and test scripts.

These figures *do not* show outstanding rates of convergence, relative to textbook partial
differential equation examples. For the errors in tests B and G, see the discussion of
free margin shape in [26]. For the errors in test I, the exact continuum solution is
not very smooth at the free boundary [38].

Previous | Up | Next |