zoomy_core.model.models.legacy.pde_generator module

zoomy_core.model.models.legacy.pde_generator module#

Automated PDE Generator for Multi-Layer SME/VAM Models.

Derives layer-averaged / moment-projected shallow water equations from the 3D Incompressible Navier-Stokes equations using Galerkin projection with arbitrary vertical basis functions.

Phase 1: Algebraic Formulation & Moment Projection

1.1 — Non-integrated INS base (continuous SymPy functions) 1.2 — Abstract Galerkin projection with weighted test functions 1.3 — Application of kinematic & dynamic boundary conditions

Phase 2: Basis Injection & Delayed Substitution

2.1 — Coordinate mapping to reference element 2.2 — Heaviside windowing for multi-layer ansatz 2.3 — Substitution into abstract projected equations

Phase 3: Custom SymPy Integration Engine

3.1 — AST parsing for Heaviside/DiracDelta 3.2 — Static piecewise collapse per layer 3.3 — DiracDelta extraction via sifting property

class zoomy_core.model.models.legacy.pde_generator.INSBase(dimension=2)#

Bases: object

3D Incompressible Navier-Stokes equations in continuous form.

State fields: u(t,x,y,z), v(t,x,y,z), w(t,x,y,z), p(t,x,y,z) Parameters: rho (density), nu (kinematic viscosity), g (gravity)

The equations are:

Continuity: du/dx + dv/dy + dw/dz = 0 X-momentum: du/dt + d(uu)/dx + d(uv)/dy + d(uw)/dz = -1/rho dp/dx + nu Lap(u) Y-momentum: dv/dt + d(vu)/dx + d(vv)/dy + d(vw)/dz = -1/rho dp/dy + nu Lap(v) Z-momentum: dw/dt + d(wu)/dx + d(wv)/dy + d(ww)/dz = -1/rho dp/dz + nu Lap(w) - g

For shallow water derivation we work in the hydrostatic limit (drop z-momentum PDE, use dp/dz = -rho*g directly).

full_args()#
continuity()#
x_momentum()#
y_momentum()#
hydrostatic_pressure()#
momentum_equations()#
class zoomy_core.model.models.legacy.pde_generator.GalerkinProjection(ins)#

Bases: object

Galerkin projection of 3D INS equations onto vertical basis functions.

Given the PDE L(u,v,w,p) = 0 over z in [b, b+H], the projected equation for test function phi_i(zeta) with weight c(zeta) is:

integral_0^1 L(…) * c(zeta) * phi_i(zeta) * H d_zeta = 0

where zeta = (z - b) / H maps z in [b, b+H] to [0, 1].

This class works with ABSTRACT symbols for:

  • phi_i(zeta): test function (not yet specified)

  • u, v, w, p: continuous fields

  • b, H: bathymetry and total depth

Parameters:

ins (INSBase) –

z_of_zeta()#
project_continuity()#

Integrate continuity equation over depth with test function.

Starting from: du/dx + dv/dy + dw/dz = 0

Multiply by H * c(zeta) * phi_i(zeta) and integrate over zeta in [0,1]:

integral_0^1 [du/dx + dv/dy] * H * c * phi_i d_zeta + integral_0^1 dw/dz * H * c * phi_i d_zeta = 0

For the dw/dz term, use dw/dz = (1/H) dw/d_zeta, so:

integral_0^1 dw/d_zeta * c * phi_i d_zeta

Integration by parts on this last term gives:

[w * c * phi_i]_0^1 - integral_0^1 w * d(c * phi_i)/d_zeta d_zeta

The boundary terms [w * c * phi_i] at zeta=0 and zeta=1 are where kinematic BCs enter (Phase 1.3).

project_momentum(component_index=0)#

Project the horizontal momentum equation for component component_index.

component_index=0 -> x-momentum component_index=1 -> y-momentum (only if dimension==2)

Starting from (x-momentum shown):

du/dt + d(uu)/dx + d(uv)/dy + d(uw)/dz + (1/rho) dp/dx = 0

Multiply by H * c(zeta) * phi_i(zeta), integrate over zeta in [0,1].

The d(uw)/dz term gets integration by parts (same as continuity):

integral dw*u/dz * H * c * phi_i d_zeta = [u*w * c * phi_i]_0^1 - integral u*w * d(c*phi_i)/d_zeta d_zeta

The pressure gradient term becomes:

integral (1/rho) dp/dx * H * c * phi_i d_zeta

Under hydrostatic assumption p(z) = p_s + rho*g*(b+H-z), the pressure gradient introduces the free surface gradient and atmospheric pressure — handled in Phase 1.3.

project_all()#

Return all projected equations (continuity + momentum components).

class zoomy_core.model.models.legacy.pde_generator.BoundaryConditions(projection)#

Bases: object

Physical boundary conditions for the depth-averaged system.

Kinematic BCs (zero relative mass flux):

At z = b: w_b = db/dt + u_b * db/dx + v_b * db/dy At z = b + H: w_s = d(b+H)/dt + u_s * d(b+H)/dx + v_s * d(b+H)/dy

Dynamic BCs:

Surface: tau_s (wind stress, atmospheric pressure gradient) Bottom: tau_b (bed friction)

Hydrostatic pressure:

p(x,y,z,t) = p_atm + rho * g * (b + H - z) => dp/dx = dp_atm/dx + rho*g * d(b+H)/dx - rho*g*dz/dx

= dp_atm/dx + rho*g * d(b+H)/dx (for depth-averaged)

Parameters:

projection (GalerkinProjection) –

kinematic_bottom()#

w at z = b: w_b = db/dt + u_b*db/dx [+ v_b*db/dy]

kinematic_surface()#

w at z = b+H: w_s = d(b+H)/dt + u_s*d(b+H)/dx [+ v_s*d(b+H)/dy]

hydrostatic_pressure_gradient(coord)#
Under hydrostatic assumption:

p = p_atm + rho * g * (b + H - z) dp/dx = dp_atm/dx + rho * g * d(b+H)/dx

Projected:

(1/rho) integral dp/dx * H * c * phi_i d_zeta = integral [1/rho * dp_atm/dx + g * d(b+H)/dx] * H * c * phi_i d_zeta = [1/rho * dp_atm/dx + g * d(b+H)/dx] * H * integral c * phi_i d_zeta

The last step holds because the pressure gradient is z-independent under hydrostatic assumption.

apply_to_continuity(projected_eq)#

Replace abstract w_s, w_b in the projected continuity equation with kinematic BC expressions.

apply_to_momentum(projected_eq, component_index=0)#

Apply BCs to projected momentum equation: 1. Replace w_s, w_b with kinematic BCs in boundary terms 2. Replace pressure gradient with hydrostatic expression 3. Add surface/bottom stress contributions

apply_all(projected_equations)#

Apply BCs to a list of projected equations from GalerkinProjection.project_all().

class zoomy_core.model.models.legacy.pde_generator.ProjectedEquation(volume, boundary_top, boundary_bottom, name='', pressure_gradient=None)#

Bases: object

Container for a single projected (depth-integrated) equation.

volume#

The volume integral terms (integrands integrated over [0,1])

boundary_top#

Boundary contribution at zeta=1 (surface)

boundary_bottom#

Boundary contribution at zeta=0 (bottom)

pressure_gradient#

Hydrostatic pressure gradient (if applicable)

name#

Human-readable identifier

full_equation()#
pprint()#
class zoomy_core.model.models.legacy.pde_generator.PiecewiseIntegrator(var, interfaces)#

Bases: object

Custom integration engine that safely handles Heaviside-windowed expressions.

SymPy’s native Piecewise integration causes AST explosion and hangs for multi-layer expressions. This engine:

  1. Scans the expression tree for Heaviside and DiracDelta terms.

  2. For smooth volume integration: divides the domain into chunks at layer interfaces, statically collapses Heaviside to 0/1 per chunk, then integrates the resulting smooth polynomial.

  3. For DiracDelta terms: extracts coefficients and applies the sifting property f(z_k) / |g'(z_k)| at each interface.

integrate(expr)#

Integrate expr over the full domain [interfaces[0], interfaces[-1]].

Splits expr into smooth (Heaviside) and singular (DiracDelta) parts, handles each with the appropriate strategy, and returns the sum.

class zoomy_core.model.models.legacy.pde_generator.LayeredAnsatz(n_layers, basis, dimension=1)#

Bases: object

Constructs the multi-layer piecewise vertical ansatz using Heaviside windowing and substitutes it into the abstract projected equations.

The ansatz for a horizontal velocity component is:

u(t, x, zeta) = sum_k sum_j u_{k,j}(t,x) * phi_j(zeta_local)

  • [H(zeta - zeta_{k-1}) - H(zeta - zeta_k)]

where:

  • k indexes the layer (k = 0, …, n_layers-1)

  • j indexes the basis function within a layer (j = 0, …, level)

  • zeta_local = (zeta - zeta_{k-1}) / (zeta_k - zeta_{k-1})

  • phi_j is the j-th basis polynomial on [0, 1]

  • H is the Heaviside step function

local_coordinate(k)#

Local coordinate within layer k: zeta_local = (zeta - zeta_k) / delta_k.

layer_thickness(k)#

Thickness of layer k in zeta-space.

window_function(k)#

Heaviside window for layer k: H(zeta - zeta_k) - H(zeta - zeta_{k+1}).

ansatz_in_layer(component, k)#
Expansion within a single layer k for given component:

sum_j u_{k,j} * phi_j(zeta_local)

full_ansatz(component)#
Full piecewise ansatz with Heaviside windowing:

u(zeta) = sum_k [sum_j u_{k,j} phi_j(zeta_local)] * W_k(zeta)

derivative_ansatz(component)#

d/d_zeta of the full ansatz. Produces: - Smooth terms: derivatives of polynomials within each window - DiracDelta terms at interfaces from Heaviside derivatives

substitute_into_integrand(integrand, component_map=None)#

Replace abstract Function calls in an integrand expression with the concrete piecewise ansatz.

Parameters:

  • integrand (Expr) – SymPy expression potentially containing Function(“u”)(t, x, zeta), etc.

  • component_map (dict, optional) – Maps function names to component keys, e.g. {“u”: “u”, “v”: “v”}. Defaults to identity mapping.

integrate_projected(integrand, test_func_index, weight_func=None)#

Full pipeline: substitute ansatz, multiply by test function and weight, then integrate using the PiecewiseIntegrator.

Parameters:

  • integrand (Expr) – Abstract integrand (from Phase 1 projection).

  • test_func_index (int) – Which test function phi_i to use (layer_k, basis_j encoded).

  • weight_func (Expr, optional) – Integration weight c(zeta). Defaults to 1.

get_all_dof_symbols(component=None)#

Return a flat list of all DOF symbols, optionally for one component.

class zoomy_core.model.models.legacy.pde_generator.InterfaceRouter(interfaces, var)#

Bases: object

Separates internal numerical interfaces from external physical boundaries.

Given a list of DiracDelta terms from the PiecewiseIntegrator, this class: 1. Checks the location (root) of each delta 2. If at z=0 or z=1 (physical boundary): discards the symbolic jump and

replaces it with the user-defined physical condition

  1. If at an internal interface (0 < z_k < 1): tags it as an unresolved numerical flux for the Riemann solver

classify_delta_terms(expr)#

Split an expression into three parts: - smooth: terms without DiracDelta - boundary_deltas: delta terms at z=0 or z=1 - internal_deltas: delta terms at internal interfaces

Returns:

ClassifiedTerms with attributes

Return type:

smooth, boundary, internal

apply_physical_bcs(classified, bc_bottom=None, bc_top=None)#

Replace boundary delta terms with user-defined physical conditions.

Parameters:

  • classified (ClassifiedTerms) –

  • bc_bottom (Expr or None) – Physical boundary condition at z=0 (e.g., tau_bottom/rho for diffusion, 0 for advection)

  • bc_top (Expr or None) – Physical boundary condition at z=1

  • applied. (Returns the modified smooth expression with BCs) –

class zoomy_core.model.models.legacy.pde_generator.InterfaceJump(expression, location, jump_type)#

Bases: object

A single interface jump term extracted from the integration.

class zoomy_core.model.models.legacy.pde_generator.ClassifiedTerms(smooth, boundary, internal)#

Bases: object

Container for terms classified by the InterfaceRouter.

property n_internal#
property n_boundary#