# 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 {file}`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 {py:meth}`~geoswe.Solver2D.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. {py:meth}`~geoswe.Solver2D.run` takes these steps for you, and `run(t_end, dt_max=...)` caps them. To control the loop yourself, call `step(dt)`:

```python
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:

1. the CFL time step over the wet cells;
2. ghost-cell fill and, under MPI, the halo exchange;
3. the residual $\mathcal{R}$: HLLC fluxes plus the well-balanced bed-slope source;
4. the explicit update $\mathbf{q}^{*} = \mathbf{q}^{n} + \Delta t\,(\mathcal{R} + \mathbf{S}_r)$, with rainfall included;
5. point-implicit **friction** on $\mathbf{q}^{*}$, under the wet/dry floor;
6. relaxation and imposition of boundary values (sponge, then the coastal stage ring or `StageBoundary`);
7. 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.