Functions
np.ndarray	_frozen_A1 (np.ndarray C, float flame_gain, float tau, np.ndarray g, np.ndarray b, complex lam)
Tuple[np.ndarray, np.ndarray]	_companion_solve (np.ndarray A0, np.ndarray A1, np.ndarray A2)
Tuple[complex, np.ndarray, int]	fixed_point_solve (np.ndarray K, np.ndarray C, np.ndarray M2, float flame_gain, float tau, np.ndarray g, np.ndarray b, complex lambda_seed, *, float tol=1e-9, int max_iter=50)
np.ndarray	residual_N (np.ndarray K, np.ndarray C, np.ndarray M2, float flame_gain, float tau, np.ndarray g, np.ndarray b, complex lam, np.ndarray p)
Tuple[complex, np.ndarray, int]	newton_polish (np.ndarray K, np.ndarray C, np.ndarray M2, float flame_gain, float tau, np.ndarray g, np.ndarray b, complex lambda0, np.ndarray p0, *, float tol=1e-12, int max_iter=30)

Detailed Description

Solvers for the active-flame nonlinear eigenvalue problem (NLEVP), verification
case V3. This is the first place ``lambda`` enters nonlinearly (through the flame
delay ``exp(-lambda*tau)``, Eq. 7/11), so it is the first place the paper's Section IV
solver hierarchy is actually needed rather than the plain linear GEVP of ``passive.py``.

The discrete operator being solved (Eq. 11, specialized to this package's cases: a
single azimuthal wavenumber m=0, a constant boundary-admittance matrix C from a compact
nozzle, and a single compact flame reference, so F(lambda) is rank 1):

    N(lambda) p = [ K + lambda*C + lambda^2*M2 - lambda*flame_gain*exp(-lambda*tau)*g@b.T ] p = 0

``K``, ``C``, ``M2`` come from ``acoustics/assembly.py``. ``g`` and ``b`` are TWO
DIFFERENT vectors, not the same one reused (see the long comment in
``acoustics/assembly.py`` above ``disk_load_vector`` for why conflating them is a real,
previously-caught bug): ``b`` (``assembly.point_sampling_vector``) is a plain point
sample of the pressure at the flame's reference location (Eq. 7a); ``g``
(``assembly.disk_load_vector``) is the r-weighted energy-injection weight, carrying the
same ``r dOmega`` measure as every other matrix here. ``flame_gain`` and ``tau`` are the
lumped Crocco interaction strength and time lag (see ``validation/v3_ntau_duct.py`` for
how they relate to the physical n, tau of Eq. 2.1 and to this module's own dispersion
relation).

Scale note: everything here is DENSE (``numpy``/``scipy.linalg``, not sparse/Krylov).
That is a deliberate scope choice for P0's verification meshes (a few hundred nodes):
the sparse/shift-invert machinery of Section IV.E is already exercised (for the
LINEAR-in-lambda case) by ``passive.py``; this module's job is to validate the
NONLINEAR delay handling, a separate concern, at a problem size where a dense solve is
simpler to write and just as fast. A production-scale (P1+) NLEVP solver would still
need the sparse Krylov kernel underneath the same two algorithms below.

Function Documentation

◆ _companion_solve()

Tuple[np.ndarray, np.ndarray] engine.stability_hifi.eigen.nlevp._companion_solve	(	np.ndarray	A0,
		np.ndarray	A1,
		np.ndarray	A2 )

protected

All eigenpairs of the quadratic pencil ``A0 + lambda*A1 + lambda^2*A2 = 0`` via
the standard companion linearization to a size-2N generalized eigenvalue problem
(Algorithm 1, step 5 — "companion-linearize to GEVP of size 2N_h").

Writing ``z = [p; lambda*p]``, one checks by direct substitution that

    [[0, I], [-A0, -A1]] z = lambda * [[I, 0], [0, A2]] z

reproduces exactly ``A0 p + lambda A1 p + lambda^2 A2 p = 0`` (the top block row
gives the trivial identity ``lambda*p = lambda*p``; the bottom block row gives the
QEP once that identity is substituted in). Solved directly and densely
(``scipy.linalg.eig``) since P0's verification meshes are small — see module
docstring for why that is an adequate, deliberate scope choice here.

Definition at line 52 of file nlevp.py.

Referenced by fixed_point_solve().

◆ _frozen_A1()

np.ndarray engine.stability_hifi.eigen.nlevp._frozen_A1	(	np.ndarray	C,
		float	flame_gain,
		float	tau,
		np.ndarray	g,
		np.ndarray	b,
		complex	lam )

protected

The coefficient of ``lambda`` in N(lambda), with the flame's delay factor
``exp(-lambda*tau)`` FROZEN at the current iterate ``lam`` (Algorithm 1, step 3).

N(lambda) = K + lambda*[C - flame_gain*exp(-lam*tau)*g@b.T] + lambda^2*M2 once
frozen — a genuine QUADRATIC eigenvalue problem in lambda, since freezing the delay
factor at a fixed number removes the only source of nonlinearity.

Definition at line 40 of file nlevp.py.

Referenced by fixed_point_solve().

◆ fixed_point_solve()

Tuple[complex, np.ndarray, int] engine.stability_hifi.eigen.nlevp.fixed_point_solve	(	np.ndarray	K,
		np.ndarray	C,
		np.ndarray	M2,
		float	flame_gain,
		float	tau,
		np.ndarray	g,
		np.ndarray	b,
		complex	lambda_seed,
		*	,
		float	tol = 1e-9,
		int	max_iter = 50 )

Algorithm 1 (Section IV.B): frozen-delay fixed point, specialized to this
package's single-flame-reference, dense-small-scale case.

At each outer iterate, freeze ``exp(-lambda*tau)`` at the current guess, solve the
resulting *linear* (in the sense of ordinary QEP, not NLEVP) quadratic eigenvalue
problem exactly via companion linearization, and take whichever of its 2N
eigenvalues is nearest the current guess as the next iterate. Converges linearly;
the paper reports 3-8 outer iterations for realistic parameters, matched here.

Returns ``(lambda, p, n_iterations)``.

Definition at line 77 of file nlevp.py.

References _companion_solve(), _frozen_A1(), and max.

◆ newton_polish()

Tuple[complex, np.ndarray, int] engine.stability_hifi.eigen.nlevp.newton_polish	(	np.ndarray	K,
		np.ndarray	C,
		np.ndarray	M2,
		float	flame_gain,
		float	tau,
		np.ndarray	g,
		np.ndarray	b,
		complex	lambda0,
		np.ndarray	p0,
		*	,
		float	tol = 1e-12,
		int	max_iter = 30 )

Bordered Newton polish (Section IV.C), solving the *true* NLEVP (delay left in,
not frozen) to quadratic convergence from the fixed-point solver's output.

The bordered system (paper's own equations, Section IV.C):

    [[N(lambda),      N'(lambda) p],   [dp     ]     [ N(lambda) p    ]
     [c^H,             0          ]] @ [dlambda] = - [ c^H p - 1       ]

with the normalization row ``c^H p = 1`` fixing scale/phase (paper uses a generic
``c``; here ``c = e_k``, a unit vector at whichever component of ``p`` has the
largest magnitude — a simple, standard pivoting choice so the constraint row is
never close to singular).

N'(lambda) = C + 2*lambda*M2 - flame_gain*exp(-lambda*tau)*(1 - lambda*tau)*outer(g,b)
(derivative of the "-lambda*flame_gain*exp(-lambda*tau)" term in N(lambda) w.r.t.
lambda; C and M2 contribute their usual constant/2*lambda derivatives).

Definition at line 113 of file nlevp.py.

References max.

◆ residual_N()

np.ndarray engine.stability_hifi.eigen.nlevp.residual_N	(	np.ndarray	K,
		np.ndarray	C,
		np.ndarray	M2,
		float	flame_gain,
		float	tau,
		np.ndarray	g,
		np.ndarray	b,
		complex	lam,
		np.ndarray	p )

N(lambda) @ p — the raw NLEVP residual, used to check solver agreement/quality.

Definition at line 106 of file nlevp.py.