Use complex frequency in nonlinear BC in eigenmode simulations#778
Draft
simlapointe wants to merge 16 commits into
Draft
Use complex frequency in nonlinear BC in eigenmode simulations#778simlapointe wants to merge 16 commits into
simlapointe wants to merge 16 commits into
Conversation
Add WavePortOperator::AddBoundaryMassBdrCoefficients and GetWavePortKn
helpers, plus SpaceOperator::GetWavePortBoundaryMassMatrix, exposing
the per-port μ⁻¹ boundary mass M^(p) and the scalar kₙ,p(ω) separately.
The existing AddExtraSystemBdrCoefficients now calls the new helper
internally, so behaviour and assembled operators are unchanged.
These helpers are foundational plumbing for upcoming changes to the
RomOperator: storing per-port projected M^(p)_r once and assembling
the wave-port contribution online as Σ_p kₙ,p(ω)·M^(p)_r recovers an
O(n²) online phase for adaptive sweeps with wave ports, and unblocks
circuit synthesis with wave ports.
Adds a unit test that compares Im{A2(ω)·v} to Σ_p kₙ,p(ω)·M^(p)·v on
the existing CPW wave-port example at three frequencies in the
configured sweep band, in serial and under MPI.
(cherry picked from commit 9248d8f)
The reduced-order driven sweep currently re-assembles the full HDM
A2(ω) and projects it onto the basis at every online frequency, which
defeats the asymptotic cost model of the reduced model: each online
solve was paying for one HDM-scale boundary form assembly plus
dim_V HDM-size mat-vecs.
Wave ports are the dominant A2 contributor in practice. Their
contribution factors as A_wp(ω) = i·Σ_p kₙ,p(ω)·M^(p)_{μ⁻¹}, with
kₙ,p the only ω-dependent factor (a scalar from the per-port
cross-section EVP). Store the per-port HDM operators in RomOperator,
project them once when the basis grows, and apply the wave-port term
online as Σ_p kₙ,p(ω)·M^(p)_r. Online cost per ω drops from
HDM-scale to O(n²).
To keep correctness for any non-wave-port A2 contributors
(second-order farfield, surface conductivity), add an
include_wave_ports flag to GetExtraSystemMatrix and call the
"excluding wave ports" path on the slow fallback. The default
behaviour is unchanged for all existing call sites.
Validated against the parent branch on cpw_wave_adaptive: port-S
agrees to ~1e-11 relative, domain-E to ~1e-13.
(cherry picked from commit 9938828)
(cherry picked from commit b5cdbd08bcf90ea3ed24c7c4da6c6e4eca35fab5)
Make the wave-port (and the co-located 2nd-order farfield ABC and surface conductivity) boundary conditions accept a genuinely complex frequency in eigenmode simulations, so the NLEPS HYBRID and SLP nonlinear eigensolvers evaluate A2(λ) as the exact analytic continuation at the complex eigenvalue (ω = -i·λ) rather than the real-axis projection A2(|Im λ|). Cross-section modal EVP (modeeigensolver / boundarymodeoperator): - AssembleAtt/AssembleAnn/Solve/AssembleFrequencyDependent take a complex ω with a `complex_omega = (Im ω != 0)` gate; for real ω the assembly reduces bit-for-bit to the previous real-ω form. The multigrid preconditioner is built at Re(ω) (the full complex operator is carried by the monolithic opA). - WavePortData::SolveKnExact(complex ω) solves the cross-section EVP at the genuinely complex frequency and returns the exact complex kₙ, side-effect-free (does not disturb the cached real-ω modal state). Boundary operators: - WavePortOperator, FarfieldBoundaryOperator and SurfaceConductivityOperator gain complex-ω AddExtraSystemBdrCoefficients overloads. The shared per-term formula (boundary mass / curl-curl / skin-depth surface impedance) is factored out so the real and complex overloads stay in sync and the real overload is unchanged. The driven / boundary-mode wave-port path keeps stamping only the real (propagating) part of kₙ; only the eigenmode path uses the full complex i·kₙ·M (line attenuation -Im(kₙ)·M on the real slot). - SpaceOperator: complex GetExtraSystemMatrix(λ) dispatcher and a complex BC frequency a3 in GetPreconditionerMatrix / AssemblePreconditioner so the preconditioner matches the exact complex system matrix. Nonlinear eigensolvers (nleps / slepc): - The QuasiNewtonSolver and SlepcNEPSolver A2 callback is now complex-only (SetExtraSystemMatrix takes std::function<...(complex)>); the production system matrix, residual, and finite-difference Jacobian all evaluate A2(λ) at the complex eigenvalue (holomorphic FD: perturb λ directly). The HYBRID seed pencil (NewtonInterpolationOperator) keeps a separate real-ω A2 for the initial guess. Whether a wave port is real or complex is decided purely by solver type (eigenmode → complex, driven/boundary-mode → real); there is no config knob. Excludes the NLEIGS split-form solver and all kₙ analytic-fit / continuation machinery (FitKnSq, polynomial pencils, FN objects); those land with NLEIGS. Adds a unit test (test-waveportimpedance) checking SolveKnExact at a complex ω against the closed-form TE10 dispersion √(εμω²−k_c²), including the near-cutoff decaying branch. The A2/kₙ·M_p factorisation invariant test is unchanged and still passes (driven path stays real-kₙ).
The HYBRID/SLP linear seed stage (the has_A2 probe, the
NewtonInterpolationOperator polynomial pencil, and the seed PEP
linear-system matrix) was assembled from the real-ω
GetExtraSystemMatrix<ComplexOperator>(double) path, which stamps only
Re(kₙ)·M. The Newton stage and the matching preconditioner (funcP at
a3 = ω, routed through the complex BC overload) use the full complex
i·kₙ·M, which for a lossy cross-section also carries -Im(kₙ)·M on the
real slot. The seed operator and its preconditioner therefore differed
by -Im(kₙ)·M, with two consequences on cpw_wave_eigen:
- the (otherwise exact, MGMaxLevels=1 + ComplexCoarseSolve) seed
preconditioner no longer matched the system matrix, so each GMRES
solve took 2 iterations instead of 1 (88 iters / 44 solves), and
- the seed eigenvalues were built from the lossless-port operator,
shifting them by ~1e-4 relative from the complex-operator result.
Make funcA2 a real-ω-signature shim over funcA2_complex, evaluating the
exact complex operator at λ = i·ω. The seed probe, pencil, and linear
system matrix now use the same wave-port stamping as the preconditioner
and the Newton stage. On cpw_wave_eigen this restores 44 iters / 44
solves (exact preconditioner) and matches the nleps-waveport eigenvalues
to machine precision.
Regenerate the cpw/wave_eigen, adapter/hybrid, and cpw/lumped_eigen
reference outputs from the complex-waveport build, where the eigenmode
nonlinear solve evaluates the wave-port / 2nd-order ABC boundary
conditions at the genuinely complex eigenvalue (ω = -i·λ) rather than
the real-axis projection ω = |Im λ|. Eigenvalues shift at the O(1/Q)
level relative to the old real-BC baseline.
Config changes:
- cpw/cpw_wave_eigen.json: N = 1 -> 2 (compute two modes).
- adapter/hybrid.json: Target 6.6 -> 8.0 GHz (move further from the
wave-port cutoff).
- cpw/cpw_lumped_eigen.json: unchanged; its reference moves only
because of the complex-ABC physics (the sole frequency-dependent
A2 term here is the 2nd-order absorbing BC).
All three converge with backward error ~1e-11 and exact (single-GMRES-
iteration) preconditioning.
The SolveKnExact-vs-TE10-dispersion checks used loose 1e-2–3e-2 tolerances. The measured error against the closed-form √(εμω²−k_c²) is ~1.5e-6 relative away from cutoff and ~9e-6 near cutoff (order-2 FEM discretization of the cross-section mode on the 8x8x4 tet mesh; eig/ksp tol are 1e-10 so solver error is negligible). Tighten all six checks to 1e-4, which keeps ~10–65x headroom above the discretization floor while actually catching regressions — in particular a real-axis-projection regression that would miss the imaginary part by O(Im ω/Re ω) ~ 9% at the Q≈10 probe point. Also set wave.max_size (mirroring IoData::CheckConfiguration) so the ARPACK cross-section eigensolve has a large enough Krylov subspace.
The HYBRID seed pencil's Interpolate() already works entirely in the complex λ plane (complex sample points, coefficients, divided differences), but its funcA2 callback was typed std::function<...(double)> and the k=0 sample extracted points[j].imag() before calling it. On the complex-waveport eigenmode path this forced a real-ω-signature shim in eigensolver.cpp that immediately reconstructed the complex point (funcA2(ω) -> funcA2_complex(i·ω)) — a pure round-trip and a needless second callback. Type NewtonInterpolationOperator's funcA2 as std::function<...(std::complex<double>)> and sample at points[j] directly. eigensolver.cpp now has a single complex A2 callback (funcA2) feeding the has_A2 probe, the seed pencil, the seed PEP system matrix, and the Newton/NEP refinement, removing the funcA2/funcA2_complex duplication. The .imag() extraction was also lossy for any off-imaginary-axis sample point; passing the complex node through is strictly more correct. No-op numerically: the removed round-trip was the identity on the imaginary axis. cpw_wave_eigen eigenvalues match the reference to ~1e-11 (iterative-solve nondeterminism) and the preconditioner stays exact (26 solves = 26 iterations).
The "Exact" suffix was meaningful on the nleps-waveport branch, where it distinguished the exact cross-section solve from the polynomial kₙ² fit. There is no fit on this branch, so "exact as opposed to what?" is confusing. Rename to SolveKnComplex / GetWavePortKnComplex, which names the actual distinguishing feature: a complex-frequency solve returning complex kₙ, pairing naturally with the real-ω GetWavePortKn (which truncates to Re(kₙ)). Comments reworded to match; pure rename, no behavior change.
ModeEigenSolver::BuildSystemMatrixA segfaulted (inside hypre's MergeDiagAndOffd, via HypreParMatrixFromBlocks) whenever the imaginary block system needed a zero-diagonal placeholder — triggered e.g. by a domain Conductivity with no LossTan, where the transverse imaginary block (Atti) is non-null but the normal imaginary block (Anni) is null. Loss-tangent did not crash because it populates Anni, so no placeholder is built. Root cause: the placeholder was built with the 4-arg HypreParMatrix(comm, glob, row_starts, &diag) constructor, which does NOT deep-copy diag — it aliases the SparseMatrix's CSR arrays (CopyCSR with mem_owner=false). The backing Vector and SparseMatrix were declared inside the `if` blocks, so they were destroyed at the block's closing brace, before HypreParMatrixFromBlocks read the placeholder's (now freed) CSR. In an optimized build the freed stack memory is reused, producing a hard segfault. Hoist the backing Vector and SparseMatrix (as unique_ptr<SparseMatrix>) to the enclosing scope so they outlive both the placeholder HypreParMatrix and the HypreParMatrixFromBlocks call. No behavioral change otherwise; the offd block is a valid empty block and is correct in parallel. Verified: adapter/hybrid.json with Conductivity=1e-6 (no LossTan), which previously segfaulted, now converges at np=1 and np=2 with matching eigenvalues and backward error ~1e-11. Pre-existing on main (not caused by the complex-waveport work); the fix applies cleanly there too.
simlapointe
changed the title
simlapointe
added
the
no-long-tests
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
In eigenmode simulations, the (complex) eigenmode of the 3D system is an input to the nonlinear frequency-dependent boundary conditions (wave port, 2nd order absorbing BC, surface conductivity). Previously, these boundary conditions took only a real frequency as input. This PR changes that to allow the freq-dependent BCs to take a complex frequency as input so that they are consistent with the eigenmode simulation. Results of eigenmode simulations with wave ports (or other nonlinear BCs) change by roughly 1/Q.
This PR is currently based on top of
simlapointe/waveport-A2-split. We should merge #772 and rebase before merging this PR.