Solvers and configuration

The three classes a run is built from. All three are defined in geoswe.solver and re-exported at the top level, which is the form this page documents them in.

class geoswe.Config(pde='baseline', flux='hllc', recon='first', time='euler', cfl=0.5, bc_x='extrapolate', bc_y='extrapolate', bc_x_left=None, bc_x_right=None, g=9.81, h_min=1e-10, h_min_cfl=0.0, wb_method='srm', friction=None, manning_n=0.0, friction_velocity_cap_ms=15.0, friction_quadratic_alpha=True, storage_courant=0.0, storage_dt_ref=0.0, rainfall=0.0, rainfall_forcing=None, manning_field=None, stage_boundary=None, well_balanced=True, rk_storage='low_storage', dtype=None)[source]

Bases: object

Solver configuration.

Config() with no arguments is the production flood scheme of the GeoSWE paper: the shallow-water equations with first-order HLLC fluxes and the SRM well-balanced bed treatment (flux='hllc', recon='first', well_balanced=True, wb_method='srm'), forward Euler at cfl=0.5, and open boundaries. A flood run usually adds three things: friction (friction='manning'), rain (rainfall_forcing), and boundaries that let water out (bc_x='fall', bc_y='fall').

Fields a first script uses (the configuration reference lists them all):

  • friction – None (off) or 'manning'. The roughness is manning_n (one value), or a map given with Solver2D.set_manning().

  • rainfall_forcing – a RainfallForcing (mm/h over time, uniform or gridded); rainfall is the constant-rate short form in m/s.

  • stage_boundary – a StageBoundary: a water level over time on chosen cells (tide, surge, a river stage).

  • bc_x, bc_y – what happens at the edges of the rectangle: 'extrapolate' (zero gradient), 'wall', 'fall' (free outflow), 'periodic'; 1D also has 'dirichlet' with bc_x_left and bc_x_right.

  • dtype – 'float32' or 'float64'; left unset it follows the backend (float32 on the GPU, float64 on the CPU).

  • h_min – the wet/dry depth threshold; g – gravity; cfl – the Courant number.

Scheme choices, for method studies:

  • flux – 'hllc' or 'lf'; time – 'euler' or 'ssprk3'.

  • recon – 'first', 'muscl', 'linear2', 'linear3', 'linear5', or 'weno5'. Use a higher-order recon together with well_balanced=False when measuring reconstruction order: the well-balanced face states are first-order by design.

  • well_balanced, wb_method – 'srm' (default) or 'audusse'.

Parameters:
pde: str = 'baseline'
flux: str = 'hllc'
recon: str = 'first'
time: str = 'euler'
cfl: float = 0.5
bc_x: str = 'extrapolate'
bc_y: str = 'extrapolate'
bc_x_left: ndarray | None = None
bc_x_right: ndarray | None = None
g: float = 9.81
h_min: float = 1e-10
h_min_cfl: float = 0.0
wb_method: str = 'srm'
friction: str | None = None
manning_n: float = 0.0
friction_velocity_cap_ms: float = 15.0
friction_quadratic_alpha: bool = True
storage_courant: float = 0.0
storage_dt_ref: float = 0.0
rainfall: float = 0.0
rainfall_forcing: object | None = None
manning_field: ndarray | None = None
stage_boundary: object | None = None
well_balanced: bool = True
rk_storage: str = 'low_storage'
dtype: str | None = None
class geoswe.Solver2D(mesh, cfg, q0, b, comm=None, dims=None)[source]

Bases: object

Two-dimensional shallow-water solver on a structured (dense) grid.

This is the reference dense path: it allocates the full padded rectangle (3, nx + 2*ngh, ny + 2*ngh) for (h, hu, hv), runs the SRM-HLLC residual (fused CUDA kernel on GPU, NumPy on CPU), the point-implicit Manning friction, rainfall and the optional depth sinks, and supports per-face boundary conditions, the gauge-driven coastal ring, and MPI domain decomposition with a two-cell halo. The compressed active-cell path (geoswe.compressed_solver) reuses its kernels and reproduces its residual bitwise.

Parameters:
set_inside_mask(mask, cfl_robust_pct=None)[source]

Choose the cells that are updated: the active set.

mask is a boolean array of the unpadded shape (nx, ny); None switches the mask off. Cells where it is True are updated as usual. Cells where it is False get a zero residual and keep the state they have, which is how open ocean or terrain outside the study area is left out. The time step still looks at every wet cell (set_cfl_ghost_mask() excludes cells from it). This is also the mask that geoswe.CompressedSolver.from_dense() packs.

cfl_robust_pct (for example 99.99) replaces the largest wave speed in the time-step limit by that percentile over this rank’s cells, which drops a few outlier cells at wet/dry fronts; None keeps the strict maximum. It is single-rank only and refused under MPI, where the maximum of the per-rank percentiles is not the global percentile; use set_cfl_ghost_mask() there. Setting it also drops the run off the fused fp32 time-step kernel and off the fused dense step.

set_manning_table(cls_padded, table)[source]

Give Manning’s n as a class index per cell plus a small table.

cls_padded is a uint8 array of the padded shape (nx + 2*ngh, ny + 2*ngh) holding each cell’s class, and table lists Manning’s n of each class. The friction kernel then reads n = table[cls]: the same result as a full field with those values, at one byte per cell instead of four. Leave Config.manning_field unset when you use the table. set_manning() builds both from an ordinary (nx, ny) array and is the simpler call.

set_manning(n)[source]

Set Manning’s roughness from a scalar or an unpadded (nx, ny) array.

The convenient form of Config.manning_field and set_manning_table(): the ghost cells are filled for you, and a float32 GPU run stores a field with at most 256 distinct values as a one-byte class table, the form the single-kernel time step reads. Friction is switched on if the configuration left it off.

depth()[source]

Water depth h as a NumPy array of shape (nx, ny), ghost cells removed.

Under MPI this is the calling rank’s block.

max_depth()[source]

Largest depth each cell has held after any step so far, as a NumPy (nx, ny) array.

This is the running maximum behind an inundation map. Before the first step it equals depth().

set_cfl_ghost_mask(mask)[source]

Mark interior cells as ghost cells for the CFL reduction.

Cells flagged here (mask True) are skipped by the fused lam_max kernel in cfl_dt. Use this for BC cells (e.g. Dirichlet / open-h ring cells) whose h, hu, hv are imposed externally each step, including them in the dt reduction artificially shrinks dt whenever the imposed state has a large lam (deep h + nonzero u from interior extrapolation).

Parameters:

mask (ndarray of bool, or None to clear.) – Interior shape (nx, ny) (no ghost). Stored flattened in C-order (matches the kernel’s per-cell index over q_int = q[:, ngh:-ngh, ngh:-ngh]).

set_storage_fraction(sigma)[source]

Sub-grid channel storage scaling for narrow-flowline cells.

sigma is an array of interior shape (nx, ny) with values in (0, 1]. For cells where σ < 1, the Euler update for h is scaled: dh/dt = rhs_h / σ. Physically the cell’s wet area is only σ·dx² (a narrow channel within the cell), so the same volumetric flux raises h by 1/σ.

Mass is conserved across the channel↔floodplain interface: per face, ΔV_channel = σ·dx²·(rhs/σ)·dt = rhs·dx²·dt matches ΔV_floodplain. (An earlier docstring asserted mass cons was violated; that was wrong; volume balance works out exactly.)

STABILITY: the 1/σ scaling on ∂h/∂t effectively requires dt <= sigma * dx / (abs(u) + sqrt(g*h)) in each σ<1 cell. Pass this 1/σ field to the CFL kernel via self._storage_inv_sigma so the per-cell lam is multiplied by 1/σ, tightening dt only where storage is narrow. Without this per-cell CFL the simple heuristic blows up (h overshoots in one step, c_sound explodes, NaN). Verified empirically on Pinellas v29: CFL=0.5 NaNs; CFL≤σ_min=0.167 runs stably for 72h.

CAVEAT: this scales mass storage only, not the face flux. For uniform-σ regions the scheme is consistent (the same 1/σ falls out of both flux and storage); at σ-jumps (channel↔floodplain) the momentum equation is off by σ relative to a full porosity-based SWE (Casulli 2009 / Sanders 2008). Empirically this manifests as a small velocity bias at C↔F interfaces, likely smaller than the gauge calibration error it’s trying to fix, but worth checking with the bench validation before relying on it.

Pass None to disable.

add_inflow(idx_i, idx_j, normal, ds, t_series, q_series, dry_frac=0.1)[source]

Add one discharge (hydrograph) inlet. Call once per river.

normal is the INWARD unit normal (nx, ny); ds the cell width across the inlet face; t_series/q_series the hydrograph knots (s, m^3/s), linearly interpolated and held flat outside their range. Indices are into the PADDED arrays. See bc.apply_inflow_discharge for the distribution rule and the dry-inlet convention.

Each inlet keeps its OWN normal, width and hydrograph, and its discharge is weighted only against its own cross-section – lumping several rivers into one call would share a single Q and a single normal between them.

set_inflow(*a, **k)[source]

Replace all inlets with a single one (convenience for one river).

set_step_forcings(spec)[source]

Hand the driver’s post-step forcings to the dense fused step (opt-in).

spec is None (off) or a dict with the band sponge (sponge: keep, amb, x0, y0, w, do_x, do_y; or None) and Green-Ampt/drain (gd: mode 1 GA, 2 drain, 3 both; cls, Ks, psi, dth, F, Fmax or None, inv_tau; or None), plus ring_mask, an interior uint8 mask of the cells whose infiltration/drain the caller applies after its ring update. When the fused step applies them it sets self._step_forcings_done for that step; otherwise the caller runs its own kernels. Results are bit-identical to the separate kernels.

fold_cells_into_next_cfl(cells)[source]

Add the current lambda of cells (padded linear indices, int32) to the next step’s CFL reduction when the fused step with forcings reduced it (set_step_forcings cfl=True).

cfl_dt()[source]

Return the global CFL-limited time step.

The per-cell wave speed uses the configured velocity norm and the CFL depth floor; ring/boundary cells flagged by set_cfl_ghost_mask() are skipped, dry partitions fall back to sqrt(g*h_min), and under MPI the maximum is reduced across ranks so every rank advances with the same dt.

step(dt=None)[source]

Advance one time step of size dt (CFL step if None).

Order of operations (Algorithm 1 of the paper): fill ghost cells and exchange halos, evaluate the residual, integrate with rainfall, apply the point-implicit friction, impose boundary/ring values, then the optional depth sinks and the running depth maximum.

Parameters:

dt (float | None)

run(t_end, max_steps=10000000, callback=None, dt_max=None)[source]

Step to t_end with CFL-sized steps, clipping the last step onto t_end.

callback(solver, step) runs after each step, and dt_max caps the step size. While it rains, the step is also kept below the CFL step of the film that the rain lays down during the step, (cfl*dx)**(2/3) / (g*R)**(1/3) for a rain rate R, so rain on a dry bed, where no wave speed limits the step, is not deposited minutes at a time. On wet ground that bound is far above the CFL step and changes nothing. Returns the number of steps taken in this call; self.nsteps accumulates across calls.

Parameters:
property q_interior

Conserved state without the ghost padding, shape (3, nx, ny).

class geoswe.Solver1D(mesh, cfg, q0, b)[source]

Bases: object

One-dimensional shallow-water solver (CPU, development and verification).

Holds the padded conserved state q of shape (2, nx + 2*ngh) (depth and unit discharge) and the bed b, and advances it with the flux, reconstruction, time integrator and friction selected in cfg. The 1D path exists for convergence tests and the Ritter/dam-break checks; production runs use Solver2D or the compressed mesh.

Parameters:
cfl_dt()[source]

Return the CFL-limited time step cfg.cfl * dx / max(|u| + c) over the interior.

step(dt=None)[source]

Advance one time step of size dt (CFL step if None): fill ghosts, evaluate the residual, integrate, apply friction.

Parameters:

dt (float | None)

run(t_end, max_steps=10000000, callback=None, dt_max=None)[source]

Step to t_end with CFL-sized steps (the last one clipped to land exactly on t_end); returns the padded state.

dt_max caps the step size, and rain keeps the step below the CFL step of the film it lays down (see Solver2D.run()).

Parameters:
property q_interior

Conserved state without the ghost padding, shape (2, nx).