Numerical methods¶

GeoSWE is a cell-centered finite-volume scheme on a uniform Cartesian grid. Each cell average is updated from numerical fluxes across its four faces plus the split source terms.

Riemann fluxes¶

The face flux comes from an approximate Riemann solver, selected by Config.flux:

"hllc": the HLLC solver. It resolves the left/right acoustic waves and the middle contact/shear wave, giving sharp shocks and good shear resolution. Recommended for flood applications.
"lf": Local Lax–Friedrichs (Rusanov). More diffusive but very robust; useful as a fallback.

Reconstruction¶

Config.recon sets how cell averages are reconstructed to face values, trading accuracy for cost and robustness:

`recon`	Order	Notes
`"first"`	1st	most robust; pairs with the well-balanced source for a consistent 1st-order scheme
`"muscl"`	2nd	MUSCL with slope limiting
`"linear2"`, `"linear3"`	2nd/3rd	linear reconstructions
`"linear5"`	5th	high-order linear (needs `ngh>=3` in 2D, `ngh>=4` in 1D)
`"weno5"`	5th	WENO, shock-capturing

Warning

recon is used only when well_balanced=False. With the default well_balanced=True the surface-reconstruction method builds the face states itself, and they are first order by construction, so muscl and weno5 give the same answer to the last bit as first. A reconstruction study therefore sets well_balanced=False, as examples/ex07_convergence_order.py does; a Config that asks for both warns at construction.

Note

Mesh2D.ngh and Mesh1D.ngh default to 4, which covers every scheme. Below that, both solvers check the halo at construction against what will actually run, not against the reconstruction radius alone, and refuse one that is too narrow with the width it needs. Three bounds feed that check:

the reconstruction stencil, 1 for first, muscl, linear2 and linear3 and 3 for linear5 and weno5, which is all the NumPy path reads;
on the GPU, the guard inside the fused kernel the configuration dispatches to, which skips every cell outside its own halo. The SRM well-balanced kernels (the default, with either flux) need 2, the non-well-balanced Lax-Friedrichs family needs 3 with every reconstruction, "first" included, and first-order HLLC without the well-balanced source needs 1. A halo one layer too narrow used to run and leave the outermost interior rows and columns frozen, with water never leaving a "fall" boundary;
in 1D, the window the residual is written into, which is asymmetric and wider than the stencil: well_balanced=True needs 1 whatever recon says, and with well_balanced=False muscl, linear2 and linear3 need 2 while linear5 and weno5 need 4.

Tip

For real-terrain flood runs the production configuration is recon="first" with the well-balanced SRM source, which is robust on noisy DEMs and wet/dry fronts. Use muscl/weno5 for smooth academic test cases.

Time integration and the CFL condition¶

Config.time selects the integrator:

"euler": forward Euler (first order in time). Cheapest; the default for large production runs.
"ssprk3": three-stage strong-stability-preserving Runge–Kutta (third order). Best for smooth, accuracy-sensitive problems.

The stable time step is set by the CFL condition,

\[ \Delta t = \mathrm{CFL}\,\frac{\min(\Delta x, \Delta y)}{\max_{h \ge h_\min}\big(V + \sqrt{g h}\big)}, \]

computed by cfl_dt() over the wet cells (a global all-reduce under MPI). The velocity norm \(V\) is \(\max(|u|,|v|)\) in the dense solver and the more conservative \(\sqrt{u^2+v^2}\) in the compressed solver (cfl_linf=True or SWE_CFL_LINF=1 selects the former there). Config.cfl is the Courant number; every run in the paper uses 0.5. run() takes these steps for you, and run(t_end, dt_max=...) caps them. To control the loop yourself, call step(dt):

s.step(dt=min(s.cfl_dt(), dt_max))

Note

On a fully dry domain no wave speed limits the step: the wave speed falls back to \(\sqrt{g h_\min}\) and cfl_dt() returns minutes to hours. While it rains, run therefore also keeps the step below the CFL step of the film that the rain lays down during the step, \((\mathrm{CFL}\,\Delta x)^{2/3}/(g R)^{1/3}\) for a rain rate \(R\). On wet ground this bound is far above the CFL step and changes nothing. Two loops do not get it: a hand-written step(cfl_dt()) loop, and run on more than one rank, where the ranks see different rain and must keep a common step. Cap the first steps yourself there.

Operator splitting¶

One forward-Euler step, in order:

the CFL time step over the wet cells;
ghost-cell fill and, under MPI, the halo exchange;
the residual \(\mathcal{R}\): HLLC fluxes plus the well-balanced bed-slope source;
the explicit update \(\mathbf{q}^{*} = \mathbf{q}^{n} + \Delta t\,(\mathcal{R} + \mathbf{S}_r)\), with rainfall included;
point-implicit friction on \(\mathbf{q}^{*}\), under the wet/dry floor;
relaxation and imposition of boundary values (sponge, then the coastal stage ring or StageBoundary);
the optional depth sinks (Green-Ampt infiltration, uniform recession, karst cap).

Friction is unconditionally stable in its point-implicit form, and the sinks are algebraic updates of \(h\) that scale the momentum by the remaining-depth ratio. On the GPU, steps 3 to 5 run as one fused kernel.

Robustness on real terrain¶

Production DEMs are noisy and create extreme states at pits, curbs, bridge decks, and bathymetry seams. GeoSWE guards against them with the dry-state limits of the wet/dry floor, the dry-bed wave speeds in the HLLC solver, the \(\sqrt{g h_\min}\) fallback in the CFL reduction, and the velocity cap in the friction step (friction_velocity_cap_ms, 15 m/s in every reported run). These are what let the same solver run a clean dam break and a continental DEM without retuning.