diff --git a/src/gsim/common/modes.py b/src/gsim/common/modes.py index 37ef8d28..c15af58c 100644 --- a/src/gsim/common/modes.py +++ b/src/gsim/common/modes.py @@ -292,13 +292,15 @@ def propagating_modes[ModeT]( min_index: float = 1.0, max_loss_ratio: float = 1.0, max_gain_ratio: float = MAX_GAIN_RATIO, + max_index: float | None = None, ) -> list[ModeT]: """Keep the Modes that propagate, dropping evanescent and spurious ones. A Mode is kept when it is guided (``Re(n_eff)`` above *min_index*), loses no more than *max_loss_ratio* per unit length than it advances in phase, and does not grow: ``Im(n_eff)`` positive beyond - *max_gain_ratio* is gain, which a passive line cannot have. + *max_gain_ratio* is gain, which a passive line cannot have. With + *max_index* set, a Mode whose ``|n_eff|`` exceeds it is dropped too. Args: modes: Every solved Mode. @@ -311,6 +313,13 @@ def propagating_modes[ModeT]( max_gain_ratio: Upper bound on ``+Im(n_eff) / Re(n_eff)``, the room left for a lossless Mode's numerical noise (default: :data:`MAX_GAIN_RATIO`). + max_index: Upper bound on ``|n_eff|``, or ``None`` for no bound + (default). A solve asked for many Modes around a lossy + conductor also returns Modes living inside the metal, with + ``|n_eff|`` near ``|sqrt(eps_metal)|``, some of them inside + the loss bound; the line Mode sits orders of magnitude below. + ``gsim.femwell.conductor_index_bound`` derives a bound from + the permittivity the solve was given. Returns: The propagating Modes, in the order they were solved. @@ -321,6 +330,7 @@ def propagating_modes[ModeT]( if ( n_eff.real > min_index and -max_loss_ratio * n_eff.real < n_eff.imag <= max_gain_ratio * n_eff.real + and (max_index is None or abs(n_eff) <= max_index) ): kept.append(mode) return kept @@ -340,6 +350,7 @@ def select_line_mode[ModeT]( rule: LineModeRule | None = None, min_index: float = 1.0, max_loss_ratio: float = 1.0, + max_index: float | None = None, degeneracy_rtol: float = 0.03, ) -> ModeT: """Select the physical line Mode from a set of solved Modes. @@ -365,6 +376,14 @@ def select_line_mode[ModeT]( min_index: Lower bound on ``Re(n_eff)`` for the default rule. max_loss_ratio: Upper bound on ``|Im(n_eff)| / Re(n_eff)`` for the default rule. + max_index: Upper bound on ``|n_eff|`` for the default rule, or + ``None`` for no bound (default). The default rule picks the + highest ``Re(n_eff)``, so a spurious Mode inside a lossy + conductor, which carries an index near + ``|sqrt(eps_metal)|``, wins over the line Mode when the solve + returns one; set this when solving more than a few Modes + around such a conductor (``gsim.femwell.conductor_index_bound`` + gives a bound from the solve's permittivity). degeneracy_rtol: Relative spread in ``Re(n_eff)`` within which two candidates count as ambiguous and a warning is issued. @@ -381,7 +400,12 @@ def select_line_mode[ModeT]( choice between them is not physically meaningful. """ candidates: Sequence[ModeT] = ( - propagating_modes(modes, min_index=min_index, max_loss_ratio=max_loss_ratio) + propagating_modes( + modes, + min_index=min_index, + max_loss_ratio=max_loss_ratio, + max_index=max_index, + ) if rule is None else list(rule(modes)) ) diff --git a/src/gsim/femwell/__init__.py b/src/gsim/femwell/__init__.py index 63376269..334dfc09 100644 --- a/src/gsim/femwell/__init__.py +++ b/src/gsim/femwell/__init__.py @@ -13,6 +13,7 @@ from gsim.femwell.adapter import ( boundary_field_ratio, + conductor_index_bound, electrode_current, elementwise_epsilon, epsilon_by_region, @@ -25,6 +26,7 @@ __all__ = [ "boundary_field_ratio", + "conductor_index_bound", "electrode_current", "elementwise_epsilon", "epsilon_by_region", diff --git a/src/gsim/femwell/adapter.py b/src/gsim/femwell/adapter.py index 4029fbd4..2b67bbf0 100644 --- a/src/gsim/femwell/adapter.py +++ b/src/gsim/femwell/adapter.py @@ -24,7 +24,9 @@ Marks-Williams power-current impedance it implies (:func:`z0_power_current`), and — for an RF line — the one reading of a selected Mode an RF line solve asks for: index, impedance and whether it is -the wall Mode (:func:`line_reading`). A conductor is named to every +the wall Mode (:func:`line_reading`). :func:`conductor_index_bound` reads, from the +permittivity a solve was given, the index above which a Mode lives inside a +conductor rather than on the line. A conductor is named to every current integral the same way, as a :class:`~gsim.common.modes.Conductor`. @@ -34,6 +36,7 @@ from __future__ import annotations +import warnings from pathlib import Path from typing import TYPE_CHECKING, Any, cast @@ -211,6 +214,93 @@ def elementwise_epsilon( return np.asarray(sampled, dtype=np.complex128) +def conductor_index_bound( + epsilon: dict[str, complex] | ArrayLike, + *, + fraction: float = 0.1, + margin: float = 10.0, +) -> float | None: + """Upper bound on ``|n_eff|`` that rejects Modes living inside a conductor. + + A mode solve asked for more than a few Modes around a lossy + conductor also returns Modes confined to the metal, with ``|n_eff|`` + close to ``|sqrt(eps_metal)|`` (about 2.6e4 for aluminium at 1 GHz). + Some lose less per unit length than they advance, so the loss bound + keeps them, and the slowest-Mode rule of + :func:`~gsim.common.modes.select_line_mode` picks one over the line + Mode. The line Mode's index is orders of magnitude lower: over a good + conductor ``n_line**2 / eps_metal`` is about + ``eps_d R' eps0 / (L' sigma)``, independent of frequency and of + order 1e-7 for an aluminium strip in oxide. Pass the bound returned + here as ``max_index`` to ``select_line_mode``. + + A region or element is a conductor when conduction dominates its + permittivity, ``|Im(eps)| > |Re(eps)|``. That alone does not make + it a safe source for the bound: a weak conductor, such as moderately + doped silicon, has ``|sqrt(eps)|`` close to the indices a line Mode + can reach. A conductor sets the bound only when + ``fraction * |sqrt(eps)|`` is at least *margin* times the largest + ``|sqrt(eps)|`` of the other regions (and above the light line); + the others are skipped with a warning. A slow-wave line Mode can sit + several times above the dielectric index (3.5 times on an aluminium + strip in oxide at 1 GHz), hence the default margin of 10. + + Args: + epsilon: The permittivity given to :func:`solve_modes`: either + ``{region_name: eps}`` or per-element values. + fraction: Bound as a fraction of a conductor's ``|sqrt(eps)|`` + (default 0.1; spurious metal Modes have been seen from 0.9 + times that value up). + margin: Smallest ratio between a conductor's bound and the + largest non-conductor index for the conductor to count + (default 10). + + Returns: + The smallest bound over the conductors that count, or ``None`` + when none does (no bound is then given from this source). + + Raises: + ValueError: If *fraction* or *margin* is not positive. + + Warns: + UserWarning: When a conductor is skipped because its bound would + not clear the dielectric index by *margin*. + """ + if fraction <= 0.0: + raise ValueError(f"fraction must be positive, got {fraction}.") + if margin <= 0.0: + raise ValueError(f"margin must be positive, got {margin}.") + values = ( + np.fromiter(epsilon.values(), dtype=np.complex128) + if isinstance(epsilon, dict) + else np.asarray(epsilon, dtype=np.complex128).ravel() + ) + is_conductor = np.abs(values.imag) > np.abs(values.real) + if not is_conductor.any(): + return None + dielectric = values[~is_conductor] + dielectric_index = ( + float(np.sqrt(np.abs(dielectric)).max()) if dielectric.size else 1.0 + ) + floor = max(margin * dielectric_index, 1.0) + bounds = fraction * np.sqrt(np.abs(values[is_conductor])) + usable = bounds >= floor + if not usable.all(): + skipped = np.unique(values[is_conductor][~usable]) + shown = ", ".join(f"{eps:.4g}" for eps in skipped[:3]) + more = f" and {skipped.size - 3} more" if skipped.size > 3 else "" + warnings.warn( + f"conductor_index_bound skipped conductor permittivity {shown}{more}: " + f"{fraction:g} * |sqrt(eps)| does not reach {margin:g} times the " + f"largest dielectric index ({dielectric_index:.4g}), so it cannot " + "tell a Mode inside it from the line Mode.", + stacklevel=2, + ) + if not usable.any(): + return None + return float(bounds[usable].min()) + + def solve_modes( msh_path: str | Path, *, @@ -252,6 +342,9 @@ def solve_modes( for RF materials with big conductive ``|Im(eps)|`` can put the shift far from the physical quasi-TEM mode; pass an explicit guess (e.g. the expected slow-wave index) there. + Asked for more than a few Modes around a lossy conductor, + the solve also returns Modes confined to the metal; select + the line Mode with ``max_index=conductor_index_bound(epsilon)``. Returns: The femwell ``Modes`` result (each mode carries ``n_eff``). diff --git a/tests/common/test_mode_selection.py b/tests/common/test_mode_selection.py index 99bdb01b..2488358e 100644 --- a/tests/common/test_mode_selection.py +++ b/tests/common/test_mode_selection.py @@ -130,6 +130,45 @@ def test_a_bound_nothing_survives_is_reported(self): select_line_mode([3.0 - 2.0j], max_loss_ratio=0.5) +class TestIndexBound: + """A Mode living inside a lossy conductor is not the line Mode (#498). + + Asked for 8 to 16 Modes on an aluminium strip in oxide, the solve also + returns Modes with ``|n_eff|`` near ``|sqrt(eps_metal)|`` (about 2.6e4 + at 1 GHz), some inside the loss bound. The slowest-Mode rule then picks + one of them over the line Mode, 5.13-4.68j on that geometry. + """ + + LINE = 5.13 - 4.68j + METAL = 23174.0 - 12017.0j # -Im/Re = 0.52, inside the loss bound + + def test_without_a_bound_the_metal_mode_wins(self): + """Unchanged default: the bound is the caller's to give.""" + assert select_line_mode([mode(self.LINE), mode(self.METAL)]).n_eff == ( + self.METAL + ) + + def test_the_bound_drops_the_metal_mode(self): + line = mode(self.LINE) + picked = select_line_mode([line, mode(self.METAL)], max_index=2.6e3) + assert picked is line + + def test_the_bound_reaches_the_candidates(self): + assert propagating_modes([self.METAL, self.LINE], max_index=2.6e3) == [ + self.LINE + ] + + def test_the_bound_is_on_the_magnitude(self): + """A Mode at the bound stays; one just above it on |n_eff| goes.""" + at_bound = 4.0 - 3.0j # |n_eff| = 5, inside the loss bound + assert propagating_modes([at_bound], max_index=5.0) == [at_bound] + assert propagating_modes([at_bound], max_index=4.99) == [] + + def test_a_bound_nothing_survives_is_reported(self): + with pytest.raises(NoLineModeError, match="No propagating line mode"): + select_line_mode([self.METAL], max_index=2.6e3) + + class TestGainBound: """A Mode that grows along the line is a numerical artefact, not a candidate. diff --git a/tests/femwell/test_adapter.py b/tests/femwell/test_adapter.py index c4920986..6217c76e 100644 --- a/tests/femwell/test_adapter.py +++ b/tests/femwell/test_adapter.py @@ -6,9 +6,11 @@ import numpy as np import pytest +from gsim.common.modes import select_line_mode from gsim.common.stack.extractor import Layer, LayerStack from gsim.common.stack.materials import MaterialProperties from gsim.femwell.adapter import ( + conductor_index_bound, elementwise_epsilon, epsilon_by_region, ) @@ -144,3 +146,66 @@ def test_user_override_wins_over_database(self, tmp_path): overrides={"si": MaterialProperties(permittivity=9.0)}, ) assert eps["core"].real == pytest.approx(9.0) + + +class TestConductorIndexBound: + """The |n_eff| above which a Mode lives inside a conductor (#498).""" + + def test_a_tenth_of_the_metal_index(self): + eps_metal = 1.0 - 6.8e8j # aluminium at 1 GHz + bound = conductor_index_bound({"metal": eps_metal, "oxide": 3.9 + 0j}) + assert bound == pytest.approx(0.1 * abs(np.sqrt(eps_metal))) + + def test_the_least_conductive_strong_conductor_sets_it(self): + bound = conductor_index_bound( + {"al": 1.0 - 6.8e8j, "w": 1.0 - 3.3e8j, "oxide": 3.9 + 0j} + ) + assert bound == pytest.approx(0.1 * abs(np.sqrt(1.0 - 3.3e8j))) + + def test_a_weak_conductor_alone_gives_no_bound(self): + """11.7-12j conducts, but 0.1|sqrt(eps)| = 0.41 is below the light line.""" + with pytest.warns(UserWarning, match="skipped conductor"): + assert conductor_index_bound({"si": 11.7 - 12j, "oxide": 3.9 + 0j}) is None + + def test_a_weak_conductor_does_not_lower_the_metal_bound(self): + """Doped silicon next to aluminium leaves the aluminium bound in place.""" + eps = {"al": 1.0 - 6.8e8j, "si": 11.7 - 1000j, "oxide": 3.9 + 0j} + with pytest.warns(UserWarning, match="skipped conductor"): + bound = conductor_index_bound(eps) + assert bound == pytest.approx(0.1 * abs(np.sqrt(1.0 - 6.8e8j))) + line = 5.13 - 4.68j + assert select_line_mode([line, 2.5 - 0.01j], max_index=bound) == line + + def test_the_margin_is_against_the_largest_dielectric_index(self): + """A conductor clears oxide by 10x but not silicon is skipped.""" + eps_cond = 1.0 - 9e4j # 0.1|sqrt(eps)| = 30 + assert conductor_index_bound({"c": eps_cond, "ox": 3.9 + 0j}) == ( + pytest.approx(30.0, rel=1e-4) + ) + with pytest.warns(UserWarning, match="skipped conductor"): + assert ( + conductor_index_bound({"c": eps_cond, "si": 11.7 + 0j, "ox": 3.9 + 0j}) + is None + ) + + def test_per_element_values(self): + eps = np.array([3.9, 3.9, 1.0 - 6.8e8j, 3.9], dtype=complex) + assert conductor_index_bound(eps) == pytest.approx( + 0.1 * abs(np.sqrt(1.0 - 6.8e8j)) + ) + + def test_a_lossy_dielectric_is_not_a_conductor(self): + """Loss below the real permittivity leaves the Mode free to sit there.""" + assert conductor_index_bound({"si": 11.7 - 5.0j, "oxide": 3.9 + 0j}) is None + + def test_the_fraction_is_a_parameter(self): + eps = {"metal": 1.0 - 6.8e8j} + default = conductor_index_bound(eps) + assert default is not None + assert conductor_index_bound(eps, fraction=0.5) == pytest.approx(5 * default) + + def test_the_fraction_and_margin_must_be_positive(self): + with pytest.raises(ValueError, match="fraction"): + conductor_index_bound({"metal": 1.0 - 6.8e8j}, fraction=0.0) + with pytest.raises(ValueError, match="margin"): + conductor_index_bound({"metal": 1.0 - 6.8e8j}, margin=0.0) diff --git a/tests/femwell/test_series_resistance.py b/tests/femwell/test_series_resistance.py index e7ce1765..fefa7fe4 100644 --- a/tests/femwell/test_series_resistance.py +++ b/tests/femwell/test_series_resistance.py @@ -28,7 +28,7 @@ from gsim.common.modes import Conductor, select_line_mode from gsim.common.transmission_line import line_params_from_neff -from gsim.femwell import solve_modes, z0_power_current +from gsim.femwell import conductor_index_bound, solve_modes, z0_power_current pytest.importorskip("femwell") pytest.importorskip("skfem") @@ -49,11 +49,9 @@ ) -@pytest.fixture(scope="module") -def strip_mesh(tmp_path_factory): +def _write_strip_mesh(path, *, near_um: float, far_um: float): import gmsh - path = tmp_path_factory.mktemp("femwell-series") / "strip.msh" gmsh.initialize() try: gmsh.option.setNumber("General.Terminal", 0) @@ -76,12 +74,12 @@ def strip_mesh(tmp_path_factory): # developer ~/.gmsh-options; gsim's own mesh paths pin these. gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0) gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0) - gmsh.option.setNumber("Mesh.MeshSizeMax", 2.0) - # 50 nm elements in and around the strip, 2 um in the far oxide. + gmsh.option.setNumber("Mesh.MeshSizeMax", far_um) + # near_um elements in and around the strip, far_um in the oxide. field = gmsh.model.mesh.field near = field.add("Box") - field.setNumber(near, "VIn", 0.05) - field.setNumber(near, "VOut", 2.0) + field.setNumber(near, "VIn", near_um) + field.setNumber(near, "VOut", far_um) field.setNumber(near, "XMin", -WIDTH_UM / 2 - 0.1) field.setNumber(near, "XMax", WIDTH_UM / 2 + 0.1) field.setNumber(near, "YMin", -THICKNESS_UM / 2 - 0.1) @@ -96,21 +94,41 @@ def strip_mesh(tmp_path_factory): return path -def _series_resistance(mesh_path, **solve_kwargs) -> float: +@pytest.fixture(scope="module") +def strip_mesh(tmp_path_factory): + """50 nm at the strip, 2 um in the far oxide.""" + path = tmp_path_factory.mktemp("femwell-series") / "strip.msh" + return _write_strip_mesh(path, near_um=0.05, far_um=2.0) + + +@pytest.fixture(scope="module") +def capped_strip_mesh(tmp_path_factory): + """25 nm at the strip, oxide capped at 0.5 um: the mesh of #498.""" + path = tmp_path_factory.mktemp("femwell-series-capped") / "strip.msh" + return _write_strip_mesh(path, near_um=0.025, far_um=0.5) + + +def _epsilon() -> dict[str, complex]: omega = 2 * np.pi * FREQ_HZ + return { + "metal": 1.0 - 1j * SIGMA_AL / (omega * EPS0), + "oxide": EPS_OXIDE + 0j, + } + + +def _series_resistance( + mesh_path, *, num_modes: int = 4, max_index=None, **solve_kwargs +) -> float: modes = solve_modes( mesh_path, - epsilon={ - "metal": 1.0 - 1j * SIGMA_AL / (omega * EPS0), - "oxide": EPS_OXIDE + 0j, - }, + epsilon=_epsilon(), wavelength_um=C0 / FREQ_HZ * 1e6, - num_modes=4, + num_modes=num_modes, metallic_boundaries=True, n_guess=3.0, **solve_kwargs, ) - mode = select_line_mode(modes) + mode = select_line_mode(modes, max_index=max_index) z0 = z0_power_current(mode, frequency_hz=FREQ_HZ, conductor=STRIP, mesh=mesh_path) params = line_params_from_neff( np.array([FREQ_HZ]), [complex(mode.n_eff)], z0_ohm=[z0] @@ -127,3 +145,21 @@ def test_order_one_is_still_available_and_reads_low(self, strip_mesh): """The first-order solve stays callable; on this mesh it is percent low.""" r_series = _series_resistance(strip_mesh, order=1) assert r_series < (1 - 3e-2) * R_DC_OHM_PER_M + + +class TestManyModes: + """More Modes than the line needs also bring Modes inside the metal (#498). + + On the capped mesh at order 2, a 12-Mode solve returns Modes with + ``|n_eff|`` of 1e4 to 5e4 confined to the strip, some inside the loss + bound; the slowest-Mode rule picked one and read R' about 1e9 off. The + conductor index bound rejects them and the line Mode is read again. + """ + + def test_the_bound_selects_the_line_mode(self, capped_strip_mesh): + r_series = _series_resistance( + capped_strip_mesh, + num_modes=12, + max_index=conductor_index_bound(_epsilon()), + ) + assert r_series == pytest.approx(R_DC_OHM_PER_M, rel=3e-3)