Add more robust gauge choice for high-order modes

CHANGELOG.md

@@ -31,6 +31,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - `MethodRandom` and `MethodRandomCustom` have been removed from the Design plugin, and `DesignSpace.run_batch` has been superseded by `.run`.
 - Design plugin has been significantly reworked to improve ease of use and allow for new optimization methods.
 - Behavior of `FieldProjector` now matches the server-side computation, which does not truncate the integration surface when it extends into PML regions.
+- Mode solver fields are more consistently normalized with respect to grid-dependent sign inversions in high order modes.
 
 ### Fixed
 - Significant speedup for field projection computations.

tests/test_plugins/test_mode_solver.py

@@ -1034,3 +1034,63 @@ def test_modes_eme_sim(mock_remote_api, local):
         _ = msweb.run(solver.to_fdtd_mode_solver())
 
     _ = solver.reduced_simulation_copy
+
+
+def make_high_order_mode_solver(sign, dim=3):
+    waveguide = td.Structure(
+        geometry=td.Box(size=(td.inf, 0.6, 0.2)),
+        medium=td.Medium(permittivity=3.47**2),
+    )
+
+    refine_box = td.MeshOverrideStructure(
+        geometry=td.Box(center=(0, sign * 0.3, 0), size=(td.inf, 0.1, 0.3)),
+        dl=[None, 0.02, 0.02],
+    )
+
+    pml = td.Boundary(plus=td.PML(), minus=td.PML())
+    periodic = td.Boundary(plus=td.Periodic(), minus=td.Periodic())
+
+    sim = td.Simulation(
+        size=(10, 2.5, 1.5 if dim == 3 else 0),
+        grid_spec=td.GridSpec.auto(
+            min_steps_per_wvl=20, wavelength=1.55, override_structures=[refine_box]
+        ),
+        structures=[waveguide],
+        medium=td.Medium(permittivity=1.44**2),
+        boundary_spec=td.BoundarySpec(x=pml, y=pml, z=pml if dim == 3 else periodic),
+        run_time=1e-12,
+    )
+
+    plane = td.Box(center=(0, 0, 0), size=(0, 2.5, 1.5))
+    freq0 = td.C_0 / 1.55
+    num_modes = 3
+
+    mode_spec = td.ModeSpec(
+        num_modes=num_modes,
+        target_neff=3.47,
+        group_index_step=False,
+    )
+
+    return ModeSolver(
+        simulation=sim,
+        plane=plane,
+        mode_spec=mode_spec,
+        freqs=[freq0],
+    )
+
+
+def test_high_order_mode_normalization():
+    # 3D simulation
+    ms1 = make_high_order_mode_solver(1)
+    ms2 = make_high_order_mode_solver(-1)
+    overlap = ms1.data.outer_dot(ms2.data).isel(mode_index_0=2, mode_index_1=2).values.item()
+    assert abs(1 - overlap) < 1e-3
+
+    # 2D simulation
+    ms1 = make_high_order_mode_solver(1, 2)
+    values = ms1.data.Ez.isel(mode_index=2).values.squeeze().real
+    assert (values[: values.size // 3] > 0).all()
+
+    ms2 = make_high_order_mode_solver(-1, 2)
+    values = ms2.data.Ez.isel(mode_index=2).values.squeeze().real
+    assert (values[: values.size // 3] > 0).all()

tidy3d/plugins/mode/mode_solver.py

@@ -37,6 +37,7 @@
     ArrayComplex3D,
     ArrayComplex4D,
     ArrayFloat1D,
+    ArrayFloat2D,
     Ax,
     Axis,
     Direction,
@@ -679,12 +680,16 @@ def _solve_all_freqs_relative(
 
         return n_complex, fields, eps_spec
 
-    def _postprocess_solver_fields(self, solver_fields):
+    def _postprocess_solver_fields(self, solver_fields, coords):
         """Postprocess `solver_fields` from `compute_modes` to proper coordinate"""
         fields = {key: [] for key in ("Ex", "Ey", "Ez", "Hx", "Hy", "Hz")}
+        diff_coords = (np.diff(coords[0]), np.diff(coords[1]))
+
         for mode_index in range(self.mode_spec.num_modes):
             # Get E and H fields at the current mode_index
-            ((Ex, Ey, Ez), (Hx, Hy, Hz)) = self._process_fields(solver_fields, mode_index)
+            ((Ex, Ey, Ez), (Hx, Hy, Hz)) = self._process_fields(
+                solver_fields, mode_index, diff_coords
+            )
 
             # Note: back in original coordinates
             fields_mode = {"Ex": Ex, "Ey": Ey, "Ez": Ez, "Hx": Hx, "Hy": Hy, "Hz": Hz}
@@ -718,7 +723,7 @@ def _solve_single_freq(
             direction=self.direction,
         )
 
-        fields = self._postprocess_solver_fields(solver_fields)
+        fields = self._postprocess_solver_fields(solver_fields, coords)
         return n_complex, fields, eps_spec
 
     def _rotate_field_coords_inverse(self, field: FIELD) -> FIELD:
@@ -770,16 +775,59 @@ def _solve_single_freq_relative(
             solver_basis_fields=solver_basis_fields,
         )
 
-        fields = self._postprocess_solver_fields(solver_fields)
+        fields = self._postprocess_solver_fields(solver_fields, coords)
         return n_complex, fields, eps_spec
 
     def _rotate_field_coords(self, field: FIELD) -> FIELD:
         """Move the propagation axis=z to the proper order in the array."""
         f_x, f_y, f_z = np.moveaxis(field, source=3, destination=1 + self.normal_axis)
         return np.stack(self.plane.unpop_axis(f_z, (f_x, f_y), axis=self.normal_axis), axis=0)
 
+    @staticmethod
+    def _weighted_coord_max(
+        array: ArrayFloat2D, u: ArrayFloat1D, v: ArrayFloat1D
+    ) -> Tuple[int, int]:
+        m = array * u.reshape(-1, 1)
+        i = np.arange(array.shape[0])
+        i = int(0.5 + (m * i.reshape(-1, 1)).sum() / m.sum())
+        m = array * v
+        j = np.arange(array.shape[1])
+        j = int(0.5 + (m * j).sum() / m.sum())
+        return i, j
+
+    @staticmethod
+    def _inverted_gauge(e_field: FIELD, diff_coords: Tuple[ArrayFloat1D, ArrayFloat1D]) -> bool:
+        """Check if the lower xy region of the mode has a negative sign."""
+        dx, dy = diff_coords
+        e_x, e_y = e_field[:2, :, :, 0]
+        e_x = e_x.real
+        e_y = e_y.real
+        e = e_x if np.abs(e_x).max() > np.abs(e_y).max() else e_y
+        abs_e = np.abs(e)
+        e_2 = abs_e**2
+        i, j = e.shape
+        while i > 0 and j > 0:
+            if (e[:i, :j] > 0).all():
+                return False
+            elif (e[:i, :j] < 0).all():
+                return True
+            else:
+                threshold = abs_e[:i, :j].max() * 0.5
+                i, j = ModeSolver._weighted_coord_max(e_2[:i, :j], dx[:i], dy[:j])
+                if abs(e[i, j]) >= threshold:
+                    return e[i, j] < 0
+                # Do not close the window for 1D mode solvers
+                if e.shape[0] == 1:
+                    i = 1
+                elif e.shape[1] == 1:
+                    j = 1
+        return False
+
     def _process_fields(
-        self, mode_fields: ArrayComplex4D, mode_index: pydantic.NonNegativeInt
+        self,
+        mode_fields: ArrayComplex4D,
+        mode_index: pydantic.NonNegativeInt,
+        diff_coords: Tuple[ArrayFloat1D, ArrayFloat1D],
     ) -> Tuple[FIELD, FIELD]:
         """Transform solver fields to simulation axes and set gauge."""
 
@@ -792,6 +840,13 @@ def _process_fields(
         E *= np.exp(-1j * phi)
         H *= np.exp(-1j * phi)
 
+        # High-order modes with opposite-sign lobes will show inconsistent sign, heavily dependent
+        # on the exact local grid to choose the previous gauge. This method flips the gauge sign
+        # when necessary to make it more consistent.
+        if ModeSolver._inverted_gauge(E, diff_coords):
+            E *= -1
+            H *= -1
+
         # Rotate back to original coordinates
         (Ex, Ey, Ez) = self._rotate_field_coords(E)
         (Hx, Hy, Hz) = self._rotate_field_coords(H)