zoomy_core.fvm.riemann_solvers module

zoomy_core.fvm.riemann_solvers module#

Symbolic Riemann solvers: Rusanov, positive, nonconservative variants.

Every index-dependent class (PositiveRusanov, NonconservativeRusanov, PositiveNonconservativeRusanov, …) auto-locates the named fields h and b (and hinv when present) in model.state / model.aux_state via FieldHandle so the same Riemann code path works whether bathymetry is part of the conservative state (legacy SWE) or lives in Qaux (chain-DAE convention).

class zoomy_core.fvm.riemann_solvers.FieldHandle(name, container, index, minus, plus, state)#

Bases: object

A named field resolved to either Q (state) or Qaux.

Looks up the location at construction time via Numerics.find_field(), then exposes the symbolic Symbol on the minus / plus / state side of a face so callers can write reconstruction code without knowing where the field lives.

name#

Field name (e.g. "h", "b", "hinv").

Type:: str

container#

Which array carries this field.

Type:: {“q”, “qaux”}

index#

Index into that array.

Type:: int

minus, plus, state

Direct references — minus and plus are the per-face symbolic state on the L / R side of a Rusanov face; state is the cell-centre reference. Using these in symbolic code produces the same lambdified output as Q[index] / Qaux[index] — a true placeholder.

Type:: sympy.Symbol

name#

container#

index#

minus#

plus#

state#

access(q_array, qaux_array)#: Return q_array[index] or qaux_array[index] depending on container. Works for symbolic ZArray and numeric numpy arrays alike.

assign(q_array, qaux_array, value)#: In-place write into the appropriate container.

class zoomy_core.fvm.riemann_solvers.Numerics(model, **params)#

Bases: Parameterized, SymbolicRegistrar

Symbolic numerics over a SystemModel.

The numerics consumes a SystemModel — the frozen operator-form snapshot of a derivation. A Model passed to the constructor is normalised once via SystemModel.from_model(); everything internal reads the SystemModel’s stored operators (flux, hydrostatic_pressure, eigenvalues, nonconservative_matrix, quasilinear_matrix) and its state / aux_state / parameters / normal.

Subclasses that depend on the location of named fields (h, b, hinv, …) call find_field() to obtain a FieldHandle — the search walks state then aux automatically, so the same numerics works whether bathymetry is in the state or in aux_state.

name = 'NumericsV2'#

model = None#

scaled_q_indices = None#

find_field(name, *, required=True)#

Return a FieldHandle for the named field.

Searches self.model.state (Q) first, then self.model.aux_state (Qaux). Caches the result in self._field_handles[name] so repeat lookups are free.

Parameters:

name (str) –
required (bool) – If True (default), raise KeyError when the field is in neither container. Otherwise return None.

has_field(name)#

local_max_eigenvalue_definition()#: Returns the opaque max_wavespeed function. The actual implementation is provided by the backend at runtime.

local_max_abs_eigenvalue(Q=None, Qaux=None, p=None, n=None)#: Called during symbolic Rusanov construction. Returns opaque max_wavespeed with the given state.

to_runtime(backend='numpy')#: To runtime.

to_runtime_numpy()#: To runtime numpy.

to_runtime_ufl()#: Lambdify all registered symbolic functions through the UFL module dict — for use by Firedrake-based backends.

numerical_flux()#: Numerical flux.

numerical_fluctuations()#: Numerical fluctuations.

class zoomy_core.fvm.riemann_solvers.Rusanov(model, **params)#

Bases: Numerics

Rusanov. (class).

name = 'RusanovV2'#

get_viscosity_identity_flux()#: Get viscosity identity flux.

get_viscosity_identity_fluctuations()#: Get viscosity identity fluctuations.

numerical_flux()#: Numerical flux.

class zoomy_core.fvm.riemann_solvers.HLL(model, **params)#

Bases: Numerics

HLL (Harten-Lax-van-Leer) approximate Riemann solver.

Model-agnostic: needs only the SystemModel’s flux, hydrostatic_pressure and eigenvalues operators. When the SystemModel carries a symbolic spectrum the wave-speed bounds are the Davis estimates — min / max over the eigenvalues of both face states. When eigenvalues is None (the model skipped the spectral derivation) it falls back to ± local_max_abs_eigenvalue, i.e. HLL collapses to local Lax-Friedrichs (a valid, more diffusive HLL).

The numerical flux is a single closed-form (branch-free) SymPy expression — clamping the wave speeds with Min(s_L, 0) / Max(s_R, 0) recovers the upwind branches without Piecewise, so it codegens cleanly to every backend.

name = 'HLLV2'#

numerical_flux()#: Numerical flux.

wave_speed_bounds(qL, qR, auxL, auxR, p, n)#: Return (s_L, s_R) — slowest / fastest signal speeds at the face.

class zoomy_core.fvm.riemann_solvers.HLLC(model, **params)#

Bases: HLL

HLLC approximate Riemann solver for the free-surface family.

Restores the contact / shear wave that HLL smears. Requires a depth field h (resolved via Numerics.find_field()); the momentum block is the first model.dimension depth-scaled Q rows. Any further depth-scaled rows (higher moments) are advected by the contact wave; non-scaled rows (e.g. bed in conservative state) pass through unchanged. Models without an h field should use HLL instead.

Region selection (F_L | F_L* | F_R* | F_R) uses the opaque conditional primitive, so it codegens to np.where / ternary expressions on every backend.

name = 'HLLCV2'#

property h_field#: Depth FieldHandle — raises KeyError if the model has no h field (such models should use HLL).

class zoomy_core.fvm.riemann_solvers.PositiveRusanov(model, **params)#

Bases: Rusanov

PositiveRusanov. (class).

Hydrostatic reconstruction follows Audusse-Bristeau-Klein: h_L* = max(0, h_L + b_L − b*), b* = max(b_L, b_R). The h, b and (optional) hinv fields are resolved via Numerics.find_field() so the same logic works whether bathymetry is part of conservative state or lives in Qaux.

name = 'PositiveRusanovV2'#

property h_field#: Depth FieldHandle.

property b_field#: Bathymetry FieldHandle.

property hinv_field#: Optional 1/h FieldHandle (None if absent).

hydrostatic_reconstruction(qL, qR, auxL, auxR)#: Hydrostatic reconstruction.

numerical_flux()#: Numerical flux.

numerical_fluctuations()#: Numerical fluctuations.

class zoomy_core.fvm.riemann_solvers.NonconservativeRusanov(model, **params)#

Bases: Rusanov

NonconservativeRusanov. (class).

name = 'NonconservativeRusanovV2'#

integration_order = 3#

get_path_integral_states()#: Get path integral states.

get_viscosity_identity_flux()#: Get viscosity identity flux.

get_viscosity_identity_fluctuations()#: Get viscosity identity fluctuations.

numerical_fluctuations()#: Numerical fluctuations.

class zoomy_core.fvm.riemann_solvers.WellBalancedNonconservativeRusanov(model, **params)#

Bases: NonconservativeRusanov

Path-conservative Rusanov with equilibrium-variable fluctuation viscosity for free-surface (lake-at-rest) well-balancing.

Both the bed-slope g·h·∂_x b and the hydrostatic pressure g·h·∂_x h live in the nonconservative product (the “Malaga” formulation; hydrostatic_pressure is empty, b is a trivial conserved state with ∂_t b = 0). The base NonconservativeRusanov.get_viscosity_identity_fluctuations() already zeros the stationary-bed row. This subclass additionally couples the depth-continuity row to the bed column so the Rusanov dissipation acts on the free-surface jump Δη = Δh + Δb instead of Δh. At lake-at-rest Δη = 0 while Δh = −Δb ≠ 0, so only the coupled form vanishes — giving exact well-balancing for both Rusanov and HLL-flavoured path integrals.

For Q = [b, h, hu] the fluctuation viscosity becomes:

[[0, 0, 0],
 [1, 1, 0],     # continuity dissipates on Δη = Δh + Δb
 [0, 0, 1]]

Model-derived — the coupling is added only when both h and b resolve to conservative-state fields via Numerics.find_field(). A plain SWE model with no bed gets the unmodified identity (no bed-indexed term at all), as required.

name = 'WellBalancedNonconservativeRusanovV2'#

get_viscosity_identity_fluctuations()#: Get viscosity identity fluctuations (equilibrium-coupled).

class zoomy_core.fvm.riemann_solvers.PositiveHLL(model, **params)#

Bases: HLL

HLL with Audusse-Bristeau-Klein hydrostatic reconstruction.

Mirrors PositiveRusanov but uses the sharper HLL two-wave numerical flux underneath instead of LF/Rusanov dissipation. Recommended for free-surface dam-break / wet-dry simulations: the hydrostatic reconstruction enforces h_face ≥ 0 (positivity), and HLL captures the rarefaction / shock fronts more accurately than Rusanov on the same mesh.

Reconstruction (same as PositiveRusanov, Audusse-Bristeau-Klein): b* = max(b_L, b_R), h_L* = max(0, h_L + b_L − b*), momentum scaled by h_L* / max(h_L, eps).

The depth field h, bathymetry b and (optional) 1/h inverse are resolved through Numerics.find_field() — the same flux code works whether bathymetry is part of the conservative state or carried in Qaux.

name = 'PositiveHLLV2'#

property h_field#

property b_field#

property hinv_field#

hydrostatic_reconstruction(qL, qR, auxL, auxR)#: Hydrostatic reconstruction.

numerical_flux()#

HLL flux evaluated on the hydrostatically-reconstructed face states — positivity-preserving and well-balanced under the lake-at-rest steady state.

The bed row is exactly zero by construction: HLL’s _state_jump zeros dq[b] when b lives in the conservative state, and a well-posed SWE has F[b, :] = 0 symbolically.

numerical_fluctuations()#

Audusse 2004 well-balancing consistency source S̃.

Audusse-Bouchut-Bristeau-Klein-Perthame 2004, SIAM J. Sci. Comput. 25(6):2050-2065, eq. (2.17)-(2.18):

F_{i+1/2,L} = F_num(U*_L, U*_R) - S̃_L F_{i+1/2,R} = F_num(U*_L, U*_R) + S̃_R

S̃_L = (0, ½g h_L² - ½g h*_L²)^T = (0, P_raw_L - P*_L)^T S̃_R = (0, ½g h*_R² - ½g h_R²)^T = (0, P*_R - P_raw_R)^T

The HR’d numerical flux alone is not well-balanced; the per- cell-per-face S̃ correction restores consistency with the original SWE (paper, Theorem 2.5). S̃ depends only on hydrostatic_pressure and is independent of the underlying Riemann solver — same formula already in PositiveRusanov.numerical_fluctuations().

super().numerical_fluctuations() chains the next entry on the MRO. For plain PositiveHLL it lands at Numerics.numerical_fluctuations() (returns zero). For PositiveNonconservativeHLL it lands at NonconservativeRusanov.numerical_fluctuations() (the DLM path-integral over the NCP) — so the same override gives well-balanced HLL whether the model carries an NCP or not.

class zoomy_core.fvm.riemann_solvers.PositiveNonconservativeRusanov(model, **params)#

Bases: PositiveRusanov, NonconservativeRusanov

PositiveNonconservativeRusanov. (class).

name = 'PositiveNonconservativeRusanovV2'#

get_path_integral_states()#: Get path integral states.

class zoomy_core.fvm.riemann_solvers.PositiveNonconservativeHLL(model, **params)#

Bases: PositiveHLL, NonconservativeRusanov

HLL conservative flux + Audusse-Bristeau-Klein hydrostatic reconstruction + path-integral NCP fluctuations.

Combines the sharper HLL two-wave numerical flux (vs. Rusanov / LF) with the same well-balanced reconstruction and path-integral NCP fluctuations as PositiveNonconservativeRusanov. The bed row is automatically excluded from the LF-style fluctuation dissipation via the bed mask in NonconservativeRusanov.get_viscosity_identity_fluctuations().

The Python MRO (PositiveNonconservativeHLL → PositiveHLL → HLL → NonconservativeRusanov → Rusanov → Numerics) resolves numerical_flux to PositiveHLL.numerical_flux (hydrostatic reconstruction → HLL combine → bed-row mask) and numerical_fluctuations to NonconservativeRusanov. numerical_fluctuations (path-integral NCP + LF identity, bed masked) — exactly the split we want.

name = 'PositiveNonconservativeHLLV2'#

get_path_integral_states()#: NCP path-integral evaluated on the hydrostatically- reconstructed face states.

class zoomy_core.fvm.riemann_solvers.QuasilinearRusanov(model, **params)#

Bases: NonconservativeRusanov

QuasilinearRusanov. (class).

name = 'QuasilinearRusanovV2'#

numerical_flux()#: Numerical flux.

class zoomy_core.fvm.riemann_solvers.PositiveQuasilinearRusanov(model, **params)#

Bases: PositiveRusanov, QuasilinearRusanov

PositiveQuasilinearRusanov. (class).

name = 'PositiveQuasilinearRusanovV2'#

numerical_flux()#: Numerical flux.

get_path_integral_states()#: Get path integral states.