Fix mode simulations with a bend when the grid is not symmetric

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
 

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(

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(

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)
 

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