From 06500963256141d4c3044453ade5f75b2eebe7d0 Mon Sep 17 00:00:00 2001 From: George-Guryev-flxcmp Date: Tue, 17 Feb 2026 16:49:47 -0500 Subject: [PATCH] feat(FXC-5013): Lumped element support for arbitrary RLC circuit --- changelog.d/3336.added.md | 1 + changelog.d/3336.planned_deprecation.md | 1 + schemas/EMESimulation.json | 37 + schemas/ModeSimulation.json | 37 + schemas/Simulation.json | 37 + schemas/TerminalComponentModeler.json | 37 + tests/test_components/test_lumped_element.py | 614 +++++++++++++- tidy3d/__init__.py | 2 + tidy3d/components/lumped_element.py | 819 ++++++++++++++++++- tidy3d/rf.py | 2 + 10 files changed, 1583 insertions(+), 4 deletions(-) create mode 100644 changelog.d/3336.added.md create mode 100644 changelog.d/3336.planned_deprecation.md diff --git a/changelog.d/3336.added.md b/changelog.d/3336.added.md new file mode 100644 index 0000000000..292355bb8a --- /dev/null +++ b/changelog.d/3336.added.md @@ -0,0 +1 @@ +Added `CircuitImpedanceModel` to support arbitrary RLC circuits in `LinearLumpedElement`, with class methods `from_spice_file()` and `from_component_list()` for instantiation from a SPICE netlist or a component list. diff --git a/changelog.d/3336.planned_deprecation.md b/changelog.d/3336.planned_deprecation.md new file mode 100644 index 0000000000..6c34248fbc --- /dev/null +++ b/changelog.d/3336.planned_deprecation.md @@ -0,0 +1 @@ +`RLCNetwork` (a deprecation warning is now issued). Use `CircuitImpedanceModel.from_component_list()` or `CircuitImpedanceModel.from_spice_file()` instead. `AdmittanceNetwork` may be renamed to `AdmittanceModel` in a future release. diff --git a/schemas/EMESimulation.json b/schemas/EMESimulation.json index 0e3f0f2937..93b22e3faf 100644 --- a/schemas/EMESimulation.json +++ b/schemas/EMESimulation.json @@ -1539,6 +1539,39 @@ ], "type": "object" }, + "CircuitImpedanceModel": { + "additionalProperties": false, + "properties": { + "a": { + "items": { + "minimum": 0, + "type": "number" + }, + "type": "array" + }, + "attrs": { + "additionalProperties": true, + "type": "object" + }, + "b": { + "items": { + "minimum": 0, + "type": "number" + }, + "type": "array" + }, + "type": { + "const": "CircuitImpedanceModel", + "default": "CircuitImpedanceModel", + "type": "string" + } + }, + "required": [ + "a", + "b" + ], + "type": "object" + }, "ClipOperation": { "additionalProperties": false, "properties": { @@ -6797,6 +6830,7 @@ "discriminator": { "mapping": { "AdmittanceNetwork": "#/$defs/AdmittanceNetwork", + "CircuitImpedanceModel": "#/$defs/CircuitImpedanceModel", "RLCNetwork": "#/$defs/RLCNetwork" }, "propertyName": "type" @@ -6805,6 +6839,9 @@ { "$ref": "#/$defs/AdmittanceNetwork" }, + { + "$ref": "#/$defs/CircuitImpedanceModel" + }, { "$ref": "#/$defs/RLCNetwork" } diff --git a/schemas/ModeSimulation.json b/schemas/ModeSimulation.json index 48b0368a98..c9a4a1f2ba 100644 --- a/schemas/ModeSimulation.json +++ b/schemas/ModeSimulation.json @@ -1741,6 +1741,39 @@ ], "type": "object" }, + "CircuitImpedanceModel": { + "additionalProperties": false, + "properties": { + "a": { + "items": { + "minimum": 0, + "type": "number" + }, + "type": "array" + }, + "attrs": { + "additionalProperties": true, + "type": "object" + }, + "b": { + "items": { + "minimum": 0, + "type": "number" + }, + "type": "array" + }, + "type": { + "const": "CircuitImpedanceModel", + "default": "CircuitImpedanceModel", + "type": "string" + } + }, + "required": [ + "a", + "b" + ], + "type": "object" + }, "ClipOperation": { "additionalProperties": false, "properties": { @@ -6284,6 +6317,7 @@ "discriminator": { "mapping": { "AdmittanceNetwork": "#/$defs/AdmittanceNetwork", + "CircuitImpedanceModel": "#/$defs/CircuitImpedanceModel", "RLCNetwork": "#/$defs/RLCNetwork" }, "propertyName": "type" @@ -6292,6 +6326,9 @@ { "$ref": "#/$defs/AdmittanceNetwork" }, + { + "$ref": "#/$defs/CircuitImpedanceModel" + }, { "$ref": "#/$defs/RLCNetwork" } diff --git a/schemas/Simulation.json b/schemas/Simulation.json index 4348ea614e..66e70fe306 100644 --- a/schemas/Simulation.json +++ b/schemas/Simulation.json @@ -2181,6 +2181,39 @@ ], "type": "object" }, + "CircuitImpedanceModel": { + "additionalProperties": false, + "properties": { + "a": { + "items": { + "minimum": 0, + "type": "number" + }, + "type": "array" + }, + "attrs": { + "additionalProperties": true, + "type": "object" + }, + "b": { + "items": { + "minimum": 0, + "type": "number" + }, + "type": "array" + }, + "type": { + "const": "CircuitImpedanceModel", + "default": "CircuitImpedanceModel", + "type": "string" + } + }, + "required": [ + "a", + "b" + ], + "type": "object" + }, "ClipOperation": { "additionalProperties": false, "properties": { @@ -9145,6 +9178,7 @@ "discriminator": { "mapping": { "AdmittanceNetwork": "#/$defs/AdmittanceNetwork", + "CircuitImpedanceModel": "#/$defs/CircuitImpedanceModel", "RLCNetwork": "#/$defs/RLCNetwork" }, "propertyName": "type" @@ -9153,6 +9187,9 @@ { "$ref": "#/$defs/AdmittanceNetwork" }, + { + "$ref": "#/$defs/CircuitImpedanceModel" + }, { "$ref": "#/$defs/RLCNetwork" } diff --git a/schemas/TerminalComponentModeler.json b/schemas/TerminalComponentModeler.json index 1773246fbd..edead74c16 100644 --- a/schemas/TerminalComponentModeler.json +++ b/schemas/TerminalComponentModeler.json @@ -2181,6 +2181,39 @@ ], "type": "object" }, + "CircuitImpedanceModel": { + "additionalProperties": false, + "properties": { + "a": { + "items": { + "minimum": 0, + "type": "number" + }, + "type": "array" + }, + "attrs": { + "additionalProperties": true, + "type": "object" + }, + "b": { + "items": { + "minimum": 0, + "type": "number" + }, + "type": "array" + }, + "type": { + "const": "CircuitImpedanceModel", + "default": "CircuitImpedanceModel", + "type": "string" + } + }, + "required": [ + "a", + "b" + ], + "type": "object" + }, "ClipOperation": { "additionalProperties": false, "properties": { @@ -9319,6 +9352,7 @@ "discriminator": { "mapping": { "AdmittanceNetwork": "#/$defs/AdmittanceNetwork", + "CircuitImpedanceModel": "#/$defs/CircuitImpedanceModel", "RLCNetwork": "#/$defs/RLCNetwork" }, "propertyName": "type" @@ -9327,6 +9361,9 @@ { "$ref": "#/$defs/AdmittanceNetwork" }, + { + "$ref": "#/$defs/CircuitImpedanceModel" + }, { "$ref": "#/$defs/RLCNetwork" } diff --git a/tests/test_components/test_lumped_element.py b/tests/test_components/test_lumped_element.py index ec2795e63f..e1b38674f0 100644 --- a/tests/test_components/test_lumped_element.py +++ b/tests/test_components/test_lumped_element.py @@ -2,13 +2,19 @@ from __future__ import annotations +import warnings + import numpy as np import pytest from pydantic import ValidationError import tidy3d as td from tidy3d.components.geometry.utils import SnapBehavior, SnapLocation -from tidy3d.components.lumped_element import network_complex_permittivity +from tidy3d.components.lumped_element import ( + Component, + NodeMapper, + network_complex_permittivity, +) def test_lumped_resistor(): @@ -480,3 +486,609 @@ def test_very_small_lateral_width_snapping(): # Should not raise division by zero error structure = element.to_structure(grid) + + +# --- SPICE parser tests --- + + +def test_parse_spice_value_no_suffix(): + """Parse numeric values without scale suffix.""" + parse = td.CircuitImpedanceModel._parse_spice_value + assert parse("50") == 50.0 + assert parse("1.5") == 1.5 + assert parse("1e-12") == 1e-12 + assert parse("+2.5") == 2.5 + assert parse("-3") == -3.0 + + +@pytest.mark.parametrize( + "s,expected", + [ + ("1K", 1e3), + ("1k", 1e3), + ("2.5K", 2.5e3), + ("1MEG", 1e6), + ("1meg", 1e6), + ("1M", 1e-3), + ("10m", 0.01), + ("10M", 0.01), + ("1u", 1e-6), + ("1U", 1e-6), + ("1n", 1e-9), + ("1p", 1e-12), + ("1f", 1e-15), + ("1G", 1e9), + ("1g", 1e9), + ("100T", 100e12), + ("1T", 1e12), + ], +) +def test_parse_spice_value_scale_suffixes(s, expected): + """Parse values per SPICE convention: T (tera), G (giga), K (kilo), MEG (mega), m/M (milli), u, n, p, f.""" + assert td.CircuitImpedanceModel._parse_spice_value(s) == pytest.approx(expected) + + +def test_parse_spice_value_invalid(): + """Invalid value or suffix raises ``ValueError``.""" + parse = td.CircuitImpedanceModel._parse_spice_value + with pytest.raises(ValueError, match="Empty value"): + parse("") + with pytest.raises(ValueError, match="Cannot parse"): + parse("abc") + with pytest.raises(ValueError, match="Unknown scale suffix"): + parse("1x") + + +def test_parse_spice_file_no_source_port_from_first_element(tmp_path): + """With no voltage source, port is taken from the first element's nodes (``_parse_spice_file``).""" + netlist = """ + Title: simple RC + * simple RC + R1 1 0 50 + C1 1 0 1p + """ + path = tmp_path / "circuit.cir" + path.write_text(netlist) + comp_list, port_plus, port_minus = td.CircuitImpedanceModel._parse_spice_file(path) + assert len(comp_list) == 2 + assert comp_list[0].element_type == "R" + assert comp_list[0].node_plus == "1" + assert comp_list[0].node_minus == "0" + assert comp_list[0].value == 50.0 + assert comp_list[0].name == "R1" + assert comp_list[1].element_type == "C" + assert comp_list[1].value == 1e-12 + assert port_plus == "1" + assert port_minus == "0" + + +def test_parse_spice_file_single_voltage_source_defines_port(tmp_path): + """Single voltage source defines port nodes (``_parse_spice_file``).""" + netlist = """ + Title: RC with V source + V1 1 0 DC 0 AC 1 + R1 1 0 50 + C1 1 0 1p + """ + path = tmp_path / "circuit.cir" + path.write_text(netlist) + comp_list, port_plus, port_minus = td.CircuitImpedanceModel._parse_spice_file(path) + assert len(comp_list) == 2 + assert port_plus == "1" + assert port_minus == "0" + + +def test_parse_spice_file_comments_and_continuation(tmp_path): + """Comment lines and + continuation are handled.""" + netlist = """ + Title + $ comment + * another comment + R1 1 0 50 + C1 1 0 + + 1p + """ + path = tmp_path / "circuit.cir" + path.write_text(netlist) + comp_list, _, _ = td.CircuitImpedanceModel._parse_spice_file(path) + assert len(comp_list) == 2 + assert comp_list[1].value == 1e-12 + + +def test_parse_spice_file_scale_suffixes_in_netlist(tmp_path): + """Scale suffixes ``K``, ``M``, ``m``, ``p`` in netlist are parsed.""" + netlist = """ + Title + R1 1 0 2K + C1 1 0 0.5p + L1 2 1 10n + """ + path = tmp_path / "circuit.cir" + path.write_text(netlist) + comp_list, _, _ = td.CircuitImpedanceModel._parse_spice_file(path) + assert comp_list[0].value == 2000.0 + assert comp_list[1].value == 0.5e-12 + assert comp_list[2].value == 10e-9 + + +def test_parse_spice_file_two_voltage_sources_raises(tmp_path): + """At most one voltage source is allowed.""" + netlist = """ + Title + V1 1 0 DC 0 + V2 2 0 DC 0 + R1 1 0 50 + """ + path = tmp_path / "circuit.cir" + path.write_text(netlist) + with pytest.raises(ValueError, match="at most one voltage source"): + td.CircuitImpedanceModel._parse_spice_file(path) + + +def test_parse_spice_file_no_rlc_raises(tmp_path): + """Netlist must contain at least one ``R``, ``C``, or ``L`` (``_parse_spice_file``).""" + netlist = """ + Title + V1 1 0 DC 0 + """ + path = tmp_path / "circuit.cir" + path.write_text(netlist) + with pytest.raises(ValueError, match="no R, C, or L"): + td.CircuitImpedanceModel._parse_spice_file(path) + + +def test_parse_spice_file_component_line_too_short_raises(tmp_path): + """Malformed component line raises.""" + netlist = """ + Title + R1 1 0 + """ + path = tmp_path / "circuit.cir" + path.write_text(netlist) + with pytest.raises(ValueError, match="Component line must have name"): + td.CircuitImpedanceModel._parse_spice_file(path) + + +def test_component_value_must_be_positive(): + """Component ``value`` must be positive (R, L, C); zero or negative raises ValidationError.""" + with pytest.raises(ValidationError): + Component(element_type="R", node_plus="1", node_minus="0", value=0.0, name="R1") + with pytest.raises(ValidationError): + Component(element_type="R", node_plus="1", node_minus="0", value=-50.0, name="R1") + with pytest.raises(ValidationError): + Component(element_type="C", node_plus="1", node_minus="0", value=0.0, name="C1") + Component(element_type="R", node_plus="1", node_minus="0", value=50.0, name="R1") # ok + + +def test_parse_spice_file_zero_value_raises(tmp_path): + """SPICE netlist with zero or negative component value raises when building ``Component``.""" + for bad_value in ("0", "-1"): + netlist = f""" +* circuit +R1 1 0 {bad_value} +""" + path = tmp_path / "bad_value.cir" + path.write_text(netlist) + with pytest.raises(ValidationError): + td.CircuitImpedanceModel._parse_spice_file(path) + + +def _effective_admittance_from_component_list( + component_list: list[Component], + frequencies: np.ndarray, + port_plus_node: str = "1", + port_minus_node: str = "0", +) -> np.ndarray: + """Compute one-port admittance from component list (same logic as ``from_component_list``).""" + return td.CircuitImpedanceModel._get_effective_admittance( + component_list, frequencies, port_plus_node=port_plus_node, port_minus_node=port_minus_node + ) + + +def test_effective_admittance_parallel_rc(tmp_path): + """Component list, SPICE (with and without voltage source), and analytical agree for parallel RC. + + Circuit: ``R`` and ``C`` both between node 1 and 0. Port 1-0. + Analytical: ``Y = 1/R + j*omega*C``. + """ + R, C = 50.0, 1e-12 + component_list = [ + Component(element_type="R", node_plus="1", node_minus="0", value=R, name="R1"), + Component(element_type="C", node_plus="1", node_minus="0", value=C, name="C1"), + ] + freqs = np.linspace(0.2e9, 8e9, 30) + port_plus, port_minus = "1", "0" + omega = 2 * np.pi * freqs + Y_analytical = 1.0 / R + 1j * omega * C + + Y_from_list = _effective_admittance_from_component_list( + component_list, freqs, port_plus_node=port_plus, port_minus_node=port_minus + ) + assert np.allclose(Y_from_list, Y_analytical, rtol=1e-14, atol=1e-20) + + # SPICE without voltage source: port from first element + netlist = f""" + Title + R1 1 0 {R} + C1 1 0 {C} + """ + path = tmp_path / "parallel_rc.cir" + path.write_text(netlist) + comp_list_spice, port_plus_s, port_minus_s = td.CircuitImpedanceModel._parse_spice_file(path) + assert port_plus_s == port_plus and port_minus_s == port_minus + Y_from_spice = _effective_admittance_from_component_list( + comp_list_spice, freqs, port_plus_node=port_plus_s, port_minus_node=port_minus_s + ) + assert np.allclose(Y_from_list, Y_from_spice, rtol=1e-14, atol=1e-20) + assert np.allclose(Y_from_spice, Y_analytical, rtol=1e-14, atol=1e-20) + + # SPICE with voltage source: port from V + netlist_v = f""" + Title + V1 1 0 DC 0 AC 1 + R1 1 0 {R} + C1 1 0 {C} + """ + path_v = tmp_path / "with_v.cir" + path_v.write_text(netlist_v) + comp_list_spice_v, port_plus_v, port_minus_v = td.CircuitImpedanceModel._parse_spice_file( + path_v + ) + assert port_plus_v == "1" and port_minus_v == "0" + Y_from_spice_v = _effective_admittance_from_component_list( + comp_list_spice_v, freqs, port_plus_node=port_plus_v, port_minus_node=port_minus_v + ) + assert np.allclose(Y_from_list, Y_from_spice_v, rtol=1e-14, atol=1e-20) + assert np.allclose(Y_from_spice_v, Y_analytical, rtol=1e-14, atol=1e-20) + + +def test_effective_admittance_series_rc(tmp_path): + """Component list, SPICE, and analytical agree for series RC. + + Circuit: node 1 -- ``R`` -- node 2 -- ``C`` -- node 0. Port 1-0. + Analytical: ``Y = 1/(R + 1/(j*omega*C)) = j*omega*C / (1 + j*omega*R*C)``. + """ + R, C = 50.0, 1e-12 + component_list = [ + Component(element_type="R", node_plus="1", node_minus="2", value=R, name="R1"), + Component(element_type="C", node_plus="2", node_minus="0", value=C, name="C1"), + ] + freqs = np.linspace(0.2e9, 8e9, 30) + port_plus, port_minus = "1", "0" + omega = 2 * np.pi * freqs + Y_analytical = (1j * omega * C) / (1.0 + 1j * omega * R * C) + + Y_from_list = _effective_admittance_from_component_list( + component_list, freqs, port_plus_node=port_plus, port_minus_node=port_minus + ) + assert np.allclose(Y_from_list, Y_analytical, rtol=1e-14, atol=1e-20) + + # SPICE: use voltage source to define port 1-0 (no single element has both 1 and 0 in series RC) + netlist = f""" + Title + V1 1 0 DC 0 AC 1 + R1 1 2 {R} + C1 2 0 {C} + """ + path = tmp_path / "series_rc.cir" + path.write_text(netlist) + comp_list_spice, port_plus_s, port_minus_s = td.CircuitImpedanceModel._parse_spice_file(path) + assert port_plus_s == "1" and port_minus_s == "0" + Y_from_spice = _effective_admittance_from_component_list( + comp_list_spice, freqs, port_plus_node=port_plus_s, port_minus_node=port_minus_s + ) + assert np.allclose(Y_from_list, Y_from_spice, rtol=1e-14, atol=1e-20) + assert np.allclose(Y_from_spice, Y_analytical, rtol=1e-14, atol=1e-20) + + +def test_effective_admittance_series_rcl(tmp_path): + """Component list, SPICE, and analytical agree for ``L`` in parallel with (``R`` in series ``C``). + + Circuit: ``L`` between 1-0, ``R`` between 1-2, ``C`` between 2-0. Port 1-0. + Analytical: ``Y = 1/(j*omega*L) + 1/(R + 1/(j*omega*C))``. + """ + R, C, L = 50.0, 1e-12, 2e-9 + component_list = [ + Component(element_type="R", node_plus="1", node_minus="2", value=R, name="R1"), + Component(element_type="C", node_plus="2", node_minus="0", value=C, name="C1"), + Component(element_type="L", node_plus="1", node_minus="0", value=L, name="L1"), + ] + freqs = np.linspace(0.2e9, 8e9, 30) + port_plus, port_minus = "1", "0" + omega = 2 * np.pi * freqs + Y_L = 1.0 / (1j * omega * L) + Z_RC = R + 1.0 / (1j * omega * C) + Y_analytical = Y_L + 1.0 / Z_RC + + Y_from_list = _effective_admittance_from_component_list( + component_list, freqs, port_plus_node=port_plus, port_minus_node=port_minus + ) + assert np.allclose(Y_from_list, Y_analytical, rtol=1e-14, atol=1e-20) + + netlist = f""" + Title + L1 1 0 {L} + R1 1 2 {R} + C1 2 0 {C} + """ + path = tmp_path / "series_rcl.cir" + path.write_text(netlist) + comp_list_spice, port_plus_s, port_minus_s = td.CircuitImpedanceModel._parse_spice_file(path) + assert port_plus_s == "1" and port_minus_s == "0", "port from first element (L1 1 0)" + Y_from_spice = _effective_admittance_from_component_list( + comp_list_spice, freqs, port_plus_node=port_plus_s, port_minus_node=port_minus_s + ) + assert np.allclose(Y_from_list, Y_from_spice, rtol=1e-14, atol=1e-20) + assert np.allclose(Y_from_spice, Y_analytical, rtol=1e-14, atol=1e-20) + + +def test_effective_admittance_port_minus_not_ground(): + """Driving-point admittance with port_minus_node not ground matches analytical. + + Circuit: node 1 -- R -- node 2 -- C -- node 0 (series RC). Port 1-2 (reference = node 2). + With node 2 as reference, the admittance at node 1 is 1/R (only R connects 1 to 2). + """ + R, C = 50.0, 1e-12 + component_list = [ + Component(element_type="R", node_plus="1", node_minus="2", value=R, name="R1"), + Component(element_type="C", node_plus="2", node_minus="0", value=C, name="C1"), + ] + freqs = np.linspace(0.2e9, 8e9, 30) + Y_analytical = np.full(len(freqs), 1.0 / R, dtype=complex) # real 1/R + + Y_from_list = _effective_admittance_from_component_list( + component_list, freqs, port_plus_node="1", port_minus_node="2" + ) + assert np.allclose(Y_from_list, Y_analytical, rtol=1e-14, atol=1e-20) + assert np.allclose(Y_from_list.real, 1.0 / R, rtol=1e-14, atol=1e-20) + assert np.allclose(Y_from_list.imag, 0.0, rtol=1e-14, atol=1e-20) + + +def test_effective_admittance_port_nodes_same_raises(): + """Port nodes must be distinct; same node raises ValueError.""" + component_list = [ + Component(element_type="R", node_plus="1", node_minus="0", value=50.0, name="R1"), + ] + freqs = np.linspace(1e9, 2e9, 5) + with pytest.raises(ValueError, match="distinct|cannot be the same"): + _effective_admittance_from_component_list( + component_list, freqs, port_plus_node="1", port_minus_node="1" + ) + with pytest.raises(ValueError, match="distinct|cannot be the same"): + _effective_admittance_from_component_list( + component_list, freqs, port_plus_node="0", port_minus_node="0" + ) + + +# --- Coverage: NodeMapper, Component, SPICE parser, CircuitImpedanceModel --- + + +def test_nodemapper_indices_by_first_appearance(): + """Node indices are assigned in order of first appearance; no special treatment for any label.""" + mapper = NodeMapper() + assert mapper.get_or_create_index("1") == 0 + assert mapper.get_or_create_index("2") == 1 + assert mapper.get_or_create_index("0") == 2 + assert mapper.total_nodes() == 3 + assert mapper.lookup_index("1") == 0 + assert mapper.lookup_index("2") == 1 + assert mapper.lookup_index("0") == 2 + # "0" and "GND" are distinct nodes (no special ground mapping) + mapper2 = NodeMapper() + assert mapper2.get_or_create_index("0") == 0 + assert mapper2.get_or_create_index("GND") == 1 + assert mapper2.get_or_create_index("gnd") == 2 + assert mapper2.total_nodes() == 3 + + +def test_nodemapper_lowercase_gnd_same_admittance_as_zero(): + """Circuit with node ``gnd`` yields same effective admittance as circuit with ``0`` (same topology).""" + R, C = 50.0, 1e-12 + freqs = np.linspace(0.2e9, 8e9, 20) + comp_list_zero = [ + Component(element_type="R", node_plus="1", node_minus="0", value=R, name="R1"), + Component(element_type="C", node_plus="1", node_minus="0", value=C, name="C1"), + ] + comp_list_gnd = [ + Component(element_type="R", node_plus="1", node_minus="gnd", value=R, name="R1"), + Component(element_type="C", node_plus="1", node_minus="gnd", value=C, name="C1"), + ] + Y_zero = _effective_admittance_from_component_list( + comp_list_zero, freqs, port_plus_node="1", port_minus_node="0" + ) + Y_gnd = _effective_admittance_from_component_list( + comp_list_gnd, freqs, port_plus_node="1", port_minus_node="gnd" + ) + assert np.allclose(Y_zero, Y_gnd, rtol=1e-14, atol=1e-20) + + +def test_parse_spice_file_comment_title_first_line_components_not_skipped(tmp_path): + """When the first line is a comment (e.g. ``* Title``), it is filtered; first component is not skipped.""" + netlist = """* Title +R1 1 0 50 +C1 1 0 1p +""" + path = tmp_path / "comment_title.cir" + path.write_text(netlist) + comp_list, port_plus, port_minus = td.CircuitImpedanceModel._parse_spice_file(path) + assert len(comp_list) == 2 + assert comp_list[0].element_type == "R" + assert comp_list[0].name == "R1" + assert comp_list[0].value == 50.0 + assert comp_list[1].element_type == "C" + assert comp_list[1].value == 1e-12 + assert port_plus == "1" and port_minus == "0" + + +def test_parse_spice_file_plain_title_first_line_skipped(tmp_path): + """When the first non-comment line does not start with R/C/L/V, it is treated as title and skipped.""" + netlist = """My circuit +R1 1 0 50 +""" + path = tmp_path / "plain_title.cir" + path.write_text(netlist) + comp_list, port_plus, port_minus = td.CircuitImpedanceModel._parse_spice_file(path) + assert len(comp_list) == 1 + assert comp_list[0].element_type == "R" + assert comp_list[0].name == "R1" + assert port_plus == "1" and port_minus == "0" + + +def test_parse_spice_file_lowercase_gnd_parsed(tmp_path): + """SPICE netlist with lowercase ``gnd`` parses; port and components use ``gnd``, same topology as ``0``.""" + netlist = """* parallel RC +R1 1 gnd 50 +C1 1 gnd 1p +""" + path = tmp_path / "lowercase_gnd.cir" + path.write_text(netlist) + comp_list, port_plus, port_minus = td.CircuitImpedanceModel._parse_spice_file(path) + assert len(comp_list) == 2 + assert comp_list[0].node_minus == "gnd" + assert comp_list[1].node_minus == "gnd" + assert port_plus == "1" and port_minus == "gnd" + # Same admittance as 1-0 circuit + comp_list_zero = [ + Component(element_type="R", node_plus="1", node_minus="0", value=50.0, name="R1"), + Component(element_type="C", node_plus="1", node_minus="0", value=1e-12, name="C1"), + ] + freqs = np.linspace(0.2e9, 8e9, 15) + Y_gnd = _effective_admittance_from_component_list( + comp_list, freqs, port_plus_node="1", port_minus_node="gnd" + ) + Y_zero = _effective_admittance_from_component_list( + comp_list_zero, freqs, port_plus_node="1", port_minus_node="0" + ) + assert np.allclose(Y_gnd, Y_zero, rtol=1e-14, atol=1e-20) + + +def test_parse_spice_file_voltage_source_between_1_and_2_infers_port(tmp_path): + """SPICE with V between nodes 1 and 2 (no 0/GND) infers port; admittance matches component_list.""" + netlist = """V1 1 2 0 +R1 1 2 50 +""" + path = tmp_path / "v_1_2.cir" + path.write_text(netlist) + comp_list, port_plus, port_minus = td.CircuitImpedanceModel._parse_spice_file(path) + assert port_plus == "1" and port_minus == "2" + assert len(comp_list) == 1 + assert comp_list[0].node_plus == "1" and comp_list[0].node_minus == "2" + freqs = np.linspace(0.2e9, 8e9, 15) + Y_spice = _effective_admittance_from_component_list( + comp_list, freqs, port_plus_node=port_plus, port_minus_node=port_minus + ) + comp_list_direct = [ + Component(element_type="R", node_plus="1", node_minus="2", value=50.0, name="R1"), + ] + Y_direct = _effective_admittance_from_component_list( + comp_list_direct, freqs, port_plus_node="1", port_minus_node="2" + ) + assert np.allclose(Y_spice, Y_direct, rtol=1e-14, atol=1e-20) + assert np.allclose(Y_spice.real, 1.0 / 50.0, rtol=1e-14, atol=1e-20) + + +def test_lookup_index_raises_when_node_not_in_circuit(): + """lookup_index raises ValueError when port node is not in the circuit.""" + component_list = [ + Component(element_type="R", node_plus="1", node_minus="0", value=50.0, name="R1"), + ] + freqs = np.linspace(1e9, 2e9, 5) + with pytest.raises(ValueError, match="not in the circuit"): + _effective_admittance_from_component_list( + component_list, freqs, port_plus_node="99", port_minus_node="0" + ) + with pytest.raises(ValueError, match="not in the circuit"): + _effective_admittance_from_component_list( + component_list, freqs, port_plus_node="1", port_minus_node="99" + ) + + +def test_component_default_name_set_in_post_init(): + """Component with empty name gets a default name in model_post_init.""" + comp = Component(element_type="R", node_plus="1", node_minus="0", value=50.0, name="") + assert comp.name + assert comp.name.startswith("R") + + +def test_admittance_network_warns_circuit_impedance_model_does_not(): + """AdmittanceNetwork emits FutureWarning; CircuitImpedanceModel does not.""" + a, b = (1.0, 0.0), (1.0,) + with pytest.warns(FutureWarning, match="AdmittanceNetwork may be renamed"): + _ = td.AdmittanceNetwork(a=a, b=b) + with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always", FutureWarning) + _ = td.CircuitImpedanceModel(a=a, b=b) + assert not any("AdmittanceNetwork" in str(m.message) for m in w) + + +def test_parse_spice_file_voltage_source_too_short_raises(tmp_path): + """Voltage source line with fewer than two nodes raises.""" + netlist = """ + Title + V1 1 + """ + path = tmp_path / "bad_v.cir" + path.write_text(netlist) + with pytest.raises(ValueError, match="Voltage source line needs at least two nodes"): + td.CircuitImpedanceModel._parse_spice_file(path) + + +def test_get_effective_admittance_unknown_component_type_raises(): + """Component with element_type not R, L, or C raises in branch admittance.""" + component_list = [ + Component(element_type="X", node_plus="1", node_minus="0", value=1.0, name="X1"), + ] + freqs = np.array([1e9]) + with pytest.raises(ValueError, match="Unknown component type"): + td.CircuitImpedanceModel._get_effective_admittance( + component_list, freqs, port_plus_node="1", port_minus_node="0" + ) + + +def test_from_component_list_returns_model(): + """from_component_list builds a CircuitImpedanceModel and covers fit path.""" + component_list = [ + Component(element_type="R", node_plus="1", node_minus="0", value=50.0, name="R1"), + Component(element_type="C", node_plus="1", node_minus="0", value=1e-12, name="C1"), + ] + freqs = np.linspace(0.2e9, 8e9, 20) + model = td.CircuitImpedanceModel.from_component_list(component_list, freqs, show_progress=False) + assert hasattr(model, "a") and hasattr(model, "b") + assert len(model.a) >= 1 and len(model.b) >= 1 + + +def test_from_spice_file_port_override(tmp_path): + """from_spice_file with port_plus_node/port_minus_node overrides uses them.""" + netlist = """ + Title + R1 1 0 50 + C1 1 0 1p + """ + path = tmp_path / "rc.cir" + path.write_text(netlist) + freqs = np.linspace(0.2e9, 8e9, 15) + model = td.CircuitImpedanceModel.from_spice_file( + path, freqs, port_plus_node="1", port_minus_node="0", show_progress=False + ) + assert hasattr(model, "a") and hasattr(model, "b") + + +def test_from_touchstone_file_not_implemented(): + """from_touchstone_file raises NotImplementedError.""" + with pytest.raises(NotImplementedError, match="not yet implemented"): + td.CircuitImpedanceModel.from_touchstone_file("dummy.s1p") + + +def test_parse_spice_file_blank_line_in_body_skipped(tmp_path): + """Blank lines (after strip) in the netlist are skipped; parsing still succeeds.""" + netlist = """ + Title + + R1 1 0 50 + """ + path = tmp_path / "with_blank.cir" + path.write_text(netlist) + comp_list, p_plus, p_minus = td.CircuitImpedanceModel._parse_spice_file(path) + assert len(comp_list) == 1 + assert comp_list[0].element_type == "R" + assert p_plus == "1" and p_minus == "0" diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index ea5b63c43a..196d451353 100644 --- a/tidy3d/__init__.py +++ b/tidy3d/__init__.py @@ -307,6 +307,7 @@ # lumped elements from .components.lumped_element import ( AdmittanceNetwork, + CircuitImpedanceModel, CoaxialLumpedResistor, LinearLumpedElement, LumpedElement, @@ -566,6 +567,7 @@ def set_logging_level(level: str) -> None: "ChargeInsulatorMedium", "ChargeToleranceSpec", "ChebSampling", + "CircuitImpedanceModel", "ClipOperation", "CoaxialLumpedResistor", "CompositeCurrentIntegral", diff --git a/tidy3d/components/lumped_element.py b/tidy3d/components/lumped_element.py index 53cd053e10..164da906f5 100644 --- a/tidy3d/components/lumped_element.py +++ b/tidy3d/components/lumped_element.py @@ -2,19 +2,28 @@ from __future__ import annotations +import re +import warnings from abc import ABC, abstractmethod +from pathlib import Path from typing import TYPE_CHECKING, Literal, Optional, Union +if TYPE_CHECKING: + from typing import Any, Callable + import numpy as np from pydantic import ( + ConfigDict, Field, NonNegativeFloat, PositiveFloat, PositiveInt, + PrivateAttr, field_validator, model_validator, ) +from tidy3d.components.dispersion_fitter import AdvancedFastFitterParam, fit from tidy3d.components.types.base import discriminated_union from tidy3d.components.validators import assert_plane, validate_name_str from tidy3d.constants import EPSILON_0, FARAD, HENRY, MICROMETER, OHM @@ -55,6 +64,62 @@ LOSS_FACTOR_INDUCTOR = 1e6 +class NodeMapper(MicrowaveBaseModel): + """Maps node names to indices in order of first appearance. + + No special meaning is assigned to any label (e.g. ``'0'`` or ``'GND'``); + the effective admittance between two nodes depends only on topology and + port choice, not on which node gets which index. + """ + + model_config = ConfigDict(frozen=False) # allow mutating _name_to_idx in get_or_create_index() + + _name_to_idx: dict[str, int] = PrivateAttr(default_factory=dict) + _counter: int = PrivateAttr(default=0) + + def get_or_create_index(self, name: str) -> int: + """Register the node if new and return its index. Use when building the mapper from components.""" + if name not in self._name_to_idx: + idx = self._counter + self._name_to_idx[name] = idx + object.__setattr__(self, "_counter", self._counter + 1) + return self._name_to_idx[name] + + def lookup_index(self, name: str) -> int: + """Return index for a node that must already be in the mapper. Use for port nodes.""" + if name not in self._name_to_idx: + raise ValueError( + f"Node {name!r} is not in the circuit (not an endpoint of any R, L, or C component)." + ) + return self._name_to_idx[name] + + def total_nodes(self) -> int: + """Number of distinct node indices (0 through N-1).""" + return len(self._name_to_idx) + + +class Component(MicrowaveBaseModel): + """Single R, L, or C branch between two nodes.""" + + model_config = ConfigDict( + frozen=False + ) # allow mutating node_plus_idx/node_minus_idx in build_incidence_matrix + + element_type: str = Field(description="Element type: 'R', 'L', or 'C'.") + node_plus: str = Field(description="Name of the plus node.") + node_minus: str = Field(description="Name of the minus node.") + value: PositiveFloat = Field(description="Nominal value of the element.") + name: str = Field(default="", description="Optional component name.") + + # Set by build_incidence_matrix; not part of schema/serialization. + _node_plus_idx: int = PrivateAttr(default=-1) + _node_minus_idx: int = PrivateAttr(default=-1) + + def model_post_init(self, __context: Any) -> None: + if not self.name: + object.__setattr__(self, "name", f"{self.element_type}{id(self)}") + + class LumpedElement(MicrowaveBaseModel, ABC): """Base class describing the interface all lumped elements obey.""" @@ -561,6 +626,11 @@ class RLCNetwork(MicrowaveBaseModel): """Class for representing a simple network consisting of a resistor, capacitor, and inductor. Provides additional functionality for representing the network as an equivalent medium. + .. deprecated:: + :class:`RLCNetwork` is deprecated. Prefer :class:`CircuitImpedanceModel` for general RLC + circuits (e.g. :meth:`CircuitImpedanceModel.from_component_list` or + :meth:`CircuitImpedanceModel.from_spice_file`). + Notes ----- @@ -805,6 +875,17 @@ def _validate_single_element(self) -> Self: raise ValueError("At least one element must be defined in the 'RLCNetwork'.") return self + @model_validator(mode="after") + def _warn_deprecated(self) -> Self: + """Emit deprecation warning when RLCNetwork is used.""" + warnings.warn( + "RLCNetwork is deprecated; use CircuitImpedanceModel.from_component_list or " + "CircuitImpedanceModel.from_spice_file for general RLC circuits.", + category=DeprecationWarning, + stacklevel=2, + ) + return self + class AdmittanceNetwork(MicrowaveBaseModel): """Class for representing a network consisting of an arbitrary number of resistors, @@ -812,6 +893,11 @@ class AdmittanceNetwork(MicrowaveBaseModel): as an admittance function. Provides additional functionality for representing the network as an equivalent medium. + .. warning:: + This class may be renamed to ``AdmittanceModel`` in a future release. For building + networks from SPICE files or component lists, use :class:`CircuitImpedanceModel` + (e.g. :meth:`CircuitImpedanceModel.from_spice_file`, :meth:`CircuitImpedanceModel.from_component_list`). + Notes ----- @@ -866,9 +952,28 @@ class AdmittanceNetwork(MicrowaveBaseModel): "The length of the ``tuple`` is equal to the order of the network.", ) + @model_validator(mode="after") + def _warn_future_rename(self) -> Self: + """Warn about upcoming rename to AdmittanceModel.""" + if type(self).__name__ == "CircuitImpedanceModel": + return self + warnings.warn( + "AdmittanceNetwork may be renamed to AdmittanceModel in a future release. " + "For building from SPICE or component lists, use CircuitImpedanceModel.", + category=FutureWarning, + stacklevel=2, + ) + return self + def _to_medium(self, scaling_factor: float) -> PoleResidue: - """Converts the :class:`AdmittanceNetwork` model directly into a :class:`PoleResidue` model - with proper scaling depending on the lumped element's dimensions.""" + """Convert to a :class:`PoleResidue` medium with geometric scaling applied. + + The stored ``(a, b)`` coefficients represent the unscaled admittance :math:`Y(s)`. + When used in a :class:`~tidy3d.LinearLumpedElement`, ``scaling_factor`` is + :meth:`RectangularLumpedElement._admittance_transfer_function_scaling` (e.g. + ``size_voltage / size_lateral`` of the cell box). The returned medium corresponds + to :math:`\\Delta \\cdot Y(s)`, so geometric scaling is applied only here. + """ a = np.array(self.a) * scaling_factor b = np.array(self.b) return PoleResidue.from_admittance_coeffs(a, b) @@ -881,7 +986,715 @@ def _as_admittance_function(self) -> tuple[tuple[float, ...], tuple[float, ...]] return (self.a, self.b) -NetworkType = discriminated_union(Union[RLCNetwork, AdmittanceNetwork]) +class CircuitImpedanceModel(AdmittanceNetwork): + """Circuit model represented by its effective admittance as a rational function in the Laplace domain. + + Subclass of :class:`AdmittanceNetwork` with the same ``(a, b)`` signature. Adds factory methods + :meth:`from_spice_file`, :meth:`from_component_list`, and :meth:`from_touchstone_file` to build + from external descriptions. Construct directly with ``a`` and ``b`` or use the factories. + + Notes + ----- + The effective one-port admittance is + :math:`Y(s) = (a_0 + a_1 s + \\dots + a_M s^M) / (b_0 + b_1 s + \\dots + b_N s^N)`. + Use the factory methods to build from a SPICE netlist or a list of R/L/C components; + the Y→ε→pole-residue fit path is used to obtain stable ``(a, b)`` coefficients. + + See Also + -------- + AdmittanceNetwork : Base class with the same ``(a, b)`` storage and medium conversion. + LinearLumpedElement : Lumped element that accepts a :class:`CircuitImpedanceModel` as its network. + """ + + @staticmethod + def _parse_spice_value(s: str) -> float: + """Parse a SPICE value string including scale suffixes. + + Parameters + ---------- + s : str + Value string (e.g. ``"1K"``, ``"10n"``, ``"2.5p"``). + + Returns + ------- + float + Parsed value in base SI units (Ohms, Henrys, or Farads). + + Raises + ------ + ValueError + If the string is empty or contains an unknown scale suffix. + + Notes + ----- + Scale suffixes follow common SPICE convention: ``T`` (tera), ``G`` (giga), + ``MEG`` (mega), ``K`` (kilo), ``m`` / ``M`` (milli), ``u`` / ``U`` (micro), + ``n`` / ``N`` (nano), ``p`` / ``P`` (pico), ``f`` / ``F`` (femto). Both ``m`` + and ``M`` are milli (1e-3); use ``MEG`` for mega (1e6). Scientific notation + (e.g. ``1e-12``) is also allowed. + """ + s = s.strip() + if not s: + raise ValueError("Empty value") + m = re.match(r"([+-]?\d*\.?\d+([eE][+-]?\d+)?)(.*)", s, re.IGNORECASE) + if not m: + raise ValueError(f"Cannot parse value: {s!r}") + num_str, suffix = m.group(1), m.group(3).strip() + scale = 1.0 + if suffix: + suf_upper = suffix.upper() + if suf_upper == "MEG": + scale = 1e6 + elif suf_upper == "K": + scale = 1e3 + elif suf_upper == "M": + # SPICE: both m and M are milli (1e-3); MEG is mega + scale = 1e-3 + else: + scales = {"F": 1e-15, "P": 1e-12, "N": 1e-9, "U": 1e-6, "G": 1e9, "T": 1e12} + scale = scales.get(suf_upper) + if scale is None: + raise ValueError(f"Unknown scale suffix: {suffix!r}") + return float(num_str) * scale + + @staticmethod + def _parse_spice_file( + spice_file: str | Path, + ) -> tuple[list[Component], str, str]: + """Parse a SPICE netlist and return components plus port nodes. + + Parameters + ---------- + spice_file : str or Path + Path to the SPICE netlist file. + + Returns + ------- + tuple[list[Component], str, str] + ``(component_list, port_plus_node, port_minus_node)``. R, C, and L elements + are converted to :class:`Component` instances. Port is taken from the single + voltage source (V) if present, otherwise from the first element's two nodes. + + Raises + ------ + ValueError + If the file contains no R/C/L components, more than one voltage source, + or a malformed component line. + + Notes + ----- + The first non-empty, non-comment line is treated as the SPICE title only if it does + not look like a component or voltage source (i.e. does not start with R, C, L, or V). + Comment lines (starting with ``*`` or ``$``) and continuation lines (starting with ``+``) + are handled. Value scale suffixes are parsed via :meth:`_parse_spice_value`. + """ + path = Path(spice_file) + text = path.read_text() + lines = [] + for raw in text.splitlines(): + line = raw.strip() + if not line or line.startswith(("*", "$")): + continue + if line.startswith("+"): + if lines: + lines[-1] = lines[-1] + " " + line[1:].strip() + continue + lines.append(line) + + component_list: list[Component] = [] + port_plus_node: Optional[str] = None + port_minus_node: Optional[str] = None + + for i, line in enumerate(lines): + toks = line.split() + if not toks: + continue + kind = toks[0][0].upper() + # First line that does not look like R/C/L/V is treated as SPICE title and skipped. + if i == 0 and kind not in ("R", "C", "L", "V"): + continue + comp_name = toks[0] + + if kind == "V": + if port_plus_node is not None: + raise ValueError( + "SPICE file must contain at most one voltage source for port detection." + ) + if len(toks) < 3: + raise ValueError(f"Voltage source line needs at least two nodes: {line!r}") + port_plus_node = toks[1] + port_minus_node = toks[2] + continue + + if kind in ("R", "C", "L"): + if len(toks) < 4: + raise ValueError( + f"Component line must have name, node+, node-, value: {line!r}" + ) + node_plus = toks[1] + node_minus = toks[2] + value = CircuitImpedanceModel._parse_spice_value(toks[3]) + comp = Component( + element_type=kind, + node_plus=node_plus, + node_minus=node_minus, + value=value, + name=comp_name, + ) + component_list.append(comp) + continue + + if not component_list: + raise ValueError("SPICE file contains no R, C, or L components.") + + if port_plus_node is None: + port_plus_node = component_list[0].node_plus + port_minus_node = component_list[0].node_minus + + return (component_list, port_plus_node, port_minus_node) + + @staticmethod + def _create_branch_admittance_matrix( + component_list: list[Component], frequency: float + ) -> np.ndarray: + """Diagonal matrix of branch admittances at a given frequency. + + Parameters + ---------- + component_list : list[Component] + R, L, and C components defining the branches. + frequency : float + Frequency in Hz. + + Returns + ------- + np.ndarray + Diagonal matrix of branch admittances. R → 1/R, C → jωC, L → 1/(jωL). + """ + branch_admittance_list = [] + omega = 2 * np.pi * frequency + for comp in component_list: + if comp.element_type == "R": + branch_admittance_list.append(1.0 / comp.value) + elif comp.element_type == "L": + branch_admittance_list.append(1.0 / (1j * omega * comp.value)) + elif comp.element_type == "C": + branch_admittance_list.append(1j * omega * comp.value) + else: + raise ValueError(f"Unknown component type: {comp.element_type}") + return np.diag(branch_admittance_list) + + @staticmethod + def _build_incidence_matrix_and_branch_admittance_factory( + component_list: list[Component], + ) -> tuple[np.ndarray, NodeMapper, Callable[[float], np.ndarray]]: + """Build incidence matrix, node mapper, and callable for branch admittance matrix. + + Parameters + ---------- + component_list : list[Component] + R, L, and C components. Node indices are assigned via the returned + :class:`NodeMapper` in order of first appearance (no special treatment + for any node label). + + Returns + ------- + A : np.ndarray + Full incidence matrix (N_nodes × N_branches). Each column has +1 at + node_plus and -1 at node_minus. + node_mapper : NodeMapper + Maps node names to indices. + branch_admittance_at : callable + Callable that takes a frequency in Hz and returns the diagonal branch + admittance matrix. + """ + node_mapper = NodeMapper() + for comp in component_list: + comp._node_plus_idx = node_mapper.get_or_create_index(comp.node_plus) + comp._node_minus_idx = node_mapper.get_or_create_index(comp.node_minus) + N_nodes = node_mapper.total_nodes() + B_components = len(component_list) + A = np.zeros((N_nodes, B_components)) + for idx, comp in enumerate(component_list): + A[comp._node_plus_idx, idx] = 1 + A[comp._node_minus_idx, idx] = -1 + + def branch_admittance_at(frequency: float) -> np.ndarray: + return CircuitImpedanceModel._create_branch_admittance_matrix(component_list, frequency) + + return A, node_mapper, branch_admittance_at + + @staticmethod + def _get_effective_admittance( + component_list: list[Component], + frequencies: np.ndarray | list[float], + port_plus_node: str = "1", + port_minus_node: str = "0", + ) -> np.ndarray: + """Compute driving-point admittance at each frequency for a one-port network. + + Uses the reduced nodal admittance matrix (reference node row removed) and Schur + complement to eliminate internal nodes, leaving the one-port admittance between + port_plus_node and port_minus_node. + + Parameters + ---------- + component_list : list[Component] + R, L, and C components defining the network. + frequencies : np.ndarray or list[float] + Frequencies in Hz at which to evaluate the admittance. + port_plus_node : str, optional + Name of the port's positive node (default ``"1"``). + port_minus_node : str, optional + Name of the port's reference (negative) node (default ``"0"``). + May be any node in the circuit (e.g. ``"0"``, ``"GND"``, or ``"2"``). + + Returns + ------- + np.ndarray + Complex driving-point admittance at each frequency, same length as ``frequencies``. + + Raises + ------ + ValueError + If ``port_plus_node`` and ``port_minus_node`` are the same, or if either + port node is not in the circuit (not an endpoint of any R, L, or C component). + """ + frequencies = np.atleast_1d(np.asarray(frequencies, dtype=float)) + Y_LE = np.zeros(len(frequencies), dtype=complex) + + A, node_mapper, branch_admittance_at = ( + CircuitImpedanceModel._build_incidence_matrix_and_branch_admittance_factory( + component_list + ) + ) + idx_plus = node_mapper.lookup_index(port_plus_node) + idx_minus = node_mapper.lookup_index(port_minus_node) + n_nodes = node_mapper.total_nodes() + + if idx_plus == idx_minus: + raise ValueError( + "Port nodes must be distinct: port_plus_node and port_minus_node cannot be the same." + ) + # Reduced matrix: remove the reference (port_minus) row. Reduced index for full node i + # (i != idx_minus) is i if i < idx_minus else i - 1. + n_red = n_nodes - 1 + idx_plus_red = idx_plus if idx_plus < idx_minus else idx_plus - 1 + eliminate_red = [i for i in range(n_red) if i != idx_plus_red] + A_red = np.delete(A, idx_minus, axis=0) + + for k, freq in enumerate(frequencies): + Y_branch = branch_admittance_at(float(freq)) + Y_red = A_red @ Y_branch @ A_red.T + if len(eliminate_red) == 0: + Y_LE[k] = Y_red[idx_plus_red, idx_plus_red] + else: + # Schur complement: keep only plus node, eliminate other non-reference nodes. + Y_bb = Y_red[np.ix_(eliminate_red, eliminate_red)] + Y_ab = Y_red[np.ix_([idx_plus_red], eliminate_red)] + Y_ba = Y_red[np.ix_(eliminate_red, [idx_plus_red])] + Y_aa_pp = Y_red[idx_plus_red, idx_plus_red] + Y_LE[k] = Y_aa_pp - (Y_ab @ np.linalg.solve(Y_bb, Y_ba))[0, 0] + + return Y_LE + + @staticmethod + def _admittance_to_eps_data( + frequencies: np.ndarray, + Y_complex: np.ndarray, + ) -> np.ndarray: + """Convert engineering-convention admittance Y(f) to equivalent complex permittivity. + + The conversion uses the relationship between admittance and the equivalent + dispersive medium used in FDTD (Pereda et al., IEEE TMTT 1999): + + .. math:: + + \\epsilon(\\omega) = 1 + \\frac{j \\, \\Delta \\, Y^*(\\omega)} + {\\omega \\, \\epsilon_0} + + with :math:`\\Delta = 1` here. Geometric scaling is applied later in + :meth:`AdmittanceNetwork._to_medium` when the model is used in a + :class:`~tidy3d.LinearLumpedElement` (via ``scaling_factor``). + + Parameters + ---------- + frequencies : np.ndarray + Frequencies in Hz (must be positive). + Y_complex : np.ndarray + Complex admittance at each frequency in engineering convention + (e.g. :math:`Y_C = j\\omega C`, :math:`Y_L = 1/(j\\omega L)`). + + Returns + ------- + np.ndarray + Complex permittivity array (same length as *frequencies*). + """ + frequencies = np.asarray(frequencies, dtype=float) + Y_complex = np.asarray(Y_complex, dtype=complex) + omega = 2 * np.pi * frequencies + return 1.0 + 1j * np.conj(Y_complex) / (omega * EPSILON_0) + + @staticmethod + def _fit_admittance_to_pole_residue( + frequencies: np.ndarray, + Y_complex: np.ndarray, + min_num_poles: int = 1, + max_num_poles: int = 5, + tolerance_rms: float = 1e-5, + show_progress: bool = True, + ) -> tuple[PoleResidue, float]: + """Fit admittance Y(f) to a :class:`~tidy3d.PoleResidue` medium via the dispersion fitter. + + Converts the engineering-convention admittance to equivalent permittivity + (see :meth:`_admittance_to_eps_data` with :math:`\\Delta=1`), then fits with the + standard dispersion fitter. Geometric scaling is applied in + :meth:`AdmittanceNetwork._to_medium` when the model is used in a + :class:`~tidy3d.LinearLumpedElement`. + + This approach has several advantages over fitting Y directly: + + * **Correct symmetry** -- permittivity has Hermitian symmetry + (:math:`\\epsilon(-\\omega) = \\epsilon^*(\\omega)`), matching the + conjugate-pair pole-residue model. + * **Correct passivity** -- the fitter's built-in passivity enforcement + (Im[eps] >= 0) directly ensures admittance passivity (Re[Y] >= 0). + * **No intermediate polynomial** -- bypasses the ``AdmittanceNetwork`` + ``(a, b)`` representation and its non-negative-coefficient constraint. + + Parameters + ---------- + frequencies : np.ndarray + Frequencies in Hz (must be positive). + Y_complex : np.ndarray + Complex admittance at each frequency in engineering convention. + min_num_poles : int + Minimum number of poles in the model. + max_num_poles : int + Maximum number of poles in the model. + tolerance_rms : float + Weighted RMS error below which the fit is considered successful. + show_progress : bool + Whether to show fitter progress bar. + + Returns + ------- + tuple[PoleResidue, float] + The fitted pole-residue medium and the weighted RMS error. + + Raises + ------ + ValueError + If ``frequencies`` is empty, lengths of ``frequencies`` and ``Y_complex`` differ, or any frequency is non-positive. + """ + frequencies = np.asarray(frequencies, dtype=float) + Y_complex = np.asarray(Y_complex, dtype=complex) + if frequencies.size == 0: + raise ValueError("frequencies must not be empty.") + if frequencies.size != Y_complex.size: + raise ValueError("frequencies and Y_complex must have the same length.") + if np.any(frequencies <= 0): + raise ValueError("All frequencies must be positive.") + + omega = 2 * np.pi * frequencies + eps_data = CircuitImpedanceModel._admittance_to_eps_data(frequencies, Y_complex) + + # Scale factor for numerical conditioning: normalize max(omega) to ~1 + scale_factor = 1.0 / (np.max(omega) + 1e-30) + + advanced_param = AdvancedFastFitterParam(show_progress=show_progress) + + (eps_inf, poles, residues), rms = fit( + omega_data=omega, + resp_data=eps_data, + min_num_poles=min_num_poles, + max_num_poles=max_num_poles, + resp_inf=None, + tolerance_rms=tolerance_rms, + scale_factor=scale_factor, + advanced_param=advanced_param, + ) + + # Build PoleResidue from fitter output + pole_pairs = tuple((complex(a), complex(c)) for a, c in zip(poles, residues)) + medium = PoleResidue(eps_inf=float(eps_inf), poles=pole_pairs) + + return medium, float(rms) + + @staticmethod + def _pole_residue_to_ab( + resp_inf: float, + poles: np.ndarray, + residues: np.ndarray, + ) -> tuple[np.ndarray, np.ndarray]: + """Convert pole-residue form (from dispersion fitter) to polynomial (a, b) ascending order. + + Fitter model: Y(s) = resp_inf - sum_i [ r_i/(s + p_i) + r_i*/(s + p_i*) ], s = j*omega. + The fitter stores only one pole per conjugate pair. To reconstruct the full rational + function N(s)/D(s) we must expand to ALL Laplace poles: + + * Complex pole p_i → two Laplace poles at -p_i and -conj(p_i) with residues -r_i and + -conj(r_i). + * Real pole p_i → a single Laplace pole at -p_i with residue -2·Re(r_i) (the conjugate + pair collapses into one pole with doubled weight). + """ + from scipy import signal + + poles = np.atleast_1d(poles) + residues = np.atleast_1d(residues) + if len(poles) == 0: + a = np.array([float(resp_inf)]) + b = np.array([1.0]) + return a, b + + # Expand to full Laplace partial-fraction form: H(s) = k + sum r_j/(s - p_j) + p_laplace = [] + r_laplace = [] + for p, r in zip(poles, residues): + if np.isreal(p): + # Conjugate pair collapses: -r/(s+p) - conj(r)/(s+p) = -2*Re(r)/(s+p) + p_laplace.append(-complex(p).real) + r_laplace.append(-2.0 * complex(r).real) + else: + # Two distinct Laplace poles + p_laplace.append(-complex(p)) + r_laplace.append(-complex(r)) + p_laplace.append(-np.conj(complex(p))) + r_laplace.append(-np.conj(complex(r))) + p_laplace = np.array(p_laplace, dtype=complex) + r_laplace = np.array(r_laplace, dtype=complex) + + k = np.array([resp_inf]) # direct term + # invres(r, p, k) -> (b_num, a_den) with H = b(s)/a(s) in descending order + b_desc, a_desc = signal.invres(r_laplace, p_laplace, k, tol=1e-12) + # Our convention: Y = a(s)/b(s), ascending order + a_asc = np.real(np.flip(b_desc)) + b_asc = np.real(np.flip(a_desc)) + b0 = b_asc[0] + if np.abs(b0) < 1e-20: + b0 = 1.0 + b_asc = b_asc / b0 + a_asc = a_asc / b0 + return a_asc, b_asc + + @staticmethod + def _eps_medium_to_admittance_ab( + medium: PoleResidue, + ) -> tuple[tuple[float, ...], tuple[float, ...]]: + """Extract admittance (a, b) coefficients from a PoleResidue ε-model. + + Uses the relation :math:`\\epsilon - 1 = -a(-s)/(s \\epsilon_0 b(-s))` between + the permittivity model and the admittance rational :math:`Y(s) = a(s)/b(s)`. + Coefficients are clipped to non-negative for :class:`CircuitImpedanceModel`. + + Parameters + ---------- + medium : PoleResidue + Fitted pole-residue medium representing the equivalent permittivity. + + Returns + ------- + tuple[tuple[float, ...], tuple[float, ...]] + ``(a, b)`` numerator and denominator coefficients in ascending monomial order. + """ + + poles_arr = np.array([complex(p) for p, _ in medium.poles]) + resid_arr = np.array([complex(c) for _, c in medium.poles]) + P, Q = CircuitImpedanceModel._pole_residue_to_ab( + float(medium.eps_inf) - 1.0, poles_arr, resid_arr + ) + P = np.asarray(P) + Q = np.asarray(Q) + b = np.array([((-1) ** k) * Q[k] for k in range(len(Q))], dtype=float) + a = np.zeros(len(P) + 1, dtype=float) + for k in range(1, len(a)): + a[k] = ((-1) ** (k + 1)) * EPSILON_0 * P[k - 1] + if len(b) == 0: + raise ValueError( + "Pole-residue to admittance conversion produced empty denominator; " + "try different fit frequencies or pole count." + ) + if len(a) == 0: + raise ValueError( + "Pole-residue to admittance conversion produced empty numerator; " + "try different fit frequencies or pole count." + ) + b0 = b[0] if abs(b[0]) > 1e-30 else 1.0 + a = a / b0 + b = b / b0 + # Clip to non-negative for CircuitImpedanceModel + a = np.maximum(np.real(a), 0.0) + b = np.maximum(np.real(b), 0.0) + a = np.trim_zeros(a, "b") + b = np.trim_zeros(b, "b") + if len(a) == 0 or len(b) == 0: + raise ValueError( + "Pole-residue to admittance conversion produced degenerate (all-zero) coefficients after " + "normalization; try different fit frequencies or pole count." + ) + if not np.all(np.isfinite(a)) or not np.all(np.isfinite(b)): + raise ValueError( + "Pole-residue to admittance conversion produced non-finite (NaN or inf) coefficients; " + "try different fit frequencies or pole count." + ) + return (tuple(float(x) for x in a), tuple(float(x) for x in b)) + + @classmethod + def from_spice_file( + cls, + spice_file: str | Path, + frequencies: np.ndarray | list[float], + port_plus_node: Optional[str] = None, + port_minus_node: Optional[str] = None, + min_num_poles: int = 1, + max_num_poles: int = 5, + tolerance_rms: float = 1e-5, + show_progress: bool = False, + ) -> Self: + """Build a :class:`CircuitImpedanceModel` from a SPICE netlist file. + + The netlist is parsed for R, C, and L elements; the port is taken from the single + voltage source (V) if present, otherwise from the first element's two nodes. + The driving-point admittance is fitted via the Y→ε→pole-residue path to obtain + stable ``(a, b)`` coefficients. + + Parameters + ---------- + spice_file : str or Path + Path to the SPICE netlist file. + frequencies : np.ndarray or list[float] + Frequencies in Hz at which to evaluate and fit the admittance. + port_plus_node : str, optional + Override port positive node (default: from netlist). + port_minus_node : str, optional + Override port negative node (default: from netlist). + min_num_poles : int, optional + Minimum number of poles for the dispersion fitter (default 1). + max_num_poles : int, optional + Maximum number of poles (default 5). + tolerance_rms : float, optional + Target weighted RMS error for the fit (default 1e-5). + show_progress : bool, optional + Whether to show the fitter progress bar (default False). + + Returns + ------- + CircuitImpedanceModel + Model with (a, b) coefficients for use in :class:`LinearLumpedElement`. + """ + component_list, port_p, port_m = cls._parse_spice_file(spice_file) + if port_plus_node is not None: + port_p = port_plus_node + if port_minus_node is not None: + port_m = port_minus_node + return cls.from_component_list( + component_list, + frequencies, + port_plus_node=port_p, + port_minus_node=port_m, + min_num_poles=min_num_poles, + max_num_poles=max_num_poles, + tolerance_rms=tolerance_rms, + show_progress=show_progress, + ) + + @classmethod + def from_component_list( + cls, + component_list: list[Component], + frequencies: np.ndarray | list[float], + port_plus_node: str = "1", + port_minus_node: str = "0", + min_num_poles: int = 1, + max_num_poles: int = 5, + tolerance_rms: float = 1e-5, + show_progress: bool = False, + ) -> Self: + """Build a :class:`CircuitImpedanceModel` from a list of R/L/C components. + + Computes the driving-point admittance at the given frequencies, converts to + equivalent permittivity, fits a pole-residue ε model with the dispersion fitter, + then extracts ``(a, b)`` coefficients for the admittance rational. + + Parameters + ---------- + component_list : list[Component] + R, L, and C components defining the one-port network. + frequencies : np.ndarray or list[float] + Frequencies in Hz at which to evaluate and fit the admittance. + port_plus_node : str, optional + Port positive node name (default ``"1"``). + port_minus_node : str, optional + Port reference (negative) node name (default ``"0"``). May be any node in the circuit. + min_num_poles : int, optional + Minimum number of poles for the dispersion fitter (default 1). + max_num_poles : int, optional + Maximum number of poles (default 5). + tolerance_rms : float, optional + Target weighted RMS error for the fit (default 1e-5). + show_progress : bool, optional + Whether to show the fitter progress bar (default False). + + Returns + ------- + CircuitImpedanceModel + Model with (a, b) coefficients for use in :class:`LinearLumpedElement`. + """ + frequencies = np.asarray(frequencies, dtype=float) + Y_complex = cls._get_effective_admittance( + component_list, + frequencies, + port_plus_node=port_plus_node, + port_minus_node=port_minus_node, + ) + medium, _ = cls._fit_admittance_to_pole_residue( + frequencies=frequencies, + Y_complex=Y_complex, + min_num_poles=min_num_poles, + max_num_poles=max_num_poles, + tolerance_rms=tolerance_rms, + show_progress=show_progress, + ) + a, b = cls._eps_medium_to_admittance_ab(medium) + return cls(a=a, b=b) + + @classmethod + def from_touchstone_file( + cls, + touchstone_file: str, + num_order: int = 2, + denom_order: int = 2, + ) -> Self: + """Build a :class:`CircuitImpedanceModel` from a Touchstone file. + + Not yet implemented. Use :meth:`from_spice_file` or :meth:`from_component_list` instead. + + Parameters + ---------- + touchstone_file : str + Path to the Touchstone file (e.g. .s1p). + num_order : int, optional + Numerator order for the rational fit (reserved for future use). + denom_order : int, optional + Denominator order for the rational fit (reserved for future use). + + Returns + ------- + CircuitImpedanceModel + Model with ``(a, b)`` coefficients (when implemented). + + Raises + ------ + NotImplementedError + Touchstone file support is not yet implemented. + """ + raise NotImplementedError( + "CircuitImpedanceModel.from_touchstone_file is not yet implemented. " + "Use CircuitImpedanceModel.from_spice_file or from_component_list." + ) + + +NetworkType = discriminated_union(Union[RLCNetwork, AdmittanceNetwork, CircuitImpedanceModel]) class LinearLumpedElement(RectangularLumpedElement): diff --git a/tidy3d/rf.py b/tidy3d/rf.py index 25230a7158..6aa4203765 100644 --- a/tidy3d/rf.py +++ b/tidy3d/rf.py @@ -17,6 +17,7 @@ # Lumped elements from tidy3d.components.lumped_element import ( AdmittanceNetwork, + CircuitImpedanceModel, CoaxialLumpedResistor, LinearLumpedElement, LumpedResistor, @@ -164,6 +165,7 @@ "BlackmanHarrisWindow", "BlackmanWindow", "ChebWindow", + "CircuitImpedanceModel", "CoaxialLumpedPort", "CoaxialLumpedResistor", "ComponentModelerDataType",