add mesh refinement option for WavePorts and validator to check if the grid is refined enough around WavePorts

Author: dmarek-flex

CHANGELOG.md

@@ -20,6 +20,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added `TriangleMesh.from_height_expression` class method to create a mesh from an analytical height function defined on a 2D grid and `TriangleMesh.from_height_grid` class method to create a mesh from height values sampled on a 2D grid.
 - Added heat sources with custom spatial dependence. It is now possible to add a `SpatialDataArray` as the `rate` in a `HeatSource`.
 - Added Transient Heat simulations. It is now possible to run transient Heat simulations. This can be done by specifying `analysis_spec` of `HeatChargeSimulation` object as `UnsteadyHeatAnalysis`. 
+- A `num_grid_cells` field to `WavePort`, which ensures that there are 5 grid cells across the port. Grid refinement can be disabled by passing `None` to `num_grid_cells`.
 
 ### Fixed
 - Fixed bug in broadband adjoint source creation when forward simulation had a pulse amplitude greater than 1 or a nonzero pulse phase.

tests/test_plugins/smatrix/terminal_component_modeler_def.py

@@ -304,7 +304,9 @@ def make_port(center, direction, type, name) -> Union[CoaxialLumpedPort, WavePor
                     normal_axis=2,
                     clockwise=direction != "+",
                 )
-
+            port_cells = None
+            if port_refinement:
+                port_cells = 5
             port = WavePort(
                 center=center,
                 size=[2 * Router, 2 * Router, 0],
@@ -314,6 +316,7 @@ def make_port(center, direction, type, name) -> Union[CoaxialLumpedPort, WavePor
                 mode_index=0,
                 voltage_integral=voltage_integral,
                 current_integral=current_integral,
+                num_grid_cells=port_cells,
             )
         return port
 

tests/test_plugins/smatrix/test_terminal_component_modeler.py

@@ -9,7 +9,11 @@
 import tidy3d as td
 from tidy3d.components.data.data_array import FreqDataArray
 from tidy3d.exceptions import SetupError, Tidy3dError, Tidy3dKeyError
-from tidy3d.plugins.microwave import CustomCurrentIntegral2D, VoltageIntegralAxisAligned
+from tidy3d.plugins.microwave import (
+    CurrentIntegralAxisAligned,
+    CustomCurrentIntegral2D,
+    VoltageIntegralAxisAligned,
+)
 from tidy3d.plugins.smatrix import (
     AbstractComponentModeler,
     CoaxialLumpedPort,
@@ -604,6 +608,69 @@ def test_wave_port_path_integral_validation():
     assert wave_port.bounds[0][2] > wave_port.voltage_integral.bounds[0][2]
 
 
+def test_wave_port_grid_validation(tmp_path):
+    """Ensure that 'num_grid_cells' is validated and works to ensure that the grid is refined around wave ports."""
+    size_port = [2, 2, 0]
+    center_port = [0, 0, -10]
+
+    voltage_path = VoltageIntegralAxisAligned(
+        center=(0.5, 0, -10),
+        size=(1.0, 0, 0),
+        extrapolate_to_endpoints=True,
+        snap_path_to_grid=True,
+        sign="+",
+    )
+
+    current_path = CurrentIntegralAxisAligned(
+        center=(0.5, 0, -10),
+        size=(0.25, 0.5, 0),
+        snap_contour_to_grid=True,
+        sign="+",
+    )
+
+    mode_spec = td.ModeSpec(num_modes=1, target_neff=1.8)
+
+    _ = WavePort(
+        center=center_port,
+        size=size_port,
+        name="wave_port_1",
+        mode_spec=mode_spec,
+        direction="+",
+        voltage_integral=voltage_path,
+        current_integral=current_path,
+        num_grid_cells=None,
+    )
+
+    with pytest.raises(pd.ValidationError):
+        _ = WavePort(
+            center=center_port,
+            size=size_port,
+            name="wave_port_1",
+            mode_spec=mode_spec,
+            direction="+",
+            voltage_integral=voltage_path,
+            current_integral=current_path,
+            num_grid_cells=2,
+        )
+
+    modeler = make_coaxial_component_modeler(
+        grid_spec=td.GridSpec.auto(wavelength=10e3),
+        port_refinement=True,
+        path_dir=str(tmp_path),
+        port_types=(WavePort, WavePort),
+    )
+    _ = modeler.sim_dict
+
+    modeler = make_coaxial_component_modeler(
+        grid_spec=td.GridSpec.auto(wavelength=10e3),
+        port_refinement=False,
+        path_dir=str(tmp_path),
+        port_types=(WavePort, WavePort),
+    )
+    with pytest.raises(SetupError):
+        _ = modeler.sim_dict
+
+
 def test_wave_port_to_mode_solver(tmp_path):
     """Checks that wave port can be converted to a mode solver."""
     modeler = make_coaxial_component_modeler(

tidy3d/plugins/smatrix/component_modelers/terminal.py

@@ -19,7 +19,7 @@
 from tidy3d.components.types import Ax
 from tidy3d.components.viz import add_ax_if_none, equal_aspect
 from tidy3d.constants import C_0, OHM
-from tidy3d.exceptions import Tidy3dError, Tidy3dKeyError, ValidationError
+from tidy3d.exceptions import SetupError, Tidy3dError, Tidy3dKeyError, ValidationError
 from tidy3d.log import log
 from tidy3d.plugins.smatrix.data.terminal import PortDataArray, TerminalPortDataArray
 from tidy3d.plugins.smatrix.ports.base_lumped import AbstractLumpedPort
@@ -142,9 +142,17 @@ def sim_dict(self) -> dict[str, Simulation]:
             port.to_load(snap_center=snap_centers[port.name]) for port in self._lumped_ports
         ]
 
+        # Add mesh overrides for any wave ports present
+        mesh_overrides = list(sim_wo_source.grid_spec.override_structures)
+        for wave_port in self._wave_ports:
+            if wave_port.num_grid_cells is not None:
+                mesh_overrides.extend(wave_port.to_mesh_overrides())
+        new_grid_spec = sim_wo_source.grid_spec.updated_copy(override_structures=mesh_overrides)
+
         update_dict = {
             "monitors": new_mnts,
             "lumped_elements": new_lumped_elements,
+            "grid_spec": new_grid_spec,
         }
 
         # This is the new default simulation will all shared components added
@@ -158,10 +166,6 @@ def sim_dict(self) -> dict[str, Simulation]:
             task_name = self._task_name(port=port)
             sim_dict[task_name] = sim_wo_source.updated_copy(sources=[port_source])
 
-        # Check final simulations for grid size at lumped ports
-        for _, sim in sim_dict.items():
-            TerminalComponentModeler._check_grid_size_at_ports(sim, self._lumped_ports)
-
         # Now, create simulations with wave port sources and mode solver monitors for computing port modes
         for wave_port in self._wave_ports:
             # Source is placed just before the field monitor of the port
@@ -174,6 +178,12 @@ def sim_dict(self) -> dict[str, Simulation]:
 
             task_name = self._task_name(port=wave_port)
             sim_dict[task_name] = sim_wo_source.copy(update=update_dict)
+
+        # Check final simulations for grid size at ports
+        for _, sim in sim_dict.items():
+            TerminalComponentModeler._check_grid_size_at_ports(sim, self._lumped_ports)
+            TerminalComponentModeler._check_grid_size_at_wave_ports(sim, self._wave_ports)
+
         return sim_dict
 
     @cached_property
@@ -242,12 +252,34 @@ def _validate_radiation_monitors(cls, val, values):
         return val
 
     @staticmethod
-    def _check_grid_size_at_ports(simulation: Simulation, ports: list[Union[AbstractLumpedPort]]):
+    def _check_grid_size_at_ports(
+        simulation: Simulation, ports: list[Union[LumpedPort, CoaxialLumpedPort]]
+    ):
         """Raises :class:`.SetupError` if the grid is too coarse at port locations"""
         yee_grid = simulation.grid.yee
         for port in ports:
             port._check_grid_size(yee_grid)
 
+    @staticmethod
+    def _check_grid_size_at_wave_ports(simulation: Simulation, ports: list[WavePort]):
+        """Raises :class:`.SetupError` if the grid is too coarse at port locations"""
+        for port in ports:
+            disc_grid = simulation.discretize(port)
+            check_axes = port.transverse_axes
+            msg_header = f"'WavePort' '{port.name}' "
+            for axis in check_axes:
+                sim_size = simulation.size[axis]
+                dim_cells = disc_grid.num_cells[axis]
+                if sim_size > 0 and dim_cells <= 2:
+                    small_dim = "xyz"[axis]
+                    raise SetupError(
+                        msg_header + f"is too small along the "
+                        f"'{small_dim}' axis. Less than '3' grid cells were detected. "
+                        "Please ensure that the port's 'num_grid_cells' is not 'None'. "
+                        "You also may need to use an 'AutoGrid' or `QuasiUniformGrid` "
+                        "for the simulation passed to the 'TerminalComponentModeler'."
+                    )
+
     def compute_power_wave_amplitudes_at_each_port(
         self, port_reference_impedances: PortDataArray, sim_data: SimulationData
     ) -> tuple[PortDataArray, PortDataArray]:

tidy3d/plugins/smatrix/ports/wave.py

@@ -18,14 +18,18 @@
 from tidy3d.components.simulation import Simulation
 from tidy3d.components.source.field import ModeSource, ModeSpec
 from tidy3d.components.source.time import GaussianPulse
-from tidy3d.components.types import Direction, FreqArray
+from tidy3d.components.structure import MeshOverrideStructure
+from tidy3d.components.types import Axis, Direction, FreqArray
 from tidy3d.constants import fp_eps
 from tidy3d.exceptions import ValidationError
 from tidy3d.plugins.microwave import CurrentIntegralTypes, ImpedanceCalculator, VoltageIntegralTypes
 from tidy3d.plugins.mode import ModeSolver
 
 from .base_terminal import AbstractTerminalPort
 
+DEFAULT_WAVE_PORT_NUM_CELLS = 5
+MIN_WAVE_PORT_NUM_CELLS = 3
+
 
 class WavePort(AbstractTerminalPort, Box):
     """Class representing a single wave port"""
@@ -63,6 +67,15 @@ class WavePort(AbstractTerminalPort, Box):
         description="Definition of current integral used to compute current and the characteristic impedance.",
     )
 
+    num_grid_cells: Optional[int] = pd.Field(
+        DEFAULT_WAVE_PORT_NUM_CELLS,
+        ge=MIN_WAVE_PORT_NUM_CELLS,
+        title="Number of Grid Cells",
+        description="Number of mesh grid cells in the transverse plane of the `WavePort`. "
+        "Used in generating the suggested list of :class:`.MeshOverrideStructure` objects. "
+        "Must be greater than or equal to 3. When set to `None`, no grid refinement is performed.",
+    )
+
     def _mode_voltage_coefficients(self, mode_data: ModeData) -> FreqModeDataArray:
         """Calculates scaling coefficients to convert mode amplitudes
         to the total port voltage.
@@ -88,10 +101,16 @@ def _mode_current_coefficients(self, mode_data: ModeData) -> FreqModeDataArray:
         return current_coeffs.squeeze()
 
     @cached_property
-    def injection_axis(self):
+    def injection_axis(self) -> Axis:
         """Injection axis of the port."""
         return self.size.index(0.0)
 
+    @cached_property
+    def transverse_axes(self) -> tuple[Axis, Axis]:
+        """Transverse axes of the port."""
+        _, trans_axes = Box.pop_axis([0, 1, 2], self.injection_axis)
+        return trans_axes
+
     @cached_property
     def _mode_monitor_name(self) -> str:
         """Return the name of the :class:`.ModeMonitor` associated with this port."""
@@ -186,6 +205,25 @@ def compute_port_impedance(
         impedance_array = impedance_calc.compute_impedance(mode_data)
         return impedance_array
 
+    def to_mesh_overrides(self) -> list[MeshOverrideStructure]:
+        """Creates a list of :class:`.MeshOverrideStructure` for mesh refinement in the transverse
+        plane of the port. The mode source requires at least 3 grid cells in the transverse
+        dimensions, so these mesh overrides will be added to the simulation to ensure that this
+        requirement is satisfied.
+        """
+        dl = [None] * 3
+        for trans_axis in self.transverse_axes:
+            dl[trans_axis] = self.size[trans_axis] / self.num_grid_cells
+
+        return [
+            MeshOverrideStructure(
+                geometry=Box(center=self.center, size=self.size),
+                dl=dl,
+                shadow=False,
+                priority=-1,
+            )
+        ]
+
     @pd.validator("voltage_integral", "current_integral")
     def _validate_path_integrals_within_port(cls, val, values):
         """Raise ``ValidationError`` when the supplied path integrals are not within the port bounds."""
@@ -213,7 +251,7 @@ def _check_voltage_or_current(cls, val, values):
         return val
 
     @pd.validator("current_integral", always=True)
-    def validate_current_integral_sign(cls, val, values):
+    def _validate_current_integral_sign(cls, val, values):
         """
         Validate that the sign of ``current_integral`` matches the port direction.
         """