build_model(m::AbstractAssembledModel; kwargs...) -> (; model, spectral, meta)

Construct the full-order model and spectral data consumed by a reduction.

This is the single contract every MORFEFerrite physics module implements. The first two fields are exactly what MORFE.parametrise takes:

(; model, spectral) = build_model(case; master = [1], forcing = …)
W, R = parametrise(model, spectral, expansion_order)

Keyword arguments are backend-specific: they may select master modes, add forcing, configure diagnostics, or request/reuse an eigensolve. The return shape is shared, so everything downstream of build_model can remain physics-independent.

The return is a NamedTuple, not a positional tuple, so that a backend can grow what it reports without breaking callers who only want model and spectral.

• model — an NthOrderModel

• spectral — a SpectralData

• meta — backend metadata such as timings, the raw spectrum, forcing records, or DOF maps. It is available to callers and reporting code but is not consumed by MORFE.

Implementations must:

• return an NthOrderModel whose linear_terms and nonlinear terms are complete, including any external system the forcing introduces — not "mostly built, the caller adds forcing";

• return a SpectralData reconciled against that model's order. Use SpectralData(model, spectrum; master = …) so MORFE owns any order reconciliation;

• apply conditioning tweaks (mode scaling, unit changes) to the raw arrays before constructing the bundle; SpectralData deliberately has no scale field, so such tweaks stay visible at the call site;

• provide the spectrum's actual conjugate pairing to SpectralData; if a full reduced variable permutation is needed, extend the master pairing through the model's external system with full_conjugate_permutation rather than a fixed literal.

Implementations must not build a MultiindexSet, a ResonanceSet or a resonance policy, validate the monomial set, derive conjugate closure, warn about resonances, or solve the cohomological equations. Those are parametrise's responsibilities and are chosen by the caller. A backend may solve its linear eigenproblem here, or reuse spectral data supplied by the caller.

free_dofs_at_nodes(dh, free_to_local, node_ids, directions) -> Vector{Int}

Free-DOF indices (the row indices of K, M and of the parametrisation W) for the given mesh nodes and directions. free_to_local is the map built by mechanical_model and available as model.info.free_to_local.

Throws if a requested node/direction is constrained, since it then has no row.

node_dof(dh, node_id, direction) -> Int

Global DOF index of direction (1 = x, 2 = y, 3 = z) at mesh node node_id. Assumes a single vector field with one DOF per direction per node, which is what mechanical_model builds.

stage_timings(info) -> Vector{Pair{String, Float64}}

Every *_time_s entry of an info NamedTuple, in declaration order, with the suffix stripped.

This is the convention that lets one writer time every backend without knowing any of them: a stage that records eig_time_s in its info is reported as "eig", and a backend adds a timed stage simply by recording it.

summary_entries(case, rom) -> Vector{Pair{String, Any}}

The physics-specific rows a backend contributes to its run summary, in the order they should be written.

Add a method for your assembled-model type to extend the summary; the shared skeleton (sizes, order, stage timings) is written by write_summary and does not need restating. Returning an empty vector — the default — is fine.

function Common.summary_entries(m::AssembledFluidModel, meta::NamedTuple)
    return ["Re0" => m.Re₀, "mode_scale" => meta.scale]
end

A row whose key repeats one the skeleton already wrote overrides it, so a physics can also correct a shared row rather than only append to it.

write_summary(io, path, case; rom = nothing, title, metadata = Pair[],
			  append = true) -> Nothing

Print a run banner and timing table to io, and write key: value rows to path.

The skeleton is shared: title, problem size, stage timings gathered by stage_timings from both case.info and the ROM meta, then whatever summary_entries the physics contributes, then metadata from the caller (which wins over both — it is the most specific).

append = true is the default because MORFE.save_rom writes this same file first, with the run's provenance (julia version, commit, timestamp). Opening it with "w" here would silently discard that.

load_comsol_grid(path, dirichlet_entity_ids; scale = 1.0) -> (grid, constrained_nodes)

Read a COMSOL .mphtxt file and return the Ferrite Grid (P18 quadratic prism cells) together with the Set{Int} of node indices lying on the Dirichlet boundaries. dirichlet_entity_ids is the set of COMSOL geometric entity IDs (1-indexed, i.e. raw file ID + 1) whose surface elements define the clamped region.

scale multiplies every node coordinate on load — use it when the mesh and the material constants are in different unit systems (e.g. scale = 1e-3 for a mesh drawn in millimetres against SI material data). Scaling the geometry here rather than afterwards keeps the DOF numbering and the entity sets untouched.

abaqus_to_gmsh(inp_file::String, gmsh_file::String)

Read an Abaqus .inp mesh and write a Gmsh .msh file.

Supported element types: C3D4/10/10M, C3D8/8R/8I, C3D20/20R, C3D6, C3D15, S3/3R, S4/4R, S6, S8/8R, CPS3/4/4R/6/8, CPE3/4/4R/6/8, CAX3/4/4R/6/8/8R, B31, B32, T3D2, T3D3.

Node IDs may be non-contiguous; they are remapped to a compact 1-based Gmsh tag space. Unknown element types are silently skipped.

abaqus_to_gmsh_linear(inp_file::String, gmsh_file::String)

Like abaqus_to_gmsh but downgrades quadratic elements to their linear (corner-nodes-only) counterparts before writing.

_read_mesh( mesh_file::String)

Reads a Comsol .mphtxt file and extracts

• n2c: node coordinates

• e2nT6: node list of elements of type T6 (triangle)

• e2qT6: label of elements of type T6 (triangle)

• e2nQ9: node list of elements of type Q9 (rectangle)

• e2qQ9: label of elements of type Q9 (rectangle)

• e2nT10: node list of elements of type T10 (triangle)

• e2qT10: label of elements of type T10 (triangle)

• e2nP18: node list of elements of type P18 (prism)

• e2qP18: label of elements of type P18 (prism)

• e2nH27: node list of elements of type H27 (hex)

• e2qH27: label of elements of type H27 (hex)

Based on the mesh.jl file from Morfe2.0.

comsol_to_gmsh(comsol_file::String, gmsh_file::String)

Reads comsol_file (.mphtxt) and returns a Gmsh (.msh) file. Warning: The Node Re-Ordering is done by hand comparing with the old comsol files from Morfe2.0. So be aware that Comsol may have changed the structure or numbering of nodes. Here is the permutation only implemented for T6, Q9 and P18

comsol_to_gmsh_linear(comsol_file::String, gmsh_file::String)

Like comsol_to_gmsh but downgrades quadratic elements to their linear counterparts before writing the Gmsh .msh file.

gmsh_to_comsol(gmsh_file::String, comsol_file::String)

Reads gmsh_file (.msh) and writes a COMSOL (.mphtxt) mesh file. This is the inverse of comsol_to_gmsh: it undoes the node-reordering permutations and converts 1-based Gmsh indexing back to 0-based COMSOL indexing.

Assumptions:

• The .msh file was produced by comsol_to_gmsh (node tags 1:nn).

• Element tags in the Gmsh file carry the COMSOL geometric-entity indices (as written by comsol_to_gmsh).

adj_tensor_type(cache) -> Type

The tensor type of the adj J coefficients, i.e. Tensor{2,dim,Float64,dim²}. The driver sizes its gradient buffers from this rather than from a hardcoded 3D alias, which is what makes the assembly loops dimension-general.

adjugate_series(J, det_ser, basis[, supp]) -> Vector{<:Tensor{2,dim}}

Exact multivariate coefficients of adj(J(θ)), truncated to the box.

Theory App. A.2 (Adjugate). Matching θ^σ in the defining identity J · adj J = det J · I gives a recurrence that is strictly triangular in graded-lex order (β ≤ σ componentwise with β ≠ σ forces |β| < |σ|):

J₀ A_σ + Σ_{0<α≤σ} J_α A_{σ−α} = c_σ I

The appendix solves this by assuming J₀ = I, which this code cannot. A curved reference configuration has J₀ ≠ I — example 07's arch is exactly that, J₀ = I + ∇w ⊗ e₁ — so the leading term must be inverted explicitly:

A_σ = J₀⁻¹ ( c_σ I − Σ_{0<α≤σ, α ∈ supp J} J_α A_{σ−α} ),   A_0 = adj(J₀)

which reduces to the appendix's form when J₀ = I. The same correction applies to Method 3, where the leading matrix is the c₀-weighted mass matrix, not M.

The sum runs only over the support of J — Nθ+1 terms for an affine map, not L — so this costs O(L·|supp J|) per quadrature point. basis.diff supplies the position of σ − α without hashing.

Exactly polynomial with deg adj J ≤ (d−1)·deg J. Truncating det J to the box does not corrupt it: A_σ reads only c_σ, which is in the box whenever A_σ is.

assemble_linear_series!(A_arr, B_arr, pd, kernel)

Fill A_arr[m], B_arr[m] with the θ^α coefficient matrices of the two linear operators kernel defines (for a structure: stiffness and mass), at multiindex basis.mset.exponents[m].

A_arr and B_arr must each be nterms(basis)-long vectors of sparse matrices sharing the DOF handler's pattern (allocate_matrix(dh)). The driver owns the cell loop and the assembly; the kernel supplies only the per-entry θ-series through linear_qp_series!.

build_inverse_determinant(strategy, dh, cv, basis, det) -> Vector{Vector{Vector{Float64}}}

The θ-series of 1/det J, indexed [cell][qp], from the already-computed det J series in the same layout. The single seam between Methods 3 and 4.

build_linear_corrections(A_arr, basis, arity) -> Vector{MultilinearMap}

−θ^α · A_α · x for every α ≠ 0, where A_arr[αidx] is the operator's θ^α coefficient matrix and arity names the derivative slot x occupies.

The base coefficient (α = 0) is skipped: it belongs in the model's linear_terms, not among its nonlinear terms. Structurally empty coefficients are skipped too — an all-zero correction is a term the solve would evaluate for nothing.

For a second-order structure on the augmented ORD = 3 model: (1,0,0) is a stiffness correction, (0,1,0) damping, (0,0,1) mass. For a first-order physics, (1,).

build_model(m::AssembledParametricModel; master, spectrum, …)
	-> (; model, spectral, meta)

Assemble the augmented full-order model and its spectral data for a parametric problem.

The θ parameters enter as frozen external states: N_EXT = n_geometry_parameters external variables with eigenvalue 0 (θ̇ = 0), so a monomial z^a θ^b in the reduced dynamics is a genuine parameter dependence rather than a transient.

What this method does, and what it deliberately leaves to the caller:

• Zero-pads the base physics' linear terms to model_order, because a correction on the highest derivative (a parametric mass) only exists on a model one order higher.

• Collects every θ ≠ 0 correction of every operator, plus the nonlinear maps.

• Reconciles the spectral data against the augmented model's order with SpectralData(model, spectrum; master). The base eigenproblem is second-order while the model is third — this is exactly the reconciliation MORFE owns, and hand-rolling it (multiplying the last block by λ versus forming a fresh λ^{k-1}ψ) is silently wrong.

• Derives the conjugate permutation with full_conjugate_permutation, never a literal — the θ states are real, hence self-conjugate, so an odd n_geometry_parameters works without a special case.

• Does not build a MultiindexSet. The θ-box truncation is a modelling choice; pass your own mset to parametrise.

master lists physical mode PAIRS, as elsewhere: pair p occupies spectrum entries 2p-1, 2p.

convolve_weight_accumulate!(dst, A, w, scale, basis)

dst[k] += scale · Σ_{α+β=k} A[α] · w[β] — the determinant weighting fused into the element accumulation.

The driver's last act per basis function is to multiply the integrand series by (1/det J)^p and add it to the element residual. Doing that as poly_mul followed by a loop materialises a length-L temporary for every single (cell, quadrature point, basis function); fusing removes it entirely.

det_adj_series(J, basis) -> (det_ser, adj_ser)

det(J(θ)) and adj(J(θ)) together, sharing one scan of J's support.

Both are exactly polynomial — App. A.1 and A.2 — which is what makes the adjugate identity J⁻¹ = adj J / det J worth using: it concentrates all the non-polynomial content of J⁻¹ into the single scalar 1/det J, left to reciprocal_series (Method 4) or AuxiliaryFieldInverseDet (Method 3).

determinant_series(J, basis[, supp]) -> Vector{Float64}

Exact multivariate coefficients of det(J(θ)), truncated to the box.

Theory App. A.1 (Determinant). The determinant is multilinear in the columns of J, so expanding each column's series independently and collecting the terms with α₁+…+α_d = σ gives c_σ directly — one d×d scalar determinant per choice of support multiindex per column. Dimension-general, and no cofactor formula appears.

Exactly polynomial with deg det J ≤ d·deg J, per parameter as well as in total.

free_dof_map(ndofs, freedofs) -> Vector{Int32}

Dense global-DOF → free-DOF index map, 0 where the DOF is constrained.

A Dict here costs one hash lookup per DOF per cell per θ-multiindex — scatter_local! runs once for each of the L coefficients of every cell — which is a six-figure number of hashes per sweep for no reason. The dense vector is 4·ndofs bytes and indexes directly. Physics modules outside this one keep their Dict form (model.info.free_to_local); only the parametric assembly path uses this.

geometry_validity_at(cache, θ; points = nothing) -> NamedTuple

Measure Method 4's validity conditions at one concrete θ.

Returns (; det_min, det_max, max_deviation, reciprocal_residual, worst_cell, worst_qp, convergent, orientation_ok):

• max_deviation — max |det J(θ) − 1| over the probed points. Must be < 1.

• det_min — must be > 0; a non-positive determinant means the coordinate transform has folded the mesh and is not invertible.

• reciprocal_residual — max |det J(θ) · (1/det J)(θ) − 1|, the TRUNCATED series measured against its own defining identity. This is the honest error of the expansion at this θ, computed from what the cache actually holds.

points restricts the sweep to a (cell, qp) subset; nothing sweeps all.

geometry_validity_report(cache; max_probe_points = 5000, cap = 4.0) -> NamedTuple

The θ-range over which Method 4's expansion is valid, measured per parameter.

Returns (; radius_pos, radius_neg, unconditional, det_at_reference), where radius_pos[i] / radius_neg[i] is the largest θᵢ along ±eᵢ (all other parameters zero) for which every probed quadrature point satisfies both |det J − 1| < 1 and det J > 0. Inf means the scan reached cap without leaving the radius.

unconditional is true when det J ≡ 1 everywhere in θ — a volume-preserving transform, for which the reciprocal series is the single term 1 and no convergence question arises.

Axis-wise, because a full box scan costs 2^Nθ corners and the axes already expose the parameter responsible. For a specific θ of interest — including cross terms — call geometry_validity_at directly.

Diagnostic only: nothing here changes the model.

inv_det_power(inv_det, n, basis) -> Vector{Float64}

n-th power of the reciprocal series (n ≥ 1), truncated to the box.

inv_det_power_series(cache, p) -> Vector{Vector{Vector{Float64}}}

The precomputed (1/det J)^p series, indexed [cell][qp]. Hoist this out of an assembly loop — the Dict lookup is not meant for the hot path.

inverse_determinant_comparison(a::PullbackCache, b::PullbackCache, θ; points = nothing)

Compare two inverse-determinant strategies at one concrete θ.

a and b must be built on the same mesh and quadrature (same (cell, qp) layout); their θ-bases may differ, since each series is evaluated in its own. Typically a is a PowerSeriesInverseDet cache (Method 4) and b an AuxiliaryFieldInverseDet one (Method 3).

Returns (; max_abs, max_rel, rms, residual_a, residual_b, worst_cell, worst_qp):

• max_abs / max_rel / rms — how far apart the two 1/det J fields are at the quadrature points, which is exactly where the difference enters the assembled operators.

• residual_a / residual_b — each method measured against det J · s = 1, its own defining identity. Neither method is used as the reference for the other; they are each compared with the truth they are both approximating.

Diagnostic only: nothing here changes a model.

jacobian_series(Js, basis) -> Vector{<:Tensor{2,dim}}

Assemble the degree-1 Jacobian series aligned to basis from Js = (J₀, ∇ψ₁, …, ∇ψ_{Nθ}): J₀ at the zero multiindex, ∇ψ_i at the unit multiindex e_i. All other coefficients are zero.

This is the AFFINE map x = x₀ + Σ_i θ_i ψ_i(x₀). For the theory's general polynomial x = Σ_α x_α θ^α, pass the coefficients as multiindex ⇒ tensor pairs instead — see the AbstractVector{<:Pair} method.

jacobian_series(Js::AbstractVector{<:Pair}, basis) -> Vector{<:Tensor{2,dim}}

The theory's general polynomial map x(x₀,θ) = Σ_α x_α(x₀) θ^α, given as multiindex => J_α pairs (App. A, opening). The multiindex may be any Nθ-element integer container; repeated multiindices accumulate, matching the Σ_α they stand for.

A coefficient whose multiindex falls outside the box throws rather than being silently dropped. Truncating the series computed from the map is this module's job; truncating the map itself would change which geometry is being modelled without saying so — the exact failure mode that made example 04 wrong.

The affine NTuple method above is the common case and builds this structure.

model_order(m::AssembledParametricModel) -> Int

The ORD the augmented NthOrderModel needs: one more than the highest derivative slot any θ-correction occupies.

A second-order structure whose mass is parametric needs ORD = 3 — a mass correction has arity (0,0,1), which only exists on a third-order model, and the fourth linear block is zero. A first-order physics with only a (1,)-arity correction needs no augmentation and must not get one.

multilinear_maps(m::ParametricMap; arity)

One MORFE MultilinearMap per θ-multiindex α: the θ^α coefficient of the kernel's form, with external multiplicity |α| so the solve multiplies it by the matching product of frozen θ states.

arity is the modal arity tuple — length ORD, with DEG in the slot the form acts on. For a second-order structure with the augmented ORD = 3 model, a quadratic displacement form is (2, 0, 0).

position_of(basis, e) -> Int

Grlex position of exponent e, or 0 when it falls outside the box (including any negative component). The hash-free counterpart of basis.index[e].

reciprocal_series(p, basis) -> Vector{Float64}

Coefficients of 1/p(θ) truncated to the box, from the graded recurrence q[0] = 1/p[0], q[γ] = -(1/p[0]) Σ_{0≠β≤γ} p[β] q[γ-β]. Graded-lex order makes each q[γ-β] (lower total degree) available before q[γ].

This is Method 4 of the theory write-up (App. A.3, "Reciprocal of the Determinant"). The recurrence is the geometric series for 1/(1 − (1 − det J)) in disguise, so it inherits that series' radius: applied to det J it is valid only where |det J − 1| < 1, and it is unbounded in degree because the reciprocal of a polynomial is not a polynomial. basis.diff supplies the position of γ − β without hashing — this runs at every quadrature point.

See geometry_validity_report for the measured radius, and AuxiliaryFieldInverseDet for Method 3, the alternative.

report_geometry_validity(cache; kwargs...) -> NamedTuple

Emit the measured convergence radius of the inverse-determinant expansion, and return the underlying geometry_validity_report.

Silent for a volume-preserving transform (det J ≡ 1), which is unconditionally convergent. Errors — not warns — when the REFERENCE configuration itself is already folded (det J ≤ 0 at θ = 0), because every series in the cache is then built about an invalid point.

report_zero_coefficients(m::AssembledParametricModel; rtol = 1e-12, io = stderr)

Emit one @info block naming the θ-multiindices that contribute nothing, and return the underlying zero_coefficient_report.

Silent when every multiindex carries something — the common case should not produce output.

series_extent(A) -> Int

Position of the last non-zero coefficient of a θ-series, 0 if it is all zero.

This is the theory's degree bound, measured rather than assumed. adj J has degree ≤ (d−1)·deg J and det J degree ≤ d·deg J, and graded-lex order makes "total degree ≤ k" a contiguous prefix, so everything past this index is exactly zero and can be skipped without any truncation. It adapts on its own: a volume-preserving transform (det J ≡ 1) reports 1 and its convolutions collapse to a single term, with no degree bookkeeping threaded anywhere.

Scanning from the end costs O(L) iszero calls and usually returns at once; testing the same thing per PAIR costs O(|prod|), which a profile of example 04 showed to be 23 % of the whole run.

sweep_all!(out, m::ParametricMap{DEG}, us::NTuple{DEG})

Fill out[:, α] with the θ^α coefficient of the kernel's form evaluated at the DEG displacement inputs us, for every α in the basis.

Computing the whole series costs the same per coefficient as computing one, so this is the only sweep the per-α closures need.

zero_coefficient_report(m::AssembledParametricModel; rtol = 1e-12) -> NamedTuple

Which θ-multiindices contribute nothing to the assembled model.

Returns (; geometry, operators, maps, exponents):

• geometry[i] — true when multiindex i is negligible in every geometry series (adj J, det J) at every quadrature point. The coordinate transform simply does not reach that degree.

• operators[i] — true when multiindex i is negligible in every linear operator's θ^α coefficient matrix. The transform may reach that degree while the weak form still annihilates it.

• maps[i] — true when multiindex i is negligible in every nonlinear form, probed at a pseudo-random state. A form that is identically zero gives zero for any input; a form that is not gives a nonzero result for a random input with probability 1. Set probe = false to skip (it costs one sweep per map).

• exponents — the multiindices themselves, for reporting.

The three are reported separately because they call for different responses: a geometry gap means the θ-basis can shrink at no cost, whereas an operator or form that annihilates a degree the geometry reaches is a property of the weak form or of a symmetry, and is worth understanding before it is relied on.

rtol is applied against the largest coefficient of the same quantity, so a model in newtons and one in millinewtons give the same answer.

Diagnostic only: nothing here changes the model. Use it to see whether a θ-basis is larger than the geometry justifies, or whether a symmetry you expected is actually present in the discretisation.

CubicCrystal(; c11, c12, c44, ρ, rotation = nothing) -> AnisotropicMaterial

Cubic crystal (silicon, germanium, …) from its three independent constants, with the crystal optionally rotated into the lab frame by rotation (a 3×3 matrix, or an angle in radians about the z axis).

rotation = nothing leaves the crystal axes unchanged; a scalar is interpreted as an angle in radians about the z axis, and a matrix is passed to rotate_voigt.

The isotropic limit is c11 = λ + 2μ, c12 = λ, c44 = μ; it reproduces the constitutive response of an SVKMaterial with those Lamé constants.

_forcings(forcing) -> Vector{<:HarmonicForcing}

Normalize the public forcing forms—nothing, one HarmonicForcing, or an abstract vector of them—to a vector. The autonomous case returns an empty vector.

assemble_KM!(K, M, dh, cv, λ, μ, ρ)
assemble_KM!(K, M, dh, cv, stress, ρ)

Assemble the three-dimensional global stiffness matrix K and mass matrix M into preallocated sparse matrices using isotropic Lamé constants or an internal stress-law object. The matrices must have the sparsity pattern associated with dh; both are mutated and the function returns nothing.

K_rs = ∫ ε(φ_r) ⊡ σ(ε(φ_s)) dΩ
M_rs = ∫ ρ φ_r · φ_s dΩ

For the Lamé overload, σ(ε) = λ tr(ε)I + 2με.

base_operators(m::AssembledParametricModel) -> (K, M)

Return the θ = 0 stiffness and mass coefficient matrices of an assembled parametric SVK case. These are the base-configuration operators used to solve the eigenproblem supplied to build_model.

build_model(m::AssembledMechanicalModel; master = [1], forcing = nothing,
			nev = ..., spectrum = nothing, expansion_order = nothing,
			n_monomials = nothing) -> (; model, spectral, meta)

Construct the full-order MORFE.NthOrderModel and MORFE.SpectralData consumed by MORFE.parametrise.

master is a nonempty, sorted vector of distinct positive physical mode-pair indices. The eigensolver returns adjacent conjugate eigenvalues, so pair p occupies spectrum entries 2p-1, 2p; for example, master = [1] gives ROM = 2 internal coordinates. Selected pairs may be non-leading and non-contiguous.

forcing accepts nothing, one HarmonicForcing, or a vector of them and represents

f(t) = Σₖ aₖ M ϕ_{pₖ} cos(Ωₖ t).

Forcing k appends external coordinates ROM+2k-1, ROM+2k with eigenvalues +iΩₖ, -iΩₖ. If Ωₖ is omitted it is abs(spectrum.eigenvalues[2pₖ-1]). The shape pair pₖ need not be a master pair.

If spectrum is omitted, spectrum is called with nev, where nev counts physical modes (and the returned spectrum contains 2nev conjugate eigenvalues). A supplied spectrum is reused without another eigensolve. In either case it must contain every master pair and every pair used as a forcing shape.

expansion_order and n_monomials size the batched-column cache owned by the FEM nonlinear terms; they do not construct or select MORFE's monomial set. For the standard total-degree set, expansion_order = p reserves binomial(ROM + N_EXT + p, p) - 1 columns. For a custom set, pass n_monomials = length(mset). n_monomials takes precedence if both are supplied; if neither is supplied, one column is reserved.

The result contains:

• model: the autonomous or forced MORFE.NthOrderModel;

• spectral: spectral data restricted to master, with the spectrum-wide adjacent conjugate pairing and the external pairing supplied by the model's external system;

• meta: backend reporting data, including the resolved forcings/frequencies, selected spectrum indices, external-state count, eigensolve time, and assembled-case metadata.

eigenfrequencies(m::AssembledMechanicalModel; nev = 10, eigensolver = nothing)
	-> Vector{ComplexF64}

Return the complex damped eigenvalues from spectrum. Physical mode p occupies adjacent conjugate entries 2p-1, 2p. For an underdamped pair, abs(imag(λ[2p-1]))/(2π) is its damped oscillation frequency in Hz; the undamped natural frequency ω used by the Rayleigh construction is distinct. For an underdamped or critically damped mode abs(λ) = ω; this identity does not hold for each individual eigenvalue of an overdamped pair.

Use the result to inspect the spectrum and choose master before parametrisation.

mechanical_model(grid::Ferrite.Grid, constrained_nodes::Set{Int};
                 material, damping, fe_order = 2, quad_order = fe_order + 1)
mechanical_model(grid::Ferrite.Grid; material, damping, dirichlet,
                 fe_order = 2, quad_order = fe_order + 1)
mechanical_model(mesh_path::AbstractString; material, damping, dirichlet,
                 scale = 1.0, fe_order = 2, quad_order = fe_order + 1)

Assemble a three-dimensional AssembledMechanicalModel from a Ferrite grid or mesh file. The displacement interpolation is Lagrange{RefShape,fe_order}()^3; the returned K, M, and C = αM + βK are restricted to free DOFs, while the model retains the Ferrite handlers needed to construct its quadratic and cubic SVK terms lazily.

material is an SVKMaterial or AnisotropicMaterial, and damping is a RayleighDamping. All three displacement components are fixed on the selected nodes or facets:

• The constrained_nodes overload accepts an already computed set of node indices.

• For a Ferrite grid or Gmsh .msh path, dirichlet is the name of a facet set.

• For a COMSOL .mphtxt path, dirichlet is a Set{Int} of 1-based boundary entity IDs (the raw COMSOL IDs plus one). scale rescales COMSOL node coordinates before assembly; it is ignored for Gmsh input.

parametric_model(dh, cv, geometry; geometry_parameter_basis, material, damping,
				 free = nothing, base = nothing,
				 inverse_determinant = PowerSeriesInverseDet())
	-> AssembledParametricModel

Assemble the parameter expansion of a three-dimensional SVK structure over a parametric mesh coordinate transform.

geometry is a provider returning (J₀, ∇ψ₁, …, ∇ψ_Nθ) per quadrature point, either analytically (geom(x₀)) or from an FE field (geom(x₀, cell, cv, q)).

geometry_parameter_basis is either one GeometryParameterBasis, shared by every form, or a NamedTuple with a required linear entry and optional quadratic and cubic entries. A missing quadratic basis falls back to linear; a missing cubic basis falls back to quadratic, then linear. Separate bases avoid over-expanding low-degree forms—for example, an isochoric affine arch has stiffness degree at most two but cubic-form degree at most four.

material is an SVKMaterial or AnisotropicMaterial. damping is a RayleighDamping and defines the assembled series C(θ) = αM(θ) + βK(θ). If free is omitted, every DOF in dh is retained; otherwise it supplies the global free-DOF indices. base is stored in the returned case for physics-side bookkeeping. inverse_determinant selects how the pullback cache constructs the series for 1/det(J).

This function returns an assembled parametric case, not an NthOrderModel. Calling build_model on it produces ORD = 3: a parameter-dependent mass is a correction on the highest derivative of the original second-order system, so the augmented representation needs one additional, zero linear block.

print_mode_table(eigenvalues; master = Int[], io = stdout)

Tabulate adjacent conjugate pairs in eigenvalues, ignoring an unmatched final entry. Each row reports real(λ) as the decay rate and imag(λ) as the damped angular frequency, together with imag(λ)/(2π) in Hz. Pairs listed in master are marked. Output is written to io and the function returns nothing.

print_resonances(mset, resonance_set, master_eigenvalues; io = stdout)

Print the monomials flagged as resonant for each inner target and any diagnostic outer targets. Inner targets are labelled with master_eigenvalues; because this function is not passed the outer eigenvalues, outer target labels use 0 as a placeholder. Inner resonances are the terms retained in the reduced dynamics R; outer flags diagnose off-manifold solves and do not add rows to R.

probe_dof(m::AssembledMechanicalModel, node, direction) -> Int

Free-DOF index of direction (1 = x, 2 = y, 3 = z) at mesh node: the row of K, M and of the parametrisation W that carries that node's displacement.

This is the bridge between a physical quantity named the way a person names it ("the transverse displacement at mid-span") and the integer MORFE.observable_polynomial wants:

u = observable_polynomial(W, SVK.probe_dof(case, 289, 2))

Throws if the node is constrained, since a constrained DOF has no row. Wraps Common.free_dofs_at_nodes for the single-node case.

resonances(eigenvalues, master, order; external_eigenvalues = ComplexF64[],
		   outer_eigenvalues = ComplexF64[], resonance_tol = 0.05,
		   resonance_tol_rel = nothing, mset = nothing)
	-> (mset, resonance_set, master_eigenvalues)

Build the monomial set and the complex-normal-form ResonanceSet that a parametrise call with the same arguments would use. Call it on the output of eigenfrequencies to preview — before paying for the cohomological solve — which monomials will be kept in the reduced dynamics.

master uses the same sorted physical-pair indexing as build_model: pair p occupies eigenvalues[2p-1:2p]. external_eigenvalues appends prescribed external coordinates to the monomial variables. If mset is omitted, all monomials of total degree 1:order in the internal and external variables are created; otherwise the supplied set is returned and order does not choose its contents.

resonance_tol is the absolute detuning threshold in the eigenvalues' frequency units (rad/s for structural spectra). When resonance_tol_rel is given, the inner target λⱼ instead uses resonance_tol_rel * abs(λⱼ).

outer_eigenvalues adds off-manifold resonance targets, populating the outer_resonances block of the returned set (query it with resonant_multiindices(rset, ROM + j)). Those flags are diagnostic: the cohomological solve reads only the inner block. It cannot be combined with resonance_tol_rel, because MORFE shares one tolerance object between the inner and outer target blocks.

rotate_voigt(D, Q) -> SMatrix{6,6}

Rotate a Voigt stiffness matrix by a 3×3 matrix Q, expected to be orthogonal. Only its size is validated.

Q maps crystal axes to lab axes: if D is expressed in the crystal frame, the result is expressed in the lab frame (pass Q' for the opposite convention).

Implemented by expanding D to the 4th-order stiffness C_ijkl, applying the tensor transformation C'_ijkl = Q_ip Q_jq Q_kr Q_ls C_pqrs, and contracting back — so no Bond-matrix convention (and its factor-of-2 traps) is involved.

spectrum(m::AssembledMechanicalModel; nev = 10, eigensolver = nothing)
	-> Spectrum

Solve the assembled model's eigenproblem. nev requests physical modes; the result holds 2nev eigenvalues and eigenmodes in adjacent conjugate pairs.

The default is StructureModalDampingEigensolver built from the model's own Rayleigh coefficients m.damping. It solves the undamped problem K ϕ = ω² M ϕ once by shift-invert, mass-normalises ϕ, and builds the damped eigenvalues and the left eigenvector blocks in closed form, which Rayleigh damping makes possible. It never forms the first-order 2n × 2n pencil, so it stays fast on large meshes.

A supplied eigensolver is used as-is. In particular, callers supplying a different modal-damping solver are responsible for making its damping consistent with m.C.

Pass the returned MORFE.Spectrum to build_model as spectrum = ... to inspect or report it without paying for a second solve—and without letting a repeated iterative eigensolve choose a different basis in a clustered eigenspace.

svk_assemble_KM!(K, M, dh, cv, λ, μ, ρ)
svk_assemble_KM!(K, M, dh, cv, material)

Assemble the three-dimensional linear stiffness K and mass M matrices in place with the Ferrite SVK backend. Supply either Lamé constants and density or an SVKMaterial/AnisotropicMaterial. K and M must be preallocated with the sparsity pattern of dh; the function returns nothing.

svk_nonlinearity(degree, dh, cv, free_to_local, n_free, λ, μ;
                 max_unique_cols = degree, fully_asymmetric = false)
svk_nonlinearity(degree, dh, cv, free_to_local, n_free, material;
                 max_unique_cols = degree, fully_asymmetric = false)

Construct a Ferrite-backed St. Venant-Kirchhoff geometric nonlinearity term of the given polynomial degree (2 for the quadratic form, 3 for the cubic form) as a MORFE.FEMMultilinearMap{2}. The constitutive law may be supplied as Lamé constants λ, μ or as an SVKMaterial/AnisotropicMaterial.

free_to_local maps global Ferrite DOFs into the n_free-component state. The max_unique_cols cache must accommodate the column batches used by the reduction; fully_asymmetric is forwarded to MORFE's multilinear-term symmetry policy.

voigt_stiffness(material) -> SMatrix{6,6}

Return the SMatrix{6,6,Float64} Voigt stiffness of a supported material. AnisotropicMaterial returns its stored matrix; SVKMaterial constructs the isotropic matrix from its Lamé constants. This is useful for inspection and for cross-checking anisotropic input.

_assemble_element!(Ke, Re_e, u_e, p_e, cv_vel, cv_pres, dof_range_u, dof_range_p, Re0)

Compute element Jacobian Ke and residual Re_e for the steady NSE at Re0.

Quadrature loop fills three blocks: • vel–vel : viscous + convective tangent (two terms) • vel–pres : −G^T (pressure gradient term in momentum) • pres–vel : G (divergence constraint)

_assemble_linear_element!(Me, ALe, u₀_e, cv_vel, cv_pres,
						  dof_range_u, dof_range_p, Re0)

Compute element mass matrix block Me (velocity DOFs only) and element linearised-operator matrix ALe (full element DOF space).

ALe contributes to Alin; the caller negates it for B₀ = −Alin.

_enforce_affine_2d_fluid_structure!(cache)

Validate and enforce the exact polynomial degrees of the one-parameter, two-dimensional affine fluid map. For F(mu) = F0 + mu*G, det(F) is quadratic and adj(F) is affine. The generic recurrence used to construct a PullbackCache can leave roundoff-sized coefficients above those degrees; those values must not be interpreted as physical convection maps merely because they are bitwise nonzero.

This is a narrowly bounded roundoff cleanup. When all coefficients above the affine analytical degree are at roundoff scale, the cache is identified as affine and those coefficients are replaced by exact zeros. A genuinely polynomial geometry provider is left unchanged; the shared fluid API continues to support the general maps accepted by PullbackCache.

_estimate_period_autocorr(signal, dt, τ_min, τ_max) -> (T_est, c_peak)

Estimate the dominant period of signal (assumed mean-subtracted) as the lag of the highest NORMALISED autocorrelation inside [τ_min, τ_max], refined to sub-sample accuracy by a parabola through the peak and its neighbours.

Returns (NaN, 0.0) when the peak sits on the window edge (no interior maximum). c_peak ∈ [−1, 1] is the normalised correlation at the peak — ≈ 1 for a genuinely periodic signal, ≈ 0 for noise; gate on it.

Zero-crossing counting is NOT robust here: the lift is a pressure functional, and pressure — the algebraic variable of the descriptor system — carries fast components invisible to the velocity mass norm that pollute the crossings. Fast components only correlate at their own short lags, outside the gated window, so the autocorrelation peak stays on the shedding period.

Internal bilinear convection action with independently selected state layouts.

Internal symmetric bilinear convection action on an ordinary physical mesh.

slaved_R1(R, ρ, η) → ComplexF64

First component of the reduced dynamics at the canonical phase z₁ = z̄₁ = ρ, η′ = η, with every promoted coordinate SLAVED to its quasi-steady value. Reduces to a plain evaluation when nothing was promoted.

The promoted coordinates are slaved, never zeroed. They carry the mean-flow distortion — ẏk is driven by z₁z̄₁ — and hand it back to the oscillator through z₁·yk, which is the dominant stabilising contribution to the Landau coefficient. Zeroing them removes it and reports this supercritical Hopf as subcritical.

On the orbit they are quasi-steady, so they solve R_k(ρ, ρ, y, η) = 0. The monomial set carries at most ONE promoted coordinate to the first power, so R is exactly AFFINE in y: this is one small linear solve, not an iteration and not an approximation.

_harmonic_R1(R, ρ, η, Ω) → ComplexF64

Fundamental (harmonic s = 1) component of R₁ on the orbit z₁ = ρe^{iΩt}, with the promoted coordinates closed by HARMONIC BALANCE.

Each monomial z₁^a z̄₁^b η^c y_j sits at harmonic s = a − b, so on the orbit the drive of promoted row k splits as b_k(t) = Σ_s b_{k,s} e^{isΩt} and its response solves

(i s Ω I − A₀) y_s = b_s

rather than A y = −b. The old quasi-steady form is the s = 0 case of this and is exact only for a drive that is constant on the orbit — true for the mean-flow modes (all of the 6-mode set), false for the −6.542148 ± 18.415390i pair, which is driven at s = ±1. There the denominators differ by a factor 2.9 in magnitude plus a large phase error: |−λ| = 19.54 against |iΩ − λ| = 6.72. The small one IS the near-resonance that made the mode resonant in the first place (detuning 6.73 < 8.43), so slaving discarded exactly the amplification worth capturing.

APPROXIMATION, deliberate: only the s = 0 part of A = ∂R_k/∂y_j is kept, so harmonics do not couple through A. A's nonzero harmonics come from monomials like z₁²y_j, which are higher order in ρ than the λ_k y_k diagonal that dominates it. Lifting this would need a block-coupled solve across harmonics.

Frequency of the linear Hopf mode — the seed for the fixed point on Ω.

_master_conjugate_pairing(λ; atol = 1e-8) -> Vector{Int}

The involution σ over a master set: σ[r] = s where λ[s] = conj(λ[r]).

Two cases, and the first is the one that makes a mode on the real axis usable at all:

• λ[r] is real ⇒ it is its own conjugate, so σ[r] = r. A non-oscillatory mode carries a single real coordinate rather than half of a complex pair.

• λ[r] is complex ⇒ its partner must also be in the master set. If it is not, that is an error rather than a fallback: a manifold spanned by one half of a conjugate pair is not invariant under conjugation, and the resulting ROM has no real realisation.

This is constructed, not detected. Detection over the selected modes fails here because the shift-invert eigensolve uses a COMPLEX shift σ and therefore returns only the modes near σ — a strongly oscillatory mode's conjugate sits near conj(σ) and is simply not computed. Those conjugates are synthesised analytically by the caller (λ̄, φ̄), and this function then pairs them.

_periodic_lift(FL, K) -> FL_per

Project the one-period lift signal onto its first K temporal harmonics (the measurement orbit spans EXACTLY one period at uniform steps, so the DFT bins are exact Fourier coefficients of the periodic signal).

The converged orbit is T-periodic by construction, so any non-harmonic content in the measured signal is numerical: under Crank–Nicolson (θ = ½) the algebraic pressure mode of the descriptor system has amplification factor −1, and seed inconsistencies ring at the Nyquist frequency forever — invisible to the velocity mass norm that governs Picard convergence, but added on top of every raw pressure-lift extremum (measured ≈ +9% at Re 49.5, growing with seed error). The physical lift harmonics decay fast (j ≥ 2 are ~5 orders below the fundamental at Re 49.5), so K = 10 is generous.

_promoted_ys(R, ρ, η, Ω) → Dict{Int, Vector{ComplexF64}}

Response of the promoted coordinates on the orbit, per harmonic: y_s solving (isΩ·I − A₀) y_s = b_s, with b_s the drive of the promoted rows at harmonic s and A₀ the s = 0 part of ∂R_k/∂y_j.

Split out of _harmonic_R1 so the branch, the base-flow shift and the activity check all read the SAME closure. Three copies of the slaving algorithm drifting apart is what moved this into src in the first place; a fourth would undo that.

_real_roots(p; rtol) → Vector{Float64}

Real roots of Σ p[k+1] x^k, by companion-matrix eigenvalues. No external dependency and no bracketing, so a root cannot be missed between grid points the way a sign-change scan misses double roots.

_realify_self_conjugate(φ, ψ, r, λ; rtol = 1e-6) -> (φ_real, ψ_real)

Rotate a self-paired mode's eigenvectors onto the real axis.

A mode with σ[r] = r is its own conjugate, which asserts modes[:, r] = conj(modes[:, r]) — the eigenvector must be real-valued, not merely have a real eigenvalue. An eigensolver returns e^{iθ}·(real vector) for an arbitrary θ, and a biorthogonal gauge adds more phase, so the assertion generally fails on arrival even though the mode is perfectly good.

The phase is removed by rotating against the largest component (largest, so the reference entry is well away from a node of the mode). What remains must be real to round-off; if it is not, the eigenvalue is not truly real or the mode is not simple, and that is reported rather than absorbed — a defective mode used as a master coordinate corrupts the whole reduction silently.

_split_free_solution(s_full, fom) -> (u0_free, p0_free)

Extract velocity and pressure components from the full solution vector, restricted to free DOFs.

MORFE.accumulate_qp!(Fe, ∇W_args::NTuple{2}, mult, element, q, dΩ, t::FluidConvection)

Accumulate convective integrand at one quadrature point:

Fe[k] += mult · φᵢ · (−½)[∇u₂·u₁ + ∇u₁·u₂] · dΩ   for velocity DOF k

∇u ⋅ v computes (v·∇u) via Tensor{2,2} × Vec{2} → Vec{2} contraction: (∇u ⋅ v)[i] = Σj ∂j ui · vj = (v·∇u)_i which is exactly the (v·∇)u material-derivative direction.

amplitude_series(R, η, N) → Vector{Float64}

Coefficients of G(ρ) = Re(R₁(ρ, ρ, η)) / ρ = Σ gₙ uⁿ, u = ρ², for the order-N truncation.

The ρ-direction sibling of eta_series_report: G's zero IS the limit cycle, so these are the coefficients whose divergence caps the branch. Closed form off R₁'s z₁^{n+1} z̄₁^n η^c family — those are the only monomials that survive on the orbit at harmonic 1, since z₁^a z̄₁^b sits at harmonic a − b.

Promoted coordinates are NOT included: this is the y-free part of R₁. Use rom_po_residual when the y feedback matters; use this when the object of interest is the series itself.

assemble_K_visc(fom) -> (K_visc_free, K_visc_rect)

Assemble the raw viscosity stiffness matrix K_visc (full DOF space):

K_visc_kl = ∫ 2 ε(φ^k) : ε(φ^l) dΩ    (velocity test and trial, P2)

Pressure DOF rows and columns are zero. Two restrictions are returned:

Kviscfree — square free×free block, for the parametric coupling g₁ acting on the perturbation (which vanishes on prescribed DOFs); Kviscrect — rectangular free×ALL block, for the base-flow forcing h₀ = K·u₀, where u₀ is nonzero on the prescribed inlet DOFs (Poiseuille).

MORFE.assemble_element!(accum, Fe, element, t::FluidConvection)

Scatter element residual into the global free-DOF accumulator. Pressure-DOF entries of Fe are zero so only velocity-DOF contributions are added.

assemble_linear_operators(s0_full, fom; Re0) -> (B0_free, B1_free, A_lin_free)

Assemble: B1free — velocity mass matrix restricted to free DOFs (singular) Alinfree — linearised NSE operator restricted to free DOFs B0free = −Alinfree

All three are sparse matrices of size nfree × nfree.

s0_full is the full-DOF steady-state solution (from solve_steady_state).

assemble_steady_nse!(K, R, s_full, fom, Re0)

Fill K (tangent) and R (residual) from the current full-DOF solution s_full. Both K and R must be pre-allocated (e.g. via allocate_matrix(dh) and zeros).

assemble_velocity_mass_full(fom) -> SparseMatrixCSC

Assemble Mvel = ∫Ω φᵢ·φⱼ dΩ for velocity shape functions only (full DOF space; pressure rows/columns stay zero).

backbone_direction(W) → NamedTuple

Does the manifold's high-degree content settle onto ONE direction? Needs only W — no eigenbasis, no B₁, no solve.

If the series Σ Wₙ ρⁿ is limited by a single singularity, its coefficient VECTORS align with that singularity's direction as n → ∞, so cos(Wₙ, W_top) → 1. Promotion works by making that direction a coordinate, and it can only work if the direction exists: a cos that stalls well below 1 means the high-degree content is spread over several directions and no single promoted mode captures it.

This is the measurement modal_growth cannot make. Ranking modes by their amplitude at the top degree answers "which mode is largest there", which is not the same question and gave a misleading answer here — it picked λ = −5.135420 on a degree-9 amplitude of 3.50e-8 against 2.23e-8 and 1.32e-8 for its neighbours, i.e. no dominance at all, and that mode's pairing α/α_Hopf = 8.41e-3 makes it a poor coordinate for unrelated reasons.

Measured on the un-promoted order-9 run, cos(Wₙ, W_top) climbs monotonically — 0.27, 0.46, 0.72, 0.87 along z₁^{n+1}z̄₁^n and 0.48, 0.78, 0.92 along the mean-flow z₁^n z̄₁^n. The direction is settling, and the mean-flow family settles faster, which is consistent with the fold at ρ_c ≈ 1.81 being a mean-flow-distortion effect. Five backbone terms is few, so read the trend, not the last digit.

branch_amplitude_scale(R, η; lo, hi, n) → Float64

Characteristic limit-cycle amplitude at parameter η, found without assuming one.

ρ has no natural size. It is a master coordinate, so its scale is whatever the eigenvector gauge gives it: under SymmetricBiorthogonal the Kármán orbit happens to sit at ρ ~ O(1), under LeftBiorthogonal the SAME physical orbit sits some 400× further out, because that gauge leaves φ at its natural length instead of dividing by √α. Anything carrying units of ρ — a continuation step, a finite-difference increment, an initial guess — is therefore meaningless as a bare number, and hardcoding one silently pins the code to a gauge. That is not hypothetical: continuation constants tuned for ρ ~ O(1) traced 5001 steps to ρ = 0.002 in the left gauge and reported it as a collapsed branch.

The fix is to measure the scale first and work in ρ/ρ_ref. Bracketing the first sign change of G(ρ) = Re(R₁)/ρ on a LOG grid spanning lo…hi is scale-free by construction: below the orbit G ≈ σ > 0, above it the saturating term wins. Returns NaN if no sign change is bracketed, and callers should fall back to 1.0 rather than propagate it.

branch_by_amplitude(R, N; re0, rho_max, n_rho, re_min, re_max, method) → rows

The limit-cycle branch, parametrised by AMPLITUDE instead of by Reynolds number.

G(ρ, η) = 0 is one equation in two unknowns; which one you solve for is a free choice, and the conventional choice is the bad one here. Continuing in Re makes ρ the unknown, so the solver works in the direction whose series is DIVERGENT over most of the range (ρ ≈ 2.2 at Re 54 and ≈ 3.5 at Re 70, against ρ_conv ≈ 1.95). Fixing ρ instead makes η the unknown, and at fixed ρ, G is exactly a POLYNOMIAL in η — solved by companion matrix, no iteration — in the direction that converges comfortably (η(70) is 58 % of its radius, last term 1.9 % at order 9).

Everything the arclength continuation needed disappears with it: no ds/dsmax tuning, no ρ_ref scaling, no Newton tolerance to match a finite-difference Jacobian, no runaway (order 7 once reached ρ = 234), and no fold-chasing — a fold in ρ(Re) is just a monotone Re(ρ). Measured at order 9 against the DNS points: ρ = 0.65 → Re 49.50, 1.13 → 50.50, 1.64 → 52.00, 2.17 → 53.83, against DNS 49.5 / 50.5 / 52.0 / 54.0.

The ρ direction is still resummed — each η-coefficient p_c(ρ) = Σ_n Re(c_{n+1,n,c}) u^n is a series in u = ρ² and gets resum'd before the η-polynomial is assembled.

Branch selection marches in ρ and takes the η root nearest the previous one, seeded from the Hopf. That is continuation, but in the variable that cannot fold.

build_imex_operators(B₀, B₁, K_visc, η_prime, T, Δt; θ = 0.5, γ = 0.0)
-> (L_klu, RHS_M, Δt_exact, n_steps)

Assemble and KLU-factorise the θ-method operators for one orbit of period T, with the step count rounded so that n_steps · Δt_exact == T exactly. γ > 0 adds artificial damping γ B₁ to the implicit operator, shifting every eigenvalue left by γ — used for transient suppression during spin-up.

build_model(m::AssembledFluidModel, eig; master, outer = Int[],
			n_monomials = 1, conjugate_permutation = nothing)
	-> (; model, spectral, meta)

Turn an assembled fluid case and a solved eigenproblem into the first-order model and spectral data a reduction consumes.

It makes no choice about modes. master and outer are index vectors into eig.eigenvalues, supplied by the caller — this method does not select, filter, rank or detect anything. Which modes span the manifold and which are off-manifold targets is a modelling decision, and it belongs in the driver where it can be read and changed, not behind a default here.

eig is whatever solve_hopf_eigenproblem returned — eigenvalues and right modes for the whole computed spectrum. Left eigenvectors are computed here, for the master modes only, because each one costs its own adjoint factorisation (see left_eigenvector) and outer modes need eigenvalues alone.

scale and normalisation set the mode gauge. They change what W and R mean — coefficients from two gauges are not comparable — so they are keywords rather than constants, and are recorded in meta.

ORD = 1. NthOrderModel((B₀, B₁), …) has two linear terms, so there are no companion derivative blocks and SpectralData takes plain FOM × ROM matrices. Nothing needs reconciling against a higher model order.

The Reynolds continuation variable η′ = 1/Re − 1/Re₀ is a single frozen external state, so N_EXT = 1 — odd. η′ is real and therefore its own conjugate, a pairing the usual adjacent-pairs formula cannot express, so the permutation is built with full_conjugate_permutation.

conjugate_permutation overrides the derived master-block pairing. Give it when the master set is not a plain sequence of adjacent conjugate pairs — for instance when real (self-conjugate) modes have been included alongside a Hopf pair.

n_monomials sizes FluidConvection's batched-column buffer. It is a buffer capacity, not a reduction concept: build_model builds no MultiindexSet and no ResonanceSet. The buffer is allocated (n_monomials, n_qp) and indexed directly, so the value must be at least the monomial count the reduction will use.

expansion_order is the convenient spelling of the same thing, matching StructuralSVK's build_model: give the total-degree truncation and the count for the ROM + 1 reduced variables is derived here. Pass exactly one of the two — passing neither would silently size the buffer for a single column.

Build the augmented ORD=2 MORFE model from verified coefficient arrays.

check_linearisation(s0_full, fom, B0_free; Re0, ε = 1e-5)

Finite-difference check: compare Alin s ≈ [R(s₀ + ε·ei) − R(s₀)] / ε for a few random free DOF directions. Prints max relative error.

close_under_conjugation(eig; rtol = 1e-4)
	-> (; eigenvalues, right_modes, hopf_index, conjugate_index)

Complete a computed spectrum with the conjugates the eigensolve could not return, and report where the Hopf pair ended up.

solve_hopf_eigenproblem shifts at a complex σ, so ARPACK returns only the modes near σ. A strongly oscillatory mode's conjugate sits near σ̄ and is simply never computed — at Re₀ = 49.03 the Kármán mode λ = 0.004 + 16.859i comes back while λ̄ = 0.004 − 16.859i does not. conjugate_index as returned by the eigensolve is then the nearest available eigenvalue rather than the true partner, and passing that pair as master makes build_model throw.

The missing halves are synthesised rather than solved for. B₀ and B₁ are real, so (λ̄, φ̄) is an eigenpair exactly whenever (λ, φ) is — this is an identity, not an approximation, and it costs nothing. It is also the only way to get the partner with the phase conjugate_permutation asserts: an independent solve would pin it only up to a scalar.

Three cases per computed mode:

• |Im λ| ≤ rtol·|λ| — the mode is real, hence its own conjugate. Nothing is added; giving it a partner would duplicate the coordinate.

• the conjugate is already in the set — nothing is added.

• otherwise conj(λ) and conj.(φ) are appended.

rtol is relative. ARPACK's numerical zero is ~1e-7 of a mode's magnitude, not machine epsilon, so an absolute threshold reads a real mode as complex and then demands a conjugate that does not exist.

hopf_index is carried through unchanged — the closure only appends — and conjugate_index is recomputed afterwards, so it now names the true partner.

compute_drag_lift(s_full, fom; Re0) -> (Cd, Cl)

Compute drag and lift coefficients by integrating the stress tensor over the cylinder boundary. Reference values (Turek–Schäfer benchmark, Re = 20): Cd ≈ 5.57, Cl ≈ 0.011.

compute_pressure_lift_weights(fom) -> Vector{Float64}

Assemble l ∈ ℝ^N such that F_L^pres = l ⋅ u_full = ∫_Γ_cyl (-p n_y) dΓ.

Sign convention matches compute_drag_lift: n is the outward normal from the fluid, so the pressure traction is -p·n and l[i] = -∫ n_y ψ_i^p dΓ for each pressure DOF i on the cylinder boundary (zero elsewhere).

Note: this captures the PRESSURE contribution only. The viscous shear traction (2ν ε(u))·n that compute_drag_lift integrates is deliberately omitted from the lift polynomial L(z); at Re ≈ 50 it contributes a few percent of the total lift, so L(z) slightly underestimates the physical lift amplitude.

compute_transformed_drag_lift(s_full, fom, geometry; μ, Re)

Full pressure-plus-viscous traction on the reference obstacle. Nanson's formula is applied before integration, n dΓ = adj(F)' N dΓ_ref; the velocity gradient uses the same composition pullback as the volume forms.

domain_area(fom) -> Float64

Quadrature measure of the fluid domain |Ω|.

domb_sykes(W) → NamedTuple

Locate and classify the singularity that limits the amplitude expansion.

For a series Σ aₙ ρⁿ whose nearest singularity is at ρ_c and behaves like (1 − ρ/ρ_c)^{-γ}, the coefficient ratios obey asymptotically

rₙ = aₙ / aₙ₋₁ ≈ (1/ρ_c) · (1 + (γ − 1)/n)

so plotting rₙ against 1/n gives a straight line whose intercept is 1/ρ_c and whose slope is (γ − 1)/ρ_c. That separates the two cases that matter here:

· γ = 1 (slope ≈ 0) — a simple POLE. This is what "the radius is set by the nearest pole, which signals an outer mode" predicts, and it is the case where carrying that mode as a coordinate should push the singularity out. · γ ∉ ℤ (slope ≠ 0) — a branch point, which is what a quadratic convolution generically produces and which no single mode is responsible for.

Only the backbone z₁^{n+1}z̄₁^n is used, so this is the ρ direction. Five degrees (1…9) is few for an extrapolation — the fit residual is returned so the reader can judge, and a two-point Richardson estimate is given alongside the least-squares one.

eta_series_report(R; re0, re_max, orders, gate_radius) → NamedTuple

Convergence of the PARAMETER expansion, read off R directly. Run this before believing any plot.

manifold_ratio_test measures W, which is coordinate-dependent and therefore not comparable between a promoted and an un-promoted run. R's coefficient FAMILIES are comparable: for every run of the same problem, R₁'s z₁η^c is the eigenvalue series and z₁²z̄₁η^c is the Landau series, whatever else was promoted. Two numbers per family:

· radius — the ratio test |a_c| / |a_{c+1}|, taken at the largest c where both are nonzero, i.e. where the sequence has settled. · truncation ratio — |a_c η^c| / |a_0| at η(re_max) with c the highest η power the order-N truncation retains. This is how much the LAST kept term still contributes at the top of the sweep, and it is the honest statement of where the Re range stops being trustworthy. It is not a promoted-run concern: un-promoted, η(70) = 6.1e-3 against a radius of 1.05e-2 is 58% of the way out, so the tail is not negligible there either.

This table is what distinguishes the three ways promotion has failed here, and it costs a deserialise and a dictionary lookup:

· inert — promoted rows carry no forcing at all, R₁ identical to the un-promoted run. · divergent — the promoted rows' pure-η forcing has its own radius (1.5e-4 measured), far inside the Hopf block's, and it drags R₁ down with it. · contaminated — radii look normal but a low-degree coefficient is wrong by orders of magnitude and every higher one inherits the offset, which is why families[2].coeffs[1] (the η-independent Landau coefficient) is reported separately: it is the one entry that must agree with the un-promoted run to full precision, and it did even in the run whose R₁ reached 1e36.

eval_perturbation_convection!(accum, s_free, fom)

Assemble the quadratic perturbation convection into accum (length nfreedpim):

accum[k] += −∫_Ω  φᵢ · (∇u · u)  dΩ     (velocity DOFs k only)

where u is the velocity part of s_free extracted element-by-element via fom.free_to_local_dpim. accum is zeroed on entry.

export_lift_csv(out, data_dir, mset, L0, L_coeffs)

Write L_coefficients.csv (constant base-flow row + lift polynomial). R_coefficients.csv is written by MORFE.save_rom in the driver.

export_lift_polynomial(out, data_dir, W, l_free, L0) -> L_coeffs

Project the lift weight vector onto the parametrisation, Lα = lᵀWα (bilinear — the adjoint would conjugate), serialise lift_polynomial.jls, and return the coefficient vector for the CSV export.

export_reduced_dynamics(out, data_dir, R, master_eigenvalues; re0, ord, nvar)

Write reduced_dynamics.txt with the first row ż₁ = R₁(z₁, z̄₁, η′) in complex form (equation 2 is the conjugate, equation 3 the parameter η̇′ = 0) and echo the nonzero monomials to out.

export_vtk_bundle(data_dir, fom, s0_full, master_eigenvalues, all_eigenvalues, all_modes)

Serialise a plain-array bundle (mesh, DOF maps, base flow, eigenmodes — no Ferrite types) for external ParaView/VTU export tooling. For a direct quadratic VTU export, load WriteVTK and call write_paraview_p2p1 with the same FOM.

find_periodic_orbit(s_guess, η_prime, T, fom, B₀, B₁, K_visc, h₀_vec; kwargs...)
→ (s_po, T_fom, E_norm, X_max, n_orbits, converged)

Find the FOM periodic orbit near the ROM prediction (s_guess, T) by Picard iteration (integrate-and-map) with an undamped spin-up and a lift-based period probe.

Keyword arguments

• Δt : time-step size (default 5e-4)

• θ : implicit weight (0.5 = Crank-Nicolson, 1 = backward Euler)

• n_spinup : undamped spin-up orbits before the probe (default 3)

• γ_spinup : OPTIONAL artificial damping −γ B₁ s during spin-up, for

			  seeds that overshoot badly (default 0.0).  WARNING: any
			  γ ≫ σ (the Hopf growth rate, ~0.1 s⁻¹ near onset) kills the
			  shedding oscillation itself within a couple of orbits and
			  collapses the trajectory onto the damped system's forced
			  equilibrium — leave at 0 unless you know the seed is bad.

• max_picard : maximum undamped Picard orbits before giving up (default 400)

• tol : convergence tolerance on ‖E‖_M / X_max (default 1e-4)

• lift_weights: if provided (Vector{Float64} of length nfreedpim), enables

			  the autocorrelation period probe on
			  `F_L(t) = dot(lift_weights, s)`; keeps the ROM period T when
			  the gated estimate is rejected.

• verbose : print residual every orbit if true

Returns

• s_po : periodic-orbit initial condition

• T_fom : period at convergence (probe estimate + tangent refinement)

• E_norm : ‖Φ(s_po,T)−s_po‖_{B₁} at convergence

• X_max : max_t ‖s(t)‖_{B₁} over the last orbit

• μ̂ : dominant transversal Floquet multiplier estimated from the

			ratio of consecutive Picard residuals at fixed T (NaN when no
			valid pair occurred, e.g. first-orbit convergence). The anchor
			residual under-reports the distance to the cycle by 1/(1−μ̂):
			one orbit maps a transversal deviation d to μ̂·d, so the
			measured closure is only (1−μ̂)·d.

• n_orbits : total number of periods integrated (spin-up + probe + Picard)

• converged : true if tolerance was met before max_picard was reached

fluid_model(meshfile; Re, newton_tol = 1e-10, newton_max_iter = 30,
			s_init = nothing, verbose = true) -> AssembledFluidModel

Assemble the incompressible Navier-Stokes problem about its steady base flow at Reynolds number Re.

Four stages, each already implemented in this module:

1. setup_fem — P2/P1 Taylor-Hood spaces, boundary conditions, DOF maps.

2. solve_steady_state — Newton solve for the base flow.

3. assemble_linear_operators — the linearised B₀ and the singular B₁.

4. assemble_K_visc — the viscous operator, restricted two ways.

Then the Reynolds-continuation scaling, which is the part that used to live in the example driver:

K_visc .*= -D
h₀_vec = -D .* (K_visc_rect * s₀_full)

D is the cylinder diameter — the reference length in ν = D/Re, so the external coordinate is η′ = 1/Re − 1/Re₀ and these two carry its coefficient. The rectangular block is required for h₀ because the base flow is nonzero on the prescribed inlet DOFs (Poiseuille), which the square free × free block drops.

meshfile is a .msh read by FerriteGmsh. Mesh generation stays with the caller: it is problem geometry, not physics. The geometry constants this module assumes (Turek–Schäfer channel, Ø 0.1 m cylinder) must match whatever produced meshfile — see fem_setup.jl.

fold_overlap(W, Φ, λ; master) → NamedTuple

WHICH outer modes span the direction the manifold is straining in. Use this to choose what to promote.

backbone_direction establishes THAT the high-degree content settles onto a single direction; this decomposes that direction, as the cosine |⟨φ_k, W_d⟩| / (‖φ_k‖ ‖W_d‖) over the top backbone degrees d.

Right eigenvectors only, and that is the point. modal_growth projects with ψ, and in this descriptor pencil the modes carrying the fastest-growing manifold content are precisely the ones whose pairing α = ψᵀB₁φ is degenerate (~1e-13) — left_eigenvector warns on every one of them that its left and right vectors are not the same mode. Ranking on that measures the adjoint solve, not the manifold. No adjoint enters here, so those modes cannot corrupt the answer, and it costs nothing extra because Φ is already in hand.

Measured on the un-promoted Kármán order-9 run, the degree-9 backbone is led by λ = −6.542148 + 18.415390i at 0.184 and −11.812255 + 19.322862i at 0.115, with every real mode below 0.038 — including all three that had been promoted on modal_growth's advice. Promoting a mode nearly orthogonal to this direction is a change of coordinates that changes nothing, which is exactly what was measured: ρ_conv 1.95 either way and the order-9 fold unmoved at Re 55.1.

homological_denominators(λ, master; max_ord, tol, pairings) → Vector{NamedTuple}

For every monomial z₁^a z̄₁^b up to max_ord, the eigenvalue combination μ = a·λ₁ + b·λ̄₁ is what the homological solve inverts against on the outer block: (L − μ B₁). A small |μ − λ_k| is a NEAR-RESONANCE — the outer mode k is nearly excited by that monomial and its manifold component is amplified by 1/|μ − λ_k|.

Returns one row per outer mode: its smallest denominator over the whole monomial set, the monomial achieving it, and (when pairings is supplied) the α ratio that says whether the mode could be promoted out of the outer block at all.

A denominator is only a CANDIDATE ranking — it is the amplification factor, not the amplitude. A large denominator with a large numerator beats a small one with a negligible numerator. Use modal_growth to find which is actually excited.

integrate_orbit!(s, n_steps, η_prime, fom, L_klu, RHS_M, h₀_vec, f2, rhs;
                 on_step = nothing) -> Bool

Advance s in place by n_steps IMEX steps:

(B₁/Δt + θ A) s⁺ = (B₁/Δt − (1−θ) A) s + f₂(s,s) + η′ h₀_vec

f2 and rhs are pre-allocated work vectors (length of s). If given, on_step(step, s) is called after every accepted step (lift recording, state storage, norm tracking). Returns false on blow-up (non-finite state).

left_eigenvector(A_lin, B_mass, λ, φ; normalisation, scale)
	-> (φ_gauged, ψ_gauged, α)

The left eigenvector for ONE mode, by adjoint shift-invert at that mode's own λ, biorthonormalised against its right partner φ.

One factorisation per mode, and there is no cheaper correct way. A single adjoint solve at the common shift σ returns its own subset of the spectrum in its own order; matching those to the right modes by eigenvalue pairs most of them with the wrong partner, which shows up as ψᵀBφ ≈ 0 — the pairing is degenerate and the mode is unusable as a master. Shifting at λ makes Aᵀ − λBᵀ singular exactly there, so ARPACK returns the partner that belongs to φ.

For the same reason the shift is λ and not conj(λ): the latter finds the conjugate mode's left vector instead, and the pairing collapses again.

Returns the gauged pair and the raw bilinear pairing α = ψᵀBφ before scaling, so a caller can see how well-conditioned the mode was.

lift_functional(m::AssembledFluidModel) -> (l_free, L0)

The lift functional restricted to the DPIM free DOFs, and the base flow's own lift.

l_free is the weight vector such that l_freeᵀ · u′ is the perturbation lift; L0 = l_freeᵀ · s₀ is the steady contribution the perturbation adds to. Together they give L(t) = L0 + l_freeᵀ·u′(t).

Pressure traction only (−p·n_y on the cylinder), matching compute_pressure_lift_weights — the viscous shear contribution is deliberately omitted there, so it is omitted here too.

lift_polynomial(W, l_free) -> (L_coeffs, mset)

Project the lift weights onto the parametrisation: L_α = lᵀ W_α, one complex coefficient per monomial, plus the monomial set they are indexed by.

The product is bilinear, not sesquilinear — transpose, not adjoint. W's coefficients are complex because the reduced coordinates are, but l_free is a real functional of the flow; conjugating it would silently give the lift of the conjugate manifold.

Evaluating the resulting polynomial at (z, z̄, η′) reproduces the perturbation lift without touching a FOM-sized vector, which is what lets the Python post-processing work from CSV alone.

The projection itself is MORFE.observable_polynomial; this wrapper only names it for the flow and keeps the (coefficients, mset) pair the exporters here expect.

load_rom_poly(csv) → ROMPoly

Read an R_coefficients.csv (exp_1…exp_NVAR, R1_re, R1_im, R2_re, …). NVAR comes from the header, so this reads a promoted and an un-promoted run alike.

make_base_forcing(h₀_vec_free) -> MultilinearMap{1}

Create a MultilinearMap{1} for the base-flow parametric forcing

h₀(η′) = −D·η′ · K_raw[free, ALL] · s₀

with multiindex=(0,) and multiplicity_external=1 (no state input, one η′ input).

h₀vecfree = −D · Kraw[free, ALL] · s₀full is a fixed vector encoding the viscous forcing direction when Re departs from Re₀: the rectangular block acts on the FULL base-flow vector, whose prescribed inlet DOFs carry the Poiseuille profile. It determines the external direction Φ_ext via the cohomological equation for monomial (0,0,1).

MORFE passes the external argument as a basis direction — the unit SVector{1,Int}([1]) here, since this external system is diagonal — so f!(accum, rext) adds rext[1] · h₀_vec to accum (η′ scaling is handled by the polynomial factorisation framework).

Do not assume the argument is integer-valued: a re-based external system (one whose linear matrix was not upper triangular) is passed the columns of its change of basis, which are complex in general. The eltype(accum) conversion below is what keeps this term correct in that case.

make_param_coupling(K_visc_free) -> MultilinearMap{1}

Create a MultilinearMap{1} for the parametric viscous coupling

g₁(s, η′) = η′ · K_visc · u′        (main.jl passes K_visc = −D·K_raw)

with multiindex=(1,) and multiplicity_external=1 (one state input, one η′ input).

Kviscfree is the viscosity stiffness restricted to free DOFs, pre-scaled by −CYLD in main.jl so that the term equals −D·η′·Kraw·u′ = −(ν−ν₀)·Kraw·u′. Its pressure rows and columns are zero, so the multiplication correctly targets only velocity DOFs in both input and output.

MORFE passes the external argument as a basis direction — the unit SVector{1,Int}([1]) here, since this external system is diagonal. The actual η′ scaling is handled by the polynomial factorisation in the cohomological solver. Hence f!(accum, s, rext) adds rext[1] · Kvisc · s = Kvisc · s to accum.

Do not assume the argument is integer-valued: when an external system's linear matrix is not upper triangular, ExternalSystem re-bases it and MORFE then passes the columns of the change of basis, which are complex in general. The eltype(accum) conversion below is what keeps this term correct in that case.

manifold_ratio_test(W) → NamedTuple

Ratio test on ‖W‖ per degree, along the amplitude backbone z₁^{n+1}z̄₁^n and along η′^n. A ratio that SETTLES to a constant is the signature of a genuine analytic singularity, and its reciprocal gives the radius; a ratio that keeps growing would instead indicate round-off contamination.

⚠ Valid WITHIN one run only. Promoted coordinates each carry a 1/√α scaling, and α spans orders of magnitude across modes, so ‖W‖ along the backbone means a different thing in each coordinate system. Compare runs on observables (lift vs DNS, the branch), never on this.

measure_orbit(s_po, η_prime, T, fom, ops, l_free, M_vel, vel_rows, area)
-> (; lift_max, lift_min, lift_mean, lift_ringing, FL, FL_per, Δt_meas,
	  avg_tke, periodicity_max, lift_amp_drift, μ_meas, ok)

Integrate TWO periods from the converged periodic-orbit point s_po and evaluate the observables plus a definition-level periodicity verification.

Period 1 stores the velocity history and the lift signal; observables come from it: lift statistics from the harmonic projection FL_per (see _periodic_lift; lift_ringing is the peak raw-minus-periodic residual — the numerical pollution amplitude), and the TKE subtracts the PERIOD MEAN (fluctuation TKE — matches tkefromgram's ζ̄ subtraction; velocity is constraint-consistent, so it needs no filtering).

Period 2 verifies periodicity BY ITS DEFINITION — the state at t must equal the state at t+T for ALL t, not just at the Picard anchor (with a finite anchor residual the deviation propagates and can transiently grow within a period; the flow is non-normal):

• periodicity_max = max_t ‖v(t+T) − v(t)‖_{M_vel} / max_t ‖v(t)‖_{M_vel} — the uniform-in-time periodicity error (velocity mass norm; pressure is excluded — its θ = ½ Nyquist ringing is handled by the harmonic projection).

• lift_amp_drift — relative change of the harmonic-projected lift amplitude between the two periods (observable stationarity).

• μ_meas = ‖v(2T)−v(T)‖ / ‖v(T)−v(0)‖ — the dominant transversal Floquet multiplier measured from the two anchor mismatches (deviation δ contracts to μδ per period, so consecutive closures scale by μ).

modal_growth(W, B₁, ψ, λ_outer) → Vector{NamedTuple}

Decompose the manifold's amplitude backbone into outer modes and rank them by how fast each grows with degree. ψ holds the left eigenvectors of the outer modes (columns), so the modal amplitude of mode k at backbone degree d is c_k = ψ_kᵀ B₁ W[:, m_d].

⚠ Do not select promotion candidates with this — use fold_overlap. The ranking needs ψ, and here the top eight rows are all modes whose pairing is degenerate at ~1e-13, where left_eigenvector warns that the left and right vectors are not the same mode. Acting on it sent three separate runs to real modes whose actual overlap with the fold direction is below 0.038, and all three changed nothing. Kept because the growth ratio is still the honest answer to a different question — how fast a given mode's component grows with degree, given a trustworthy ψ.

pade(c, L, M) → Function

[L/M] Padé approximant of Σ c[k+1] u^k, returned as a callable in u.

Requires L + M + 1 coefficients. Returns nothing when the Toeplitz system is singular — that happens for a genuinely degenerate series and must not be papered over, because a silently-wrong approximant is far worse than a missing one.

parametric_model(case::AssembledFluidModel, geometry; ...)

Assemble the fixed-domain composition pullback. geometry returns (I, ∇ψ) at each quadrature point. The external ordering is always μ first, then ξ; both are frozen and real.

prepare_energy_gram(fom) -> (M_vel, vel_rows, area)

Order-independent part: velocity mass restricted to the velocity rows of the free_dpim DOF set, the row indices of velocity DOFs within the free vector, and |Ω|. Call ONCE per FEM setup.

promoted_amplitude(R, ρ, η) → Float64

Largest promoted-coordinate response on the orbit at (ρ, η), max_s max_k |y_{k,s}|.

The activity check. A promoted run whose branch has this at 0 is a change of coordinates that changed nothing: {y = 0} is invariant, R₁ never sees a y-monomial, and the ROM is the un-promoted one in disguise. That is exactly what a normal-form tolerance tight enough to keep η out of R_k produces, and it went unnoticed for two runs because every physical quantity agreed — of course it did, it was the same ROM.

promoted_equilibrium(R, η) → Vector{ComplexF64}

Equilibrium of the promoted coordinates on the trivial (z₁ = 0) branch: y*(η) solving R_k(0, 0, y, η) = 0, which on that slice is the linear system A₀ y = −b₀.

y* is not zero, and requiring it to be is what broke two redesigns. Once φk is a master direction, the base flow at Re ≠ Re₀ has a component along it, and y*(η) is exactly that component expressed in the new coordinate — the base-flow shift, not an error. The s = 0 harmonic of `promoted_ysalready solves for it, so the branch and the observables are taken abouty*(η)` and always were.

What is worth checking is its SIZE. A shift comparable to the orbit amplitude means the promoted coordinate has stopped representing the base flow, and in practice means the pure-η forcing has diverged — which eta_series_report measures directly.

resum(c; method) → Function

Callable summation of Σ c[k+1] u^k. method:

· :taylor — plain partial sum. Reproduces the pre-resummation behaviour exactly, so it is the control every comparison should carry. · :pade — diagonal (or near-diagonal) Padé from all available coefficients.

Falls back to :taylor rather than returning nothing, so a caller never silently loses a curve; the fallback is visible because the resummed and Taylor curves then coincide. Fewer than three coefficients is such a case — there is nothing to fit.

rom_hopf_eta(R; ε, η0, tol, max_iter) → Float64

η′ at which the linearised ROM growth rate vanishes — the true Hopf point. Re₀ is only the FOM's expansion point, not necessarily the critical Reynolds number, so this root is generally nonzero.

rom_invariants(p; η_step) → NamedTuple

· σ, ω — the Hopf eigenvalue λ = σ + iω at η′ = 0. Fully invariant. · c101 — ∂λ/∂η′. Fully invariant. Computed by a central difference, because in a PROMOTED run the mean-flow coordinates respond to η′ directly and feed back through z₁·yk, so the [1,0,1] coefficient alone would miss part of it. · `c210eff— the effective Landau coefficient after slaving. Scales as |c|² under z → cz, so comparable only at a fixedMODESCALE. ·c210ratio— Im/Re of c₂₁₀. The |c|² cancels, so this is the gauge-free fingerprint — and the quantity whose SIGN the conjugate-pairing bug flipped. ·criticality` — sign(Re c₂₁₀).

rom_palc_step(ρ, η, τ, Δs, R; tol, max_iter)
→ (ρ, η, T, τ, n_iter, converged)

One pseudo-arclength step: predictor Δs along τ, then a 2×2 Newton corrector on [F(ρ,η); τ·(p − last) − Δs].

Unit tangent to the branch F(ρ, η′) = 0, oriented consistently with τ_prev.

rom_po_R1(R, ρ, η; Ω0, tol, maxit) → (R₁, Ω, converged)

Fundamental of R₁ and the orbit frequency, solved together. Ω enters the harmonic closure and is itself read off R₁, so the two are found by fixed point — seeded from the linear Hopf frequency, and reported as non-converged rather than silently accepted.

Angular frequency Ω = Im(R₁)/ρ at the periodic orbit.

Periodic-orbit residual F(ρ, η′) = Re(R₁); vanishes on the limit-cycle branch.

roots_at_re(R, re; re0, rho_lo, rho_hi, n) → Vector{Float64}

Every ρ > 0 with F(ρ, η(re)) = 0, by sign change on a log grid then bisection. Used to continue a branch where the arclength corrector cannot, and to cross-check a traced one.

MORFE.scatter_qp!(∇W_col, W_free, element, t::FluidConvection)

Extract velocity DOFs from the free-DOF state vector W_free, compute velocity value and gradient at each QP, and store as FluidVelQP in ∇W_col.

Pressure DOFs in W_free are ignored — only velocity DOFs feed into the convective nonlinearity.

seed_from_branch(branch, re_target) -> (ρ, T) or nothing

Linear interpolation of the branch amplitude ρ and period T at re_target (columns eta,Re,rho,omega,T; rows ordered along the branch).

setupfem(meshfileorgrid; obstacletag = "Cylinder", referencelength = 0.1, quadratureorder = 6, channelheight = 0.41, boundaryconditions = nothing) -> NamedTuple

Load meshfile, or use an already-loaded two-dimensional Ferrite grid, and build all FEM objects for the P2/P1 Taylor-Hood cylinder-flow problem. The grid overload permits topology-preserving geometry continuation without serialising one mesh file per parameter value.

First component of the reduced dynamics on the orbit, promoted coordinates closed.

solve_hopf_eigenproblem(A_lin, B_mass; nev, sigma_re, sigma_im,
						target_freq = nothing,
						normalisation = SymmetricBiorthogonal(),
						scale = 1.0, tol = 0.0, maxiter = 3000,
						ncv = nothing, close_conjugates = true,
						conjugate_rtol = 1e-4, verbose = true)
solve_hopf_eigenproblem(B::Tuple; kwargs...)
	-> (; eigenvalues, right_modes, hopf_index, conjugate_index)

Compute eigenvalues of A_lin y = λ B_mass y by shift-invert ARPACK, close the result under conjugation, and point at the Hopf pair. The second form takes an assembled model's B = (B₀, B₁) and applies the descriptor system's own sign; prefer it.

sigma_re offsets the shift from the imaginary axis; sigma_im targets a frequency band. Neither affects which mode is selected — only the factorisation. tol, maxiter, and ncv are forwarded to ARPACK. Their defaults preserve the historical solver configuration.

The Hopf mode is the eigenvalue with the smallest |Re λ| among those with Im λ > 0. That heuristic is reliable near Re_c, where the shedding mode IS the least damped; away from it another oscillatory mode can sit closer to the imaginary axis and be picked silently, so pass target_freq (rad/s) to pin the frequency instead.

The two gauge choices, both explicit

• normalisation — see AbstractModeNormalisation. Defaults to the historical SymmetricBiorthogonal.

• scale — a further uniform factor applied to both sides, purely for conditioning. The Kármán case uses 1e-2, which makes φᴴBψ = 1e-4 rather than 1. SpectralData deliberately has no scale field so that such a tweak stays visible where it is made; this keyword is that visibility.

The spectrum comes back closed under conjugation

The shift σ is complex, so ARPACK returns only the modes near σ — a strongly oscillatory mode's conjugate sits near σ̄ and is never computed. The result is therefore passed through close_under_conjugation before it is returned, which appends the missing halves analytically (exact, because A_lin and B_mass are real) and makes conjugate_index name the true partner of hopf_index.

That is on by default because the raw conjugate_index is a footgun: it is an argmin over what ARPACK happened to return, so it names the nearest available mode, and handing that pair to build_model as master throws. Two consequences worth stating:

• nev is a lower bound on length(eigenvalues), not the count. Closure appends.

• Appended conjugates go at the end, not next to their partners, so the result is not a sequence of adjacent pairs.

close_conjugates = false returns exactly what ARPACK produced. Pass it when a caller must post-process the raw modes before closing — for instance one that phase-aligns a tracked eigenvector and needs the synthesised partner to inherit that phase. conjugate_rtol is forwarded as close_under_conjugation's rtol; it is relative because ARPACK's numerical zero is ~1e-7 of a mode's magnitude, not machine epsilon.

Returns a NamedTuple so the fields are named at every call site: eigenvalues is a Vector{ComplexF64} in ARPACK order with any synthesised conjugates appended, right_modes is the matching n × length(eigenvalues) matrix, and hopf_index / conjugate_index locate the Hopf pair within them.

solve_hopf_eigenproblem(B::Tuple; kwargs...)

Solve the eigenproblem of an assembled model's linear operators, B = (B₀, B₁) as AssembledFluidModel stores them.

The descriptor system is B₁ẋ + B₀x = F, so its eigenproblem is -B₀y = λB₁y: the sign belongs to the equation, not to the caller. Spelling it out at every call site (solve_hopf_eigenproblem(-case.B[1], case.B[2]; ...)) is an easy sign to get wrong and gives no error when it is, only a spectrum reflected about the imaginary axis, so this method takes the tuple whole. Every keyword is forwarded unchanged.

solve_steady_state(fom; Re0, tol = 1e-10, max_iter = 30, s_init = nothing)
-> (u0_free, p0_free, s_full)

Newton–Raphson solver for the steady NSE at Reynolds number Re0.

Returns: u0free — free-DOF velocity perturbation base flow (length nufree) p0free — free-DOF pressure base flow (length npfree) s0_full — full-DOF solution vector (including prescribed BCs)

The initial guess applies the selected inhomogeneous velocity boundary policy with zero interior values, or uses s_init (a full-DOF vector) when given — used to warm-start a continuation in Re, e.g. following the (unstable) steady branch above Re_c. Convergence is monitored by the free-DOF residual ℓ²-norm.

sweep_branch_in_re(R, ρ_start, re_start; re0, re_max, dre, rho_max)
→ Vector{NTuple{6,Float64}}

Continue a branch by stepping Re and solving F(ρ, η) = 0 for ρ at each step, following the root nearest the previous one.

This exists because a turning point in ρ is NOT a turning point in Re. Pseudo-arclength continuation parametrises by arclength in (ρ, η) and has to negotiate the fold; where its corrector cannot — the 6-mode order-7 branch oscillated across one turning point 73 times without resolving it — the branch is still perfectly single-valued as ρ(Re) and a plain sweep walks straight through. A root scan confirms the solutions are there: the 6-mode order-9 branch runs ρ = 1.53 at Re 52 to 9.04 at Re 69.5 while PALC gave up at Re 52.6.

A jump limit rejects hops onto a disconnected root: this ROM's truncated polynomial also carries spurious high-amplitude roots (ρ ≈ 9–17), which are not the physical branch.

tke_series(G, A, η, N; consistent = false) → Vector{Float64}

Coefficients of the period-averaged fluctuation TKE, T(ρ) = Σ t_k u^k, u = ρ².

Closed form — no orbit sampling. On z₁ = ρe^{iθ} the monomial z₁^a z̄₁^b η^c sits at harmonic s = a − b, so:

· a = b is constant on the orbit and is removed by the period-mean subtraction (it IS the mean-flow distortion — correctly excluded from a FLUCTUATION energy); · a product ζ_m ζ_n survives the period average only when s_m + s_n = 0.

so t_k = ½ Σ_{s_m + s_n = 0} Re(G_mn η^{c_m + c_n}) collected at k = (a_m+b_m+a_n+b_n)/2. Equivalent to figures.py::tke_from_state with NS → ∞, and agrees with it to round-off; having it as a series is what allows the sum to be resummed instead of merely evaluated.

G and A are the Gram and exponent table written by write_energy_gram. N truncates on the CORE degree, matching truncate_dynamics.

consistent decides which coefficients you get, and the right choice depends on what you do with them.

· false (default) — every product of two kept monomials, so the series runs to ρ^{2N}. This is what figures.py evaluates and it is the more ACCURATE Taylor sum: ‖u'_N‖² has error O(ρ^{N+2}) against O(ρ^{N+1}) for a degree-N truncation. But its coefficients above total degree N are INCOMPLETE — they are products of known terms missing the contributions of the unknown ones. · true — only pairs with core(m) + core(n) ≤ N, every coefficient exact.

Resummation must use consistent = true. Padé fits the coefficient SEQUENCE, so feeding it the incomplete tail makes it model an artefact: at Re 54 the full sequence gives −106 % against DNS while the exact one gives −0.0 %. A Taylor sum merely adds the tail on and is barely harmed; a rational fit is led by it.

trace_limit_cycle_branch(R; re0, re_max, re_min, rho_max, ds0, max_steps)
→ Vector{NTuple{6,Float64}}   # (η, Re, ρ, Ω, T, fold)

PALC continuation of F(ρ, η′) = 0 from the Hopf point. Returns one row per branch point, with fold counting how many times the branch has turned in Re — downstream code segments the sheets on that counter instead of re-deriving it from Re decreasing, which mis-segments whenever consecutive rows repeat.

Three defects fixed relative to the original solve_rom.jl:trace_branch:

· it stopped at re < re_c - 0.5, abandoning any branch that folds and would come back up — folds are exactly what PALC exists to traverse; · a stalled step was still pushed as a row, so a collapsing Δs produced 82 identical rows at one point (order-7 of the 6-mode run) before the arclength floor hit; · nothing bounded ρ, so past a fold it would chase a disconnected root out to ρ ≈ 17.

truncate_dynamics(R, N) → ReducedDynamics

Zero every monomial whose CORE degree — (z₁, z̄₁, η′), excluding promoted coordinates — exceeds N. Because the cohomological solve is graded this IS the order-N reduced dynamics, bit-exact.

Truncating on the TOTAL degree instead would drop the mean-flow coupling z₁·y_k at low N and silently return a ROM with no saturation mechanism — the Landau coefficient would flip sign. Promoted coordinates are first order by construction and are not part of the order hierarchy.

velocity_dof_mask(fom) -> BitVector

Mark every global DOF that belongs to the velocity field :u.

write_energy_gram(data_dir, W, M_vel, vel_rows, area)

Per-order part: G = (Wvelᵀ Mvel Wvel)/|Ω| and the exponent table, written as tkegramre.csv / tkegramim.csv / tkeavector.csv into data_dir.