zoomy_core.model.models.ml_vam_wd#
ML-VAM-WD — ML-VAM with the wet/dry (discharge-form) divergence
constraint of 86BQYCRP Eq. 19. EXPERIMENTAL sibling of ml_vam.py.
Byte-identical to ml_vam.py except for ONE block — the per-layer
divergence-constraint emission in MLVAMWD.derive_model (search
86BQYCRP Eq. 19). The stock model keeps the σ-mass-flux compound
(1/h)·∂_x(h·q_k) folded, which leaves a 1/h INSIDE the emitted elliptic
operator and makes it SINGULAR as h→0 at a thinning front. Here the row is
distributed and multiplied by l_ℓ·h² so the 1/h cancels analytically and
the pressure block stays invertible (P→0 automatically at wet/dry fronts,
no threshold). Same constraint manifold wherever h>0.
Kept separate so the validated MLVAM is untouched; fold back
once the 2-layer Escalante bump survives the front (cid 182).
Second departure: the LDNH constraint structure (Fernández-Nieto/Parisot/
Penel/Sainte-Marie 2018; Escalante–Morales de Luna 2021 eq. 8/10; the
multilayer version 32K4N2P5 §3.1). ml_vam.py SPENDS the bottom/interface
kinematic condition to eliminate the top ŵ mode; here ŵ_top stays a
STATE (with its own momentum_z row) and the kinematic residual is re-added
as an INDEPENDENT constraint row, paired with the extra pressure unknown that
the N^q = N^w + 1 degree pairing prescribes. Spaces per layer:
u ∈ P_Nu (Nu+1 free modes, q_ℓ_k)
w ∈ P_{Nu+1} (Nu+2 free modes, r_ℓ_k — NO closed top mode)
p ∈ P_{Nu+2} (Nu+2 free modes, P_ℓ_j — top mode closed by the trace)
constraints per layer = mass[1 … Nu+1] + R_kbc = Nu+2
free pressure modes = P_ℓ_0 … P_ℓ_{Nu+1} = Nu+2 ✔ square
MEASURED (STEP-1 rank study, symbolic Jacobian on atomised modal columns, at
(N,Nu) = (1,1), (2,1), (1,2)): the per-layer relation set
{mass[0…Nu+1], R_kbc, R_top} has Nu+4 rows of rank Nu+3 — EXACTLY ONE
dependency, mass[Nu+1] = ∓R_kbc − R_top, i.e. the top mass mode is pure
boundary data. Consequences, both honoured below:
R_top(surface/interface KBC) is REDUNDANT — do not add it;mass[0]is NOT free — it is the only relation supplying∂_t h_ℓ/Gand stays the layer depth equation (measured gain from freeing it: 0);freeing
ŵ_topand addingR_kbcbuys +1 independent constraint — which is precisely the one extra pressure degree ofN^q = N^w + 1.
⚠ HISTORY: an earlier pass widened p to P_{Nu+2} and matched it with the
mode-Nu+2 LAYER-MASS projection. That row is non-zero but DEPENDENT (it
differs from mode Nu+1 by twice the bottom KBC, which was then substituted
away), so the pressure block was rank Nu+1 for Nu+2 unknowns. The
extra pressure degree needs the KINEMATIC row, never another mass mode.
Layers are coupled by KinematicBC(mass_flux=G_ℓ) with the Hörnschemeyer
fraction-multiplier closure (h_ℓ = l_ℓ·h).
The N pressure top-mode closures cascade DOWNWARD through the stack: surface
p_N(1) = 0; every lower layer’s top trace equals the layer-above’s bottom
trace (downward convention — the upward one makes the elliptic block singular,
see the ml_vam thesis notebook). The traces are FULL modal traces including
the closed top modes, so the chain resolves from the surface to the bed.
The N kinematic conditions are NO LONGER closures. w_ℓ(0) = ∂t z_{ℓ-½} + u_ℓ(0)·∂x z_{ℓ-½} + G_{ℓ-½}/ρ (bottom KBC w_1(0) = u_1(0)·∂x b on the
bed) is emitted as the constraint row constraint_kbc_{ℓ}, whose Lagrange
multiplier is the extra per-layer pressure unknown — the bottom/interface
pressure q_b of Escalante 2021 eq. (10).
Per layer the advective σ-mass flux ω̃ is resolved to its CLOSED form
ω̃(ζ) = G_{ℓ-½}/ρ − ∫₀^ζ (∂t h_ℓ + ∂x(h_ℓ ũ_ℓ)) (Escalante’s
eq(3)→(4) step with the interface mass-flux offset) — the definition form
differs by {layer-mass, constraint} combinations, equivalent ON the
constraint manifold but not in the off-manifold predictor.
Interface momentum transfer: u* is the surgical mean-trace swap on the
x-momentum rows (u-traces are genuinely one-sided). The w-rows get NO
swap — the kinematic closures make the w-traces of adjacent layers exactly
continuous, so a w* transfer correction is identically zero.
SystemModel.from_model(MLVAMWD(...)) is the square DAE (state
[b, h, q_ℓk, r_ℓk, P_ℓj]; the per-layer divergence constraints are the
zero-mass-matrix rows); MLVAMWD.chorin_split(dt) returns the structural
predictor / pressure / corrector sub-systems for
ChorinSplitVAMSolver.
Module Contents#
Classes#
ML-VAM in the LDNH constraint structure with the wet/dry (discharge-form) constraint scaling of 86BQYCRP Eq. 19. |
Data#
API#
- zoomy_core.model.models.ml_vam_wd.zeta#
‘Symbol(…)’
- class zoomy_core.model.models.ml_vam_wd.MLVAMWD(init_functions=True, **params)#
Bases:
zoomy_core.model.basemodel.ModelML-VAM in the LDNH constraint structure with the wet/dry (discharge-form) constraint scaling of 86BQYCRP Eq. 19.
Two departures from
MLVAM:the per-layer divergence constraints are emitted multiplied by
l_ℓ·h²(the kinematic row byl_ℓ·h) with the1/hcancelled ANALYTICALLY, so the elliptic block stays invertible ash→0;ŵ_topis a STATE, not a closure.w ∈ P_{Nu+1}keeps allNu+2modes (each with its ownmomentum_zrow) and the bottom/interface kinematic condition returns as the constraintconstraint_kbc_{ℓ}. That buys the +1 independent constraint (measured, seederive_model) which pairs withp ∈ P_{Nu+2}’s extra unknown:Nu+2constraints ↔Nu+2free pressure modes per layer.
EXPERIMENTAL sibling, not a replacement: kept as its own model so the validated
MLVAMis untouched while this is proven. Fold back intoml_vam.pyonly once the 2-layer Escalante bump survives the thinning front (cid 182).Initialization
Declare variables and parameters, run the derivation, and (unless
_finalize_lazy) finalize for the SystemModel.Args
init_functions: register the operator methods as symbolic functions during finalization.
- n_layers#
‘Integer(…)’
- level#
‘Integer(…)’
- dimension#
‘Integer(…)’
- closures#
‘List(…)’
- derive_model()#
- chorin_split(dt=None, *, system_model=None)#
Structural Chorin split (predictor / pressure / corrector).