Add rational surface impedance boundary condition#770
Conversation
Add a SurfaceRationalImpedanceOperator implementing a Robin boundary whose per-square surface impedance is a user-specified rational function of frequency, Zs(s) = N(s)/D(s) with s = iω. Numerator and denominator polynomial coefficients (zeros and poles) are given as lists in the new boundaries["RationalImpedance"] config block, allowing any passive lumped-network response. Being a general function of frequency, the boundary contributes iω/Zs(iω) to the frequency-dependent system matrix A2(ω) (like SurfaceConductivityOperator), and is available for the driven, eigenmode, and boundary-mode problem types. Coefficients are nondimensionalized as a_k/(Z0*tc^k) and b_k/tc^k. Wiring spans configfile (RationalImpedanceData and parser), iodata (nondimensionalization, problem-type warnings, nonlinear-eigenmode sparsity-pattern handling), spaceoperator (assembly, auxiliary-space marker, multiple-BC check), the JSON schema, and the docs. Add a verification example (examples/rational_impedance) comparing, on a parallel-plate TEM line, a parallel-RLC sheet realized as a RationalImpedance against the same RLC realized as a LumpedPort; the two S11 responses agree to solver tolerance.
…undary
Validate user-specified rational surface impedances and warn on physically
questionable inputs, following the Palace Mpi::Warning convention:
- Warn when the numerator/denominator degree difference exceeds one, since a
positive-real (passive) impedance has |deg(N) - deg(D)| <= 1.
- Warn, once per boundary, when Re{Zs(iw)} < 0 at an evaluated frequency, with
a relative tolerance so lossless reactive terminations do not trigger it.
Real-valued coefficients (enforced by the JSON schema and parser) guarantee a
Hermitian Zs(iw), and hence a real time-domain response.
Add the RationalImpedance boundary attributes to the aggregated boundary attribute list built in BoundaryData. Without this, external rational-impedance boundaries raised a spurious "no associated boundary condition, PMC assumed" warning, and internal rational-impedance sheets were skipped by the interface-element (mesh cracking) logic, so the boundary term was not applied on them.
Document the boundaries["RationalImpedance"] configuration block in the config
reference: the Numerator/Denominator polynomial coefficients of the per-square
impedance Zs(s) = N(s)/D(s) with s = iw, the real-coefficient Hermitian
guarantee, the passivity warnings (degree excess and Re{Zs} < 0), and the
frequency-domain-only restriction (driven, eigenmode, boundary mode).
laylagi
left a comment
laylagi
left a comment
There was a problem hiding this comment.
Thanks for putting this together. This is a useful addition, and the driven RLC comparison is a nice way to validate the basic behavior. I left a few comments, mostly around edge cases and solver integration.
| (dbc_marker[i] || farfield_marker[i] || surf_sigma_marker[i] || | ||
| surf_z_Rs_marker[i] || surf_z_Ls_marker[i] || lumped_port_Rs_marker[i] || | ||
| lumped_port_Ls_marker[i] || wave_port_marker[i] || floquet_port_marker[i]); | ||
| surf_z_Rs_marker[i] || surf_z_Rs_marker[i] || surf_z_Ls_marker[i] || |
There was a problem hiding this comment.
| surf_z_Rs_marker[i] || surf_z_Rs_marker[i] || surf_z_Ls_marker[i] || | |
| surf_z_Rs_marker[i] || surf_rz_marker[i] || surf_z_Ls_marker[i] || |
Was this intentional? surf_z_Rs_marker appears twice, while surf_rz_marker is defined above and used later for mutual-exclusion checking. Should one of those entries be surf_rz_marker[i]?
| { | ||
| for (auto attr : data.attributes) | ||
| { | ||
| MFEM_VERIFY(!impedance_marker[attr - 1], |
There was a problem hiding this comment.
Could we reorder this validation so we check that attr is in range and exists on the mesh before indexing impedance_marker[attr - 1]? For example:
if (attr <= 0 || attr > bdr_attr_max || !bdr_attr_marker[attr - 1])
{
bdr_warn_list.insert(attr);
continue;
}
MFEM_VERIFY(!impedance_marker[attr - 1],
"Multiple definitions of rational impedance boundary properties for "
"boundary attribute " << attr << "!");
impedance_marker[attr - 1] = 1;The duplicate-definition check is useful, but it should only run for valid attributes. Otherwise invalid input like boundary attribute 0 or an attribute larger than bdr_attr_max can index out of bounds before the warning/ignore path runs.
I noticed a couple of existing surface-BC operators have the same ordering issue, so those probably deserve a follow-up cleanup, but I don't think we should copy that pattern into new code.
| const std::complex<double> Z = N / D; | ||
| // Passivity necessary condition on the imaginary axis: Re{Zs(iω)} >= 0. Warn once per | ||
| // boundary (relative tolerance avoids false alarms on lossless reactive terminations). | ||
| if (!bdr.warned_passivity && Z.real() < -1.0e-9 * std::abs(Z)) | ||
| { | ||
| const double f_ghz = omega * freq_scale / (2.0 * M_PI); | ||
| Mpi::Warning("Rational impedance boundary (attribute {:d}) is not passive at " | ||
| "f = {:.4f} GHz: Re(Zs) = {:.3e} < 0!\n", | ||
| bdr.attr_list.Size() ? bdr.attr_list[0] : -1, f_ghz, Z.real()); | ||
| bdr.warned_passivity = true; | ||
| } | ||
| for (auto attr : bdr.attr_list) | ||
| { | ||
| const double sc = bdr.attr_scaling.at(attr); | ||
| const std::complex<double> coef = s / (Z * sc); // iω·Ys per square |
There was a problem hiding this comment.
Could we compute the coefficient directly as s * D / (N * sc) instead of forming Z = N / D and then computing s / Z? The two forms are algebraically equivalent when N(s) and D(s) are nonzero, but the direct admittance form is more robust. At an impedance pole where D(s) = 0 and N(s) != 0, Zs = N / D is infinite, while the Robin contribution s / Zs = s * D / N should go to zero. Forming Z first can introduce inf/nan intermediates unnecessarily.
| // A rational surface impedance contributes only to the frequency-dependent system matrix | ||
| // A2(ω); it has no time-domain (transient) or static realization here. | ||
| MFEM_VERIFY(impedance.empty() || problem_type == ProblemType::DRIVEN || | ||
| problem_type == ProblemType::EIGENMODE || |
There was a problem hiding this comment.
Could you clarify and validate eigenmode support for RationalImpedance?
As implemented, this BC contributes through A2(omega), so eigenmode solves trigger the nonlinear eigenmode path. The nonlinear paths I checked evaluate A2 using abs(Im(lambda)), rather than the full complex eigenvalue.
I tried adapting the PR’s parallel RLC example to an eigenmode config using the same equivalent rational coefficients:
Numerator = [2.2588181359893346e-06, 0.0]
Denominator = [9.542690323697736e-20, 5.995849163282243e-09, 376.730313668]
This case should reduce to s / Zs(s) = 1/Ls + (1/Rs)s + Cs s^2, but the rational version still triggered nonlinear refinement through A2(omega) and did not converge reliably.
Could cases where s / Zs(s) reduces to a quadratic polynomial be assembled into K/C/M instead of A2(omega)? For genuinely nonquadratic rational impedances, I think we need an eigenmode validation test or documented limitations.
There was a problem hiding this comment.
Heads up that an upcoming change affects this. Once #778 is merged, A2 will be evaluated using the full complex eigenvalue,
| // A2(ω); it has no time-domain (transient) or static realization here. | ||
| MFEM_VERIFY(impedance.empty() || problem_type == ProblemType::DRIVEN || | ||
| problem_type == ProblemType::EIGENMODE || | ||
| problem_type == ProblemType::BOUNDARYMODE, |
There was a problem hiding this comment.
Could you double-check whether ProblemType::BOUNDARYMODE should be allowed here? I may be missing it, but I don’t see SurfaceRationalImpedanceOperator constructed or used BoundaryModeOperator.
If rational impedance is not applied in boundary-mode solves yet, maybe we should remove ProblemType::BOUNDARYMODE from this allowed list for now, or add the missing wiring if support is intended.
3 participants
This PR adds a
RationalImpedanceboundary condition: a surface (Robin) impedance boundary whose per-square impedance is an arbitrary rational function of frequency,Zs(s) = N(s)/D(s)withs = iω, given by numerator and denominator polynomial coefficients. It generalizes the parallel-RLCImpedanceboundary to any passive lumped-network response. It closes #553 and #642.