Bringing back symmetries into the mode solver

CHANGELOG.md

@@ -7,6 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
  - Clarified license terms to not include scripts written using the tidy3d python API.
+ - Simulation symmetries are now enabled but currently only affect the mode solver, if the mode plane lies on the simulation center and there's a symmetry.
+ - Validator that mode objects with symmetries are either entirely in the main quadrant, or lie on the symmetry axis.
 
 ### Changed
 - Fixed a bug in python 3.6 where polyslab vertices loaded as List of List.

requirements.txt

@@ -12,7 +12,6 @@ PyYAML
 boto3
 requests
 plotly==5.5.0
-chart-studio==1.1.0
 
 # required to get xarray to not complain
 dask

tests/test_components.py

@@ -819,3 +819,57 @@ def surface_monitor_helper(center, size, monitor_test):
     # z+ surface
     assert monitor_surfaces[5].center == (center[0], center[1], center[2] + size[2] / 2.0)
     assert monitor_surfaces[5].size == (size[0], size[1], 0.0)
+
+
+""" modes """
+
+
+def test_mode_object_syms():
+    """Test that errors are raised if a mode object is not placed right in the presence of syms."""
+    g = GaussianPulse(freq0=1, fwidth=0.1)
+
+    # wrong mode source
+    with pytest.raises(SetupError) as e_info:
+        sim = Simulation(
+            center=(1.0, -1.0, 0.5),
+            size=(2.0, 2.0, 2.0),
+            grid_size=(0.01, 0.01, 0.01),
+            run_time=1e-12,
+            symmetry=(1, -1, 0),
+            sources=[ModeSource(size=(2, 2, 0), direction="+", source_time=g)],
+        )
+
+    # wrong mode monitor
+    with pytest.raises(SetupError) as e_info:
+        sim = Simulation(
+            center=(1.0, -1.0, 0.5),
+            size=(2.0, 2.0, 2.0),
+            grid_size=(0.01, 0.01, 0.01),
+            run_time=1e-12,
+            symmetry=(1, -1, 0),
+            monitors=[ModeMonitor(size=(2, 2, 0), name="mnt", freqs=[2], mode_spec=ModeSpec())],
+        )
+
+    # right mode source (centered on the symmetry)
+    sim = Simulation(
+        center=(1.0, -1.0, 0.5),
+        size=(2.0, 2.0, 2.0),
+        grid_size=(0.01, 0.01, 0.01),
+        run_time=1e-12,
+        symmetry=(1, -1, 0),
+        sources=[ModeSource(center=(1, -1, 1), size=(2, 2, 0), direction="+", source_time=g)],
+    )
+
+    # right mode monitor (entirely in the main quadrant)
+    sim = Simulation(
+        center=(1.0, -1.0, 0.5),
+        size=(2.0, 2.0, 2.0),
+        grid_size=(0.01, 0.01, 0.01),
+        run_time=1e-12,
+        symmetry=(1, -1, 0),
+        monitors=[
+            ModeMonitor(
+                center=(2, 0, 1), size=(2, 2, 0), name="mnt", freqs=[2], mode_spec=ModeSpec()
+            )
+        ],
+    )

tests/test_grid.py

@@ -102,8 +102,8 @@ def test_sim_nonuniform_large():
     bound_coords = sim.grid.boundaries.x
     dls = np.diff(bound_coords)
 
-    dl_min = grid_size_x[0]
-    dl_max = grid_size_x[-1]
+    dl_min = dls[0]
+    dl_max = dls[-1]
 
     # checks the bounds were adjusted correctly
     # (smaller than sim size as is, but larger than sim size with one dl added on each edge)
@@ -115,20 +115,13 @@ def test_sim_nonuniform_large():
 
     # tests that PMLs were added correctly
     for i in range(num_layers_pml_x):
-        assert np.diff(bound_coords[i : i + 2]) == grid_size_x[0]
-        assert np.diff(bound_coords[-2 - i : len(bound_coords) - i]) == grid_size_x[-1]
-
-    # tests that all the grid sizes are in there
-    for size in grid_size_x:
-        assert size in dls
+        assert np.diff(bound_coords[i : i + 2]) == dls[0]
+        assert np.diff(bound_coords[-2 - i : len(bound_coords) - i]) == dls[-1]
 
     # tests that nothing but the grid sizes are in there
     for dl in dls:
         assert dl in grid_size_x
 
-    # tests that it gives exactly what we expect
-    # assert np.all(bound_coords == np.array([-12,-10,-8,-6,-4,-2,0,1,4,7,10,13,16]))
-
 
 def test_sim_grid():
 
@@ -143,6 +136,36 @@ def test_sim_grid():
         assert np.all(b == np.array([-2, -1, 0, 1, 2]))
 
 
+def test_sim_symmetry_grid():
+    """tests that a grid symmetric w.r.t. the simulation center is created in presence of
+    symmetries."""
+
+    grid_size = [2, 1, 3, 2]
+    sim = td.Simulation(
+        center=(1, 1, 1),
+        size=(10, 10, 10),
+        grid_size=(grid_size, grid_size, grid_size),
+        pml_layers=[
+            td.PML(num_layers=2),
+        ]
+        * 3,
+        symmetry=(0, 1, -1),
+    )
+
+    coords_x, coords_y, coords_z = sim.grid.boundaries.to_list
+
+    # Assert coords size is even for the non-symmetry axis but odd otherwise
+    assert coords_x.size % 2 == 0
+    assert coords_y.size % 2 != 0
+    assert coords_z.size % 2 != 0
+
+    # Assert the dls along the symmetric axes are symmetric
+    dls_y = np.diff(coords_y)
+    dls_z = np.diff(coords_z)
+    assert np.all(dls_y[dls_y.size // 2 - 1 :: -1] == dls_y[dls_y.size // 2 :])
+    assert np.all(dls_z[dls_z.size // 2 - 1 :: -1] == dls_z[dls_z.size // 2 :])
+
+
 def test_sim_pml_grid():
 
     sim = td.Simulation(

tidy3d/components/data.py

@@ -9,9 +9,10 @@
 import numpy as np
 import h5py
 
-from .types import Numpy, Direction, Array, numpy_encoding, Literal, Ax
+from .types import Numpy, Direction, Array, numpy_encoding, Literal, Ax, Coordinate, Symmetry
 from .base import Tidy3dBaseModel
 from .simulation import Simulation
+from .grid import YeeGrid
 from .mode import ModeSpec
 from .viz import add_ax_if_none, equal_aspect
 from ..log import log, DataError
@@ -397,6 +398,53 @@ def colocate(self, x, y, z) -> xr.Dataset:
         # import pdb; pdb.set_trace()
         return xr.Dataset(centered_data_dict)
 
+    # pylint:disable=too-many-locals
+    def apply_syms(self, new_grid: YeeGrid, sym_center: Coordinate, symmetry: Symmetry):
+        """Create a new AbstractFieldData subclass by interpolating on the supplied ``new_grid``,
+        using symmetries as defined by ``sym_center`` and ``symmetry``."""
+
+        new_data_dict = {}
+        yee_grid_dict = new_grid.yee.grid_dict
+        # Defines how field components are affected by a positive symmetry along each of the axes
+        component_sym_dict = {
+            "Ex": [-1, 1, 1],
+            "Ey": [1, -1, 1],
+            "Ez": [1, 1, -1],
+            "Hx": [1, -1, -1],
+            "Hy": [-1, 1, -1],
+            "Hz": [-1, -1, 1],
+        }
+
+        for field, scalar_data in self.data_dict.items():
+            new_data = scalar_data.data
+
+            # Get new grid locations
+            yee_coords = yee_grid_dict[field].to_list
+
+            # Apply symmetries
+            zipped = zip("xyz", yee_coords, sym_center, symmetry)
+            for dim, (dim_name, coords, center, sym) in enumerate(zipped):
+                # There shouldn't be anything to do if there's no symmetry on this axis
+                if sym == 0:
+                    continue
+
+                # Get indexes of coords that lie on the left of the symmetry center
+                flip_inds = np.where(coords < center)[0]
+
+                # Get the symmetric coordinates on the right
+                coords[flip_inds] = 2 * center - coords[flip_inds]
+
+                # Interpolate. There generally shouldn't be values out of bounds except potentially
+                # when handling modes, in which case we set such values to zero.
+                new_data = new_data.interp({dim_name: coords}, kwargs={"fill_value": 0.0})
+
+                # Apply the correct +/-1 for the field component
+                new_data[{dim_name: flip_inds}] *= sym * component_sym_dict[field][dim]
+
+            new_data_dict[field] = type(scalar_data)(values=new_data.values, **new_data.coords)
+
+        return type(self)(data_dict=new_data_dict)
+
 
 class FreqData(MonitorData, ABC):
     """Stores frequency-domain data using an ``f`` dimension for frequency in Hz."""

tidy3d/components/grid.py

@@ -98,6 +98,19 @@ class YeeGrid(Tidy3dBaseModel):
         description="Coordinates of the locations of all three components of the magnetic field.",
     )
 
+    @property
+    def grid_dict(self):
+        """The Yee grid coordinates associated to various field components as a dictionary."""
+        yee_grid_dict = {
+            "Ex": self.E.x,
+            "Ey": self.E.y,
+            "Ez": self.E.z,
+            "Hx": self.H.x,
+            "Hy": self.H.y,
+            "Hz": self.H.z,
+        }
+        return yee_grid_dict
+
 
 class Grid(Tidy3dBaseModel):
     """Contains all information about the spatial positions of the FDTD grid.
@@ -227,7 +240,7 @@ def _dual_steps(self) -> Coords:
         primal_steps = self._primal_steps.dict(exclude={TYPE_TAG_STR})
         dsteps = {}
         for (key, psteps) in primal_steps.items():
-            dsteps[key] = (psteps + np.roll(psteps, -1)) / 2
+            dsteps[key] = (psteps + np.roll(psteps, 1)) / 2
 
         return Coords(**dsteps)
 
@@ -292,14 +305,14 @@ def _yee_e(self, axis: Axis):
         return Coords(**yee_coords)
 
     def _yee_h(self, axis: Axis):
-        """E field yee lattice sites for axis."""
+        """H field yee lattice sites for axis."""
 
         boundary_coords = self.boundaries.dict(exclude={TYPE_TAG_STR})
 
-        # initially set all to the minus bounds
+        # initially set all to centers
         yee_coords = {key: self._avg(val) for key, val in boundary_coords.items()}
 
-        # average the axis index between the cell boundaries
+        # set the axis index to the minus bounds
         key = "xyz"[axis]
         yee_coords[key] = self._min(boundary_coords[key])
 

tidy3d/components/mode.py

@@ -6,7 +6,7 @@
 
 from ..constants import MICROMETER
 from .base import Tidy3dBaseModel
-from .types import Symmetry, Axis2D
+from .types import Axis2D
 from ..log import SetupError
 
 
@@ -17,7 +17,7 @@ class ModeSpec(Tidy3dBaseModel):
 
     Example
     -------
-    >>> mode_spec = ModeSpec(num_modes=3, target_neff=1.5, symmetries=(1, -1))
+    >>> mode_spec = ModeSpec(num_modes=3, target_neff=1.5)
     """
 
     num_modes: pd.PositiveInt = pd.Field(
@@ -28,13 +28,6 @@ class ModeSpec(Tidy3dBaseModel):
         None, title="Target effective index", description="Guess for effective index of the mode."
     )
 
-    symmetries: Tuple[Symmetry, Symmetry] = pd.Field(
-        (0, 0),
-        title="Tangential symmetries",
-        description="Symmetries to apply to the modes in the two tangential axes. "
-        "Values of (0, 1,-1) correspond to (none, even, odd) symmetries, respectvely.",
-    )
-
     num_pml: Tuple[pd.NonNegativeInt, pd.NonNegativeInt] = pd.Field(
         (0, 0),
         title="Number of PML layers",