tidal.solver package #

classmethod from_flat(layout, grid_shape, flat)[source]#

Wrap an existing flat array (no copy).

The caller must ensure the array is not unexpectedly mutated elsewhere.

Parameters:

layout (StateLayout)
grid_shape (tuple[int, ...])
flat (np.ndarray)

Return type:

property grid_shape: tuple[int, ...]#: Spatial grid shape.

property layout: StateLayout#: The state layout this FieldSet wraps.

max_norm()[source]#

Maximum absolute value across all slots.

Return type:: float

momenta_dict()[source]#

Alias for velocities_dict (transition aid).

Return type:: dict[str, ndarray]

property momentum_names: tuple[str, ...]#: Alias for velocity_names (transition aid).

rebind(flat)[source]#

Replace the underlying flat array reference (zero-copy, zero-alloc).

Used by solver loops to reuse a single FieldSet instead of allocating a new one each timestep. Auxiliary fields are cleared.

Parameters:: flat (ndarray)
Return type:: None

set_aux(name, data)[source]#

Used to inject constraint velocities from IDA’s yp vector.

Parameters:

name (str)
data (ndarray)

Return type:

None

property slot_names: tuple[str, ...]#: All slot names in order.

velocities_dict()[source]#

Return dict of velocity name → grid-shaped view (zero-copy).

Return type:: dict[str, ndarray]

property velocity_names: tuple[str, ...]#

Ordered velocity slot names (e.g. ("v_phi_0",)).

Only present for second-order equations.

classmethod zeros(layout, grid_shape)[source]#

Create a FieldSet initialized to zero.

Parameters:

layout (StateLayout)
grid_shape (tuple[int, ...])

Return type:

class tidal.solver.GridInfo(bounds, shape, periodic, bc=None, axis_bcs=None)[source]#

Bases: object

Immutable Cartesian grid descriptor.

Parameters:

bounds (tuple of (lo, hi) pairs) – Domain bounds per spatial axis, e.g. ((0, 10), (0, 10)).
shape (tuple of int) – Number of grid cells per axis, e.g. (64, 64).
periodic (tuple of bool) – Whether each axis is periodic, e.g. (True, True).
bc (tuple of str, optional) – Legacy per-axis BC type strings (e.g. ("neumann", "periodic")).
axis_bcs (tuple of AxisBCSpec, optional) – Structured per-axis BC specs with per-side values and Robin support. Takes precedence over bc when set.

Examples

>>> g = GridInfo(bounds=((0, 10),), shape=(64,), periodic=(True,))
>>> g.ndim
1
>>> g.dx
(0.15625,)
>>> g.num_points
64

axes_coords(axis)[source]#

1-D array of cell centres along axis.

Returns shape (shape[axis],) — a single 1-D vector, not a broadcasted grid. Useful for building per-axis coordinate arrays without the overhead of full meshgrid.

Raises:: IndexError – If axis is out of range.
Parameters:: axis (int)
Return type:: ndarray

axis_bcs: tuple[AxisBCSpec, ...] | None = None#

bc: tuple[str, ...] | None = None#

property bc_types: tuple[str, ...]#

Per-axis BC type as simple strings.

Always returns tuple[str, ...] — one of "periodic", "neumann", "dirichlet", or "robin" per axis.

For structured AxisBCSpec entries, uses the low-side kind as the representative type (matching the solver’s ghost cell convention). This is the canonical form stored in metadata and used by the energy module for BC-aware gradient computation.

property cell_coords: ndarray#

Cell-centre coordinates, shape (*shape, ndim).

For a 1D grid with bounds (0, 10) and shape (4,), the cell centres are at [1.25, 3.75, 6.25, 8.75] (centred in each cell).

Unlike py-pde’s CartesianGrid.cell_coords (which returns a list of arrays with inconsistent shapes), this always returns a single ndarray with the coordinate dimension as the last axis.

Cached on first access (GridInfo is immutable).

coord_arrays()[source]#

Broadcasted coordinate arrays, one per axis.

Returns ndim arrays, each of shape self.shape, containing the coordinate value along that axis at every grid point. Equivalent to np.meshgrid(*1d_axes, indexing="ij").

Cached on first access (GridInfo is immutable).

Return type:: tuple[ndarray, …]

property dx: tuple[float, ...]#

dx = (hi - lo) / N).

Type:: Grid spacing per axis (cell-centred

property effective_bc: tuple[str, ...] | tuple[AxisBCSpec, ...]#

Per-axis BC specs.

Returns AxisBCSpec tuple if axis_bcs is set (structured BCs with per-side values), otherwise falls back to string tuple from bc or infers from periodic.

property ndim: int#: Number of spatial dimensions.

property num_points: int#: Total number of grid cells (product of shape).

bounds: tuple[tuple[float, float], ...]#

shape: tuple[int, ...]#

periodic: tuple[bool, ...]#

class tidal.solver.RHSEvaluator(spec, grid, coeff_eval, bc=None)[source]#

Bases: object

Evaluate RHS of field equations with resolved coefficients.

Applies spatial operators to field data and multiplies by resolved coefficients (constant, parameter-overridden, position-dependent, or time-dependent).

Parameters:

spec (EquationSystem) – Parsed JSON equation specification.
grid (GridInfo) – Spatial grid.
coeff_eval (CoefficientEvaluator) – Coefficient resolver with caching.
bc (str or tuple of str, optional) – Boundary conditions for spatial operators.

begin_timestep(t)[source]#

Notify the coefficient evaluator of a new timestep.

Parameters:: t (float)
Return type:: None

evaluate(eq_idx, fields, t=0.0)[source]#

Compute RHS for a single equation.

Parameters:

eq_idx (int) – Index of the equation in spec.equations.
fields (FieldSet) – Current field state.
t (float) – Current simulation time.

Returns:

Grid-shaped result array. Warning: the returned array is an internal buffer and may be overwritten by the next call to evaluate(). Callers must copy if they need to persist it.

Return type:

np.ndarray

evaluate_by_field(field_name, fields, t=0.0)[source]#

Compute RHS for the equation governing field_name.

Raises:

KeyError – If field_name has no associated equation.

Parameters:

field_name (str)
fields (FieldSet)
t (float)

Return type:

np.ndarray

exception tidal.solver.SimulationDivergedError[source]#

Bases: RuntimeError

Raised when simulation fields become non-finite or exceed a norm threshold.

This indicates a physical or numerical instability — e.g. a coupled mass matrix with a negative eigenvalue causing exponentially growing modes.

class tidal.solver.SolverResult[source]#

Bases: TypedDict

Typed result dict returned by all TIDAL solver functions.

All four solvers (CVODE, IDA, scipy, leapfrog) return dicts with these exact keys. Using a TypedDict rather than dict[str, Any] gives pyright full key/value type checking at every call site.

t: NDArray[np.float64]#

y: NDArray[np.float64]#

success: bool#

message: str#

class tidal.solver.StabilityResult[source]#

Bases: object

Result from check_pointwise_mass_stability().

Separates fatal stability errors (negative eigenvalues) from informational notes (e.g. asymmetric matrix detected).

errors: list[str]#

notes: list[str]#

class tidal.solver.StateLayout(slots, num_points, field_slot_map, velocity_slot_map, dynamical_fields)[source]#

Bases: object

Describes the mapping between flat state vector and named fields.

Built from an EquationSystem spec. The flat vector is partitioned into contiguous blocks of num_points elements, one per slot.

Parameters:

slots (tuple[SlotInfo, ...])
num_points (int)
field_slot_map (dict[str, int])
velocity_slot_map (dict[str, int])
dynamical_fields (tuple[str, ...])

slots#

Ordered slot descriptors.

Type:: tuple[SlotInfo, …]

num_points#

Grid points per slot (grid.num_points).

Type:: int

field_slot_map#

Maps field name → slot index.

Type:: dict[str, int]

velocity_slot_map#

Maps field name → velocity slot index (second-order fields only).

Type:: dict[str, int]

dynamical_fields#

Names of dynamical fields (time_order >= 2), in order.

Type:: tuple[str, …]

property algebraic_indices: list[int]#

Flat indices of algebraic (constraint) variables for IDA.

IDA needs to know which entries in the state vector are algebraic (not differential) so it can handle them appropriately.

property constraint_slot_groups: tuple[tuple[int, slice, str], ...]#: Pre-computed (slot_idx, flat_slice, field_name) for constraint slots.

property drift_slot_pairs: tuple[tuple[slice, slice], ...]#

Pre-computed (field_slice, vel_slice) for zero-copy drift.

Allows y[field_slice] += dt * y[vel_slice] without copying velocity data into a separate buffer.

property dynamical_field_slot_groups: tuple[tuple[int, slice, int], ...]#: Pre-computed (slot_idx, flat_slice, vel_slot_idx) for 2nd-order fields.

property first_order_slot_groups: tuple[tuple[int, slice, str], ...]#: Pre-computed (slot_idx, flat_slice, field_name) for 1st-order fields.

classmethod from_spec(spec, num_points)[source]#

Build layout from an equation system specification.

Follows the same ordering as py-pde FieldCollection: for each equation, emit a field slot; if second-order, also emit a velocity slot immediately after.

Parameters:

spec (EquationSystem)
num_points (int)

Return type:

StateLayout

property momentum_slot_map: dict[str, int]#: Alias for velocity_slot_map (transition aid).

property num_slots: int#: Number of slots in the state vector.

property slot_name_to_idx: dict[str, int]#: Map from slot name to slot index. Cached on frozen dataclass.

slot_slice(slot_idx)[source]#

Return the flat-array slice for a given slot index.

Parameters:: slot_idx (int)
Return type:: slice

property total_size: int#: Total flat vector length (num_slots * num_points).

property velocity_slot_groups: tuple[tuple[int, slice, str], ...]#: Pre-computed (slot_idx, flat_slice, field_name) for velocity slots.

slots: tuple[SlotInfo, ...]#

num_points: int#

field_slot_map: dict[str, int]#

velocity_slot_map: dict[str, int]#

dynamical_fields: tuple[str, ...]#

tidal.solver.apply_operator(name, data, grid, bc=None)[source]#

Look up name in the registry and apply it to data.

Raises:

ValueError – If name is not a known operator.

Parameters:

name (str)
data (np.ndarray)
grid (GridInfo)
bc (BCSpec | None)

Return type:

np.ndarray

tidal.solver.check_cfl_stability(spec, grid, dt)[source]#

Check CFL stability condition for explicit time-steppers.

Returns a list of warning strings (empty if all clear). The CFL condition for the wave equation is dt <= dx / c where c is the maximum wave speed (estimated from the laplacian coefficient).

Parameters:

spec (EquationSystem)
grid (GridInfo)
dt (float)

Return type:

tidal.solver.check_mass_sign(coeff_eval, spec)[source]#

Check for sign-changing position-dependent mass terms.

Returns a list of warning strings for tachyonic diagnostics.

Parameters:

coeff_eval (CoefficientEvaluator)
spec (EquationSystem)

Return type:

tidal.solver.check_pointwise_mass_stability(coeff_eval, spec, grid)[source]#

Check eigenvalues of the pointwise mass/coupling matrix M[i,j](x,y).

Builds M[i,j](x,y) from identity-operator terms in dynamical equations (time_derivative_order > 0), then verifies that all eigenvalues are positive at every grid point. A negative eigenvalue indicates an exponentially growing (tachyonic) mode.

Constraint equations (time_derivative_order == 0) are excluded because their algebraic form 0 = +m²φ + ... produces mass signs opposite to the dynamical form d²t φ = -m²φ + .... Including both in one matrix creates a block-indefinite matrix with spurious negative eigenvalues (false positives for vector field systems).

This check runs once pre-simulation using the pre-computed spatial cache in CoefficientEvaluator — zero runtime cost during the actual simulation.

Returns a StabilityResult with errors (instability) and notes (informational diagnostics like asymmetry detection).

Notes

The stability condition is that the potential matrix M (defined as the negative of the identity-operator coefficient matrix) is positive-semidefinite at every grid point. For coupled scalars with Gaussian coupling G(x,y) = g0*exp(-r^2/2R^2), this reduces to mPhi2 * mChi2 > G(x,y)^2 everywhere, i.e. mPhi2 * mChi2 > g0^2 at the coupling peak.

Parameters:

coeff_eval (CoefficientEvaluator)
spec (EquationSystem)
grid (GridInfo)

Return type:

StabilityResult

tidal.solver.operator_min_dim(name)[source]#

Return the minimum spatial dimension required by name.

Raises:: ValueError – If name is not a recognized operator.
Parameters:: name (str)
Return type:: int

tidal.solver.validate_field_references(spec)[source]#

Check that all term field references point to valid fields.

Raises:: ValueError – If a field reference is invalid.
Parameters:: spec (EquationSystem)
Return type:: None

tidal.solver.validate_operator_dimensions(spec)[source]#

Check that all operators are compatible with the spec’s spatial dimension.

Raises:: ValueError – If any operator requires more dimensions than spec.spatial_dimension.
Parameters:: spec (EquationSystem)
Return type:: None

Submodules#

tidal.solver.analytical_jacobian module#

Analytical Jacobian for time-independent linear systems.

For time-independent linear systems (most TIDAL examples), the IDA Jacobian J = dF/dy + cj * dF/dyp consists of two constant sparse matrices. Instead of finite-difference approximation (which requires O(n_colors) residual evaluations per Newton step), we precompute dF/dy and dF/dyp once and supply them analytically.

Three delivery modes depending on system size:

Dense tier (N <= DENSE_THRESHOLD): jacfn fills a 2D numpy array with dF_dy + cj * dF_dyp. ~5.3x speedup vs FD.
Sparse tier (DENSE_THRESHOLD < N <= SPARSE_THRESHOLD): jacfn fills a 1D CSC data array with dF_dy.data + cj * dF_dyp.data, paired with a sparsity pattern for SuperLU_MT direct factorisation. Requires sksundae >= 1.1.2. IDA: ~2.5x, CVODE: ~1.2-1.4x.
GMRES tier (N > SPARSE_THRESHOLD): jactimes provides an analytical Jacobian-vector product Jv = dF_dy @ v + cj * dF_dyp @ v using sparse matrix-vector products, eliminating finite-difference residual evaluations per GMRES iteration.

Performance optimizations:

COO accumulation: build_jacobian_matrices() uses _COOAccumulator to append block triples and do a single CSC conversion, avoiding O(N²)-per-block LIL slice assignment.
Circulant operators: For all-periodic BCs, build_operator_matrix() probes a single grid point and tiles the stencil via modular arithmetic — O(nnz) instead of O(N²) probing.

Position-dependent (but time-independent) coefficients are supported: the spatial grid is fixed, so the Jacobian is still constant.

tidal.solver.analytical_jacobian.build_operator_matrix(operator, grid, bc)[source]#

Build the N x N sparse matrix for a single spatial operator.

For all-periodic BCs, uses a fast O(nnz) circulant construction (single-probe + tiling). Otherwise falls back to O(N²) column-by-column probing.

Parameters:

operator (str)
grid (GridInfo)
bc (BCSpec | None)

Return type:

SparseMatrix

tidal.solver.analytical_jacobian.build_jacobian_matrices(spec, layout, grid, bc, parameters)[source]#

Build the analytical dF/dy and dF/dyp sparse matrices.

Mirrors the residual handler structure in ida.py exactly: constraint (3 sub-cases), velocity, dynamical_field, first_order.

Parameters:

spec (EquationSystem) – Parsed equation specification.
layout (StateLayout) – State vector layout.
grid (GridInfo) – Spatial grid.
bc (BCSpec or None) – Boundary conditions.
parameters (dict[str, float]) – Resolved parameter values.

Returns:

(dF_dy, dF_dyp) – Sparse Jacobian component matrices.

Return type:

tuple[SparseMatrix, SparseMatrix]

tidal.solver.analytical_jacobian.try_analytical_jacobian(options, spec, layout, grid, bc, parameters, *, solver='ida')[source]#

Try to configure analytical Jacobian for time-independent systems.

Mutates options in-place if successful. Returns True on success, False if the system has time-dependent coefficients (caller should fall back to the finite-difference tier system).

Three delivery modes by system size:

Dense (N <= DENSE_THRESHOLD): 2D jacfn fills dense array.
Sparse (DENSE_THRESHOLD < N <= SPARSE_THRESHOLD): 1D CSC jacfn + SuperLU_MT direct factorisation. Falls through to GMRES if nnz exceeds SUPERLU_NNZ_LIMIT.
GMRES (N > SPARSE_THRESHOLD): jactimes provides analytical Jacobian-vector product for iterative GMRES.

Parameters:

solver (str) – "ida" or "cvode". Controls the callback signature.
options (dict[str, Any])
spec (EquationSystem)
layout (StateLayout)
grid (GridInfo)
bc (BCSpec | None)
parameters (dict[str, float])

Return type:

bool

tidal.solver.coefficients module#

Coefficient evaluation with multi-level caching for TIDAL solvers.

CoefficientEvaluator resolves operator term coefficients, handling: - Constant numeric coefficients (no symbolic expression) - Parameter-overridable symbolic coefficients (e.g. "m2", "-kappa") - Position-dependent coefficients (evaluated on the spatial grid) - Time-dependent coefficients (re-evaluated each timestep)

Caching levels:: L0: Constants pre-resolved at init (no coordinate/time dependence) L1: Mathematica → Python string cache (computed once) L2: Spatial-only arrays (position-dependent, NOT time-dependent) L3: Per-timestep deduplication (time+position dependent)

Reuses the evaluation pipeline from tidal.symbolic._eval_utils — no reimplementation of Mathematica→Python conversion or eval() logic.

class tidal.solver.coefficients.CoefficientEvaluator(spec, grid, parameters=None)[source]#

Bases: object

Resolve operator term coefficients with multi-level caching.

Parameters:

spec (EquationSystem) – Parsed JSON equation specification.
grid (GridInfo) – Spatial grid (for position-dependent coefficient evaluation).
parameters (dict[str, float]) – Runtime parameter overrides (e.g. {"m2": 1.0}).

Raises:

ValueError – If a symbolic coefficient references an unknown parameter and no default numeric value is available.

resolve(term, t=0.0, *, eq_idx=-1, term_idx=-1)[source]#

Resolve effective coefficient for an operator term.

Parameters:

term (OperatorTerm) – The operator term whose coefficient to resolve.
t (float) – Current simulation time.
eq_idx (int) – Equation index (for cache lookup). -1 disables caching.
term_idx (int) – Term index within equation (for cache lookup). -1 disables caching.

Returns:

Scalar for constant coefficients, grid-shaped array for position-dependent ones.

Return type:

float | NDArray[np.float64]

begin_timestep(t)[source]#

Clear per-timestep cache (L3).

Call at the start of each timestep to ensure time-dependent coefficients are re-evaluated. Skipped when all coefficients are time-independent (common case — avoids empty dict.clear()).

Parameters:: t (float)
Return type:: None

all_constant()[source]#

Check if every RHS term has a constant (scalar) coefficient.

Returns True when all coefficients are pre-resolved at L0 (no position or time dependence). This enables the analytical Jacobian optimization for constant-coefficient linear systems.

Return type:: bool

all_time_independent()[source]#

Check if every RHS term has a time-independent coefficient.

Return type:: bool

check_periodic_coefficient_continuity(periodic, *, rtol=None)[source]#

Warn if position-dependent coefficients are discontinuous at periodic boundaries.

The check uses a leak metric that estimates the IBP energy leak:

leak = (jump / scale) * (boundary_magnitude / scale)

Parameters:

periodic (tuple[bool, ...]) – Per-axis periodicity flags.
rtol (float, optional) – Solver relative tolerance. When provided, the warn and error thresholds scale as sqrt(rtol) — stricter for machine-precision solvers (modal, leapfrog) and more lenient for coarse exploratory runs. When None, the default thresholds (warn=0.01, error=0.25) are used.

Raises:

ValueError – If the leak metric exceeds the error threshold.

Return type:

None

tidal.solver.constraint_solve module#

Pre-solve algebraic constraints to produce consistent initial conditions.

Before IDA starts, constraint equations (time_order=0) must be satisfied. For example, Gauss’s law laplacian(A_0) = -div(pi) requires A_0 to be non-trivially solved when the source (pi fields) is nonzero.

Three-tier solver hierarchy (automatically selected):

Tier 1 — FFT (O(N log N)): All axes periodic, all self-coefficients constant, all operators have known Fourier multipliers. Uses modified wavenumbers k_mod = (2/dx)*sin(k*dx/2) for exact FD-consistency.

Tier 2 — Operator probing → sparse matrix (O(N²) build, O(N) solve): Universal fallback. Applies apply_operator() to unit vectors e_j to build the operator matrix column by column. Automatically handles position-dependent coefficients, arbitrary BCs, and unknown/future operators.

Tier 3 — Iterative (future): CG/GMRES for very large grids where direct factorization is impractical.

Design choices:

Tier 2 probing reuses the exact same apply_operator() code path as the simulation, guaranteeing mathematical consistency.
Position-dependent coefficients (NDArrays from CoefficientEvaluator) multiply element-wise in probing, yielding correct spatially-varying matrix entries. This is critical for curved-spacetime or background-field constraints where translational symmetry is broken.
FFT eligibility is checked conservatively: any position-dependent self-coeff, non-periodic axis, or unknown multiplier → fall back to Tier 2.
Zero-mode handling for singular operators (pure Poisson/periodic): enforce zero mean on the solution by setting u_hat[0,...,0] = 0, with a compatibility check that the source has zero mean.

References

Modified wavenumbers for finite differences: Lele, J. Comp. Phys., 1992.
Spectral constraint projection: Dedalus (Burns et al., PRR 2020).

tidal.solver.constraint_solve.pre_solve_constraints(spec, grid, y0, *, bc=None, parameters=None, t=0.0)[source]#

Solve algebraic constraints to produce consistent initial conditions.

Only processes equations where time_derivative_order == 0 AND constraint_solver.enabled == True. Returns a copy of y0 with constraint field slots overwritten by the solution.

Parameters:

spec (EquationSystem) – Parsed JSON equation specification.
grid (GridInfo) – Spatial grid.
y0 (NDArray[np.float64]) – Initial state vector (flat, from StateLayout).
bc (str or tuple of str, optional) – Boundary conditions for spatial operators.
parameters (dict[str, float], optional) – Runtime parameter overrides for symbolic coefficients.
t (float) – Time at which to evaluate coefficients (usually t_span[0]).

Returns:

Updated y0 with constraint fields solved.

Return type:

NDArray[np.float64]

Warns:

UserWarning – When the FFT solver encounters singular modes (null space of the operator) and regularizes by setting u_hat = 0 at those modes. This is a numerical gauge choice (zero-mean). To disable automatic constraint solving for a field, set constraint_solver.enabled = false in the JSON spec.

tidal.solver.constraint_solve.ensure_consistent_ic(spec, grid, y0, *, bc=None, parameters=None, t=0.0, strict=True)[source]#

Unified constraint IC solver.

Given user-supplied initial conditions, iteratively solves ALL constraint equations (time_derivative_order == 0) to produce consistent initial data. Handles three cases uniformly:

Standard constraints (field has self-terms): solve for the constraint field. Applies to all constraint equations, not just constraint_solver.enabled=True.
Subsidiary constraints (no self-terms, one free source field): solve for the free dynamical field. E.g., transverse gauge constraint gradient_x(h_4) + gradient_x(h_7) = 0 determines h_7 = -h_4 when h_4 is user-initialized.
Verification (all referenced fields determined): check that the constraint is satisfied; error or warn if not.

The algorithm iterates until no more fields can be determined. Each iteration may unlock new solvable constraints (cascade).

Parameters:

spec (EquationSystem) – Parsed JSON equation specification.
grid (GridInfo) – Spatial grid.
y0 (NDArray[np.float64]) – Initial state vector (flat, from StateLayout).
bc (str or tuple of str, optional) – Boundary conditions for spatial operators.
parameters (dict[str, float], optional) – Runtime parameter overrides for symbolic coefficients.
t (float) – Time at which to evaluate coefficients.
strict (bool) – If True (default), raise ValueError when a constraint is violated and cannot be solved. If False, issue a warning.

Returns:

Updated y0 with all solvable constraint fields determined.

Return type:

NDArray[np.float64]

Raises:

ValueError – If strict=True and any constraint is violated after solving.

tidal.solver.cvode module#

SUNDIALS/CVODE integration for TIDAL — adaptive ODE solver for wave systems.

Pure ODE solver using the CVODE module from SUNDIALS via scikit-sundae. Supports both BDF (stiff, order 1-5) and Adams (non-stiff, order 1-12) methods with tolerance-controlled adaptive time stepping.

For wave (second-order) equations, the system is reduced to first-order ODE form using E-L velocity formulation:

dq/dt = v (trivial kinematic) dv/dt = E-L RHS (from spatial operators)

Constraint fields (time_order=0) are frozen at initial values — correct for gauge-fixed systems (e.g. Coulomb gauge A_0 = 0).

Reference: Hindmarsh et al., “SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers”, ACM TOMS 31(3), 2005. scikit-sundae: NREL, BSD-3 license.

tidal.solver.cvode.solve_cvode(spec, grid, y0, t_span, *, bc=None, parameters=None, method='BDF', rtol=1e-08, atol=1e-10, max_step=0.0, max_num_steps=50000, num_snapshots=101, snapshot_callback=None, progress=None)[source]#

Solve a TIDAL equation system using SUNDIALS/CVODE.

Parameters:

spec (EquationSystem) – Parsed equation specification.
grid (GridInfo) – Spatial grid.
y0 (np.ndarray) – Initial state vector (flat).
t_span (tuple[float, float]) – (t_start, t_end).
bc (str or tuple, optional) – Boundary conditions.
parameters (dict[str, float], optional) – Runtime parameter overrides for symbolic coefficients.
method (str) – Integration method: 'BDF' (stiff, default) or 'Adams' (non-stiff).
rtol (float) – Relative and absolute tolerances.
atol (float) – Relative and absolute tolerances.
max_step (float) – Maximum step size. 0 (default) = unbounded (SUNDIALS default).
max_num_steps (int) – Maximum solver steps per output interval.
num_snapshots (int) – Number of output time points.
snapshot_callback (callable, optional) – Called as callback(t, y) at each output time.
progress (SimulationProgress | None)

Returns:

Result dictionary with keys: t, y, success, message.

Return type:

Warns:

UserWarning – If constraint fields (time_order=0) are present — they remain frozen at initial values.

tidal.solver.fields module#

Typed field container for TIDAL solvers.

FieldSet owns named field data on a grid, backed by a single contiguous flat numpy array. Named access returns zero-copy views (numpy slices).

This consolidates the state packing/unpacking logic previously duplicated across ida.py, leapfrog.py, _simulate.py, _io.py, and _writer.py.

class tidal.solver.fields.FieldSet(layout, grid_shape, data=None)[source]#

Bases: object

Typed container owning named field data on a grid.

Backed by a single contiguous flat np.ndarray. Named access returns zero-copy views (numpy slices) into this array.

Parameters:

layout (StateLayout) – Describes the slot structure (field/velocity names and ordering).
grid_shape (tuple[int, ...]) – Shape of the spatial grid (e.g. (64,) or (32, 32)).
data (np.ndarray, optional) – Flat data array of length layout.total_size. If None, initializes to zeros.

set_aux(name, data)[source]#

Used to inject constraint velocities from IDA’s yp vector.

Parameters:

name (str)
data (ndarray)

Return type:

None

property flat: ndarray#: Underlying flat vector (no copy). Used by IDA/leapfrog.

property layout: StateLayout#: The state layout this FieldSet wraps.

property grid_shape: tuple[int, ...]#: Spatial grid shape.

property field_names: tuple[str, ...]#

Ordered field slot names (e.g. ("phi_0", "A_0")).

Excludes velocity slots.

property velocity_names: tuple[str, ...]#

Ordered velocity slot names (e.g. ("v_phi_0",)).

Only present for second-order equations.

property momentum_names: tuple[str, ...]#: Alias for velocity_names (transition aid).

property slot_names: tuple[str, ...]#: All slot names in order.

fields_dict()[source]#

Return dict of field name → grid-shaped view (zero-copy).

Return type:: dict[str, ndarray]

velocities_dict()[source]#

Return dict of velocity name → grid-shaped view (zero-copy).

Return type:: dict[str, ndarray]

momenta_dict()[source]#

Alias for velocities_dict (transition aid).

Return type:: dict[str, ndarray]

as_dict()[source]#

Return dict of all slot names + aux → grid-shaped views.

Return type:: dict[str, ndarray]

copy()[source]#

Deep copy (new flat array, copied aux).

Return type:: FieldSet

classmethod zeros(layout, grid_shape)[source]#

Create a FieldSet initialized to zero.

Parameters:

layout (StateLayout)
grid_shape (tuple[int, ...])

Return type:

classmethod from_flat(layout, grid_shape, flat)[source]#

Wrap an existing flat array (no copy).

The caller must ensure the array is not unexpectedly mutated elsewhere.

Parameters:

layout (StateLayout)
grid_shape (tuple[int, ...])
flat (np.ndarray)

Return type:

rebind(flat)[source]#

Replace the underlying flat array reference (zero-copy, zero-alloc).

Used by solver loops to reuse a single FieldSet instead of allocating a new one each timestep. Auxiliary fields are cleared.

Parameters:: flat (ndarray)
Return type:: None

classmethod from_dict(layout, grid_shape, slot_data)[source]#

Pack named arrays into a FieldSet.

Missing slots default to zero. Extra keys not in the layout are silently ignored.

Parameters:

layout (StateLayout) – State layout descriptor.
grid_shape (tuple[int, ...]) – Spatial grid shape.
slot_data (dict[str, np.ndarray]) – Mapping from slot name → grid-shaped array.

Return type:

max_norm()[source]#

Maximum absolute value across all slots.

Return type:: float

check_finite()[source]#

Return True iff all values are finite (no NaN or Inf).

Return type:: bool

tidal.solver.grid module#

GridInfo — lightweight Cartesian grid descriptor for TIDAL.

Replaces py-pde’s CartesianGrid with a minimal frozen dataclass that provides only what TIDAL needs: grid spacing for FD stencils, cell coordinates for position-dependent coefficient evaluation, and periodicity flags for BC dispatch.

class tidal.solver.grid.GridInfo(bounds, shape, periodic, bc=None, axis_bcs=None)[source]#

Bases: object

Immutable Cartesian grid descriptor.

Parameters:

bounds (tuple of (lo, hi) pairs) – Domain bounds per spatial axis, e.g. ((0, 10), (0, 10)).
shape (tuple of int) – Number of grid cells per axis, e.g. (64, 64).
periodic (tuple of bool) – Whether each axis is periodic, e.g. (True, True).
bc (tuple of str, optional) – Legacy per-axis BC type strings (e.g. ("neumann", "periodic")).
axis_bcs (tuple of AxisBCSpec, optional) – Structured per-axis BC specs with per-side values and Robin support. Takes precedence over bc when set.

Examples

>>> g = GridInfo(bounds=((0, 10),), shape=(64,), periodic=(True,))
>>> g.ndim
1
>>> g.dx
(0.15625,)
>>> g.num_points
64

bounds: tuple[tuple[float, float], ...]#

shape: tuple[int, ...]#

periodic: tuple[bool, ...]#

bc: tuple[str, ...] | None = None#

axis_bcs: tuple[AxisBCSpec, ...] | None = None#

property ndim: int#: Number of spatial dimensions.

property num_points: int#: Total number of grid cells (product of shape).

property effective_bc: tuple[str, ...] | tuple[AxisBCSpec, ...]#

Per-axis BC specs.

Returns AxisBCSpec tuple if axis_bcs is set (structured BCs with per-side values), otherwise falls back to string tuple from bc or infers from periodic.

property bc_types: tuple[str, ...]#

Per-axis BC type as simple strings.

Always returns tuple[str, ...] — one of "periodic", "neumann", "dirichlet", or "robin" per axis.

property dx: tuple[float, ...]#

dx = (hi - lo) / N).

Type:: Grid spacing per axis (cell-centred

axes_coords(axis)[source]#

1-D array of cell centres along axis.

Returns shape (shape[axis],) — a single 1-D vector, not a broadcasted grid. Useful for building per-axis coordinate arrays without the overhead of full meshgrid.

Raises:: IndexError – If axis is out of range.
Parameters:: axis (int)
Return type:: ndarray

property cell_coords: ndarray#

Cell-centre coordinates, shape (*shape, ndim).

For a 1D grid with bounds (0, 10) and shape (4,), the cell centres are at [1.25, 3.75, 6.25, 8.75] (centred in each cell).

Unlike py-pde’s CartesianGrid.cell_coords (which returns a list of arrays with inconsistent shapes), this always returns a single ndarray with the coordinate dimension as the last axis.

Cached on first access (GridInfo is immutable).

coord_arrays()[source]#

Broadcasted coordinate arrays, one per axis.

Returns ndim arrays, each of shape self.shape, containing the coordinate value along that axis at every grid point. Equivalent to np.meshgrid(*1d_axes, indexing="ij").

Cached on first access (GridInfo is immutable).

Return type:: tuple[ndarray, …]

tidal.solver.ida module#

SUNDIALS/IDA integration for TIDAL — DAE solver for mixed systems.

Builds IDA-compatible residual functions from TIDAL equation specs. Handles arbitrary mixes of: - Second-order (wave) equations via E-L velocity form - First-order (diffusion/transport) equations - Algebraic (constraint) equations

Euler-Lagrange velocity form: - Velocity slot: dv/dt = E-L RHS (second-order equation) - Field slot: dq/dt = v (trivial kinematic)

Reference: Hindmarsh et al., “SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers”, ACM TOMS 31(3), 2005. scikit-sundae: NREL, BSD-3 license.

tidal.solver.ida.build_residual_fn(spec, layout, grid, bc=None, *, parameters=None, rtol=None)[source]#

Build an IDA-compatible residual function from a TIDAL equation spec.

The returned function has signature resfn(t, y, yp, res) where res is written in-place.

For each slot in the state vector:

Constraint (time_order=0): res = RHS(y, t) (algebraic, = 0)
First-order field (time_order=1): res = yp - RHS(y, t)
Second-order field slot: res = yp - v (trivial kinematic)
Second-order velocity slot: res = yp - RHS(y, t) (E-L equation)

Parameters:

spec (EquationSystem) – Parsed JSON equation specification.
layout (StateLayout) – State vector layout descriptor.
grid (GridInfo) – Spatial grid.
bc (str or tuple of str, optional) – Boundary conditions for spatial operators.
parameters (dict[str, float], optional) – Runtime parameter overrides for symbolic coefficients. When provided, enables position-dependent and time-dependent coefficient evaluation via CoefficientEvaluator.
rtol (float | None)

Returns:

IDA residual function with signature (t, y, yp, res) -> None.

Return type:

Callable

tidal.solver.ida.solve_ida(spec, grid, y0, t_span, *, bc=None, parameters=None, num_snapshots=101, rtol=1e-08, atol=1e-10, max_steps=50000, snapshot_callback=None, calc_initcond=None, allow_inconsistent_ic=False, progress=None)[source]#

Solve a TIDAL equation system using SUNDIALS/IDA.

Parameters:

spec (EquationSystem) – Parsed equation specification.
grid (GridInfo) – Spatial grid.
y0 (np.ndarray) – Initial state vector (flat).
t_span (tuple[float, float]) – (t_start, t_end).
bc (str or tuple, optional) – Boundary conditions.
parameters (dict[str, float], optional) – Runtime parameter overrides for symbolic coefficients.
num_snapshots (int) – Number of output time points.
rtol (float) – Relative and absolute tolerances.
atol (float) – Relative and absolute tolerances.
max_steps (int) – Maximum solver steps.
snapshot_callback (callable, optional) – Called as callback(t, y) at each output time.
calc_initcond (str, optional) – IDA initial condition calculation mode. "yp0" (default for mixed DAE) corrects derivatives given y0. "y0" corrects algebraic variables given yp0 — use this for constraint solving where the algebraic field values are unknown.
allow_inconsistent_ic (bool) – If False (default), raise ValueError when constraint equations are violated and cannot be solved. If True, issue a warning instead.
progress (SimulationProgress | None)

Returns:

Result dictionary with keys: t, y, success, message.

Return type:

Warns:

UserWarning – When constraint pre-solve encounters singular modes (FFT path), or when a constraint field is detected as needing gauge regularization (pure Laplacian + periodic BCs → field[0] = 0 pinning). These are numerical choices to resolve the null space of the operator. To disable for a specific field, set constraint_solver.enabled = false in the JSON spec.

tidal.solver.leapfrog module#

Symplectic integrators for TIDAL systems (leapfrog / Yoshida).

Symplectic integrators for second-order E-L equations using velocity form:: dq/dt = v (trivial kinematic) dv/dt = E-L RHS (from equations[] array)

Two integration orders are available via --leapfrog-order:

Order 2 (default): Stormer-Verlet KDK scheme. Shadow Hamiltonian error O(dt^2), one force evaluation per step (with caching).
Order 4: Yoshida triple-composition S4(dt) = S2(w1*dt) o S2(w2*dt) o S2(w3*dt). Shadow Hamiltonian error O(dt^4), three force evaluations per step (fused-kick scheme: adjacent half-kicks merge between sub-steps).

Why Yoshida is faster at equal accuracy: each individual Yoshida step is ~3x more expensive (3 vs 1 force evaluations via fused-kick caching). However, the O(dt^4) accuracy means far fewer steps are needed to reach a given error target. For example, to reach ~4e-6 accuracy on coupled scalars (N=256): - Verlet requires dt=0.005 (4000 steps, 0.21s) - Yoshida requires dt=0.040 (500 steps, 0.07s) – 3x faster The speedup comes from both the accuracy-to-cost ratio and the fused-kick optimization (Forest & Ruth, 1990; Hairer et al., 2006 §II.4).

CFL note: the effective CFL limit is reduced by max(|wi|) ~ 1.70 because the middle sub-step w2 < 0 has |w2| > 1. The CLI auto-adjusts dt by this factor when --leapfrog-order 4 is specified.

Time-dependent forces: sub-step time is tracked correctly through all three Yoshida sub-steps, so time-dependent coefficients (background fields, time-varying sources) are handled correctly for general systems.

Both orders preserve a shadow Hamiltonian to machine precision, giving zero secular energy drift for conservative systems.

Not appropriate for: dissipative systems, absorbing BCs, constraint damping, or energy outflow – use IDA instead.

References

Hairer, Lubich, Wanner, “Geometric Numerical Integration”, Springer, 2006. Chapter VI: Symplectic Integration of Hamiltonian Systems.

Yoshida, H. (1990). “Construction of higher order symplectic integrators”. Physics Letters A, 150(5-7), pp. 262-268.

Forest, E. & Ruth, R.D. (1990). “Fourth-order symplectic integration”. Physica D, 43(1), pp. 105-117. (Fused-kick composition method.)

tidal.solver.leapfrog.compute_force(spec, layout, grid, bc, y, t=0.0, rhs_eval=None, out=None, fieldset=None)[source]#

Compute dv/dt = E-L RHS for all velocity slots.

Parameters:

out (np.ndarray, optional) – Pre-allocated output buffer of length layout.total_size. If provided, filled in-place (avoids allocation). If None, a fresh array is allocated.
fieldset (FieldSet, optional) – Reusable FieldSet. If provided, rebound to y in-place (avoids per-call object allocation).
spec (EquationSystem)
layout (StateLayout)
grid (GridInfo)
bc (BCSpec | None)
y (np.ndarray)
t (float)
rhs_eval (RHSEvaluator | None)

Return type:

np.ndarray

tidal.solver.leapfrog.compute_velocity(layout, y, out=None)[source]#

Read dq/dt = v directly from velocity slots in the state vector.

Parameters:

out (np.ndarray, optional) – Pre-allocated output buffer. If provided, filled in-place.
layout (StateLayout)
y (ndarray)

Return type:

ndarray

tidal.solver.leapfrog.solve_leapfrog(spec, grid, y0, t_span, dt, *, bc=None, parameters=None, snapshot_interval=None, snapshot_callback=None, progress=None, order=2)[source]#

Solve a TIDAL Hamiltonian system using symplectic integration.

Works for second-order (wave) equations with optional constraint fields. Constraint fields (time_order=0) are frozen at their initial values — correct for gauge-fixed systems (e.g. Coulomb gauge A_0 = 0).

Parameters:

spec (EquationSystem) – Parsed equation specification.
grid (GridInfo) – Spatial grid.
y0 (np.ndarray) – Initial state vector (flat: [q_0, v_0, q_1, v_1, …]).
t_span (tuple[float, float]) – (t_start, t_end).
dt (float) – Fixed time step.
bc (str or tuple, optional) – Boundary conditions.
parameters (dict[str, float], optional) – Runtime parameter overrides for symbolic coefficients.
snapshot_interval (float, optional) – Time between snapshots. If None, only initial and final states saved.
snapshot_callback (callable, optional) – Called as callback(t, y) at each snapshot time.
order (int, optional) – Integration order: 2 (Störmer-Verlet, default) or 4 (Yoshida). Order 4 uses triple composition of order 2 sub-steps with weights from Yoshida (1990), achieving O(dt^4) accuracy at 3x force cost per step. The fused-kick scheme (Forest & Ruth, 1990) merges adjacent half-kicks to reduce 6→3 force evaluations per step.
progress (SimulationProgress | None)

Returns:

Result dictionary with keys: t, y, success, message.

Return type:

Raises:

ValueError – If the system contains first-order (diffusion/transport) equations, or if order is not 2 or 4.

Warns:

UserWarning – If constraint fields (time_order=0) are present — they remain frozen at initial values.

tidal.solver.modal module#

Fourier modal solver — exact spectral time evolution for linear PDEs.

Transforms the spatial grid to Fourier space, builds a per-mode evolution matrix, and eigendecomposes to obtain the exact solution y(t) = exp(A·t)·y₀. Eliminates spatial discretization error entirely and provides machine-precision solutions for time-independent linear systems.

Applicable to any linear PDE system with: - Flat (Minkowski) metric - All-periodic boundary conditions - Time-independent coefficients (position-dependent OK via convolution) - Operators with known exact Fourier multipliers

Two algorithm paths are used depending on coefficient structure: - Constant coefficients: per-mode eigendecomposition with block-aware independent

blocks (machine-precision, ~14x faster).

Position-dependent coefficients: Krylov matrix exponential (expm_multiply) which is backward-stable for non-normal convolution matrices where eigendecomposition gives incorrect results due to pseudospectral overflow.

References

Moler & Van Loan (2003), SIAM Review 45(1):3-49 — matrix exponential. Hairer, Lubich & Wanner (2006), Geometric Numerical Integration, §4. Burns et al. (2020), Phys. Rev. Research 2:023068 — pseudo-spectral. Raffelt & Stodolsky (1988), PRD 37:1237 — mixing-matrix formalism.

tidal.solver.modal.can_use_modal(spec, grid, bc)[source]#

Check whether the modal solver is applicable to this system.

Requirements (checked in order): 1. Flat metric (volume_element is None) 2. Constraints (time_order=0) must be Fourier-eliminable via Schur complement 3. All boundary conditions periodic 4. All RHS operators have exact Fourier multipliers 5. No time-dependent coefficients

Parameters:

spec (EquationSystem)
grid (GridInfo)
bc (BCSpec | None)

Return type:

bool

tidal.solver.modal.solve_modal(spec, grid, y0, t_span, *, bc=None, parameters=None, rtol=1e-08, atol=1e-10, num_snapshots=101, snapshot_callback=None, progress=None)[source]#

Solve a TIDAL equation system using Fourier modal decomposition.

Transforms the spatial grid to Fourier space, builds the per-mode or full evolution matrix, eigendecomposes, and evaluates the exact solution at each output time.

Parameters:

spec (EquationSystem) – Parsed equation specification (from JSON).
grid (GridInfo) – Spatial grid (must be all-periodic).
y0 (np.ndarray) – Initial state vector (flat).
t_span (tuple[float, float]) – (t_start, t_end).
bc (str or tuple, optional) – Boundary conditions (must be all-periodic).
parameters (dict[str, float], optional) – Runtime parameter overrides for symbolic coefficients.
rtol (float) – Tolerances (unused for eigendecomposition; reserved for solve_ivp fallback with time-dependent coefficients).
atol (float) – Tolerances (unused for eigendecomposition; reserved for solve_ivp fallback with time-dependent coefficients).
num_snapshots (int) – Number of output time points.
snapshot_callback (callable, optional) – Called as callback(t, y) at each output time.
progress (SimulationProgress, optional) – Progress tracker for tqdm display.

Returns:

Dict with keys: t, y, success, message.

Return type:

SolverResult

tidal.solver.operators module#

Spatial operators on plain numpy arrays (finite-difference and spectral).

All operators accept an np.ndarray (field data) and a GridInfo descriptor, returning an np.ndarray of the same shape. No py-pde ScalarField / VectorField wrappers — pure numpy with explicit ghost-cell boundary handling for FD, or FFT-based spectral differentiation for periodic domains.

Operator modes#

Two modes are available, selectable at runtime:

Finite-difference (default): Ghost-cell stencils with 2nd, 4th, or 6th-order accuracy (set_fd_order). Works with all BC types.
Spectral (set_spectral(True)): FFT-based pseudo-spectral operators via np.fft.rfft/np.fft.irfft. Requires all-periodic BCs. Achieves exponential convergence for smooth problems — typically N=64-128 instead of N=512-1024 for equivalent accuracy.

The pseudo-spectral approach computes derivatives in Fourier space (multiply by ik for gradient, -k^2 for Laplacian) and returns to physical space. Position-dependent coefficients are applied in physical space after the operator, following the standard Dedalus/py-pde pseudo-spectral pattern.

Ref: Burns et al. (2020), “Dedalus: A flexible framework for numerical simulations with spectral methods”, Physical Review Research 2:023068.

Finite-difference accuracy orders#

Supports 2nd, 4th, and 6th-order central-difference stencils, selectable at runtime via set_fd_order(). Higher orders use wider stencils (more ghost cells) for dramatically improved spatial accuracy:

Order 2 (default): 3-point stencil, 1 ghost cell/side.
Order 4: 5-point stencil, 2 ghost cells/side.
Order 6: 7-point stencil, 3 ghost cells/side.

Stencil coefficients from Fornberg (1988), “Generation of Finite Difference Formulas on Arbitrarily Spaced Grids”, Mathematics of Computation 51(184), pp. 699-706, Table 1 (central differences).

Boundary condition types#

"periodic" — wraparound via np.roll
"neumann" — zero normal derivative (mirror ghost cell)
"dirichlet"— zero value (antisymmetric ghost cell)
"robin" – mixed: d_n f + gamma*f = beta

All non-periodic BCs are unified via the ghost-cell formula ghost = const + factor * interior_cell where (const, factor) are determined by BC type, value, and grid spacing. See SideBCSpec for the derivation. For higher-order stencils (order > 2), ghost cells are filled layer-by-layer from interior outward (recursive ghost fill; LeVeque, “Finite Difference Methods for ODEs and PDEs”, SIAM, 2007, ch. 2).

The bc parameter may be: - A single string applied to all axes (e.g. "periodic") - A tuple of strings, one per axis (e.g. ("neumann", "periodic")) - A tuple of AxisBCSpec objects for per-side / non-zero / Robin BCs

References

Burns et al. (2020), Phys. Rev. Research 2:023068 — pseudo-spectral methods.
Fornberg (1988), Math. Comp. 51(184), pp. 699-706 — stencil coefficients.
LeVeque (2007), “Finite Difference Methods for ODEs and PDEs”, SIAM — convergence theory.
Zwicker (2020), J. Open Source Software 5(48) — ghost-cell virtual-point approach.

tidal.solver.operators.set_fd_order(order)[source]#

Set the finite-difference accuracy order for all spatial operators.

Must be called before any operator evaluation (typically at simulation startup). Clears the ghost-cell padding cache since buffer shapes change.

Parameters:: order (int) – Accuracy order: 2 (3-point stencil), 4 (5-point), or 6 (7-point).
Raises:: ValueError – If order is not 2, 4, or 6.
Return type:: None

tidal.solver.operators.get_fd_order()[source]#

Return the current FD accuracy order.

Return type:: int

tidal.solver.operators.get_n_ghost()[source]#

Return the current number of ghost cells per side.

Return type:: int

tidal.solver.operators.set_spectral(enable)[source]#

Enable or disable FFT-based spectral operators.

When enabled, all spatial operators (gradient, Laplacian, cross-derivative, biharmonic) use np.fft.rfft/np.fft.irfft instead of FD stencils. Requires all boundary conditions to be periodic.

Must be called before any operator evaluation (typically at simulation startup). Clears the wavenumber cache when toggled.

Parameters:: enable (bool) – True to use spectral operators, False for finite-difference.
Return type:: None

tidal.solver.operators.get_spectral()[source]#

Return whether spectral (FFT) operators are active.

Return type:: bool

tidal.solver.operators.get_wavenumbers(n, dx)[source]#

Return cached wavenumber array for rfft of length n.

Wavenumbers: k = 2*pi*rfftfreq(n, d=dx) — the angular frequencies corresponding to the real-FFT output bins.

Parameters:

n (int) – Number of grid points along the axis.
dx (float) – Grid spacing.

Returns:

Wavenumber array of length n // 2 + 1.

Return type:

np.ndarray

class tidal.solver.operators.SideBCSpec(kind, value=0.0, derivative=0.0, gamma=0.0)[source]#

Bases: object

Boundary condition for one side (low or high) of one axis.

All first-order BCs map to a ghost-cell formula:

ghost = const + factor * interior_cell

where (const, factor) depend on the BC kind:

Neumann (df/dn = D): const = dx*D, factor = +1
Dirichlet (f = V): const = 2*V, factor = -1
Robin (d_n f + gamma*f = beta): const = 2*dx*beta / (gamma*dx + 2), factor = (2 - gamma*dx) / (gamma*dx + 2)

Setting D=0, V=0 recovers the homogeneous cases.

Parameters:

kind (str)
value (float)
derivative (float)
gamma (float)

kind: str#

"dirichlet", "neumann", or "robin".

Type:: BC type

value: float = 0.0#: Dirichlet boundary value V, or Robin inhomogeneity β.

derivative: float = 0.0#

Neumann ghost-cell offset D.

Ghost cell is placed at interior + dx * D. For homogeneous Neumann (zero normal flux), set D=0. The sign is the same at both low and high boundaries (py-pde convention; see Zwicker, JOSS 2020).

gamma: float = 0.0#: Robin coefficient gamma in d_n f + gamma*f = beta. Must be >= 0.

ghost_params(dx)[source]#

Return (const, factor) for ghost = const + factor * interior.

Parameters:: dx (float) – Grid spacing in the normal direction.
Return type:: tuple[float, float]

class tidal.solver.operators.AxisBCSpec(periodic=False, low=None, high=None)[source]#

Bases: object

Boundary condition for one axis: periodic or distinct low/high sides.

For periodic axes, low and high must be None. For non-periodic axes, both must be SideBCSpec instances.

Parameters:

periodic (bool)
low (SideBCSpec | None)
high (SideBCSpec | None)

periodic: bool = False#

low: SideBCSpec | None = None#

high: SideBCSpec | None = None#

tidal.solver.operators.is_periodic_bc(bc_entry)[source]#

Check whether a single axis BC entry is periodic.

Handles both legacy string BCs and structured AxisBCSpec objects, so callers don’t need to type-check manually.

Parameters:: bc_entry (str | AxisBCSpec)
Return type:: bool

tidal.solver.operators.BCSpec = str | tuple[str, ...] | tuple[tidal.solver.operators.AxisBCSpec, ...]#: Single BC string, per-axis strings, or per-axis AxisBCSpec objects.

tidal.solver.operators.gradient(data, axis, grid, bc=None)[source]#

Central-difference gradient ∂f/∂x_i.

Stencil width adapts to the current FD order (set via set_fd_order):

Order 2 (3-point): (f[i+1] - f[i-1]) / (2·dx)
Order 4 (5-point): (f[i-2]/12 - 2f[i-1]/3 + 2f[i+1]/3 - f[i+2]/12) / dx
Order 6 (7-point): (-f[i-3]/60 + 3f[i-2]/20 - 3f[i-1]/4 + 3f[i+1]/4 - 3f[i+2]/20 + f[i+3]/60) / dx

Reference: Fornberg (1988), Mathematics of Computation 51(184), Table 1.

Parameters:

data (np.ndarray) – Field values, shape must match grid.shape.
axis (int) – Spatial axis along which to differentiate.
grid (GridInfo) – Grid descriptor.
bc (str or tuple of str, optional) – Boundary conditions. Defaults to grid periodicity inference.

Returns:

Gradient array, same shape as data.

Return type:

np.ndarray

tidal.solver.operators.directional_laplacian(data, axis, grid, bc=None)[source]#

Second derivative ∂²f/∂x_i² (FD stencil or spectral FFT).

Parameters:

data (np.ndarray)
axis (int)
grid (GridInfo)
bc (BCSpec | None)

Return type:

np.ndarray

tidal.solver.operators.laplacian(data, grid, bc=None)[source]#

Full Laplacian nabla^2 f = sum_i d^2f/dx_i^2 (FD or spectral).

Parameters:

data (np.ndarray)
grid (GridInfo)
bc (BCSpec | None)

Return type:

np.ndarray

tidal.solver.operators.cross_derivative(data, axis1, axis2, grid, bc=None)[source]#

Mixed second derivative d^2f/(dx_i dx_j) (FD or spectral).

FD mode: successive application of the 1D gradient operator along each axis, inheriting the current FD order from set_fd_order().

Spectral mode: successive FFT-based gradients (exact for band-limited functions).

Reference: LeVeque (2007), “Finite Difference Methods”, SIAM, S1.4.

Parameters:

data (np.ndarray)
axis1 (int)
axis2 (int)
grid (GridInfo)
bc (BCSpec | None)

Return type:

np.ndarray

tidal.solver.operators.biharmonic(data, grid, bc=None)[source]#

Biharmonic operator nabla^4 f = nabla^2(nabla^2 f) (FD or spectral).

Parameters:

data (np.ndarray)
grid (GridInfo)
bc (BCSpec | None)

Return type:

np.ndarray

tidal.solver.operators.identity(data, grid, bc=None)[source]#

Identity operator — returns data unchanged (no copy).

Parameters:

data (np.ndarray)
grid (GridInfo)
bc (BCSpec | None)

Return type:

np.ndarray

tidal.solver.operators.operator_min_dim(name)[source]#

Return the minimum spatial dimension required by name.

Raises:: ValueError – If name is not a recognized operator.
Parameters:: name (str)
Return type:: int

tidal.solver.operators.apply_operator(name, data, grid, bc=None)[source]#

Look up name in the registry and apply it to data.

Raises:

ValueError – If name is not a known operator.

Parameters:

name (str)
data (np.ndarray)
grid (GridInfo)
bc (BCSpec | None)

Return type:

np.ndarray

tidal.solver.progress module#

Simulation progress bar for TIDAL solvers (tqdm-based).

Provides a thin wrapper around tqdm that tracks simulation time t and displays elapsed time, ETA, and percentage completion. Designed to be called from solver step loops or RHS functions with negligible overhead.

Key features:

disable=True makes tqdm a complete no-op (zero overhead when quiet).
mininterval=0.5 throttles terminal writes to ~2 Hz.
Monotonic t tracking prevents backwards progress from IDA Jacobian finite-difference evaluations.
Output goes to sys.stderr so it doesn’t interfere with piped stdout.

class tidal.solver.progress.SimulationProgress(t_start, t_end, *, solver_name='', disable=False)[source]#

Bases: object

tqdm-based progress bar driven by simulation time.

Parameters:

t_start (float) – Simulation start time.
t_end (float) – Simulation end time.
solver_name (str) – Label for the progress bar (e.g. "IDA", "CVODE").
disable (bool) – If True, completely disables output (zero overhead).

set_phase(phase)[source]#

Update the progress bar description to show current phase.

Useful for indicating solver initialization before the step loop begins, so users know the solver hasn’t hung on large systems.

Parameters:: phase (str)
Return type:: None

update(t)[source]#

Report current simulation time.

Monotonic: ignores t values less than the previous maximum, which can occur during IDA Jacobian finite-difference evaluations. Clamps to total to prevent tqdm warnings when leapfrog overshoots.

Parameters:: t (float)
Return type:: None

finish()[source]#

Set progress to 100% and close the bar.

Return type:: None

tidal.solver.rhs module#

Unified RHS evaluation for TIDAL solvers.

RHSEvaluator applies spatial operators with resolved coefficients, replacing the duplicated inner loops in ida.py and leapfrog.py.

Depends on FieldSet (Phase 1) and CoefficientEvaluator (Phase 2).

class tidal.solver.rhs.RHSEvaluator(spec, grid, coeff_eval, bc=None)[source]#

Bases: object

Evaluate RHS of field equations with resolved coefficients.

Applies spatial operators to field data and multiplies by resolved coefficients (constant, parameter-overridden, position-dependent, or time-dependent).

Parameters:

spec (EquationSystem) – Parsed JSON equation specification.
grid (GridInfo) – Spatial grid.
coeff_eval (CoefficientEvaluator) – Coefficient resolver with caching.
bc (str or tuple of str, optional) – Boundary conditions for spatial operators.

begin_timestep(t)[source]#

Notify the coefficient evaluator of a new timestep.

Parameters:: t (float)
Return type:: None

evaluate(eq_idx, fields, t=0.0)[source]#

Compute RHS for a single equation.

Parameters:

eq_idx (int) – Index of the equation in spec.equations.
fields (FieldSet) – Current field state.
t (float) – Current simulation time.

Returns:

Grid-shaped result array. Warning: the returned array is an internal buffer and may be overwritten by the next call to evaluate(). Callers must copy if they need to persist it.

Return type:

np.ndarray

evaluate_by_field(field_name, fields, t=0.0)[source]#

Compute RHS for the equation governing field_name.

Raises:

KeyError – If field_name has no associated equation.

Parameters:

field_name (str)
fields (FieldSet)
t (float)

Return type:

np.ndarray

tidal.solver.scipy_solver module#

scipy.integrate.solve_ivp wrapper for TIDAL — explicit adaptive ODE solver.

Provides DOP853 (8th-order Dormand-Prince) as default method, which excels for smooth non-stiff wave problems with fewer function evaluations than implicit BDF for the same accuracy.

Also supports implicit methods (Radau, BDF) with Jacobian sparsity from the existing TIDAL sparsity builder.

Reference: Dormand & Prince, “A family of embedded Runge-Kutta formulae”, J. Comp. Appl. Math. 6, 1980.

tidal.solver.scipy_solver.solve_scipy(spec, grid, y0, t_span, *, bc=None, parameters=None, method='DOP853', rtol=1e-08, atol=1e-10, max_step=inf, num_snapshots=101, snapshot_callback=None, progress=None)[source]#

Solve a TIDAL equation system using scipy.integrate.solve_ivp.

Parameters:

spec (EquationSystem) – Parsed equation specification.
grid (GridInfo) – Spatial grid.
y0 (np.ndarray) – Initial state vector (flat).
t_span (tuple[float, float]) – (t_start, t_end).
bc (str or tuple, optional) – Boundary conditions.
parameters (dict[str, float], optional) – Runtime parameter overrides for symbolic coefficients.
method (str) – Integration method: 'DOP853' (default), 'RK45', 'Radau', 'BDF', 'RK23', 'LSODA'.
rtol (float) – Relative and absolute tolerances.
atol (float) – Relative and absolute tolerances.
max_step (float) – Maximum step size. np.inf (default) = unbounded.
num_snapshots (int) – Number of output time points.
snapshot_callback (callable, optional) – Called as callback(t, y) at each output time.
progress (SimulationProgress | None)

Returns:

Result dictionary with keys: t, y, success, message.

Return type:

Warns:

UserWarning – If constraint fields (time_order=0) are present — they remain frozen at initial values.

tidal.solver.sparsity module#

Analytical Jacobian sparsity pattern for TIDAL IDA solver.

Builds the sparsity structure of the IDA Jacobian J = dF/dy + cj * dF/dyp directly from the equation spec and operator stencils — no numerical probing. This avoids the O(N^2) cost of sksundae.jacband.j_pattern() (which allocates a dense N x N array) and gives exact structure in milliseconds.

The sparsity pattern is returned as a scipy.sparse.csc_matrix suitable for passing to sksundae.ida.IDA(linsolver='sparse', sparsity=pattern), which uses SuperLU_MT for direct factorisation.

Reference: SuperLU_MT, Li & Demmel, ACM TOMS 29(2), 2003.

tidal.solver.sparsity.operator_stencil_offsets(name, ndim)[source]#

Return the grid-coordinate offsets for operator name.

When spectral mode is active (get_spectral() == True), all spatial operators return dense coupling (all possible offsets) because FFT-based differentiation couples every grid point to every other. Identity and first_derivative_t remain local.

Parameters:

name (str) – Operator name (e.g. "laplacian", "gradient_x").
ndim (int) – Number of spatial dimensions.

Returns:

Offset tuples, each of length ndim.

Return type:

list[tuple[int, …]]

tidal.solver.sparsity.build_jacobian_sparsity(spec, layout, grid, _bc)[source]#

Build the analytical Jacobian sparsity pattern for IDA.

Constructs the nonzero structure of J = dF/dy + cj * dF/dyp from the equation system’s operator stencils and slot-coupling structure. No numerical probing – the pattern is exact and built in O(nnz) time.

Parameters:

spec (EquationSystem) – Parsed equation specification.
layout (StateLayout) – State vector layout.
grid (GridInfo) – Spatial grid (periodic flags used for wrapping).
_bc (str or tuple of str, optional) – Boundary conditions string (kept for API compatibility with ida.py call site; periodicity is read from grid.periodic).

Returns:

Binary CSC sparsity pattern of shape (N_total, N_total).

Return type:

SparseMatrix

tidal.solver.state module#

State vector layout and flat ↔ field-dict conversions for TIDAL solvers.

IDA and leapfrog operate on flat numpy arrays. This module provides:

StateLayout: describes which slots map to which fields/velocities
state_to_flat / flat_to_fields: convert between flat vectors and named field dictionaries

class tidal.solver.state.SlotInfo(name, field_name, kind, time_order, dynamical_index=None)[source]#

Bases: object

Metadata for a single slot in the flat state vector.

Each slot holds grid.num_points contiguous values in C-order.

Parameters:

name (str)
field_name (str)
kind (str)
time_order (int)
dynamical_index (int | None)

name#

Field or velocity name (e.g. "phi_0" or "v_phi_0").

Type:: str

field_name#

The physical field this slot belongs to (always the field name, even for velocity slots).

Type:: str

kind#

One of "field", "velocity", "constraint".

Type:: str

time_order#

Original LHS time-derivative order (0, 1, or 2).

Type:: int

dynamical_index#

Index into the dynamical-field subset (for K matrix lookup). None for constraints and first-order fields.

Type:: int | None

name: str#

field_name: str#

kind: str#

time_order: int#

dynamical_index: int | None = None#

class tidal.solver.state.StateLayout(slots, num_points, field_slot_map, velocity_slot_map, dynamical_fields)[source]#

Bases: object

Describes the mapping between flat state vector and named fields.

Built from an EquationSystem spec. The flat vector is partitioned into contiguous blocks of num_points elements, one per slot.

Parameters:

slots (tuple[SlotInfo, ...])
num_points (int)
field_slot_map (dict[str, int])
velocity_slot_map (dict[str, int])
dynamical_fields (tuple[str, ...])

slots#

Ordered slot descriptors.

Type:: tuple[SlotInfo, …]

num_points#

Grid points per slot (grid.num_points).

Type:: int

field_slot_map#

Maps field name → slot index.

Type:: dict[str, int]

velocity_slot_map#

Maps field name → velocity slot index (second-order fields only).

Type:: dict[str, int]

dynamical_fields#

Names of dynamical fields (time_order >= 2), in order.

Type:: tuple[str, …]

slots: tuple[SlotInfo, ...]#

num_points: int#

field_slot_map: dict[str, int]#

velocity_slot_map: dict[str, int]#

dynamical_fields: tuple[str, ...]#

property momentum_slot_map: dict[str, int]#: Alias for velocity_slot_map (transition aid).

classmethod from_spec(spec, num_points)[source]#

Build layout from an equation system specification.

Follows the same ordering as py-pde FieldCollection: for each equation, emit a field slot; if second-order, also emit a velocity slot immediately after.

Parameters:

spec (EquationSystem)
num_points (int)

Return type:

StateLayout

property slot_name_to_idx: dict[str, int]#: Map from slot name to slot index. Cached on frozen dataclass.

property total_size: int#: Total flat vector length (num_slots * num_points).

property num_slots: int#: Number of slots in the state vector.

slot_slice(slot_idx)[source]#

Return the flat-array slice for a given slot index.

Parameters:: slot_idx (int)
Return type:: slice

property velocity_slot_groups: tuple[tuple[int, slice, str], ...]#: Pre-computed (slot_idx, flat_slice, field_name) for velocity slots.

property dynamical_field_slot_groups: tuple[tuple[int, slice, int], ...]#: Pre-computed (slot_idx, flat_slice, vel_slot_idx) for 2nd-order fields.

property drift_slot_pairs: tuple[tuple[slice, slice], ...]#

Pre-computed (field_slice, vel_slice) for zero-copy drift.

Allows y[field_slice] += dt * y[vel_slice] without copying velocity data into a separate buffer.

property first_order_slot_groups: tuple[tuple[int, slice, str], ...]#: Pre-computed (slot_idx, flat_slice, field_name) for 1st-order fields.

property constraint_slot_groups: tuple[tuple[int, slice, str], ...]#: Pre-computed (slot_idx, flat_slice, field_name) for constraint slots.

property algebraic_indices: list[int]#

Flat indices of algebraic (constraint) variables for IDA.

IDA needs to know which entries in the state vector are algebraic (not differential) so it can handle them appropriately.

tidal.solver.state.state_to_flat(fields, layout)[source]#

Pack named field arrays into a single flat vector.

Parameters:

fields (dict[str, np.ndarray]) – Mapping from slot name (e.g. "phi_0", "v_phi_0") to grid-shaped arrays.
layout (StateLayout) – State layout descriptor.

Returns:

Flat vector of length layout.total_size.

Return type:

np.ndarray

tidal.solver.state.flat_to_fields(y, layout, shape)[source]#

Unpack a flat vector into named field arrays.

Parameters:

y (np.ndarray) – Flat state vector.
layout (StateLayout) – State layout descriptor.
shape (tuple[int, ...]) – Grid shape to reshape each slot’s data into.

Returns:

Mapping from slot name to grid-shaped array.

Return type:

dict[str, np.ndarray]

tidal.solver.validation module#

Solver-agnostic validation for TIDAL equation specifications.

Extracted from pde_builder.py so that IDA, leapfrog, and any future solver can share the same validation logic without importing py-pde.

All functions are module-level (no shared state) and raise ValueError on invalid specs per the fail-fast-and-loud convention.

tidal.solver.validation.validate_operator_dimensions(spec)[source]#

Check that all operators are compatible with the spec’s spatial dimension.

Raises:: ValueError – If any operator requires more dimensions than spec.spatial_dimension.
Parameters:: spec (EquationSystem)
Return type:: None

tidal.solver.validation.validate_field_references(spec)[source]#

Check that all term field references point to valid fields.

Raises:: ValueError – If a field reference is invalid.
Parameters:: spec (EquationSystem)
Return type:: None

tidal.solver.validation.check_cfl_stability(spec, grid, dt)[source]#

Check CFL stability condition for explicit time-steppers.

Returns a list of warning strings (empty if all clear). The CFL condition for the wave equation is dt <= dx / c where c is the maximum wave speed (estimated from the laplacian coefficient).

Parameters:

spec (EquationSystem)
grid (GridInfo)
dt (float)

Return type:

tidal.solver.validation.check_mass_sign(coeff_eval, spec)[source]#

Check for sign-changing position-dependent mass terms.

Returns a list of warning strings for tachyonic diagnostics.

Parameters:

coeff_eval (CoefficientEvaluator)
spec (EquationSystem)

Return type: