NumPy backend#
The NumPy backend is the reference implementation every other backend is
verified against. It ships inside zoomy_core at library/zoomy_core/zoomy_core/fvm/
— there is no separate zoomy_numpy package. Same Model / SystemModel API.
Pure Python / NumPy / SciPy; no external solver dependency.
Cell-centred finite-volume; symbolic Riemann solver; NCP path integral.
Solvers are
param.Parameterized— every knob is class-level and auto-introspected by the GUI / CLI.
Install#
pip install zoomy_core
Solver inventory#
Class |
File |
What it adds |
|---|---|---|
|
|
Explicit RK1 / SSP-RK2 flux + symbolic Riemann + source. The reference. |
|
|
Positive (hydrostatic-reconstruction) Rusanov + η=h+b well-balanced reconstruction + wet/dry. Needs |
|
|
Free-surface with path-conservative Roe ` |
|
|
Explicit flux + implicit diffusion + implicit source (Newton / Newton-GMRES). |
|
|
IMEX with the free-surface positive Rusanov flux. |
|
|
Chorin pressure projection: hyperbolic → viscous diffusion → Poisson correction. |
|
|
Pressure-splitting for free-surface flow (velocities at |
|
|
IMEX on extruded |
|
|
Monolithic index-1 DAE-PDE via Ascher-Ruuth-Spiteri IMEX-RK (ARS232 / ARS343). Order-1 correctness reference for the VAM chain. |
|
|
Chorin split of the VAM DAE into predictor + elliptic pressure + corrector sub-systems (from |
Common base: Solver (fvm/solver_numpy.py:124, a param.Parameterized) holds
the settings parameter, state allocation, runtime construction, and the
aux-refresh (it fills the local aux formula leg plus the LSQ-gradient leg from
sm.aux_registry); every solver inherits it. Within NumPy the implicit solvers
form a subclass chain ColumnIntegratingSolver ← SplittingSolver ← IMEXSolver,
and the free-surface variants override only the numerics/reconstruction build.
Operator-splitting stages#
[ explicit hyperbolic flux + source ] -- RK1 / SSP-RK2 on F, B·∂_x Q, S_e
▼
[ implicit / explicit diffusion ] -- A:∇Q (optional, IMEX)
▼
[ pressure Poisson correction ] -- free-surface / incompressible
▼
[ state update ] -- update_q, update_qaux
HyperbolicSolver runs stage 1; IMEXSolver adds stage 2; SplittingSolver
adds the Chorin correction in stage 3; DAESolver collapses everything into
a monolithic stage Newton (singular mass_matrix allowed);
ChorinSplitVAMSolver runs predictor → elliptic → corrector on three
sub-system models from split_for_pressure.
Tunable param knobs#
A solver carries only time-stepping knobs. Numerical knobs
(reconstruction order, limiter, eigenvalue regularisation) are not solver
params — they live on the NumericalSystemModel so the same solver can run any
discretisation. Passing them to a solver constructor raises TypeError.
# HyperbolicSolver — time stepping only
time_end = param.Number(default=0.1, bounds=(0, None))
min_dt = param.Number(default=1e-6, bounds=(0, None))
compute_dt = param.Parameter(default=None) # default: timestepping.adaptive(CFL=0.9)
# Numerical knobs — on the NumericalSystemModel specs, NOT the solver:
ReconstructionSpec(order=2, limiter="venkatakrishnan") # 1=PWC, 2=MUSCL
RegularizationSpec(eigenvalue_eps=1e-8)
# IMEXSolver (plain class attrs, not yet param-typed)
source_mode = "auto"; implicit_tol = 1e-8; implicit_maxiter = 8
gmres_tol = 1e-7; gmres_maxiter = 40; jv_backend = "analytic"
# SplittingSolver
pressure_gmres_tol = param.Number(default=1e-6)
pressure_gmres_maxiter = param.Integer(default=200)
viscosity = param.Number(default=0.01)
DAESolver: method ("ars232" | "ars343"), newton_tol, newton_maxit,
h_index, nc_integration_order, well_balanced, jacobian_mode
("sparse_fd" | "dense_fd"). ChorinSplitVAMSolver: pressure_tol,
pressure_maxit, time_order, pressure_solver ("gmres" | "lu"),
riemann_solver ("hr" | "ncp"). Aux-derivative limiting (for a shock
that overshoots into the elliptic/source block) is model-declared via
zoomy_core.model.numerics.limit (a limited_derivative tag), not a solver knob.
How a solver consumes a Model / SystemModel#
Every setup_simulation(mesh, model, write_output=True) normalises once:
# inside Solver.setup_simulation (paraphrased)
mesh = ensure_lsq_mesh(mesh, model)
if not isinstance(model, SystemModel):
model = SystemModel.from_model(model) # freeze operators
self.sm = model
runtime = NumpyRuntimeModel.from_system_model(model) # lambdify operators
The runtime is the only thing the timestep loop touches: routed solver-tag
callables (flux, nonconservative_flux, hydrostatic_pressure, source /
source_explicit, diffusion_matrix, boundary_conditions,
boundary_gradients, update_variables). See
../authoring/system-model.md and
../authoring/numerics.md for tag routing.
End-to-end example: shallow water dam break#
Shallow water is SME(level=0): the level-0 member of the shallow-moment
hierarchy, derived by the same operators as every other level — no hand-written
system anywhere. A 1-D case runs the dimension=2 model on a 1-D mesh; the
state is then [b, h, q_0].
import numpy as np
import zoomy_core.model.boundary_conditions as BC
import zoomy_core.model.initial_conditions as IC
from zoomy_core.fvm import timestepping
from zoomy_core.fvm.solver_numpy import HyperbolicSolver
from zoomy_core.mesh import LSQMesh
from zoomy_core.model.models import SME
from zoomy_core.model.models.closures import ManningFriction, Newtonian, StressFree
from zoomy_core.numerics import NumericalSystemModel, ReconstructionSpec
from zoomy_core.systemmodel import SystemModel
# Shallow water is the level-0 member of the shallow-moment hierarchy.
model = SME(
level=0, dimension=2,
closures=[Newtonian(), ManningFriction(), StressFree()],
boundary_conditions=BC.BoundaryConditions(
[BC.Extrapolation(tag="left"), BC.Extrapolation(tag="right")]),
)
sm = SystemModel.from_model(model) # state [b, h, q_0]
def ic(x): # 2:1 dam break at x = 5
return np.array([0.0, 2.0 if x[0] < 5.0 else 1.0, 0.0])
sm.initial_conditions = IC.UserFunction(function=ic)
sm.aux_initial_conditions = IC.Constant(constants=lambda k: np.zeros(k))
# Numerical knobs go on the NumericalSystemModel, not the solver:
nsm = NumericalSystemModel.from_system_model(
sm, reconstruction=ReconstructionSpec(order=2, limiter="minmod"))
mesh = LSQMesh.create_1d(domain=(0.0, 10.0), n_inner_cells=300)
solver = HyperbolicSolver(time_end=0.5, compute_dt=timestepping.adaptive(CFL=0.9))
Q, Qaux = solver.solve(mesh, nsm, write_output=False)
LSQMesh.create_1d from zoomy_core.mesh is the in-package 1-D constructor;
richer unstructured 2-D / 3-D meshes come from zoomy_mesh (sibling subrepo at
library/zoomy_mesh/). The snippet runs as-is in the
conda environment or after
pip install zoomy_core.