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 centers 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-center coordinates, shape (*shape, ndim).

For a 1D grid with bounds (0, 10) and shape (4,), the cell centers are at [1.25, 3.75, 6.25, 8.75] (centered 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-centered

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, ...]#

exception tidal.solver.KineticEvaluationError[source]#

Bases: RuntimeError

Raised when a kinetic_coefficient_symbolic cannot be resolved to a number.

The modal builders need a concrete M to populate the mass matrix. When evaluation fails — an unbound symbol (a parameter missing from --param), an unsupported construct, or a NaN/Inf result — there is no correct fallback: proceeding with M = 1 silently changes the physics (GH #447). Raised instead of a bare ValueError on purpose: the conversion-stability probe catches (LinAlgError, ValueError) and converts it into a “tachyonic” verdict, which would mislabel a configuration error as a physics result.

class tidal.solver.PerturbativeResult(orders=<factory>, total=<factory>, validity=<factory>, spec=None, full_spec=None)[source]#

Bases: object

Output of PerturbativeSolver.solve().

Parameters:

orders (list[Any])
total (dict[str, Any])
validity (dict[str, Any])
spec (EquationSystem | None)
full_spec (EquationSystem | None)

orders#

Per-order SolverResult-shaped dicts. orders[0] is Pass 0 (base), orders[n] is the Pass-n correction for n ≥ 1.

Type:: list[dict]

total#

Combined solution Σ orders[n].y in physical space, aligned with orders[0].t.

Type:: dict

validity#

Diagnostic metrics from _compute_validity():

"validity_param" — max_k( ε · ω · t_end ) across all small parameters and simulated modes. The dimensionless secular-phase-error bound. < 0.1 is safe, > 1 breaks the EFT. #284 corrected this from the previous (over-counting) ε·ω²·t formula.
"omega_max" — maximum simulated frequency max|λ| (so damped/tachyonic modes register at their full magnitude).
"warn_level" — "ok", "warn", or "error". The overall level is the worse of correction_level and base_level.
"correction_level" — warn band for the secular-error side only (thresholds: 0.1, 1.0).
"base_stability_param" — max Re(λ) · t_end. A strictly-positive value means the base theory has an unstable mode (tachyon, Jeans, ghost). > 30 → overflow; > 1e-10 → warn-level #285.
"base_level" — warn band for the base-stability side.
"max_real_lambda" — max Re(λ) across all blocks.
"dominant_tachyon_field" — field/velocity slot owning the most-positive Re(λ) (or None if base is stable).
"correction_drops" — list of diagnostic records for any correction terms or equations that _build_source_matrix_k could not route (e.g. target field has no dynamical slot, or equation is for a demoted constraint field — see #272). Empty list on a correct pipeline for every shipped JSON at shipped parameters; any non-empty list is a signal to investigate before trusting Pass 1 magnitudes.

Type:: dict

spec#

The base spec (post-Gap-B demotion) whose layout matches total["y"] and each orders[n]["y"]. Callers that hand total to SimulationData.from_solver_result() must use this spec, NOT the unreduced full spec — the demoted layout has a different slot ordering and a reduced component set. The CLI previously rebound spec = solver.base_spec after the solve; the explicit attribute removes that workaround (#276).

Type:: EquationSystem | None

full_spec#

The unreduced spec as loaded from JSON (pre-Gap-B). Kept for diagnostic purposes — e.g. inspecting original order_in_eps tags or reconstructing the implicit-dynamical field’s kinetic coefficient.

Type:: EquationSystem | None

full_spec: EquationSystem | None = None#

spec: EquationSystem | None = None#

orders: list[Any]#

total: dict[str, Any]#

validity: dict[str, Any]#

class tidal.solver.PerturbativeSolver(spec)[source]#

Bases: object

Orchestrate Pass 0 + Pass n corrections for a tagged EquationSystem.

Usage#

>>> solver = PerturbativeSolver(spec)
>>> res = solver.solve(y0, grid, t_span, order=1, parameters={"b5": 0.01})
>>> q_total = res.total["y"]  # combined physical-space trajectory

The solver is a thin orchestrator over solve_modal() with return_eigendata=True and solve_modal_pass1(). It assumes the modal backend; non-modal backends are out of scope for the v6 plan (tracked in a follow-up issue).

has_corrections()[source]#

Return True when the spec contains any order_in_eps > 0 terms.

Return type:: bool

property max_order: int#: Maximum order_in_eps available in the spec.

solve(y0, grid, t_span, *, order=1, parameters=None, num_snapshots=101, small_parameters=None)[source]#

Solve iteratively to the requested order.

Semantics of the result#

result.orders[n] is a mathematical component of the perturbative expansion, not a physical trajectory. orders[0] is the order-0 base evolution; orders[n] for n >= 1 is the Duhamel-integrated correction. All physical observables (energy, conversion probability, conserved quantities) must be measured on ``result.total``, never on individual orders[n].

Scope and limitations (v6 / Phase 2)#

The iterative Parker-Simon path used here resolves the equation side cleanly for theories with declared [perturbation] small_parameters. However, the Lagrangian Perturbative Substitution (LPS, see tidal/wolfram/PerturbativeReduction.wl) that produces the canonical Hamiltonian throws on the constraint-promotion case — when the small parameter promotes an algebraic-constraint field at epsilon=0 to a higher-derivative dynamical field at epsilon!=0. The b5*Rtilde^2 PGT torsion family (graviton_torsion, torsion_gertsenshtein, torsion_gertsenshtein_combined) triggers this abort. Their equation-side dynamics computed here remain valid; only the canonical Hamiltonian is unavailable via the standard methods. See docs/tex/perturbative_reduction_constraint_barrier.tex and issue #321.

param y0:: Flat initial state for Pass 0. Pass 1+ starts at zero.
type y0:: ndarray
param grid:: Spatial grid (must satisfy modal-solver eligibility).
type grid:: GridInfo
param t_span:: (t_start, t_end).
type t_span:: (float, float)
param order:: Highest order of correction to compute. 0 returns just the base solution.
type order:: int, default 1
param parameters:: Runtime parameter overrides. Small parameters (those tracked in order_in_eps) should appear here with their physical values so the correction coefficients resolve.
type parameters:: dict, optional
param num_snapshots:: Number of output time points.
type num_snapshots:: int
param small_parameters:: Names of the small parameters for the validity monitor. When absent, read from spec.metadata.get("perturbation", {}).get("small_parameters", []).
type small_parameters:: list[str], optional
returns:: Per-order SolverResults + the combined trajectory + the validity diagnostics.
rtype:: PerturbativeResult
raises ValueError:: If order exceeds the maximum order_in_eps carried by the spec.
raises NotImplementedError:: If order >= 2 — Pass 2 against q¹(t) is not yet implemented. See #273.
raises RuntimeError:: If a Pass-n solver call returns success=False. Failing loud rather than silently producing a degraded result.

Parameters:

y0 (NDArray[np.float64])
grid (GridInfo)
t_span (tuple[float, float])
order (int)
parameters (dict[str, float] | None)
num_snapshots (int)
small_parameters (list[str] | None)

Return type:

PerturbativeResult

Parameters:: spec (EquationSystem)

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, implicit_dynamical_fields=frozenset({}))[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, ...])
implicit_dynamical_fields (frozenset[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 — or in the implicit-dynamical sector (GH #457) — also emit a velocity slot immediately after. There is ONE layout rule for all backends; backends that cannot represent the implicit-dynamical sector refuse at their entry (RHSEvaluator / modal_jax).

Parameters:

spec (EquationSystem)
num_points (int)

Return type:

StateLayout

implicit_dynamical_fields: frozenset[str] = frozenset({})#

Order-0 fields in the implicit-dynamical sector (GH #457).

These have field + velocity slots like order >= 2 fields (their inter-constraint time-derivative couplings make them dynamical), but SlotInfo.time_order keeps the true LHS fact (0). Consumers that key “frozen constraint” on time_order == 0 must exclude these.

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.

The check is on the summed self-mass of each equation, not on individual terms: several identity terms may act on the same field, and only their sum is the physical mass.

Returns a list of warning strings for tachyonic diagnostics.

Parameters:

coeff_eval (CoefficientEvaluator)
spec (EquationSystem)

Return type:

tidal.solver.find_independent_blocks(A, threshold=1e-14)[source]#

Find independent (decoupled) blocks in an evolution matrix.

Analyzes the sparsity pattern of A: slots i and j are coupled if |A[i,j]| > threshold or |A[j,i]| > threshold. Returns a list of slot-index groups (connected components).

This prevents degenerate-eigenvalue mixing when np.linalg.eig processes a block-diagonal matrix with repeated eigenvalues across independent blocks — a common situation in symmetric multi-field theories (e.g. Gertsenshtein h₅↔a₁ + h₇↔a₂).

Parameters:

A (ndarray[tuple[Any, ...], dtype[complex128]])
threshold (float)

Return type:

list[list[int]]

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.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, return_eigendata=False, _zero_nyquist=True)[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.
return_eigendata (bool, optional) – When True, attaches a "eigendata" key to the returned SolverResult containing per-block eigendecomposition data (V, D_diag, V_inv, alpha), the mode_k grid, the dynamical state_layout, and (when constraints are present) a schur_ops sub-dict with the recovery matrix and the original → reduced slot map. Used by the v6 iterative PerturbativeSolver to drive a Pass 1 closed-form Duhamel evaluation that reuses Pass 0’s eigendecomposition. Not supported for position-dependent coefficient systems (they use expm_multiply without an explicit eigendecomposition); attempting to combine the two raises NotImplementedError.
_zero_nyquist (bool)

Returns:

Dict with keys: t, y, success, message. When return_eigendata=True also has "eigendata".

Return type:

Raises:

ValueError – If spec contains time-derivative orders > 2, or incompatible options are given for the position-dependent coefficient path.

tidal.solver.solve_modal_pass1(eigendata, correction_spec, grid, t_eval, *, parameters=None)[source]#

Solve the Pass 1 correction dy⁽¹⁾/dt = A·y⁽¹⁾ + M_src·y⁰(t).

Reuses eigendata from a prior solve_modal(..., return_eigendata=True) call on the base (order-0) spec. correction_spec is the spec.filter_by_order(n) spec containing only the order-n RHS terms. The closed-form Duhamel kernel is applied per mode in the eigenbasis of the base operator — no ODE integration.

Parameters:

eigendata (dict) – Output of Pass 0 with return_eigendata=True. Must include "blocks", "mode_k", "state_layout", and (for constraint systems) "schur_ops".
correction_spec (EquationSystem) – The correction-only equation system (RHS terms filtered by order_in_eps). The LHS structure is unused — only the right-hand side operator terms matter for the source.
grid (GridInfo) – Same grid used in Pass 0 (must match the eigendata’s mode_k).
t_eval (array of floats) – Evaluation times. Pass 1’s IC is zero, so t_eval[0] can be any reference time; output at t = t_eval[0] is zero.
parameters (dict, optional) – Runtime parameter overrides for the correction coefficients.

Returns:

t, y, success, message (from SolverResult) plus:

y_hat_dyn — Fourier-space dynamical output with shape (n_snapshots, n_slots, n_modes). Consumed by PerturbativeSolver() to recover constraint fields via Schur (v6 Gap C).
correction_drops — diagnostic records for any correction terms or equations that could not be routed. Aggregated by the driver onto PerturbativeResult.validity (#272).

Return type:

PerturbativePass1Result

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 factorization. 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 factorization. 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.

property parameters: Mapping[str, float]#

Read-only view of the parameter values this evaluator resolves against.

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]

classify_bsm_separability(bsm_symbols)[source]#

Detect which position-dependent RHS terms factor as c_BSM(α) · c_geom(x).

GH #384 Phase A′ infrastructure. Re-run if bsm_symbols changes. See tidal/symbolic/_separable.py for the AST-walker contract.

Parameters:: bsm_symbols (set[str])
Return type:: None

bsm_separable_factors(eq_idx, term_idx)[source]#

Return (c_BSM_str, c_geom_str) if the term was classified separable.

Populated by a prior classify_bsm_separability call; else None.

Parameters:

eq_idx (int)
term_idx (int)

Return type:

tuple[str, str] | None

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-centered

axes_coords(axis)[source]#

1-D array of cell centers 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-center coordinates, shape (*shape, ndim).

For a 1D grid with bounds (0, 10) and shape (4,), the cell centers are at [1.25, 3.75, 6.25, 8.75] (centered 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

Kinetic-coefficient (mass-matrix) handling — four primary dispatches#

Equations with a non-trivial kinetic_coefficient_symbolic (M ẍ = K x rather than ẍ = K x) reach the solver via four distinct code paths plus the Pass 1 correction-matrix builder. All five must consume M identically, or cross-path regressions appear silently. Two such regressions to date:

GH #367 (fixed v0.42.0, 2026-05-19): path 3 (pos-dep convolution) silently dropped M⁻¹ for theories with non-trivial kinetic_coefficient_symbolic.
GH #379 (fixed 2026-05-24): path 4 (the pos-dep + constraint case) did not exist; the dispatch silently routed those theories to path 3 which has no constraint handling, causing 10²⁵² divergence at any nonzero coupling. Now handled by _build_convolution_matrix_with_constraints() (convolution analog of path 2’s Schur elimination).

Fast path — _build_per_mode_matrices() (called when needs_reduction = False and no position-dependent coefficients). Builds the first-order evolution matrix directly as A[velocity_slot, target_slot] = M⁻¹(field) · coeff · multiplier. M⁻¹ is pre-computed once from tidal.solver._kinetic.build_inverse_kinetic_diag() and folded in via the shared tidal.solver._kinetic.velocity_row_scale() helper.
Generalized-eig path — _build_evolution_matrices() (called when the spec has constraints or time-derivative RHS operators, and no position-dependent coefficients). Populates M_mat from kinetic_coefficient_symbolic on the diagonal and solves (K − λM) v = 0 via scipy.linalg.eig(A, B) with QZ decomposition. Handles rank-deficient M via null-space projection and Schur elimination for hidden algebraic constraints. Structurally distinct from paths 1, 3, and 4: writes M to the generalized eigenvalue problem rather than post-multiplying a velocity-row contribution, so it intentionally does NOT call velocity_row_scale().
Position-dependent convolution path — _build_convolution_matrix() (called when the spec has position-dependent coefficients — RHS or kinetic, GH #427 — and no algebraic constraints). Builds A[m, m'] = ĉ(m − m') · op_mult(k_{m'}) over the full Fourier-mode block, threading M⁻¹ through both the constant-coefficient and the convolution-coupling paths via velocity_row_scale(). A position-dependent kinetic yields an ndarray scale — a REAL-SPACE M⁻¹(x) profile folded into the coefficient before the FFT (_conv_block_with_kinetic()); it must never multiply a k-space block. Rejects (via assertion) systems with time_derivative_order == 0 equations to make the GH #379 bug class a static error.
Position-dependent + constraint convolution path — _build_convolution_matrix_with_constraints() (called when both pos-dep coefficients and algebraic constraints are present). Builds convolution sub-matrices A_dd, A_dc_field/vel, K_cd, K_cc and Schur-eliminates: recovery = -K_cc⁻¹ · K_cd, A_reduced = A_dd + A_dc_field · recovery (plus B_lhs pre-solve for v_constraint and Newton correction for ẍ_constraint references in dyn RHS). Mirrors path 2’s structure but with full 2D convolution matrices instead of per-mode 3-tensors. Applies M⁻¹ on dynamical velocity rows via velocity_row_scale(), same contract as paths 1 and 3.

Plus the Pass 1 correction-matrix builder _build_pass1_source_matrices() (the perturbative Duhamel correction path; see tidal.solver.perturbative_driver.PerturbativeSolver), which also emits velocity-row entries into a per-mode source matrix and likewise consumes M⁻¹ through velocity_row_scale().

Cross-builder contract (enforced by the regression test tests.test_solver_kinetic_consistency.TestAllModalPathsRespectKinetic): any matrix builder that populates A[velocity_slot, target_slot] entries for a second-order equation MUST apply tidal.solver._kinetic.velocity_row_scale(field_name, m_inv)() to the contribution. Path 2 (generalized-eig) is the documented exception — it writes M to a separate matrix that participates in the eigenvalue problem. Adding a new matrix builder without satisfying this contract was exactly the GH #367 regression; please don’t repeat it. Path 4 (pos-dep + constraint) follows the same contract.

The choice between dispatches is transparent to callers of solve_modal() but load-bearing for reviewers: all five sites must be updated together when the kinetic handling changes. The time-domain backends (cvode/ida/ leapfrog/scipy) use tidal.solver._kinetic.build_inverse_kinetic_diag() as a shared entry point, keeping the six backends (modal-fast, modal-genEig, modal-conv, modal-conv-Schur, cvode, ida, leapfrog/scipy) numerically consistent per tests.test_solver_kinetic_consistency.

Algorithm paths for 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 6. No position-dependent kinetic coefficients (GH #421) — position-

dependent RHS coefficients are fine (handled via the ĉ(k−k′) convolution paths), but a position-dependent M(x) needs the mass-side convolution M̂(k−k′), which is not implemented (GH #427). Auto-selection falls through to the time-domain backends, which evaluate such kinetics on the grid (GH #382).

Parameters:

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

Return type:

bool

tidal.solver.modal.find_independent_blocks(A, threshold=1e-14)[source]#

Find independent (decoupled) blocks in an evolution matrix.

Analyzes the sparsity pattern of A: slots i and j are coupled if |A[i,j]| > threshold or |A[j,i]| > threshold. Returns a list of slot-index groups (connected components).

Parameters:

A (ndarray[tuple[Any, ...], dtype[complex128]])
threshold (float)

Return type:

list[list[int]]

exception tidal.solver.modal.SingularPencilError[source]#

Bases: RuntimeError

The per-mode pencil λB − A is singular (det ≡ 0): true gauge freedom.

Raised by _pencil_deflate() when an (α, β) ≈ (0, 0) pair is found — some state direction is genuinely undetermined by the equations. There is no correct evolution for such a direction; a silent choice (freeze, minimum-norm, Tikhonov) is a wrong theory with exit code 0. Gauge-fix the theory ([[gauge]] in the TOML) or file the spec on GH #457’s follow-up.

class tidal.solver.modal.PencilDiagnostics[source]#

Bases: object

Pinned (gauge-quotiented) directions recorded during deflation.

GH #468 “pin + certify”: every state direction the gauge quotient sets to zero (the explicit min-norm choice) is recorded here, and the builders turn the records into pin_overlap[mode, slot] — the squared norm of the pinned subspace’s rows on that state slot. A measurement whose field support carries no pinned content is provably independent of the pin (certified); one that does is flagged with the magnitude by tidal measure instead of being reported as if it were gauge-invariant.

record(tag, pinned, tau)[source]#

Record one deflation’s pinned basis ((n, w_dim), may be empty).

Parameters:

tag (tuple[str, int])
pinned (ndarray[tuple[Any, ...], dtype[complex128]])
tau (float)

Return type:

None

to_dict()[source]#

JSON-ready summary for run metadata (solver_diagnostics).

Return type:: dict[str, Any]

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, return_eigendata=False, _zero_nyquist=True)[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.
return_eigendata (bool, optional) – When True, attaches a "eigendata" key to the returned SolverResult containing per-block eigendecomposition data (V, D_diag, V_inv, alpha), the mode_k grid, the dynamical state_layout, and (when constraints are present) a schur_ops sub-dict with the recovery matrix and the original → reduced slot map. Used by the v6 iterative PerturbativeSolver to drive a Pass 1 closed-form Duhamel evaluation that reuses Pass 0’s eigendecomposition. Not supported for position-dependent coefficient systems (they use expm_multiply without an explicit eigendecomposition); attempting to combine the two raises NotImplementedError.
_zero_nyquist (bool)

Returns:

Dict with keys: t, y, success, message. When return_eigendata=True also has "eigendata".

Return type:

Raises:

ValueError – If spec contains time-derivative orders > 2, or incompatible options are given for the position-dependent coefficient path.

tidal.solver.modal.solve_modal_pass1(eigendata, correction_spec, grid, t_eval, *, parameters=None)[source]#

Solve the Pass 1 correction dy⁽¹⁾/dt = A·y⁽¹⁾ + M_src·y⁰(t).

Parameters:

eigendata (dict) – Output of Pass 0 with return_eigendata=True. Must include "blocks", "mode_k", "state_layout", and (for constraint systems) "schur_ops".
correction_spec (EquationSystem) – The correction-only equation system (RHS terms filtered by order_in_eps). The LHS structure is unused — only the right-hand side operator terms matter for the source.
grid (GridInfo) – Same grid used in Pass 0 (must match the eigendata’s mode_k).
t_eval (array of floats) – Evaluation times. Pass 1’s IC is zero, so t_eval[0] can be any reference time; output at t = t_eval[0] is zero.
parameters (dict, optional) – Runtime parameter overrides for the correction coefficients.

Returns:

t, y, success, message (from SolverResult) plus:

y_hat_dyn — Fourier-space dynamical output with shape (n_snapshots, n_slots, n_modes). Consumed by PerturbativeSolver() to recover constraint fields via Schur (v6 Gap C).
correction_drops — diagnostic records for any correction terms or equations that could not be routed. Aggregated by the driver onto PerturbativeResult.validity (#272).

Return type:

PerturbativePass1Result

tidal.solver.modal_jax module#

JAX-accelerated modal solver — phase 1 (constant-coefficient path only).

Replaces the per-mode scipy.linalg.expm loop in tidal.solver.modal._evolve_per_mode_pade() with jax.vmap over jax.scipy.linalg.expm, and evolves snapshots via a JIT-compiled jax.lax.scan. All setup (FFT, per-mode matrix construction, block detection) reuses the existing scipy path without modification.

Installation:

uv sync --extra jax

Requires jax>=0.4.20 and jaxlib>=0.4.20. Double-precision is enabled automatically (jax_enable_x64 = True).

Phase-1 restrictions (raise NotImplementedError): - Position-dependent coefficient systems (_evolve_full_matrix path). - Systems requiring Schur constraint elimination (needs_reduction). - return_eigendata=True (eigendecomposition data not exposed by JAX path).

The TIDAL_MODAL_SPARSE_IC env var is ignored in this path. JAX computes expm for all modes simultaneously via vmap; masking out zero-IC modes is not cost-free under vmap (all modes are traced uniformly). The eigenvalue-based divergence pre-check (inherited from the scipy path) still runs outside JIT and catches tachyonic instabilities.

The runtime snapshot-by-snapshot divergence guard is not applied inside the JIT-compiled evolution loop; the pre-check is the only divergence gate.

References

Bradbury et al. (2018). JAX: composable transformations of Python+NumPy. Higham (2009). SIAM Review 51(4):747-764 — matrix exponential.

tidal.solver.modal_jax.solve_modal_jax(spec, grid, y0, t_span, *, bc=None, parameters=None, rtol=1e-08, atol=1e-10, num_snapshots=101, snapshot_callback=None, progress=None, return_eigendata=False, _zero_nyquist=True)[source]#

Solve using the JAX-accelerated modal backend (phase 1).

Matches tidal.solver.modal.solve_modal() in signature and return type. Phase-1 restrictions (see module docstring) raise NotImplementedError rather than silently falling back to scipy, so callers that request an unsupported path get a clear error.

Parameters:

spec (EquationSystem)
grid (GridInfo)
y0 (np.ndarray)
t_span (tuple[float, float])
bc (BCSpec | None)
parameters (dict[str, float] | None)
rtol (float)
atol (float)
num_snapshots (int)
snapshot_callback (Callable[[float, np.ndarray], None] | None)
progress (SimulationProgress | None)
return_eigendata (bool)
_zero_nyquist (bool)

Return type:

:param : :param snapshot_callback: Same semantics as solve_modal(). :param progress: Same semantics as solve_modal(). :param return_eigendata: Same semantics as solve_modal(). :param _zero_nyquist: Same semantics as solve_modal().

Raises:

NotImplementedError – For position-dependent coefficients, constraint systems, or return_eigendata=True.
SimulationDivergedError – When eigenvalue pre-check detects tachyonic instability.
ValueError – When time_derivative_order > 2.

Parameters:

spec (EquationSystem)
grid (GridInfo)
y0 (np.ndarray)
t_span (tuple[float, float])
bc (BCSpec | None)
parameters (dict[str, float] | None)
rtol (float)
atol (float)
num_snapshots (int)
snapshot_callback (Callable[[float, np.ndarray], None] | None)
progress (SimulationProgress | None)
return_eigendata (bool)
_zero_nyquist (bool)

Return type:

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.perturbative_driver module#

Iterative perturbative driver (v6 plan, Stage 5).

Given an EquationSystem whose terms carry an order_in_eps annotation (emitted by ExportJSON.wl when the theory has a [perturbation] section), this module composes the modal solver’s Pass 0 base evolution with one or more Pass 1 Duhamel corrections into a single PerturbativeResult.

The full architecture is documented in docs/tex/perturbative_reduction_design.tex (“Perturbative Reduction: Implementation Design”); the driver’s constant-coefficient base-spec limitation is tracked in GitHub issue #383. Key guarantees:

Ghost-free by construction: the LHS operator at every order is the base 2nd-order operator; higher-derivative terms appear only as sources on the previous order’s solution.
Theory-agnostic: driven purely by order_in_eps tags and operator strings. No per-theory classification.
Reuses Pass 0’s eigendecomposition for every Pass 1+ call, so the added cost is ~1.5x a single Pass 0 at machine precision.

class tidal.solver.perturbative_driver.PerturbativeResult(orders=<factory>, total=<factory>, validity=<factory>, spec=None, full_spec=None)[source]#

Bases: object

Output of PerturbativeSolver.solve().

Parameters:

orders (list[Any])
total (dict[str, Any])
validity (dict[str, Any])
spec (EquationSystem | None)
full_spec (EquationSystem | None)

orders#

Per-order SolverResult-shaped dicts. orders[0] is Pass 0 (base), orders[n] is the Pass-n correction for n ≥ 1.

Type:: list[dict]

total#

Combined solution Σ orders[n].y in physical space, aligned with orders[0].t.

Type:: dict

validity#

Diagnostic metrics from _compute_validity():

"validity_param" — max_k( ε · ω · t_end ) across all small parameters and simulated modes. The dimensionless secular-phase-error bound. < 0.1 is safe, > 1 breaks the EFT. #284 corrected this from the previous (over-counting) ε·ω²·t formula.
"omega_max" — maximum simulated frequency max|λ| (so damped/tachyonic modes register at their full magnitude).
"warn_level" — "ok", "warn", or "error". The overall level is the worse of correction_level and base_level.
"correction_level" — warn band for the secular-error side only (thresholds: 0.1, 1.0).
"base_stability_param" — max Re(λ) · t_end. A strictly-positive value means the base theory has an unstable mode (tachyon, Jeans, ghost). > 30 → overflow; > 1e-10 → warn-level #285.
"base_level" — warn band for the base-stability side.
"max_real_lambda" — max Re(λ) across all blocks.
"dominant_tachyon_field" — field/velocity slot owning the most-positive Re(λ) (or None if base is stable).
"correction_drops" — list of diagnostic records for any correction terms or equations that _build_source_matrix_k could not route (e.g. target field has no dynamical slot, or equation is for a demoted constraint field — see #272). Empty list on a correct pipeline for every shipped JSON at shipped parameters; any non-empty list is a signal to investigate before trusting Pass 1 magnitudes.

Type:: dict

spec#

Type:: EquationSystem | None

full_spec#

Type:: EquationSystem | None

orders: list[Any]#

total: dict[str, Any]#

validity: dict[str, Any]#

spec: EquationSystem | None = None#

full_spec: EquationSystem | None = None#

class tidal.solver.perturbative_driver.PerturbativeSolver(spec)[source]#

Bases: object

Orchestrate Pass 0 + Pass n corrections for a tagged EquationSystem.

Usage#

>>> solver = PerturbativeSolver(spec)
>>> res = solver.solve(y0, grid, t_span, order=1, parameters={"b5": 0.01})
>>> q_total = res.total["y"]  # combined physical-space trajectory

property max_order: int#: Maximum order_in_eps available in the spec.

has_corrections()[source]#

Return True when the spec contains any order_in_eps > 0 terms.

Return type:: bool

solve(y0, grid, t_span, *, order=1, parameters=None, num_snapshots=101, small_parameters=None)[source]#

Solve iteratively to the requested order.

Semantics of the result#

Scope and limitations (v6 / Phase 2)#

param y0:: Flat initial state for Pass 0. Pass 1+ starts at zero.
type y0:: ndarray
param grid:: Spatial grid (must satisfy modal-solver eligibility).
type grid:: GridInfo
param t_span:: (t_start, t_end).
type t_span:: (float, float)
param order:: Highest order of correction to compute. 0 returns just the base solution.
type order:: int, default 1
param parameters:: Runtime parameter overrides. Small parameters (those tracked in order_in_eps) should appear here with their physical values so the correction coefficients resolve.
type parameters:: dict, optional
param num_snapshots:: Number of output time points.
type num_snapshots:: int
param small_parameters:: Names of the small parameters for the validity monitor. When absent, read from spec.metadata.get("perturbation", {}).get("small_parameters", []).
type small_parameters:: list[str], optional
returns:: Per-order SolverResults + the combined trajectory + the validity diagnostics.
rtype:: PerturbativeResult
raises ValueError:: If order exceeds the maximum order_in_eps carried by the spec.
raises NotImplementedError:: If order >= 2 — Pass 2 against q¹(t) is not yet implemented. See #273.
raises RuntimeError:: If a Pass-n solver call returns success=False. Failing loud rather than silently producing a degraded result.

Parameters:

y0 (NDArray[np.float64])
grid (GridInfo)
t_span (tuple[float, float])
order (int)
parameters (dict[str, float] | None)
num_snapshots (int)
small_parameters (list[str] | None)

Return type:

PerturbativeResult

Parameters:: spec (EquationSystem)

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 factorization.

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, implicit_dynamical_fields=frozenset({}))[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, ...])
implicit_dynamical_fields (frozenset[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, ...]#

implicit_dynamical_fields: frozenset[str] = frozenset({})#

Order-0 fields in the implicit-dynamical sector (GH #457).

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.

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.

The check is on the summed self-mass of each equation, not on individual terms: several identity terms may act on the same field, and only their sum is the physical mass.

Returns a list of warning strings for tachyonic diagnostics.

Parameters:

coeff_eval (CoefficientEvaluator)
spec (EquationSystem)

Return type: