zoomy_core.fvm.solver_dae_numpy module

DAE solver (numpy backend) — index-1 DAE-PDE integrator via Ascher-Ruuth-Spiteri IMEX-Runge-Kutta time stepping.

Sits next to HyperbolicSolver and IMEXSolver in zoomy_core.fvm and reuses every shared piece of the framework:

• SystemModel is the single symbolic source-of-truth; SystemModel.from_model auto-scans all non-state Function and Derivative atoms into aux_state + a structured SystemModel.aux_registry.

• NumpyRuntimeModel lambdifies the operator matrices into (Q, Qaux, p) → ndarray callables.

• ensure_lsq_mesh() promotes the input mesh to LSQMesh, exposing compute_derivatives for the registry- driven update_qaux walker.

• Boundary-condition objects (Extrapolation, Lambda, InflowOutflow, …) supply face values via face_value.

• Solver.update_q / Solver.update_qaux (registry-aware default) fire after every step.

Time integration: Ascher-Ruuth-Spiteri IMEX-Runge-Kutta (ARS232 / ARS343) on the singular DAE M(Q)·∂_t Q = R(Q, Qaux, p). Aux is updated once per step (lagged through the Newton iteration), keeping the Jacobian small and reusing the standard Solver.update_qaux hook.

Spatial scheme: Rusanov flux + non-conservative path-integral fluctuations, mirroring NonconservativeRusanov but driven by the SystemModel runtime (so derivative-aux entries like ∂_x h flow through as ordinary Qaux components, no special-casing). Boundary faces use BC.face_value(Q_inner, Qaux_inner, normal, d_face, time, parameters).

class zoomy_core.fvm.solver_dae_numpy.DAESolver(**kwargs)

Bases: Solver

Index-1 DAE-PDE solver with IMEX-RK time stepping.

At reconstruction_order=1 this solver is the validated correctness reference for the VAM chain: lake-at-rest is preserved to machine precision and a perturbation propagates with bounded mass loss.

At reconstruction_order=2 the spatial scheme is also correct — lake-at-rest over a bump is well-balanced to ~1e-14 (the η = h+b SurfaceReconstruction + the cell-interior non-conservative integral telescope exactly), and the slope limiter is frozen through the Newton iteration so f_I is a smooth function of the stage unknown. But the time integration does not converge: the monolithic IMEX-ARK stage Jacobian is ill-conditioned (cond ~1e7), concentrated in the algebraic pressure-constraint rows. At that conditioning the finite-difference Jacobian (step ~1e-7) is unreliable, so the stage Newton degrades from quadratic to slow linear convergence and does not reach newton_tol within newton_maxit.

This is structural, not a bug: the monolithic DAE couples a well-conditioned hyperbolic evolution block and an ill-conditioned elliptic pressure-constraint block into a single FD Jacobian. The fix is the Chorin / projection split — an explicit hyperbolic predictor (where the order-2 reconstruction lives, no Jacobian needed) plus a separate linear elliptic pressure solve (which can be preconditioned properly). The split solver is the supported home for order 2; this class stays the order-1 reference.

time_end = 0.1

method = 'ars343'

newton_tol = 1e-09

newton_maxit = 30

h_index = 0

nc_integration_order = 3

reconstruction_order = 1

limiter = 'venkatakrishnan'

well_balanced = True

jacobian_mode = 'sparse_fd'

compute_dt = None

setup_simulation(mesh, model, *, write_output=False)

Build operators + state once, including BC plumbing, output infrastructure, and the registry-driven Qaux walker.

f_E(t, Y)

f_I(t, Y)

Implicit RHS: per-cell M_evol⁻¹·R_evol on evolution rows; constraint residual on algebraic rows.

Aux is held fixed at self._sim_Qaux (refreshed once per timestep) — the lagged-aux Newton convention.

J_I(t, Y)

project_to_manifold(Q, *, time=0.0, tol=1e-10, maxit=20)

Newton-project Q onto the algebraic constraint manifold f_I[alg] = 0 by adjusting algebraic state entries.

Aborts (returns the original Q) if the Newton residual does not decrease or if any iterate becomes non-finite — the caller can then proceed with the un-projected IC and let the IMEX-ARK Newton (which uses a more conservative update) do the projection.

step(dt)

One IMEX-ARK step on the stored state.

Convention: refresh self._sim_Qaux via update_qaux once per step (lagged through the Newton iteration), apply update_q to the new state after. The slope-limiter coefficients are frozen the same way — computed once from the lagged state so f_I stays smooth in the Newton unknown.

run_simulation()

solve(mesh, model, *, write_output=True)

export_vtk(*, field_names=None, aux_field_names=None, filename='dae_output', skip_aux=False)

Post-process HDF5 snapshots to a VTK time series for Paraview. Requires write_output=True at setup time.

name = 'DAESolver'