Numerics¶

What this shows: the numerical settings of SWAN: the propagation scheme, how SWAN decides that a computation has converged, and how to choose a time step.

Prerequisites: Tutorial 5: Choosing model settings and Tutorial 6: A nonstationary hindcast.

You will learn:

• which propagation schemes PROP offers, and when to use them
• how NUMERIC sets the stopping criteria and iterations
• how to check the Courant number of a nonstationary time step

Data used: none.

The SWAN manual describes every option; SWAN's defaults suit most models.

Setup¶

In [1]:

from datetime import timedelta

import numpy as np

from rompy_swan.components.numerics import NUMERIC, PROP
from rompy_swan.subcomponents.numerics import (
    ACCUR,
    BSBT,
    DIRIMPL,
    GSE,
    NONSTAT,
    SIGIMPL,
    STAT,
    STOPC,
)
from rompy_swan.subcomponents.time import Delt

from datetime import timedelta import numpy as np from rompy_swan.components.numerics import NUMERIC, PROP from rompy_swan.subcomponents.numerics import ( ACCUR, BSBT, DIRIMPL, GSE, NONSTAT, SIGIMPL, STAT, STOPC, ) from rompy_swan.subcomponents.time import Delt

1. Propagation schemes: PROP¶

Without PROP, SWAN uses a higher-order scheme: SORDUP in stationary and S&L in nonstationary computations. PROP changes it:

• BSBT: first-order upwind. Robust and unconditionally stable, but more diffusive; use it with large time steps or when SWAN asks for it;
• GSE: the S&L scheme with a correction for the garden-sprinkler effect (swell breaking into separate rays over long distances), given the swell age.

In [2]:

print(PROP(scheme=BSBT()).render())
print(PROP(scheme=GSE(waveage=Delt(delt=timedelta(hours=12)))).render())

print(PROP(scheme=BSBT()).render()) print(PROP(scheme=GSE(waveage=Delt(delt=timedelta(hours=12)))).render())

PROP BSBT
PROP GSE waveage=43200.0 SEC

2. Convergence and iterations: NUMERIC¶

SWAN solves each step iteratively. STOPC (the default criterion) stops when the wave height and its curvature change little between iterations at enough grid points (npnts, in %). mode sets the maximum number of iterations:

• STAT(mxitst) for stationary computations (default 50);
• NONSTAT(mxitns) for nonstationary time steps (the manual advises 1 iteration per step with steps of 10 minutes or less).

In [3]:

print(NUMERIC(stop=STOPC(dabs=0.005, drel=0.01, curvat=0.005, npnts=99.5, mode=STAT(mxitst=50))).render())
print(NUMERIC(stop=STOPC(dabs=0.005, drel=0.01, curvat=0.005, npnts=99.5, mode=NONSTAT(mxitns=1))).render())

print(NUMERIC(stop=STOPC(dabs=0.005, drel=0.01, curvat=0.005, npnts=99.5, mode=STAT(mxitst=50))).render()) print(NUMERIC(stop=STOPC(dabs=0.005, drel=0.01, curvat=0.005, npnts=99.5, mode=NONSTAT(mxitns=1))).render())

NUMERIC STOPC dabs=0.005 drel=0.01 curvat=0.005 npnts=99.5 STATIONARY mxitst=50
NUMERIC STOPC dabs=0.005 drel=0.01 curvat=0.005 npnts=99.5 NONSTATIONARY mxitns=1

ACCUR is the older criterion, kept for compatibility. NUMERIC ... STAT mxitst=0 is a quick way to check that SWAN reads every input without computing.

In [4]:

print(NUMERIC(stop=ACCUR(drel=0.02, dhoval=0.02, dtoval=0.02, npnts=98.0, mode=STAT(mxitst=0))).render())

print(NUMERIC(stop=ACCUR(drel=0.02, dhoval=0.02, dtoval=0.02, npnts=98.0, mode=STAT(mxitst=0))).render())

NUMERIC ACCUR drel=0.02 dhoval=0.02 dtoval=0.02 npnts=98.0 STATIONARY mxitst=0

DIRIMPL, SIGIMPL, CTHETA and CSIGMA tune the implicit schemes in direction and frequency space, which rarely needs changing:

In [5]:

print(NUMERIC(dirimpl=DIRIMPL(cdd=0.5), sigimpl=SIGIMPL(css=0.5)).render())

print(NUMERIC(dirimpl=DIRIMPL(cdd=0.5), sigimpl=SIGIMPL(css=0.5)).render())

NUMERIC DIRIMPL cdd=0.5 SIGIMPL css=0.5

3. Choosing a nonstationary time step¶

The time step is the interval of the run period. The SWAN manual advises at most 10 minutes, and a Courant number (how many cells a wave crosses in one step) below 10 with the higher-order scheme; otherwise use PROP BSBT.

The fastest waves are the longest swell: in deep water, their group velocity is g T / (4 π). For the 0.02° grid of the tutorials, at 32°S:

In [6]:

period = 25.0  # s, the longest period of the spectral grid (flow = 0.04 Hz)
group_velocity = 9.81 * period / (4 * np.pi)
dx = 0.02 * 111_000 * np.cos(np.deg2rad(32.0))  # m, the east-west cell size
for minutes in [10, 20]:
    courant = group_velocity * minutes * 60 / dx
    print(f"{minutes} min: Courant number {courant:.1f}")

period = 25.0 # s, the longest period of the spectral grid (flow = 0.04 Hz) group_velocity = 9.81 * period / (4 * np.pi) dx = 0.02 * 111_000 * np.cos(np.deg2rad(32.0)) # m, the east-west cell size for minutes in [10, 20]: courant = group_velocity * minutes * 60 / dx print(f"{minutes} min: Courant number {courant:.1f}")

10 min: Courant number 6.2
20 min: Courant number 12.4

With 10 minutes the Courant number is about 6. With 20 minutes it exceeds 10, and SWAN stops the nonstationary computation with an error in PRINT, as Tutorial 6 mentions. For small domains, a series of stationary computations (COMPUTE_STAT with nonstationary times) is an alternative to long steps; see Stationary and nonstationary computations.

Summary¶

• PROP chooses the propagation scheme; BSBT is the robust choice for large steps.
• NUMERIC sets the stopping criterion and the iterations per computation.
• Keep nonstationary steps at 10 minutes or less, with a Courant number below 10.