From b4954bcd975621dbba7c4e94dd88ed094e5ddb9d Mon Sep 17 00:00:00 2001 From: momchil Date: Mon, 21 Oct 2024 16:14:18 +0200 Subject: [PATCH] Fix mode simulations with a bend when the grid is not symmetric --- CHANGELOG.md | 1 + docs/notebooks | 2 +- tests/test_plugins/test_mode_solver.py | 52 ++++++++++++++++++++++++++ tidy3d/plugins/mode/mode_solver.py | 8 ++++ tidy3d/plugins/mode/solver.py | 14 ++++++- tidy3d/plugins/mode/transforms.py | 6 +-- 6 files changed, 78 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7864949c7..609005d9e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -45,6 +45,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - Bug in `Simulation.subsection` where lumped elements were not being correctly removed. - Bug when adding 2D structures to the `Simulation` that are composed of multiple overlapping polygons. - Fields stored in `ModeMonitor` objects were computed colocated to the grid boundaries, but the built-in `ModeMonitor.colocate=False` was resulting in wrong results in some cases, most notably if symmetries are also involved. +- Small inaccuracy when applying a mode solver `bend_radius` when the simulation grid is not symmetric w.r.t. the mode plane center. Previously, the radius was defined w.r.t. the middle grid coordinate, while now it is correctly applied w.r.t. the plane center. ## [2.7.2] - 2024-08-07 diff --git a/docs/notebooks b/docs/notebooks index a2f1abe29..3ad84c09b 160000 --- a/docs/notebooks +++ b/docs/notebooks @@ -1 +1 @@ -Subproject commit a2f1abe2980bfae56b8bdeb969a96e6bed5f4a72 +Subproject commit 3ad84c09bfa80ceb3cb4cdf894472430acf2a1fb diff --git a/tests/test_plugins/test_mode_solver.py b/tests/test_plugins/test_mode_solver.py index 7a0fa38a6..851b0ea1d 100644 --- a/tests/test_plugins/test_mode_solver.py +++ b/tests/test_plugins/test_mode_solver.py @@ -646,6 +646,58 @@ def test_mode_solver_angle_bend(): _ = ms.to_monitor(freqs=np.array([1.0, 2.0]) * 1e12, name="mode_mnt") +def test_mode_bend_radius(): + """Test that the bend radius is correctly applied to the center of the mode plane in the case + of an auto-grid that is not symmetric w.r.t. that center, and that nominally identical + waveguides produce the same modes when the bend center if shifted away from the mode plane + center, after a post-processing transformation to correct for that.""" + + simulation = td.Simulation( + size=(10, 10, 10), + grid_spec=td.GridSpec(wavelength=1.0), + # grid_spec=td.GridSpec.uniform(dl=0.04), + structures=[WAVEGUIDE], + run_time=1e-12, + ) + mode_spec1 = td.ModeSpec( + num_modes=3, + bend_radius=5, + bend_axis=1, + ) + + # plane centered on the waveguide center + plane1 = td.Box(center=(0, 0, 0), size=(3, 0, 2)) + + # plane centered away from the waveguide center + plane2 = td.Box(center=(0.5, 0, 0), size=(4, 0, 2)) + # mode spec with a radius such that the radius w.r.t. the waveguide center should be the same + mode_spec2 = mode_spec1.updated_copy(bend_radius=5.5) + + ms1 = ModeSolver( + simulation=simulation, + plane=plane1, + mode_spec=mode_spec1, + freqs=[td.C_0 / 1.0], + ) + ms2 = ms1.updated_copy(plane=plane2, mode_spec=mode_spec2) + data1 = ms1.solve() + data2 = ms2.solve() + + print(data1.n_complex) + print(data2.n_complex * 5 / 5.5) + + # Plot fields + # _, ax = plt.subplots(3, 2) + # for mode_index in range(3): + # ms1.plot_field("Ex", ax=ax[mode_index, 0], mode_index=mode_index) + # ms2.plot_field("Ex", ax=ax[mode_index, 1], mode_index=mode_index) + # plt.show() + + # The mode field dependence is E0 * exp(1j * n * R * k0 * phi) so to switch from one radius + # to another we have n' * R' = n * R -> n' = n * R / R' + assert np.allclose(data1.n_complex, data2.n_complex * 5.5 / 5.0) + + def test_mode_solver_2D(): """Run mode solver in 2D simulations.""" mode_spec = td.ModeSpec( diff --git a/tidy3d/plugins/mode/mode_solver.py b/tidy3d/plugins/mode/mode_solver.py index 60baa818e..f3bc2fa2e 100644 --- a/tidy3d/plugins/mode/mode_solver.py +++ b/tidy3d/plugins/mode/mode_solver.py @@ -179,6 +179,12 @@ def normal_axis(self) -> Axis: """Axis normal to the mode plane.""" return self.plane.size.index(0.0) + @staticmethod + def plane_center_tangential(plane) -> tuple[float, float]: + """Mode lane center in the tangential axes.""" + _, plane_center = plane.pop_axis(plane.center, plane.size.index(0.0)) + return plane_center + @cached_property def solver_symmetry(self) -> Tuple[Symmetry, Symmetry]: """Get symmetry for solver for propagation along self.normal axis.""" @@ -723,6 +729,7 @@ def _solve_single_freq( symmetry=symmetry, direction=self.direction, precision=self._precision, + plane_center=self.plane_center_tangential(self.plane), ) fields = self._postprocess_solver_fields( @@ -778,6 +785,7 @@ def _solve_single_freq_relative( direction=self.direction, solver_basis_fields=solver_basis_fields, precision=self._precision, + plane_center=self.plane_center_tangential(self.plane), ) fields = self._postprocess_solver_fields( diff --git a/tidy3d/plugins/mode/solver.py b/tidy3d/plugins/mode/solver.py index 71bb053e8..1b2d35c28 100644 --- a/tidy3d/plugins/mode/solver.py +++ b/tidy3d/plugins/mode/solver.py @@ -49,6 +49,7 @@ def compute_modes( symmetry=(0, 0), direction="+", solver_basis_fields=None, + plane_center: tuple[float, float] = None, ) -> Tuple[Numpy, Numpy, EpsSpecType]: """ Solve for the modes of a waveguide cross-section. @@ -84,6 +85,10 @@ def compute_modes( Direction of mode propagation. solver_basis_fields If provided, solve for modes in this basis. + plane_center + The center of the mode plane along the tangential axes of the global simulation. Used + in case of bend modes to offset the coordinates correctly w.r.t. the bend radius, which + is assumed to refer to the distance from the bend center to the mode plane center. Returns ------- @@ -153,7 +158,14 @@ def compute_modes( new_coords, jac_e, jac_h = angled_transform(new_coords, angle_theta, angle_phi) if bend_radius is not None: - new_coords, jac_e_tmp, jac_h_tmp = radial_transform(new_coords, bend_radius, bend_axis) + if plane_center is None: + raise ValueError( + "When using a nonzero 'bend_radius', the 'plane_center' must also " + "be provided to correctly offset the coordinates." + ) + new_coords, jac_e_tmp, jac_h_tmp = radial_transform( + new_coords, bend_radius, bend_axis, plane_center + ) jac_e = np.einsum("ij...,jp...->ip...", jac_e_tmp, jac_e) jac_h = np.einsum("ij...,jp...->ip...", jac_h_tmp, jac_h) diff --git a/tidy3d/plugins/mode/transforms.py b/tidy3d/plugins/mode/transforms.py index b6ad24eb1..53372b395 100644 --- a/tidy3d/plugins/mode/transforms.py +++ b/tidy3d/plugins/mode/transforms.py @@ -11,7 +11,7 @@ import numpy as np -def radial_transform(coords, radius, bend_axis): +def radial_transform(coords, radius, bend_axis, plane_center): """Compute the new coordinates and the Jacobian of a polar coordinate transformation. After offsetting the plane such that its center is a distance of ``radius`` away from the center of curvature, we have, e.g. for ``bend_axis=='y'``: @@ -45,8 +45,8 @@ def radial_transform(coords, radius, bend_axis): norm_axis = 0 if bend_axis == 1 else 1 # Center the new coordinates such that the radius is at the center of the plane - u = coords[0] + (norm_axis == 0) * (radius - coords[0][Nx // 2]) - v = coords[1] + (norm_axis == 1) * (radius - coords[1][Ny // 2]) + u = coords[0] + (norm_axis == 0) * (radius - plane_center[0]) + v = coords[1] + (norm_axis == 1) * (radius - plane_center[1]) new_coords = (u, v) """The only nontrivial derivative is dwdz and it only depends on the coordinate in the