Sumedh Joshi bio photo

Sumedh Joshi

Applied Mathematician

Email Twitter Facebook LinkedIn Instagram Github

A necessary (but not sufficient) condition for numerical stability is the Courant-Freidrich-Lewy condition which restricts the maximum timestep:

where is the grid spacing and the characteristic velocity of the flow. This is easy to calculate if the grid is regular, but what if your grid is a discretized slice of the ocean, and consequently looks something like this:

A 540,000 node grid with 15 GLL points per direction in each of the 200 x 12 elements (the boxes drawn represent each element).

in which each of the tiny boxes themselves represent a 15 x 15 mesh of Gauss-Lobatto-Legendre quadrature points. This is the case in a high-order element-based code like the spectral multi-domain penalty method, and the computation of the CFL condition (and namely the local deformation ) is not quite so straightforward.

Computing the grid spacing in the master element

Call the coordinates on the master element , and noting that there are GLL points in each direction, these coordinates are vectors . The mapping functions that define the physical grid are defined on each element as and . Call the collection of all of these maps and (drop the subscripts. What we want, for computing the CFL condition, is some measure of and , the spacing between grid points in a deformed domain like the one shown above.

Probably there are many ways to do this, but the one I chose was to make use of the fact that it is easy to get a measure of and on the master element. For example for the one-dimensional case, compute using a centered finite difference approximation:

for , with and . The same idea extends to two-dimensions and computing along with .

Computing the spacing everywhere in grid

Remember that we can easily compute derivatives in and ; that is the whole point of defining these coordinates, to be able to use them in constructing spectral differentiation matrices and computing numerical derivatives. Using the mapping functions and we can then use the chain rule to write the grid spacing on each element as (again dropping the in the subscript) as functions of derivatives in coordinates:

These quantities are shown in the figures below; the function is basically constant and so is kind of boring (the grid is uniformly spaced in the horizontal), but the function has some interesting character as the bathymetry slopes and shallows.

Top: the local grid spacing in x. Bottom: the local grid spacing in z.

Computing a maximum time-step

Using these two quantities, and , we can define the maximum time-step allowed in the simulation as

which ensures that we are satisfying and everywhere in the domain . All of the quantities on the right hand side of the above are computable, and are defined above.

Example: shoaling internal wave

An internal wave on the grid shown before; this initial velocity disturbance will propagate to the right and shoal (steepen and potentially break). The magnitude of the velocity vector is what is visualized.
Same plot as above but zoomed into where the velocity is non-zero.

Setting an initial velocity field that represents a propagating internal wave (as shown in the two pictures above), we can compute the quantity :

The minimum value of this function is the maximum allowable time-step as determined by the CFL condition. Of course, this function is not constant, and so parts of the grid will have too fine a time-step (i.e. if there is nothing going on in the velocity field there). But this inefficiency is inescapable since the grid has to advance in time together, and so we have to take the minimum value over the whole grid. Corresponding to the the internal wave field shown above, is shown in the two figures below.

The time-step restriction all over the grid. Of course we have to take the minimum of these since the time-step has to be the same over the whole grid.
The time-step restriction all over the grid (zoomed in to just where the velocity is non-zero). Of course we have to take the minimum of these since the time-step has to be the same over the whole grid.

For this particular field, the maximum time-step allowed, as computed by this method, turns out to be seconds, and is governed by a tiny region in the middle of wave field near an element boundary. This makes sense, near the center of the wave the velocity is the highest, and near an element boundary the grid spacing is the most fine.

Finally, what if we skipped this computation and instead used a heuristic to estimate the time-step? If instead of doing the above computation, we were to establish a minimum time-step by guessing an average grid spacing and the using the maximum velocity, we would get a very different (and too large) answer of seconds. If we attempted to be more conservative and calculated the difference between closest two grid points and divided by the maximum velocity in the grid, we’d get a far too conservative answer of seconds. The point being: don’t skip this computation.