diff --git a/CHANGELOG.md b/CHANGELOG.md index 9d6a823cfc..b9e6fc9425 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/requirements.txt b/requirements.txt index e1a8a62620..3acf27edc8 100644 --- a/requirements.txt +++ b/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 diff --git a/tests/test_components.py b/tests/test_components.py index e0487478a2..e1513a3b8e 100644 --- a/tests/test_components.py +++ b/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() + ) + ], + ) diff --git a/tests/test_grid.py b/tests/test_grid.py index 2432a563ea..49270fe99f 100644 --- a/tests/test_grid.py +++ b/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( diff --git a/tidy3d/components/data.py b/tidy3d/components/data.py index 960d0d7435..252e1234f5 100644 --- a/tidy3d/components/data.py +++ b/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.""" diff --git a/tidy3d/components/grid.py b/tidy3d/components/grid.py index 85b6b77b5f..4817fd63b5 100644 --- a/tidy3d/components/grid.py +++ b/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]) diff --git a/tidy3d/components/mode.py b/tidy3d/components/mode.py index 8529cfa0ef..7fe9068569 100644 --- a/tidy3d/components/mode.py +++ b/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", diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py index b667377ecf..7494ebe401 100644 --- a/tidy3d/components/simulation.py +++ b/tidy3d/components/simulation.py @@ -11,6 +11,7 @@ from descartes import PolygonPatch from .validators import assert_unique_names, assert_objects_in_sim_bounds +from .validators import validate_mode_objects_symmetry from .geometry import Box from .types import Symmetry, Ax, Shapely, FreqBound, GridSize from .grid import Coords1D, Grid, Coords @@ -114,6 +115,17 @@ class Simulation(Box): # pylint:disable=too-many-public-methods units=SECOND, ) + symmetry: Tuple[Symmetry, Symmetry, Symmetry] = pydantic.Field( + (0, 0, 0), + title="Symmetries", + description="Tuple of integers defining reflection symmetry across a plane " + "bisecting the simulation domain normal to the x-, y-, and z-axis, respectvely. " + "Each element can be ``0`` (no symmetry), ``1`` (even, i.e. 'PMC' symmetry) or " + "``-1`` (odd, i.e. 'PEC' symmetry). " + "Note that the vectorial nature of the fields must be taken into account to correctly " + "determine the symmetry value.", + ) + structures: List[Structure] = pydantic.Field( [], title="Structures", @@ -143,17 +155,6 @@ class Simulation(Box): # pylint:disable=too-many-public-methods "and periodic boundary conditions will be used.", ) - symmetry: Tuple[Symmetry, Symmetry, Symmetry] = pydantic.Field( - (0, 0, 0), - title="Symmetries", - description="Tuple of integers defining reflection symmetry across a plane " - "bisecting the simulation domain normal to the x-, y-, and z-axis, respectvely. " - "Each element can be ``0`` (no symmetry), ``1`` (even, i.e. 'PMC' symmetry) or " - "``-1`` (odd, i.e. 'PEC' symmetry). " - "Note that the vectorial nature of the fields must be taken into account to correctly " - "determine the symmetry value.", - ) - shutoff: pydantic.NonNegativeFloat = pydantic.Field( 1e-5, title="Shutoff Condition", @@ -188,17 +189,6 @@ class Simulation(Box): # pylint:disable=too-many-public-methods """ Validating setup """ - @pydantic.validator("symmetry", always=True) - def not_supported_yet(cls, val): - """Log an error if non 0 value supplied.""" - if any(sym_val != 0 for sym_val in val): - raise SetupError( - "Symmetry not supported in this version of tidy3d, " - "but will be included in coming release." - "For now, if symmetry is required, use the originl tidy3d." - ) - return val - @pydantic.validator("pml_layers", always=True) def set_none_to_zero_layers(cls, val): """if any PML layer is None, set it to an empty :class:`PML`.""" @@ -207,6 +197,8 @@ def set_none_to_zero_layers(cls, val): _structures_in_bounds = assert_objects_in_sim_bounds("structures") _sources_in_bounds = assert_objects_in_sim_bounds("sources") _monitors_in_bounds = assert_objects_in_sim_bounds("monitors") + _mode_sources_symmetries = validate_mode_objects_symmetry("sources") + _mode_monitors_symmetries = validate_mode_objects_symmetry("monitors") # assign names to unnamed structures, sources, and mediums # _structure_names = set_names("structures") @@ -849,13 +841,8 @@ def num_pml_layers(self) -> List[Tuple[float, float]]: List[Tuple[float, float]] List containing the number of absorber layers in - and + boundaries. """ - num_layers = [] - for pml_axis, pml_layer in enumerate(self.pml_layers): - if self.symmetry[pml_axis] != 0: - num_layers.append((0, pml_layer.num_layers)) - else: - num_layers.append((pml_layer.num_layers, pml_layer.num_layers)) - return num_layers + + return [(pml.num_layers, pml.num_layers) for pml in self.pml_layers] @property def pml_thicknesses(self) -> List[Tuple[float, float]]: @@ -1112,7 +1099,8 @@ def num_time_steps(self) -> int: return len(self.tmesh) - def _make_bound_coords_uniform(self, dl, center, size, num_layers): + @staticmethod + def _make_bound_coords_uniform(dl, center, size): """creates coordinate boundaries with uniform mesh (dl is float)""" num_cells = int(np.floor(size / dl)) @@ -1120,22 +1108,21 @@ def _make_bound_coords_uniform(self, dl, center, size, num_layers): # Make sure there's at least one cell num_cells = max(num_cells, 1) - # snap to grid, recenter, and add PML + # snap to grid, recenter size_snapped = dl * num_cells bound_coords = center + np.linspace(-size_snapped / 2, size_snapped / 2, num_cells + 1) - bound_coords = self._add_pml_to_bounds(num_layers, bound_coords) return bound_coords @staticmethod - def _make_bound_coords_nonuniform(dl, center, size, num_layers): + def _make_bound_coords_nonuniform(dl, center, size): """creates coordinate boundaries with non-uniform mesh (dl is arraylike)""" # get bounding coordinates dl = np.array(dl) bound_coords = np.array([np.sum(dl[:i]) for i in range(len(dl) + 1)]) - # shift coords to center at center of simulation along dimension - bound_coords = bound_coords - np.sum(dl) / 2 + center + # place the middle boundary at the center of the simulation along dimension + bound_coords += center - bound_coords[bound_coords.size // 2] # chop off any coords outside of simulation bounds bound_min = center - size / 2 @@ -1151,11 +1138,29 @@ def _make_bound_coords_nonuniform(dl, center, size, num_layers): while bound_coords[-1] + dl_max <= bound_max: bound_coords = np.append(bound_coords, bound_coords[-1] + dl_max) - # add PML layers in using dl on edges - for _ in range(num_layers[0]): - bound_coords = np.insert(bound_coords, 0, bound_coords[0] - dl_min) - for _ in range(num_layers[1]): - bound_coords = np.append(bound_coords, bound_coords[-1] + dl_max) + return bound_coords + + def _make_bound_coords(self, dim): + """Creates coordinate boundaries along dimension ``dim`` and handle PML and symmetries""" + + dl = self.grid_size[dim] + center = self.center[dim] + + # Make uniform or nonuniform boundaries depending on dl input + if isinstance(dl, float): + bound_coords = self._make_bound_coords_uniform(dl, center, self.size[dim]) + else: + bound_coords = self._make_bound_coords_nonuniform(dl, center, self.size[dim]) + + # Add PML layers in using dl on edges + bound_coords = self._add_pml_to_bounds(self.num_pml_layers[dim], bound_coords) + + # Enforce a symmetric grid by taking only the coordinates to the right of the center. + # Note that there's always a boundary at the simulation center. + if self.symmetry[dim] != 0: + bound_coords = bound_coords[bound_coords >= center] + bound_coords = np.append(2 * center - bound_coords[:0:-1], bound_coords) + return bound_coords @property @@ -1168,13 +1173,8 @@ def grid(self) -> Grid: :class:`Grid` storing the spatial locations relevant to the simulation. """ cell_boundary_dict = {} - zipped_vals = zip("xyz", self.grid_size, self.center, self.size, self.num_pml_layers) - for key, dl, center, size, num_layers in zipped_vals: - if isinstance(dl, float): - bound_coords = self._make_bound_coords_uniform(dl, center, size, num_layers) - else: - bound_coords = self._make_bound_coords_nonuniform(dl, center, size, num_layers) - cell_boundary_dict[key] = bound_coords + for dim, key in enumerate("xyz"): + cell_boundary_dict[key] = self._make_bound_coords(dim) boundaries = Coords(**cell_boundary_dict) return Grid(boundaries=boundaries) @@ -1233,6 +1233,40 @@ def wvl_mat_min(self) -> float: n_max, _ = AbstractMedium.eps_complex_to_nk(eps_max) return wvl_min / n_max + def min_sym_box(self, box: Box) -> Box: + """Compute the smallest Box restricted to the first quadrant in the presence of symmetries + that fully covers the original Box when symmetries are applied. + + Parameters + ---------- + box : :class:`Box` + Rectangular geometry. + + Returns + ------- + new_box : :class:`Box` + The smallest Box such that any point in ``box`` is either in ``new_box`` or can be + mapped from ``new_box`` using the simulation symmetries. + """ + + bounds_min, bounds_max = box.bounds + bmin_new, bmax_new = [], [] + + for (center, sym, bmin, bmax) in zip(self.center, self.symmetry, bounds_min, bounds_max): + if sym == 0 or center < bmin: + bmin_new.append(bmin) + bmax_new.append(bmax) + else: + if bmax < center: + bmin_new.append(2 * center - bmax) + bmax_new.append(2 * center - bmin) + else: + # bmin <= center <= bmax + bmin_new.append(center) + bmax_new.append(max(bmax, 2 * center - bmin)) + + return Box.from_bounds(bmin_new, bmax_new) + def discretize_inds(self, box: Box) -> List[Tuple[int, int]]: """Start and stopping indexes for the cells that intersect with a :class:`Box`. diff --git a/tidy3d/components/validators.py b/tidy3d/components/validators.py index 8657dc1d42..e926c2c8b2 100644 --- a/tidy3d/components/validators.py +++ b/tidy3d/components/validators.py @@ -64,6 +64,37 @@ def field_has_unique_names(cls, val): return field_has_unique_names +def validate_mode_objects_symmetry(field_name: str): + """If a Mode object, this checks that the object is fully in the main quadrant in the presence + of symmetry along a given axis, or else centered on the symmetry center.""" + + obj_type = "ModeSource" if field_name == "sources" else "ModeMonitor" + + @pydantic.validator(field_name, allow_reuse=True, always=True) + def check_symmetry(cls, val, values): + """check for intersection of each structure with simulation bounds.""" + sim_center = values.get("center") + for position_index, geometric_object in enumerate(val): + if geometric_object.type == obj_type: + bounds_min, _ = geometric_object.bounds + for dim, sym in enumerate(values.get("symmetry")): + if ( + sym != 0 + and bounds_min[dim] < sim_center[dim] + and geometric_object.center[dim] != sim_center[dim] + ): + raise SetupError( + f"Mode object '{geometric_object}' " + f"(at `simulation.{field_name}[{position_index}]`) " + "in presence of symmetries must be in the main quadrant, " + "or centered on the symmetry axis." + ) + + return val + + return check_symmetry + + def assert_unique_names(field_name: str, check_mediums=False): """makes sure all elements of a field have unique .name values""" @@ -85,19 +116,22 @@ def field_has_unique_names(cls, val, values): def assert_objects_in_sim_bounds(field_name: str): - """makes sure all objects in field are at least partially inside of simulation bounds/""" + """Makes sure all objects in field are at least partially inside of simulation bounds.""" @pydantic.validator(field_name, allow_reuse=True, always=True) def objects_in_sim_bounds(cls, val, values): """check for intersection of each structure with simulation bounds.""" - sim_bounds = Box(size=values.get("size"), center=values.get("center")) + sim_center = values.get("center") + sim_box = Box(size=values.get("size"), center=sim_center) + for position_index, geometric_object in enumerate(val): - if not sim_bounds.intersects(geometric_object.geometry): + if not sim_box.intersects(geometric_object.geometry): raise SetupError( f"'{geometric_object}' " - f"(at `simulation.{field_name}[{position_index}]`)" - "is completely outside of simulation domain" + f"(at `simulation.{field_name}[{position_index}]`) " + "is completely outside of simulation domain." ) + return val return objects_in_sim_bounds diff --git a/tidy3d/components/viz.py b/tidy3d/components/viz.py index f307dd7065..d60e3b70ce 100644 --- a/tidy3d/components/viz.py +++ b/tidy3d/components/viz.py @@ -6,7 +6,6 @@ import matplotlib.pylab as plt from pydantic import BaseModel -import chart_studio.plotly as py import plotly.graph_objects as go import numpy as np @@ -209,6 +208,7 @@ def get_plot_params(self) -> PatchParams: return PatchParams(alpha=0.4, edgecolor="black") +# pylint:disable=too-many-locals, too-many-branches, too-many-statements def plotly_sim(sim, x=None, y=None, z=None): """Make a plotly plot.""" @@ -302,10 +302,12 @@ def plotly_sim(sim, x=None, y=None, z=None): pml_thick = pml.num_layers * dl_edge pml_size = 2 * np.array(sim.size) pml_size[dir_index] = pml_thick - rmin[i] += sign * pml_thick - rmax[i] += sign * pml_thick + if sign == 1: + rmax[dir_index] += pml_thick + else: + rmin[dir_index] -= pml_thick pml_center = np.array(sim.center) + sign * np.array(sim.size) / 2 - pml_center[dir_index] + sign * pml_size / 2 + pml_center[dir_index] += sign * pml_thick / 2 pml_box = sim.geometry pml_box.center = pml_center.tolist() pml_box.size = pml_size.tolist() diff --git a/tidy3d/plugins/mode/mode_solver.py b/tidy3d/plugins/mode/mode_solver.py index f0f03ba86a..098cd6ac4b 100644 --- a/tidy3d/plugins/mode/mode_solver.py +++ b/tidy3d/plugins/mode/mode_solver.py @@ -15,6 +15,7 @@ from ...components import ModeSource, GaussianPulse from ...components.types import Direction from ...components.data import ScalarFieldData, FieldData +from ...log import SetupError from .solver import compute_modes @@ -109,8 +110,23 @@ def solve(self, mode_spec: ModeSpec) -> List[ModeInfo]: normal_axis = self.plane.size.index(0.0) + # get the in-plane grid coordinates on which eps and the mode fields live + plane_grid = self.simulation.discretize(self.plane) + + # restrict to a smaller plane if symmetries present in the simulation + plane_sym = self.simulation.min_sym_box(self.plane) + plane_grid_sym = self.simulation.discretize(plane_sym) + + # Coords and symmetry arguments to the solver (restricted to in-plane) + _, solver_coords = self.plane.pop_axis(plane_grid_sym.boundaries.to_list, axis=normal_axis) + mode_symmetry = list(self.simulation.symmetry) + for dim in range(3): + if self.simulation.center[dim] != self.plane.center[dim]: + mode_symmetry[dim] = 0 + _, solver_symmetry = self.plane.pop_axis(mode_symmetry, axis=normal_axis) + # Get diagonal epsilon components in the plane - (eps_xx, eps_yy, eps_zz) = self.get_epsilon() + (eps_xx, eps_yy, eps_zz) = self.get_epsilon(plane_sym) # get rid of normal axis eps_xx = np.squeeze(eps_xx, axis=normal_axis) @@ -122,35 +138,23 @@ def solve(self, mode_spec: ModeSpec) -> List[ModeInfo]: (eps_xx, eps_yy, eps_zz), axis=normal_axis ) - # get the in-plane grid coordinates on which eps and the mode fields live - plane_grid = self.simulation.discretize(self.plane) - _, plane_coords = self.plane.pop_axis(plane_grid.boundaries.to_list, axis=normal_axis) - - # note: from this point on, in waveguide coordinates (propagating in z) - # construct eps_cross section to feed to mode solver eps_cross = np.stack((eps_wg_xx, eps_wg_yy, eps_wg_zz)) - # Nx, Ny = eps_cross.shape[1:] - # if mode_spec.symmetries[0] != 0: - # eps_cross = np.stack(tuple(e[Nx // 2, :] for e in eps_cross)) - # if mode_spec.symmetries[1] != 0: - # eps_cross = np.stack(tuple(e[:, Ny // 2] for e in eps_cross)) - + # Compute the modes mode_fields, n_eff_complex = compute_modes( eps_cross=eps_cross, - coords=plane_coords, + coords=solver_coords, freq=self.freq, mode_spec=mode_spec, + symmetry=solver_symmetry, ) - def rotate_field_coords(e_field, h_field): + def rotate_field_coords(field): """move the propagation axis=z to the proper order in the array""" - Ex, Ey, Ez = np.moveaxis(e_field, source=3, destination=1 + normal_axis) - e_rot = np.stack(self.simulation.unpop_axis(Ez, (Ex, Ey), axis=normal_axis), axis=0) - Hx, Hy, Hz = np.moveaxis(h_field, source=3, destination=1 + normal_axis) - h_rot = np.stack(self.simulation.unpop_axis(Hz, (Hx, Hy), axis=normal_axis), axis=0) - return (e_rot, h_rot) + f_x, f_y, f_z = np.moveaxis(field, source=3, destination=1 + normal_axis) + f_rot = np.stack(self.plane.unpop_axis(f_z, (f_x, f_y), axis=normal_axis), axis=0) + return f_rot modes = [] for mode_index in range(mode_spec.num_modes): @@ -164,20 +168,9 @@ def rotate_field_coords(e_field, h_field): E *= np.exp(-1j * phi) H *= np.exp(-1j * phi) - # # Handle symmetries - # if mode.symmetries[0] != 0: - # E_half = E[:, 1:, ...] - # H_half = H[:, 1:, ...] - # E = np.concatenate((+E_half[:, ::-1, ...], E_half), axis=1) - # H = np.concatenate((-H_half[:, ::-1, ...], H_half), axis=1) - # if mode.symmetries[1] != 0: - # E_half = E[:, :, 1:, ...] - # H_half = H[:, :, 1:, ...] - # E = np.concatenate((+E_half[:, :, ::-1, ...], E_half), axis=2) - # H = np.concatenate((-H_half[:, :, ::-1, ...], H_half), axis=2) - # Rotate back to original coordinates - (Ex, Ey, Ez), (Hx, Hy, Hz) = rotate_field_coords(E, H) + (Ex, Ey, Ez) = rotate_field_coords(E) + (Hx, Hy, Hz) = rotate_field_coords(H) # apply -1 to H fields if a reflection was involved in the rotation if normal_axis == 1: @@ -191,9 +184,7 @@ def rotate_field_coords(e_field, h_field): # note: re-discretizing, need to make consistent. data_dict = {} for field_name, field in fields.items(): - plane_grid = self.simulation.discretize(self.plane) - plane_coords = plane_grid[field_name] - xyz_coords = [plane_coords.x, plane_coords.y, plane_coords.z] + xyz_coords = plane_grid_sym[field_name].to_list xyz_coords[normal_axis] = [self.plane.center[normal_axis]] data_dict[field_name] = ScalarFieldData( x=xyz_coords[0], @@ -203,8 +194,11 @@ def rotate_field_coords(e_field, h_field): values=field[..., None], ) + field_data = FieldData(data_dict=data_dict).apply_syms( + plane_grid, self.simulation.center, self.simulation.symmetry + ) mode_info = ModeInfo( - field_data=FieldData(data_dict=data_dict), + field_data=field_data, mode_spec=mode_spec, mode_index=mode_index, n_eff=n_eff_complex[mode_index].real, @@ -215,12 +209,12 @@ def rotate_field_coords(e_field, h_field): return modes - def get_epsilon(self): + def get_epsilon(self, plane): """Compute the diagonal components of the epsilon tensor in the plane.""" - eps_xx = self.simulation.epsilon(self.plane, "Ex", self.freq) - eps_yy = self.simulation.epsilon(self.plane, "Ey", self.freq) - eps_zz = self.simulation.epsilon(self.plane, "Ez", self.freq) + eps_xx = self.simulation.epsilon(plane, "Ex", self.freq) + eps_yy = self.simulation.epsilon(plane, "Ey", self.freq) + eps_zz = self.simulation.epsilon(plane, "Ez", self.freq) return np.stack((eps_xx, eps_yy, eps_zz), axis=0) diff --git a/tidy3d/plugins/mode/solver.py b/tidy3d/plugins/mode/solver.py index e4fce7c7c9..ff2becaa0f 100644 --- a/tidy3d/plugins/mode/solver.py +++ b/tidy3d/plugins/mode/solver.py @@ -17,6 +17,7 @@ def compute_modes( coords, freq, mode_spec, + symmetry=(0, 0), angle_theta=0.0, angle_phi=0.0, ) -> Tuple[Numpy, Numpy]: @@ -117,7 +118,7 @@ def compute_modes( those positions, unless a PMC symmetry is specifically requested. The PMC symmetry is imposed by modifying the backward derivative matrices.""" dmin_pmc = [False, False] - if mode_spec.symmetries[0] != 1: + if symmetry[0] != 1: # PEC at the xmin edge eps_tensor[1, 1, :Ny] = pec_val eps_tensor[2, 2, :Ny] = pec_val @@ -126,7 +127,7 @@ def compute_modes( dmin_pmc[0] = True if Ny > 1: - if mode_spec.symmetries[1] != 1: + if symmetry[1] != 1: # PEC at the ymin edge eps_tensor[0, 0, ::Ny] = pec_val eps_tensor[2, 2, ::Ny] = pec_val @@ -144,7 +145,7 @@ def compute_modes( Dmats = D_mats((Nx, Ny), dLf, dLb, dmin_pmc) # PML matrices; do not impose PML on the bottom when symmetry present - dmin_pml = np.array(mode_spec.symmetries) == 0 + dmin_pml = np.array(symmetry) == 0 Smats = S_mats(omega, (Nx, Ny), mode_spec.num_pml, dLf, dLb, dmin_pml) # Add the PML on top of the derivatives; normalize by k0 to match the EM-possible notation diff --git a/tidy3d/web/config.py b/tidy3d/web/config.py index c003a7eecf..21f9f48957 100644 --- a/tidy3d/web/config.py +++ b/tidy3d/web/config.py @@ -3,7 +3,7 @@ from dataclasses import dataclass -SOLVER_VERSION = "release-22.1.3" +SOLVER_VERSION = "release-22.1.4" @dataclass