zoomy_core.model.operations

`zoomy_core.model.operations`#

Composable projection operations on INS equations (Expression, SigmaTransform, Recombine, …).

The user drives every step; there are no hidden assumptions.

Module Contents#

Classes#

`StateSpace`	Shared symbolic state: coordinates, velocity fields, pressure, stress tensor, bathymetry, and physical parameters.
`SymbolicBase`	Base for all symbolic objects. Provides name and notebook display.
`RelationMixin`	Orientation surface shared by `Expression` and `Equation`.
`Expression`	Symbolic expression for PDE terms.
`DepthIntegralResult`	Result of depth-integrating a term over [lower, upper].
`IBPResult`	Structured result from integration by parts.
`Operation`	A structural transformation on a tree node.
`Multiply`	Multiply an Expression (or tree of Expressions) by a factor.
`DepthIntegrate`	Depth-integrate all equations over [b, b+h] w.r.t. z.
`ApplyKinematicBCs`	Apply kinematic BCs globally to combined boundary terms.
`EvaluateIntegrals`	Run `sympy.integrate` on every `Integral` node in the expression.
`Expand`	Substitute a vertical field with its polynomial ansatz expansion in unevaluated Sum form.
`ProjectBasisIntegrals`	Resolve every `Integral(..., (ζ, a, b))` via basis-aware cache lookup.
`SimplifyIntegrals`	Evaluate integrals with constant integrand and remove zero integrals.
`Integrate`	Integrate each additive term of an equation w.r.t. `var` over [lower, upper].
`IntegralTransform`	Affine change of variable applied to every `Integral` in an expression. Maps each integration interval to a reference interval (default: the unit interval `[0, 1]`) via the affine map:
`IsolateBasisIntegrand`	Rewrite each `Integral(integrand, (var, a, b))` so the integrand contains only factors that depend on `var`.
`MapBasisToReference`	Rewrite `phi_k(arg) → phi_k((arg − b) / h)` in every leaf.
`PartialIntegrate`	Rewrite a matching `Integral(f, (var, lo, hi))` to a running integral with variable upper bound: `Integral(f, (var, lo, upper))`.
`AffineProjection`	Transform vertical coordinate: z = ζ·(upper−lower) + lower, dz = (upper−lower)·dζ.
`SigmaTransform`	Change of variables `(t, x, [y,] z) -> (t, x, [y,] ζ)` with `ζ = (z - b) / h` on a PDE leaf, rewriting derivatives by the chain rule.
`Recombine`	Fold expanded expressions back into compact/conservative shapes, e.g. `h·∂_t α + α·∂_t h` -> `∂_t(h α)`.
`ProductRule`	Product rule — auto-routes per term based on outermost shape.
`Newtonian`	Newtonian fluid: tau_ij = mu * (du_i/dx_j + du_j/dx_i), mu = rho * nu.
`Inviscid`	Inviscid fluid: all tau_ij = 0.
`materials`	Material model library. Usage: materials.newtonian(state)
`LayerMeanClosure`	Specialise a σ-resolved system to L piecewise-constant layers.
`KinematicBC`	Kinematic boundary condition at a moving interface.
`NoTangentialBoundaryStress`	Zero tangential normal stress at both surface and bottom.
`HydrostaticPressure`	p = p_atm + rho * g * (eta - z)
`StressFreeSurface`	Stress-free surface: τ·n\|_{z=η} = 0.
`ZeroAtmosphericPressure`	p_atm = 0 (no atmospheric pressure).
`assumptions`	Assumptions library.
`Basis`	Vertical polynomial basis for the SME expansion `u = Σ α_k φ_k(ζ)`.
`Symmetrize`	`a → a/2 + rule(a/2)`.
`ResolveDummy`	Replace every `dummy(arg)` call across every equation in a model.

Functions#

`is_bracket_body`	`True` iff `integrand` is a Galerkin BRACKET body — a pure function of the integration variable `var`, depending on NOTHING else (integer mode indices `i`/`k` excepted, since they parametrise the bracket).
`is_bracket`	`True` iff `integral` is a Galerkin BRACKET: an `sp.Integral` over a CONSTANT interval whose body satisfies `is_bracket_body` — hence a pure NUMBER, evaluable by basis quadrature. A running `∫_0^ζ` (variable upper limit) is therefore NOT a bracket; it is a ζ-function that may sit inside one.
`register_integration_rule`	Register a Zoomy integration rule.
`integrate_cache_stats`	Return `(hits, misses)` for the cached-integrate helper.
`integrate_cache_clear`	Drop all cached integrals and reset stats (useful between runs).
`KinematicBCBottom`	`w\|_{z=b} = ∂_t b + u_b ∂_x b [+ v_b ∂_y b]` (physical z).
`KinematicBCSurface`	`w\|_{z=η} = ∂_t η + u_s ∂_x η [+ v_s ∂_y η]` (physical z).
`integrate_by_parts`	Standalone IBP: integral d(f)/dvar * g dvar = [fg]_a^b - integral f dg/dvar dvar Returns IBPResult(integrate, boundary_upper, boundary_lower).
`gauss_legendre_integrate`

Data#

InterfaceKBC

API#

class zoomy_core.model.operations.StateSpace(dimension=2)#

Shared symbolic state: coordinates, velocity fields, pressure, stress tensor, bathymetry, and physical parameters.

dimension is the physical space dimension: 2 = xz plane (1D horizontal shallow water) 3 = xyz space (2D horizontal shallow water)

The vertical coordinate is always z. Horizontal coordinates are x (and y if dimension=3). The horizontal_dim property gives the number of horizontal directions (1 or 2).

Initialization

property horizontal_dim#

property has_y#

class zoomy_core.model.operations.SymbolicBase(name='')#

Base for all symbolic objects. Provides name and notebook display.

Initialization

class zoomy_core.model.operations.RelationMixin#

Orientation surface shared by Expression and Equation.

An expression/equation OPTIONALLY carries a {lhs: rhs} relation in _as_relation — set either by Expression.solve_for (which isolates a symbol) or by constructing from sp.Eq(lhs, rhs). When present it acts as a substitution under .apply (the dispatcher reads _as_relation) and exposes a dict-like .lhs / .rhs / .subs_map surface. So a SOLVED equation is the relation — there is no separate Relation / Substitution class.

When NOT oriented, the object is a plain residual expr == 0: .lhs is the residual, .rhs is 0, .subs_map is empty.

property oriented: bool#: True when this carries a {lhs: rhs} relation.

property subs_map: dict#: The {lhs: rhs} substitution rules (empty dict if not oriented).

property lhs#: Left side — the solved-for symbol when oriented, else the residual.

property rhs#: Right side — the solution when oriented, else 0.

class zoomy_core.model.operations.Expression(expr, name='', term_groups=None, solver_groups=None, term_tags=None, tag_order=None)#

Bases: zoomy_core.model.operations.RelationMixin, zoomy_core.model.operations.SymbolicBase

Symbolic expression for PDE terms.

Supports:

Term access: .terms, [i], len(), iteration
Projection: .project(test, var, domain, …)
Integration by parts: .ibp(var, test_weight, domain) -> IBPResult
Apply conditions: .apply({old: new}, …) or .apply(relation)
SymPy: .subs(), .simplify(), .expand(), .doit()
Notebook: repr_latex

Initialization

property terms#

apply_to_term(index, *operations)#

Apply operations to a specific term and return the full expression.

Usage:

    xmom.apply_to_term(5, ProductRule())
    xmom.apply_to_term(5, {old: new})

.. note:: Public callers should prefer xmom[5].apply(...) — Expression.__getitem__ returns a live accessor whose .apply mutates the parent in place. This method is kept as the internal split-apply-merge primitive that the accessor uses; library code may continue to call it directly.

keep_groups(*group_names)#: Return a new Expression keeping only the specified tags’ terms.

drop_groups(*group_names)#: Return a new Expression dropping the specified tags’ terms.

tag(**named_groups)#

Attach physical tags to terms for display.

Each name=value kwarg is split into additive terms and each term is recorded as belonging to name in the per-term tag map.

Tags follow the terms through .apply / .subs / .simplify: when a transformation replaces a tagged term, the new term(s) inherit the tag. Terms that cancel (via simplify) silently disappear from the view — no parallel-sum math, no cross-tag contamination.

Usage:

    xmom = xmom.tag(pressure=grad_p, stress=div_tau)

property tags#

Dict of current-expr terms grouped by physical tag.

Computed on demand by walking Add.make_args(self.expr) and assigning each term to the tag recorded in _term_tags. Terms not in the main expression (cancelled by simplify, etc.) do not appear — the view stays in sync with the single source of truth.

property untagged#: Sum of current-expr terms that have no physical tag.

solver_tag(**named_groups)#

Attach/replace solver tag groups for operator routing.

Tag names are normalized via the solver-tag alias catalog, so .solver_tag(convection=...) and .solver_tag(flux=...) both write to the canonical flux slot.

Solver tags are dropped per-tag whenever apply() / simplify() changes the tagged sub-expression (structural equality check).

Usage:

    xmom = xmom.solver_tag(flux=F_x, source=g * h * dbdx)

property solver_tags#: Dict of solver tags (canonical names), or empty dict if none set.

get_solver_tag(name)#

Return the sp.Expr for a solver tag, or None if not set.

name can be any alias; it is normalized to canonical form.

untagged_remainder()#

Return self.expr - sum(solver_tags.values()), simplified.

Zero means every term has been routed to a solver tag; a non-zero remainder is what drives the model’s untagged_policy (warn/error).

auto_solver_tag(*, state_vars, time_var, coords, parameters=())#

Algorithmically assign solver tags by structural matching.

Each additive term of self.expr is classified:

coeff * Derivative(q, t) with q a state variable → time_derivative.
coeff * Derivative(F, x_i) where coeff is state-variable- free → flux (conservative form).
coeff * Derivative(q_j, x_i) with q_j a state variable and coeff containing state vars → nonconservative_flux.
coeff * Derivative(p, x_i) where p is a declared parameter and coeff contains state vars → also nonconservative_flux (the term couples a state coefficient to a parameter gradient — e.g. bathymetry-gradient forcing).
No derivative and purely algebraic in state_vars ∪ parameters space → source.
Anything else (unclosed fields like u(t, x, z), unmatched Derivative shapes) stays in the untagged remainder.

Parameters

state_vars (iterable of sympy Functions / Symbols): The evolution variables (α_k(t, x), interface heights, etc.).
time_var (sympy Symbol): The time coordinate (for time_derivative detection).
coords (iterable of sympy Symbols): Spatial coordinates (x, y, …) for flux / NCP.
parameters (iterable, optional): Known external Function/Symbol atoms (bathymetry b(t, x), prescribed fluxes, etc.). These may freely appear in any term without the term being marked “unclosed”, and they may appear as the inner of a spatial Derivative in the NCP tagging.
Returns a new Expression with solver_groups populated.
Physical tags are preserved.

project(test_func, var, domain=(0, 1), weight=S.One, scale=S.One, numerical=False, order=4)#: Galerkin projection: integral(expr * test * weight * scale, (var, a, b)).

ibp(var, test_weight, domain=(0, 1), scale=S.One)#: Integration by parts on the outermost Derivative w.r.t. var. Returns IBPResult(integrate, boundary_upper, boundary_lower).

apply(*conditions, rhs=_UNSET)#

Apply substitutions / Relations / Operations — per term.

The transformation is applied to every additive term of self.expr individually. Each output term inherits the tag of its parent term (follow-the-term), so the per-term tag dict stays in sync with the evolving main expression automatically.

Shorthand: expr.apply(lhs, rhs=0) is equivalent to expr.apply({lhs: 0}). Valid only with a single positional argument which is a sympy expression (or Expression).

After Relation/dict-style applies the result is lightly simplified via _simplify_preserve_integrals (linearity-only derivatives, cancellations across the full expr). Terms that cancel vanish from both self.expr and the tag dict.

Operations are invoked via their unified __call__ interface; when the op produces multiple output pieces per input term (DepthIntegrate etc.) all pieces inherit the parent tag.

expand()#: Return a new Expression with sympy expand applied.

copy()#: Return an independent Expression with the same expr, tags, relation.

solve_for(variable)#

Return a new Expression whose _as_relation maps variable to the solution of self.expr == 0.

Pipeline-style counterpart to _NodeProxy.solve_for: does not mutate a system tree. Downstream model.apply(result) picks up the relation and performs the substitution.

depth_integrate(lower, upper, var, method='auto')#

Depth-integrate this expression over [lower, upper] w.r.t. var.

Args

lower, upper: Integration bounds (typically b and b+H, both functions of t and x).
var: The vertical coordinate z.
method: 'auto' detects the derivative direction; 'leibniz' pulls a horizontal derivative outside; 'fundamental_theorem' handles vertical derivatives; 'direct' keeps Integral(expr, (var, lower, upper)).

Returns

DepthIntegralResult with .volume, .boundary_upper, .boundary_lower, or a plain Expression for 'direct' and 'fundamental_theorem'.

subs(*args, **kwargs)#

split_integrals()#: Distribute each Integral(Add(...), lim) into a sum of single-term Integrals, recursing through any Derivative wrapper. Returns a new Expression whose .terms exposes one logical equation term per Integral.

merge_integrals()#: Group sibling Integrals with matching (limits, deriv-wrapper) signature into a single Integral(Σ c·f, lim). Inverse of split_integrals; useful when cross-integrand cancellations need sympy’s Add canonicalisation to see them in one integrand.

simplify()#: Simplify: expand + cancel, preserving Integral and Derivative(Integral) terms.

doit()#

has(*args)#

atoms(*types)#

xreplace(rule)#

property free_symbols#

map(fn)#

Apply fn to each term, reassemble into a single Expression.

fn receives an Expression (single term) and must return either an Expression or a DepthIntegralResult. DepthIntegralResults are assembled (volume + boundaries) before summing.

Example

integrated = expr.map(lambda t: t.depth_integrate(b, eta, z))

map_with_bcs(fn, bcs)#

Like map(), but collects boundary terms and applies BCs globally.

This is the correct way to depth-integrate a full equation: boundary terms from ALL terms are combined first, then BCs are applied once (so cross-term cancellations happen properly).

Parameters

fn (callable): Applied to each term. Must return DepthIntegralResult or Expression.
bcs (list of Relation): Kinematic BCs etc. applied to the combined boundary expression.

Returns

Expression: The fully depth-integrated equation with BCs applied.

classify(t=None, x=None, z=None)#

Classify each term by its role in the PDE.

Returns a dict: {role: Expression} where role is one of: ‘temporal’, ‘convective’, ‘diffusive’, ‘source’.

Detection rules:

Has d/dt → temporal
Has d/dx (first-order) of a product → convective flux
Has d²/dz² or d/dz of d/dz → diffusive
Otherwise → source (algebraic)

property temporal#: View: only temporal (d/dt) terms.

property convective#: View: only convective flux (d/dx) terms.

project_onto_basis(basis, level, field_map, z_var, lower=None, upper=None, test_mode=None)#

Project a depth-integrated equation onto a polynomial basis.

Replaces every Integral(f(u,...), (z, b, eta)) by substituting the basis expansion u(z) = sum alpha_k phi_k(zeta) and evaluating the resulting integrals using the SymbolicIntegrator.

If test_mode is an integer, multiplies each integral by the test function phi_{test_mode}(zeta) before evaluating (Galerkin projection for a specific mode). If test_mode=None, returns the scalar integral (e.g. for the mass equation where no test function is needed).

Parameters

basis (Basisfunction class (e.g. Legendre_shifted))
level (int)
field_map (dict): Maps the original Function name to a list of SymPy Symbols for the basis coefficients. Example: {‘u’: [alpha_0, alpha_1, alpha_2]}
z_var (Symbol): The vertical coordinate (z) that appears in the integrals.
lower, upper (sympy expressions (optional)): The integration bounds (b, eta). If None, detected from the first Integral found.
test_mode (int or None): If int, project onto test function phi_{test_mode}.

Returns

Expression: With all depth integrals replaced by basis matrix products.

latex(strip_args=False, multiline=False)#

LaTeX representation.

Parameters

strip_args (bool): u(t,x,z) → u. Partial derivatives preserved.
multiline (bool): Render as \begin{aligned} with one group per line (requires term_groups).

describe(header=True, final_equation=True, parameters=False, strip_args=True)#

Composable description of this expression.

Returns a Description that renders as markdown in Jupyter.

Parameters

header (bool): Show expression name + term count.
final_equation (bool): Show the symbolic equation.
parameters (bool): List free symbols.
strip_args (bool): Display u instead of u(t, x, z).

zoomy_core.model.operations.is_bracket_body(integrand, var)#

True iff integrand is a Galerkin BRACKET body — a pure function of the integration variable var, depending on NOTHING else (integer mode indices i/k excepted, since they parametrise the bracket).

This is THE definition of “this ∫ is a bracket” and the single source of truth reused everywhere it is needed:

the LaTeX printer — render ∫_0^1 (…) dζ as ⟨…⟩ only when this holds;
ExtractBrackets — the naming gate (a t/x-bearing body is not yet a bracket);
the bracket-quadrature step that turns ⟨…⟩ into NUMBERS — pick the integrals to evaluate via the basis instead of sympy.

A factor that carries any non-index continuous symbol (a_i(t,x), h, τ̃, ω, a parameter g) disqualifies it.

zoomy_core.model.operations.is_bracket(integral)#: True iff integral is a Galerkin BRACKET: an sp.Integral over a CONSTANT interval whose body satisfies is_bracket_body — hence a pure NUMBER, evaluable by basis quadrature. A running ∫_0^ζ (variable upper limit) is therefore NOT a bracket; it is a ζ-function that may sit inside one.

class zoomy_core.model.operations.DepthIntegralResult(volume, boundary_upper, boundary_lower)#

Result of depth-integrating a term over [lower, upper].

Attributes

volume: the volume integral (Expression)
boundary_upper: boundary term at z=upper (Expression)
boundary_lower: boundary term at z=lower (Expression)

The full integral = volume + boundary_upper - boundary_lower. Kinematic BCs can be applied to the boundary terms via .apply_bcs().

Initialization

apply_bcs(bc_lower=None, bc_upper=None)#

assemble()#

Combine all terms: volume + boundary_upper + boundary_lower.

The sign convention is that each component already carries its correct sign. For Leibniz: upper = -f(eta)*d(eta)/dx, lower = +f(b)*db/dx. For fundamental theorem: upper = +f(eta), lower = -f(b).

class zoomy_core.model.operations.IBPResult(integrate, boundary_upper, boundary_lower)#

Structured result from integration by parts.

Attributes

integrate: the volume integral (Expression)
boundary_upper: boundary term at upper limit (Expression)
boundary_lower: boundary term at lower limit (Expression)

Initialization

apply_bcs(bc_lower=None, bc_upper=None)#

assemble()#

class zoomy_core.model.operations.Operation(name='', description=None, log_level=None)#

Bases: zoomy_core.model.operations.SymbolicBase

A structural transformation on a tree node.

A node is either an Expression (leaf) or a Zstruct of nodes (intermediate). The unified entry point is __call__(node):

Subclasses override _apply_leaf to transform a leaf Expression. The default __call__ dispatches over the tree, recursing into intermediate Zstructs.
Rank-changing Operations (Multiply with outer=True on a leaf, or contraction on a Zstruct) override __call__ directly and may return a node of different shape than they received.

Every Operation carries name + description so System.apply can record what was done in the derivation history.

log_level (1–5) is the operation’s VERBOSITY in the derivation graph: 1 = major (always reported — the default for every op), rising to 5 = minor bookkeeping (e.g. Simplify canonicalisation). A history view can then filter, e.g. describe(..., max_log_level=1) shows only the major steps. Set it per class (log_level = 5) or per instance (Op(..., log_level=3)).

Initialization

log_level#: 1

rank_changes_leaf#: False

whole_leaf_op#: False

system_level#: False

single_term_only#: False

class zoomy_core.model.operations.Multiply(factor, outer=False, name=None, description=None)#

Bases: zoomy_core.model.operations.Operation

Multiply an Expression (or tree of Expressions) by a factor.

Default behaviour (outer=False) — rank-preserving: Scales every term of a leaf Expression by factor. Tags follow the terms through the per-term dispatch in Expression.apply, so model.apply(Multiply(1/rho)) scales the whole system and keeps physical tags intact.

outer=True — rank-changing (leaf → Zstruct): factor must be a Zstruct of scalar factors (e.g. test functions phi_0, phi_1, ...). The leaf Expression is promoted to a Zstruct with one child Expression per factor, keyed test_<index>. Each child is leaf.expr * factor_k.

Per-term tag inheritance: each parent tag is re-keyed through the multiplication (same rule as Expression.apply), so physical tags follow into every test_<k> child.

Usage:

    model.apply(Multiply(1 / state.rho))                 # scalar scaling
    phis = Zstruct(**{f"phi_{k}": ... for k in range(N)})
    model.momentum.x.apply(Multiply(phis, outer=True))   # Galerkin test

Initialization

property rank_changes_leaf#

class zoomy_core.model.operations.DepthIntegrate(state)#

Bases: zoomy_core.model.operations.Operation

Depth-integrate all equations over [b, b+h] w.r.t. z.

Applies Leibniz rule and fundamental theorem term-by-term. Boundary values (w at z=b, u at z=b+h, etc.) remain as Subs objects. Apply ApplyKinematicBCs to see the cancellations.

Initialization

class zoomy_core.model.operations.ApplyKinematicBCs(state)#

Bases: zoomy_core.model.operations.Operation

Apply kinematic BCs globally to combined boundary terms.

Evaluates Subs boundary values, applies kinematic BCs at surface and bottom, and simplifies. The Leibniz boundary u-terms cancel with the fundamental theorem w-terms.

Must be applied immediately after DepthIntegrate.

Initialization

zoomy_core.model.operations.register_integration_rule(rule)#

Register a Zoomy integration rule.

A rule is (integrand, limits) -> Optional[sp.Expr]:

limits is the full (var, lower, upper) triple passed to _cached_integrate (matching sympy.integrate’s call shape).
Return the definite integral if the pattern matches — the caller does no further transformation.
Return None if the rule doesn’t apply. _cached_integrate tries rules in registration order; the first non-None wins.

zoomy_core.model.operations.integrate_cache_stats()#: Return (hits, misses) for the cached-integrate helper.

zoomy_core.model.operations.integrate_cache_clear()#: Drop all cached integrals and reset stats (useful between runs).

class zoomy_core.model.operations.EvaluateIntegrals(state=None, name='evaluate_integrals', description=None)#

Bases: zoomy_core.model.operations.Operation

Run sympy.integrate on every Integral node in the expression.

After AffineProjection + basis expansion the volume integrals look like Integral(polynomial_in_zeta, (zeta, 0, 1)). This op collapses each such integral to a scalar.

Walks Derivative(Integral(...), x) too: the inner Integral is evaluated first, then the outer Derivative survives symbolically.

Uses _cached_integrate so repeated integrands across leaves (every test_k clone after Multiply(basis.phi, outer=True)) cost one sympy.integrate call each, not one-per-leaf.

Leaves unevaluatable Integrals untouched (e.g. non-polynomial ones or ones whose result would need a closed form sympy can’t find).

Initialization

whole_leaf_op#: True

class zoomy_core.model.operations.Expand(field, basis, coefficients, state, name='expand', description=None)#

Bases: zoomy_core.model.operations.Operation

Substitute a vertical field with its polynomial ansatz expansion in unevaluated Sum form.

Replaces every call field(t, x, [y,] arg_z) with

`Sum(coeff_fn(k) · basis.phi_fn(k, (arg_z − b)/h), (k, 0, L))`

— a single sp.Sum atom rather than the unrolled Add. Two placeholders make this work as one symbolic chunk:

basis.phi_fn is the basis’s opaque 2-arg Function (carries _basis = basis); concrete polynomial substitution is deferred until EvaluateIntegrals.
coeff_fn is a private sympy Function class named coeff_<basis.symbol> and carrying a _coeff_table class attribute holding the pre-declared coefficient list. It stays opaque through ProductRule, AffineProjection, and Integrate (the Sum carries through as a single mathematical term). EvaluateIntegrals substitutes coeff_fn(k_concrete) → coefficients[k_concrete] after the Sum has been .doit()-ed.

Boundary evaluations of the field (arg_z = b, arg_z = b+h) are handled automatically: (arg_z − b)/h becomes 0 or 1 respectively, the Sum is Σ_k amp_k · phi_fn(k, 0|1), and EvaluateIntegrals’ final cleanup pass resolves these via basis.resolve_atoms once the Sum has been .doit()-ed.

Why Sum-form: at higher levels (L ≥ 1), the unrolled Add form has L+1 terms per expanded field — products of two expanded fields explode quadratically, and ProductRule would have to be applied to each one individually. The single-Sum representation keeps the equation compact and makes ProductRule a single apply per expansion, which is the way the math reads on paper.

Parameters

field (sympy Function call): The vertical field, e.g. state.u (which is u(t, x, z)).
basis (Basisfunction): The basis instance. basis.phi_fn is the 2-arg opaque Function used inside the Sum.
coefficients (sequence of sympy Function calls): Pre-declared horizontal coefficient functions (U_0(t, x), U_1(t, x), …) — one per basis level. Length must equal basis.level + 1. Substituted into the Sum by EvaluateIntegrals after Sum.doit() unrolls the symbolic index to concrete integers.
state (StateSpace): The state space — used to read state.b and state.h for the affine ζ-map.

Initialization

class zoomy_core.model.operations.ProjectBasisIntegrals(basis_cache, name='project_basis_integrals', description=None)#

Bases: zoomy_core.model.operations.Operation

Resolve every Integral(..., (ζ, a, b)) via basis-aware cache lookup.

Replaces the old “expand polynomial basis + run sympy.integrate” chain. Every integral whose integrand contains opaque phi_k(arg) nodes is mapped to a BasisIntegralCache query, with ζ-independent factors pulled out first so the kernel the cache sees is a clean basis-only product.

Steps, applied per leaf expression:

Canonicalize Subs — rewrite Subs(Derivative(phi_k(ξ), ξ), ξ, arg) (the sympy chain-rule leftover) back to Derivative(phi_k(arg), arg).
Distribute Adds — Integral(Add(t1, t2), …) → Σ Integral(ti, …) so each term can be factored independently.
Factor constants out + cache query — for each Integral(kernel, (var, a, b)):
- split the Mul integrand into var-independent factor C and var-dependent kernel K;
- if K contains any phi_k, hand K to the cache and multiply its result by C;
- otherwise leave the Integral alone.

The cache returns a sympy scalar / polynomial; nested running integrals inside K get evaluated via the cache too (their polynomial results are substituted into the outer integrand on the next pass, and the fixpoint loop below handles this).

No sympy.integrate is ever called on an opaque phi_k: the cache concretizes phi_k against the basis polynomial first, then does polynomial integration only.

Initialization

class zoomy_core.model.operations.SimplifyIntegrals(state)#

Bases: zoomy_core.model.operations.Operation

Evaluate integrals with constant integrand and remove zero integrals.

∫ 0 dz → 0
∫ c dz → c·h (if integrand has no z-dependence)
∂/∂x ∫ 0 dz → 0

Initialization

class zoomy_core.model.operations.Integrate(var, lower, upper, method='auto')#

Bases: zoomy_core.model.operations.Operation

Integrate each additive term of an equation w.r.t. var over [lower, upper].

Per-term strategy is driven by method:

"auto" (default) — detect the outermost derivative direction per term and pick Leibniz / fundamental theorem / direct accordingly.
"leibniz" — force the Leibniz rule (for terms containing ∂/∂x_i where x_i != var): ∫ ∂_x f dvar = ∂_x [∫ f dvar] − f|_upper · ∂_x upper + f|_lower · ∂_x lower.
"fundamental_theorem" — for terms containing ∂/∂var as the outermost derivative: ∫ ∂_var f dvar = f|_upper − f|_lower.
"direct" — keep as an unevaluated Integral(expr, (var, lower, upper)).
"analytical" — run sympy.integrate on the whole expression (what the original small Integrate did; useful for z-momentum analytic integration when the integrand has no horizontal derivatives).

Boundary terms produced by Leibniz / fundamental-theorem are kept as Subs(f, var, bound) expressions — no kinematic BCs are applied here. Resolve them later via an explicit .apply({Subs-pattern: value}) or a convenience Relation (ApplyKinematicBCs, StressFreeSurface, …).

Partial integration (upper == var) is supported by the "direct" and "analytical" methods — sympy will produce a running integral.

Usage:

    zmom.apply(Integrate(state.z, state.b, state.z))             # partial
    zmom.apply(Integrate(state.z, state.b, state.eta))           # full depth
    zmom.apply(Integrate(state.z, state.b, state.b + state.h,
                         method="analytical"))                    # sympy.integrate

Initialization

whole_leaf_op#: True

class zoomy_core.model.operations.IntegralTransform(ref_interval=(0, 1), ref_name='\\hat{\\zeta}', name='integral_transform', description=None)#

Bases: zoomy_core.model.operations.Operation

Affine change of variable applied to every Integral in an expression. Maps each integration interval to a reference interval (default: the unit interval [0, 1]) via the affine map:

`var = a + (b - a) · (ref - ref_lo) / (ref_hi - ref_lo)`

with Jacobian |dφ/dref| = (b - a) / (ref_hi - ref_lo).

Unlike substitution-based approaches (expr.subs(z, ζ·h + b)), this is strictly local to each Integral: only the integration variable is replaced, only inside the integrand of the Integral it’s bound to. Sibling occurrences of var in the expression — boundary Subs terms, the outer z of a running integral after an outer transform — are untouched, because those aren’t bound by this Integral.

Nested Integrals are handled naturally: the walk recurses into integrands first (bottom-up), so the innermost Integral is transformed before its outer host. After the outer transform replaces var throughout its own integrand, any inner Integral’s now-renamed limits (in the old var) get rewritten in terms of the outer’s reference variable — so a running integral Integral(g(ẑ), (ẑ, b, z)) sitting inside Integral(f(z), (z, b, b+h)) correctly transforms to Integral(g(…), (ẑ′, 0, 1)) · (b + h·ζ − b) after the outer z → b + h·ζ step.

Parameters

ref_interval (tuple (ref_lo, ref_hi), default (0, 1)): The target reference interval. Always the unit interval in typical use.
ref_name (str, optional): Base name for the reference-variable Dummies (rendered with a hat by default). Each transformed Integral gets its own Dummy — different Dummies never collide even with the same display name.

Initialization

whole_leaf_op#: True

class zoomy_core.model.operations.IsolateBasisIntegrand(name='isolate_basis_integrand', description=None)#

Bases: zoomy_core.model.operations.Operation

Rewrite each Integral(integrand, (var, a, b)) so the integrand contains only factors that depend on var.

Two passes per Integral:

Distribute Integral(Add(t1, t2, …), …) → Σ Integral(t_i, …), so each term can be analysed independently.
For each (now single-term) Integral, split the multiplicative factors of the integrand into a var-independent coefficient C and a var-dependent kernel K, and rewrite as C · Integral(K, (var, a, b)).

After this, every Integral has the form coeff · Integral(kernel, …) where kernel only contains var (and any opaque basis functions / derivatives evaluated at var). This is the shape a basismatrix-lookup pattern-matches against: ∫_0^1 phi_k(z)·phi_l(z) dz, ∫_0^1 z·phi_k(z)·phi_l(z) dz, etc.

The Op walks Integrals bottom-up so nested running integrals are isolated first. Each Integral’s integration variable var is read directly from its own limit tuple — there’s no ambiguity, and different Integrals can carry different vars (e.g. ζ̂_0, ζ̂_1, …).

Initialization

whole_leaf_op#: True

class zoomy_core.model.operations.MapBasisToReference(b, h, name='map_basis_to_reference', description=None)#

Bases: zoomy_core.model.operations.Operation

Rewrite phi_k(arg) → phi_k((arg − b) / h) in every leaf.

With the phi_k(state.z) convention the basis is opaque on the raw coordinate z. After IntegralTransform substitutes z → ζ·h + b inside an integrand, integrand basis evaluations read phi_k(ζ·h + b) — semantically the physical basis at point ζ·h + b. This Op states the FEM convention that the physical basis is the reference basis composed with the affine map: phi_k(z) := phi_k_ref((z − b)/h). After the rewrite, sympy auto-simplifies arguments — (b − b)/h → 0, ((b + h) − b)/h → 1, ((ζ h + b) − b)/h → ζ — so every phi_k call lands on a single argument in [0, 1]: the canonical pattern a basismatrix lookup recognises.

Also forces the chain rule on any Derivative(f(t, x, arg), v) whose third argument got rewritten to a (b, h)-dependent form by IntegralTransform — sympy’s structural subs keeps a Subs wrapper around such Derivatives instead of distributing, so we explicitly .doit() them so the chain-rule metric terms materialise inside the integrand.

Parameters

b, h (sympy expressions): The bottom topography and column height entering the affine map z = ζ·h + b. Typically state.b and state.h.

Initialization

class zoomy_core.model.operations.PartialIntegrate(var, upper_symbol=None, lo_match=None, hi_match=None, name='partial_integrate')#

Bases: zoomy_core.model.operations.Operation

Rewrite a matching Integral(f, (var, lo, hi)) to a running integral with variable upper bound: Integral(f, (var, lo, upper)).

The result is left symbolic — no evaluation. Downstream, the expression can be promoted to a Qaux variable via DerivedModel.promote_to_qaux and rewritten into a runtime primitive (e.g. column_partial_integrate).

Matching

Every Integral whose integration variable equals var is inspected. Additional matchers lo_match / hi_match restrict the rewrite to Integrals whose lower/upper bounds also match (structural equality); pass None (default) to leave them unrestricted.

The new upper bound is upper_symbol (default: var itself, producing the classic running integral ∫_lo^{var} f dvar).

Usage:

    # depth integral [b, b+h] → running integral [b, z]
    eq.apply(PartialIntegrate(state.z, lo_match=state.b, hi_match=state.eta))

    # unrestricted: rewrite every ∫ dz to ∫_·^z dz
    eq.apply(PartialIntegrate(z))

Initialization

class zoomy_core.model.operations.AffineProjection(state, lower=None, upper=None, zeta=None, name=None, description=None, *, rewrite_basis_args=True)#

Bases: zoomy_core.model.operations.Operation

Transform vertical coordinate: z = ζ·(upper−lower) + lower, dz = (upper−lower)·dζ.

Transforms z-integrals whose bounds match the configured (lower, upper):

    ∫_{lower}^{upper} f(z) dz → (upper − lower) · ∫_0^1 f(ζ·(upper−lower)+lower) dζ

With no lower / upper given, defaults to state.b / state.eta — the full free-surface transform used by the single-layer SME walkthrough. Pass layer-specific bounds (e.g. lower=state.b, upper=z_1) in the multi-layer case so each layer’s integrals get mapped to their own ζ-interval.

Only Integral nodes whose bounds structurally equal the configured pair are rewritten; other integrals (including already-transformed (zeta, 0, 1) ones) are left untouched. That makes it safe to apply multiple AffineProjections back-to-back on a branch with several integration ranges.

Initialization

class zoomy_core.model.operations.SigmaTransform(state, src=None, *, sigma=None, name='sigma_transform', description=None)#

Bases: zoomy_core.model.operations.Operation

Change of variables (t, x, [y,] z) -> (t, x, [y,] ζ) with ζ = (z - b) / h on a PDE leaf, rewriting derivatives by the chain rule.

Only first-order derivatives are supported.

References

Kowalski, Torrilhon 2019, https://doi.org/10.4208/cicp.OA-2017-0263

Initialization

class zoomy_core.model.operations.Recombine(rules=None, vars=(), aliases=None, name=None, description=None)#

Bases: zoomy_core.model.operations.Operation

Fold expanded expressions back into compact/conservative shapes, e.g. h·∂_t α + α·∂_t h -> ∂_t(h α).

rules picks from Recombine.RULES by name (anti_product_rule, apply_aliases, collapse_trivial_derivative; default all) or takes callables fn(expr, **kwargs); vars and aliases are passed on as kwargs.

Initialization

RULES#: None

whole_leaf_op#: True

class zoomy_core.model.operations.ProductRule(variables=None, direction='inverse', fields=None)#

Bases: zoomy_core.model.operations.Operation

Product rule — auto-routes per term based on outermost shape.

Default direction="both" — applied to a leaf, dispatches per additive term:

Bare Derivative(Π f_i, v) outer → forward: ∂_v(Π f_i) → Σ_i (Π_{j≠i} f_j) · ∂_v(f_i).
Mul carrying a first-order Derivative(f, v) factor → inverse: coeff · ∂_v(f) → ∂_v(coeff · f) − ∂_v(coeff).doit() · f (gets us into divergence form for the integration step). When the term carries MORE THAN ONE derivative factor — e.g. a σ-metric coefficient ∂_x(ζh) multiplying the field derivative ∂_x ũ — the rule pulls out the field-bearing outermost first-order derivative and keeps the other derivative factors inside the conserved flux coeff. (The old rule bailed whenever a second derivative factor was present, so the conservative fold no-op’d on exactly the raw-field σ-metric terms the pre-ansatz Galerkin projection needs.)
Anything else → unchanged (no-op).

Use direction="inverse" / direction="forward" to force one direction.

variables=None (default) — act on derivatives w.r.t. any Symbol named t, x, y, or z. Pass variables=[state.t, state.x, ...] for an explicit whitelist by Symbol identity.

fields=None (default) — when a term carries several allowed first-order derivatives, the field-bearing one (the derivative whose argument contains a user field application) is pulled out; pass fields=[u, w, ...] (heads or applications) to disambiguate explicitly against geometry fields (h, b).

Cross-term combinations (h·∂_t α + α·∂_t h → ∂_t(h·α)) are not the operation’s job — apply ProductRule once on the leaf and sympy’s Add canonicalisation in the next simplify pulls sibling terms together.

Initialization

single_term_only#: False

class zoomy_core.model.operations.Newtonian(state: zoomy_core.model.operations.StateSpace, nu=None)#

Bases: zoomy_core.model.operations.Expression

Newtonian fluid: tau_ij = mu * (du_i/dx_j + du_j/dx_i), mu = rho * nu.

Initialization

class zoomy_core.model.operations.Inviscid(state: zoomy_core.model.operations.StateSpace)#

Bases: zoomy_core.model.operations.Expression

Inviscid fluid: all tau_ij = 0.

Initialization

class zoomy_core.model.operations.materials#

Material model library. Usage: materials.newtonian(state)

newtonian#: None

inviscid#: None

class zoomy_core.model.operations.LayerMeanClosure(state, sigma, u_field, u_layer, *, weights=None, max_u_power=3, name='layer_mean_closure', description=None)#

Bases: zoomy_core.model.operations.Operation

Specialise a σ-resolved system to L piecewise-constant layers.

For each leaf that depends on the σ-coordinate, project onto each layer α = 1..L: integrate symbolically over the layer’s σ-bounds [σ_{α-1/2}, σ_{α+1/2}] while keeping u(t, x, σ) opaque, then resolve the surviving bulk integrals using the layer-mean closure ∫_{layer α} u(σ) dσ = l_α · u_α(t, x). Interior IBP rules handle f(σ) · u^k · ∂_σ u terms via the identity u^k · ∂_σ u = (1/(k+1)) ∂_σ u^{k+1}.

The boundary atoms u(t, x, σ_{α±1/2}) and w(t, x, σ_{α±1/2}) are deliberately left opaque in the resulting per-layer leaves. That is the user’s modeling choice — apply KinematicBC with mass_flux=G_{α+1/2}/ρ at each interface to resolve w, and apply an interface-velocity choice (donor-cell / central / sign-upwind) for the surviving u boundary atoms. AHS26 (9) is the standard donor-cell-by-G choice; other closures are equally valid and now visible as separate substitutions rather than buried in the projection.

Why not substitute a piecewise ansatz before integration:

Piecewise(u_α, σ in layer α) substituted globally would force sympy’s integration to pick a specific value at the layer boundary (dominant-branch convention), silently committing the user to a specific (asymmetric) interface velocity.
Σ u_α · 𝟙_{[σ_{α-1/2}, σ_{α+1/2}]}(σ) (Heaviside form) trips sympy’s Meijer-G fallback when products of Heavisides meet ∂_σ in the integrand.

Keeping u opaque and integrating term-by-term avoids both traps and surfaces the modeling choice explicitly.

Parameters

state (StateSpace): The state-space the operation lives on.
sigma (sympy.Symbol): σ-coordinate of the σ-mapped system (typically state.zeta_ref).
u_field (sympy.Function call): The opaque u(t, x, σ) atom whose layer-mean closure we use.
u_layer (sequence of sympy Function calls): Per-layer mean velocities u_α(t, x), one per layer α.
weights (sequence of sympy expressions, optional): Layer weights l_α summing to 1. Default: uniform l_α = 1/L.
max_u_power (int, optional): Maximum monomial power of u_field to recognise in IBP patterns f(σ) · u^k · ∂_σ u (default 3 — covers ∂_σ u, u·∂_σ u (momentum convection), u²·∂_σ u (cubic shallow-water) and u³·∂_σ u).
name, description: Standard operation labels for derivation history.

Initialization

rank_changes_leaf#: True

class zoomy_core.model.operations.KinematicBC(state=None, interface=None, src=None, *, at=None, w=None, u=None, v=None, rho=None, mass_flux=None, name=None, description=None)#

Bases: zoomy_core.model.operations.Expression

Kinematic boundary condition at a moving interface.

One explicit signature covering every kinematic-BC pattern in the codebase — solid bed, free surface, σ-domain boundary in a σ-transformed system, and every internal interface in a multi-layer derivation.

The interface interface(t, x[, y]) moves with the fluid:

.. math::

w\big|_\text{at}
  \;=\; \partial_t \,\text{interface}
  + u\big|_\text{at}\, \partial_x \,\text{interface}
  \;[\,+\, v\big|_\text{at}\, \partial_y \,\text{interface}\,]
  \;+\; \frac{\text{mass\_flux}}{\rho}\,.

Parameters

state (StateSpace)
interface (sympy expression): The moving surface, a function of (t, x[, y]). Typical choices: state.b (bottom), state.eta (free surface), a freshly-minted Function("z_k")(t, x[, y]) for an internal multi-layer interface.
at (sympy expression, optional): Vertical-coordinate value at which the velocity fields are sampled. Defaults to interface — the natural physical-z choice “evaluate at the interface itself”. Use an explicit at when the system has been mapped to a reference vertical (e.g. at=sp.S.Zero after SigmaTransform for the bottom; at=sp.S.One for the surface; per-layer σ values in ML). The substitution is purely structural — replaces state.z by at in the velocity arg lists — so it works in physical (t, x[, y], z) or post-σ (t, x[, y], ζ) form regardless.
mass_flux (sympy expression, optional): Mass flux m per unit area per unit time crossing the interface. Adds mass_flux / ρ to the RHS. Default None for impermeable interfaces (bottom / free surface). Use a fresh Function for internal layer interfaces whose flux is an unknown of the resulting multi-layer system.
name, description (str, optional): History labels.

Examples

Single-layer in physical z (default at=interface):

    sys.apply(KinematicBC(state, state.b))
    sys.apply(KinematicBC(state, state.eta))

Single-layer after SigmaTransform — evaluate at the σ endpoints rather than at the (now opaque) physical interface:

    sys.apply(KinematicBC(state, state.b,   at=sp.S.Zero))
    sys.apply(KinematicBC(state, state.eta, at=sp.S.One))

Multi-layer (interfaces z_k; internals carry a mass flux m_k; each layer’s state has its own σ but the formula is the same — only interface, at, and mass_flux change per call):

    for k, z_k in enumerate(interfaces):
        m_k = mass_fluxes[k] if 0 < k < N else None
        sys_layer_above.apply(KinematicBC(state_above, z_k, at=0, mass_flux=m_k))
        sys_layer_below.apply(KinematicBC(state_below, z_k, at=1, mass_flux=m_k))

Initialization

zoomy_core.model.operations.KinematicBCBottom(state: zoomy_core.model.operations.StateSpace)#: w|_{z=b} = ∂_t b + u_b ∂_x b [+ v_b ∂_y b] (physical z).

zoomy_core.model.operations.KinematicBCSurface(state: zoomy_core.model.operations.StateSpace)#: w|_{z=η} = ∂_t η + u_s ∂_x η [+ v_s ∂_y η] (physical z).

zoomy_core.model.operations.InterfaceKBC#: None

class zoomy_core.model.operations.NoTangentialBoundaryStress(state: zoomy_core.model.operations.StateSpace)#

Bases: zoomy_core.model.operations.Expression

Zero tangential normal stress at both surface and bottom.

Closes the Leibniz boundary terms tau_xx(z=b), tau_xx(z=eta) (and y-components in 3D) that remain after depth-integrating ∂_x tau_xx. Uses a Newtonian-style “no flow variation” closure at the vertical boundaries.

Initialization

class zoomy_core.model.operations.HydrostaticPressure(state: zoomy_core.model.operations.StateSpace)#

Bases: zoomy_core.model.operations.Expression

p = p_atm + rho * g * (eta - z)

Initialization

class zoomy_core.model.operations.StressFreeSurface(state: zoomy_core.model.operations.StateSpace)#

Bases: zoomy_core.model.operations.Expression

Stress-free surface: τ·n|_{z=η} = 0.

For a free surface with normal n ≈ (0,0,1), this gives: τ_xz|{z=η} = 0, τ_zz|{z=η} = 0 (and τ_yz|_{z=η} = 0 in 3D)

Initialization

class zoomy_core.model.operations.ZeroAtmosphericPressure(state: zoomy_core.model.operations.StateSpace)#

Bases: zoomy_core.model.operations.Expression

p_atm = 0 (no atmospheric pressure).

Initialization

class zoomy_core.model.operations.assumptions#

Assumptions library.

kinematic_bc_bottom#: None

kinematic_bc_surface#: None

hydrostatic_pressure#: None

stress_free_surface#: None

zero_atmospheric_pressure#: None

class zoomy_core.model.operations.Basis(state: zoomy_core.model.operations.StateSpace, basisfunction_cls, level=0, alpha_name='alpha', **basis_kwargs)#

Vertical polynomial basis for the SME expansion u = Σ α_k φ_k(ζ).

Wraps a Basisfunction class (Legendre_shifted etc.) and exposes symbolic pieces for the apply-based walkthrough:

basis.phi — Zstruct(phi_0, …, phi_N) of sympy expressions in state.zeta_ref (the reference-element symbol on [0, 1]). Drop-in factor for Multiply(basis.phi, outer=True) after AffineProjection has collapsed the opaque state.zeta to state.zeta_ref.
basis.alpha — Zstruct(alpha_0, …, alpha_N) of Function coefficients α_k(t, x, [y]). Using Functions (not bare Symbols) keeps downstream partial derivatives honest and makes linear-stability analysis straightforward.
basis.expand(field) — single substitution dict covering the volume evaluation (in ζ-transformed form) and both boundary evaluations at z=b / z=η. Feed it to model.apply(...).

Usage:

    from zoomy_core.model.derivation.basisfunctions import Legendre_shifted
    basis = Basis(state, Legendre_shifted, level=2)
    model.momentum.x.apply(Multiply(basis.phi, outer=True))
    model.apply(AffineProjection(state))
    model.apply(basis.expand(state.u))

The volume substitution uses the ζ-transformed key field.subs(z, ζ·h + b) — apply it after AffineProjection so it matches what’s under the integrals.

Initialization

Wrap a Basisfunction subclass for the apply-based pipeline.

basis_kwargs are forwarded to the basisfunction_cls constructor — use them for bases that need extra structure beyond the level index (e.g. PiecewiseConstant(interfaces=[state.b, z_1, state.eta])). When basis_kwargs is non-empty, level is ignored and re-read from the instantiated basisfunction.

expand(field)#

Substitution dict for every bulk / boundary evaluation of field.

Four keys:

field → Σ α_k(t,x) · φ_k((z−b)/h) (pointwise form — matches field(t,x,z) in inner running integrals produced by the w-closure pipeline).
field.subs(z, ζ·h + b) → Σ α_k(t,x) · φ_k(ζ) (ζ-transformed form occurring inside outer ζ-integrals).
field.subs(z, b) → Σ α_k · φ_k(0) (bottom eval).
field.subs(z, b + h) → Σ α_k · φ_k(1) (surface eval).

Including the pointwise key is free — it only fires where the expression structurally carries field(t, x, z) — and makes the w-closure’s inner ∫_b^z ∂_x u dz' close cleanly once u has been expanded.

layer_expand(field, layer_idx, zeta_transformed=False)#

Per-layer closure substitution for a LayeredBasis: inside layer layer_idx the field is closed against that layer’s coefficients.

Args

field: Variable being closed (e.g. state.u).
layer_idx: Index into self._bf.interfaces specifying the layer.
zeta_transformed: False: the bulk key is field itself; True (after AffineProjection): field.subs(z, ζ·h_i + z_i).

Returns

Substitution dict with three entries (volume, lower and upper interface of the current layer).

Raises

ValueError: If the wrapped basis function is not a LayeredBasis.

layer_phi(layer_idx)#

Zstruct of inner test functions in the shared ζ symbol.

Equivalent to basis.phi but for one layer of a LayeredBasis. Use after AffineProjection, inside integrands, when you need φ_inner_k(ζ) multiplied into an EvaluateIntegrals target.

Raises ValueError on a non-layered basis.

layer_phi_of_z(layer_idx)#

Zstruct of inner test functions at ζ_i(z) = (z − z_i)/h_i.

Equivalent to basis.phi_of_z but for one layer of a LayeredBasis. Use before depth-integration — Multiply(basis.layer_phi_of_z(i), outer=True) on a layer’s branch produces its Galerkin test equations, one per inner moment, each carrying the rescaled φ_inner_k((z − z_i)/h_i) factor that the Leibniz rule will dismantle.

Raises ValueError on a non-layered basis.

zoomy_core.model.operations.integrate_by_parts(f, g, var, domain=(0, 1))#: Standalone IBP: integral d(f)/dvar * g dvar = [f*g]_a^b - integral f * dg/dvar dvar Returns IBPResult(integrate, boundary_upper, boundary_lower).

zoomy_core.model.operations.gauss_legendre_integrate(expr, var, a, b, order=4)#

class zoomy_core.model.operations.Symmetrize(rule)#

Bases: zoomy_core.model.operations.Operation

a → a/2 + rule(a/2).

When rule = ProductRule() and the term has the self-pair shape c · f · ∂_v(f), the cancellation under Equation.simplify leaves a clean conservative atom c/2 · ∂_v(f²). Used for the gravity self-pair fold in SME / VAM.

Initialization

class zoomy_core.model.operations.ResolveDummy(dummy, value)#

Bases: zoomy_core.model.operations.Operation

Replace every dummy(arg) call across every equation in a model.

Promoted from the legacy notebook Model.resolve_dummy method to a full Operation. The new Model surface is minimal (add_equation / remove_equation / apply / describe) and substitutions go via model.apply(ResolveDummy(dummy, value)).

Parameters

dummy: A sympy Function placeholder (e.g. phi_u_fn, psi_u, omega) that the derivation introduced as an opaque modal-ansatz / weight stand-in.
value:
- sympy expression / scalar → single substitution (replace every dummy(arg) with value);
- callable f(*args) → sp.Expr → call-with-args substitution (dummy(a, b) → value(a, b));
- list / tuple → branch each equation containing the dummy into one sub-equation per list element (outer product on the dummy axis). Each sub-equation is named "{eq_name}_{i}".

Initialization

apply_to_model(model)#

zoomy_core.model.operations

Contents

zoomy_core.model.operations#

Module Contents#

Classes#

Functions#

Data#

API#

`zoomy_core.model.operations`#