Well-balanced schemes

A scheme is well-balanced if it preserves steady states exactly at the discrete level. The essential one for flood modeling is lake at rest: still water (\(\mathbf{u}=0\)) with a flat free surface \(\eta = h + b = \text{const}\) over arbitrary bathymetry must stay perfectly still. A naive discretization of the bed-slope source does not balance the hydrostatic pressure flux and spawns spurious currents that can swamp a real flood signal. Preserving that balance is the C-property.

Enable well-balancing with Config.well_balanced=True and choose the variant with Config.wb_method:

wb_method	Scheme
"audusse"	Audusse hydrostatic reconstruction (simpler)
"srm"	Xia (2017) Surface-Reconstruction Method, the production default for real terrain

from geoswe import Config
cfg = Config(well_balanced=True, wb_method="srm")    # both are the defaults

Warning

"audusse" is kept for the bed-discretisation comparison, not for terrain. With the default flux="hllc" no fused kernel is wired for it, so a GPU run falls back to the Python residual, which the dispatch records as about 100 times slower, and warns once:

RuntimeWarning: wb_method='audusse' with flux='hllc' takes the slow Python RHS path (the fused kernel is only wired for wb_method='srm')

The Audusse kernel that flux="lf" does reach has an unresolved wet/dry instability on real bathymetry: thin cells accumulate momentum and their velocity grows without bound as the depth falls. Use "srm" on real terrain and keep "audusse" for smooth cases.

Note

Well-balanced reconstruction is first-order in space; pair it with recon="first" for a consistent first-order well-balanced scheme. This is the combination used for the continental-scale runs, where DEM noise and wet/dry fronts make robustness paramount.

Verifying the C-property

examples/ex03_lake_at_rest_2d.py puts still water over a Gaussian bump and checks that no current develops:

eta0 = to_host(s.q_interior[0]) + bed          # the free surface before stepping
for _ in range(50):
    s.step(dt=s.cfl_dt())
hu = to_host(s.q_interior[1]); hv = to_host(s.q_interior[2])
eta1 = to_host(s.q_interior[0]) + bed
assert np.max(np.hypot(hu, hv)) < 1e-10        # no spurious momentum
assert np.max(np.abs(eta1 - eta0)) < 1e-10     # surface unchanged

With wb_method="srm" the residual currents are at machine precision: the example, which pins dtype="float64", reports a maximum momentum of 4.7e-15 and a surface drift of 1.3e-15 m on the CPU. This is also enforced by tests/test_well_balanced.py.