simulation.subsection to include PML

CustomGridBoundaries to precisely replicate grids from a simulation of identical size

CHANGELOG.md

@@ -12,6 +12,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 - Error if field projection monitors found in 2D simulations, except `FieldProjectionAngleMonitor` with `far_field_approx = True`. Support for other monitors and for exact field projection will be coming in a subsequent Tidy3D version.
+- Support for grid specifications from grid cell boundary coordinates via `CustomGridBoundaries` that subclasses from `GridSpec1d`.
+- More convenient mesh importing from another simulation through `grid_spec = GridSpec.from_grid(sim.grid)`.
 
 ### Fixed
 - Error when loading a previously run `Batch` or `ComponentModeler` containing custom data.

docs/api/discretization.rst

@@ -11,6 +11,7 @@ Discretization
    tidy3d.AutoGrid
    tidy3d.UniformGrid
    tidy3d.CustomGrid
+   tidy3d.CustomGridBoundaries
    tidy3d.Coords
    tidy3d.FieldGrid
    tidy3d.YeeGrid

tests/test_components/test_grid_spec.py

@@ -305,3 +305,40 @@ def test_zerosize_dimensions():
             ),
             run_time=1e-12,
         )
+
+
+def test_custom_grid_boundaries():
+    custom = td.CustomGridBoundaries(coords=np.linspace(-1, 1, 11))
+    grid_spec = td.GridSpec(grid_x=custom, grid_y=custom, grid_z=custom)
+    source = td.PointDipole(
+        source_time=td.GaussianPulse(freq0=3e14, fwidth=1e14), polarization="Ex"
+    )
+
+    # matches exactly
+    sim = td.Simulation(
+        size=(2, 2, 2),
+        sources=[source],
+        grid_spec=grid_spec,
+        run_time=1e-12,
+        medium=td.Medium(permittivity=4),
+        boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()),
+    )
+    assert np.allclose(sim.grid.boundaries.x, custom.coords)
+
+    # chop off
+    sim_chop = sim.updated_copy(size=(1, 1, 1))
+    assert np.allclose(sim_chop.grid.boundaries.x, np.linspace(-0.4, 0.4, 5))
+
+    sim_chop = sim.updated_copy(size=(1.2, 1, 1))
+    assert np.allclose(sim_chop.grid.boundaries.x, np.linspace(-0.6, 0.6, 7))
+
+    # expand
+    sim_expand = sim.updated_copy(size=(4, 4, 4))
+    assert np.allclose(sim_expand.grid.boundaries.x, np.linspace(-2, 2, 21))
+
+    # pml
+    num_layers = 10
+    sim_pml = sim.updated_copy(
+        boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML(num_layers=num_layers))
+    )
+    assert np.allclose(sim_pml.grid.boundaries.x, np.linspace(-3, 3, 31))

tidy3d/__init__.py

@@ -119,7 +119,13 @@
 from .components.geometry.polyslab import PolySlab
 from .components.geometry.primitives import Cylinder, Sphere
 from .components.grid.grid import Coords, Coords1D, FieldGrid, Grid, YeeGrid
-from .components.grid.grid_spec import AutoGrid, CustomGrid, GridSpec, UniformGrid
+from .components.grid.grid_spec import (
+    AutoGrid,
+    CustomGrid,
+    CustomGridBoundaries,
+    GridSpec,
+    UniformGrid,
+)
 from .components.heat.boundary import ConvectionBC, HeatBoundarySpec, HeatFluxBC, TemperatureBC
 from .components.heat.data.monitor_data import TemperatureData
 from .components.heat.data.sim_data import HeatSimulationData
@@ -287,6 +293,7 @@ def set_logging_level(level: str) -> None:
     "UniformGrid",
     "CustomGrid",
     "AutoGrid",
+    "CustomGridBoundaries",
     "Box",
     "Sphere",
     "Cylinder",

tidy3d/components/grid/grid_spec.py

@@ -8,7 +8,7 @@
 import numpy as np
 import pydantic.v1 as pd
 
-from ...constants import C_0, MICROMETER, fp_eps
+from ...constants import C_0, MICROMETER, inf
 from ...exceptions import SetupError
 from ...log import log
 from ..base import Tidy3dBaseModel
@@ -136,6 +136,78 @@ def _add_pml_to_bounds(num_layers: Tuple[int, int], bounds: Coords1D) -> Coords1
         add_right = bounds[-1] + last_step * np.arange(1, num_layers[1] + 1)
         return np.concatenate((add_left, bounds, add_right))
 
+    @staticmethod
+    def _postprocess_unaligned_grid(
+        axis: Axis,
+        simulation_box: Box,
+        machine_error_relaxation: bool,
+        bound_coords: Coords1D,
+    ) -> Coords1D:
+        """Postprocess grids whose two ends  might be aligned with simulation boundaries.
+        This is to be used in `_make_coords_initial`.
+
+        Parameters
+        ----------
+        axis : Axis
+            Axis of this direction.
+        structures : List[StructureType]
+            List of structures present in simulation, the first one being the simulation domain.
+        machine_error_relaxation : bool
+            When operations such as translation are applied to the 1d grids, fix the bounds
+            were numerically within the simulation bounds but were still chopped off.
+        bound_coords : Coord1D
+            1D grids potentially unaligned with the simulation boundary
+
+        Returns
+        -------
+        :class:`.Coords1D`:
+            1D coords to be used as grid boundaries.
+
+        """
+        center, size = simulation_box.center[axis], simulation_box.size[axis]
+        # chop off any coords outside of simulation bounds, beyond some buffer region
+        # to take numerical effects into account
+        bound_min = np.nextafter(center - size / 2, -inf, dtype=np.float32)
+        bound_max = np.nextafter(center + size / 2, inf, dtype=np.float32)
+
+        if bound_max < bound_coords[0] or bound_min > bound_coords[-1]:
+            axis_name = "xyz"[axis]
+            raise SetupError(
+                f"Simulation domain does not overlap with the provided grid in '{axis_name}' direction."
+            )
+
+        if size == 0:
+            # in case of zero-size dimension return the boundaries between which simulation falls
+            ind = np.searchsorted(bound_coords, center, side="right")
+
+            # in case when the center coincides with the right most boundary
+            if ind >= len(bound_coords):
+                ind = len(bound_coords) - 1
+
+            return bound_coords[ind - 1 : ind + 1]
+
+        else:
+            bound_coords = bound_coords[bound_coords <= bound_max]
+            bound_coords = bound_coords[bound_coords >= bound_min]
+
+            # if not extending to simulation bounds, repeat beginning and end
+            dl_min = bound_coords[1] - bound_coords[0]
+            dl_max = bound_coords[-1] - bound_coords[-2]
+            while bound_coords[0] - dl_min >= bound_min:
+                bound_coords = np.insert(bound_coords, 0, bound_coords[0] - dl_min)
+            while bound_coords[-1] + dl_max <= bound_max:
+                bound_coords = np.append(bound_coords, bound_coords[-1] + dl_max)
+
+            # in case operations are applied to coords, it's possible the bounds were numerically within
+            # the simulation bounds but were still chopped off, which is fixed here
+            if machine_error_relaxation:
+                if np.isclose(bound_coords[0] - dl_min, bound_min):
+                    bound_coords = np.insert(bound_coords, 0, bound_coords[0] - dl_min)
+                if np.isclose(bound_coords[-1] + dl_max, bound_max):
+                    bound_coords = np.append(bound_coords, bound_coords[-1] + dl_max)
+
+            return bound_coords
+
 
 class UniformGrid(GridSpec1d):
     """Uniform 1D grid. The most standard way to define a simulation is to use a constant grid size in each of the three directions.
@@ -176,8 +248,6 @@ def _make_coords_initial(
             Axis of this direction.
         structures : List[StructureType]
             List of structures present in simulation, the first one being the simulation domain.
-        **kwargs:
-            Other arguments all go here.
 
         Returns
         -------
@@ -199,6 +269,50 @@ def _make_coords_initial(
         return center - size / 2 + np.arange(num_cells + 1) * dl_snapped
 
 
+class CustomGridBoundaries(GridSpec1d):
+    """Custom 1D grid supplied as a list of grid cell boundary coordinates.
+
+    Example
+    -------
+    >>> grid_1d = CustomGridCoords(boundaries=[-0.2, 0.0, 0.2, 0.4, 0.5, 0.6, 0.7])
+    """
+
+    coords: Coords1D = pd.Field(
+        ...,
+        title="Grid Boundary Coordinates",
+        description="An array of grid boundary coordinates.",
+        units=MICROMETER,
+    )
+
+    def _make_coords_initial(
+        self,
+        axis: Axis,
+        structures: List[StructureType],
+        **kwargs,
+    ) -> Coords1D:
+        """Customized 1D coords to be used as grid boundaries.
+
+        Parameters
+        ----------
+        axis : Axis
+            Axis of this direction.
+        structures : List[StructureType]
+            List of structures present in simulation, the first one being the simulation domain.
+
+        Returns
+        -------
+        :class:`.Coords1D`:
+            1D coords to be used as grid boundaries.
+        """
+
+        return self._postprocess_unaligned_grid(
+            axis=axis,
+            simulation_box=structures[0].geometry,
+            machine_error_relaxation=False,
+            bound_coords=self.coords,
+        )
+
+
 class CustomGrid(GridSpec1d):
     """Custom 1D grid supplied as a list of grid cell sizes centered on the simulation center.
 
@@ -241,16 +355,14 @@ def _make_coords_initial(
             Axis of this direction.
         structures : List[StructureType]
             List of structures present in simulation, the first one being the simulation domain.
-        *kwargs
-            Other arguments all go here.
 
         Returns
         -------
         :class:`.Coords1D`:
             1D coords to be used as grid boundaries.
         """
 
-        center, size = structures[0].geometry.center[axis], structures[0].geometry.size[axis]
+        center = structures[0].geometry.center[axis]
 
         # get bounding coordinates
         dl = np.array(self.dl)
@@ -263,49 +375,12 @@ def _make_coords_initial(
         else:
             bound_coords += self.custom_offset
 
-        # chop off any coords outside of simulation bounds, beyond some buffer region
-        # to take numerical effects into account
-        buffer = fp_eps * size
-        bound_min = center - size / 2 - buffer
-        bound_max = center + size / 2 + buffer
-
-        if bound_max < bound_coords[0] or bound_min > bound_coords[-1]:
-            axis_name = "xyz"[axis]
-            raise SetupError(
-                f"Simulation domain does not overlap with the provided custom grid in '{axis_name}' direction."
-            )
-
-        if size == 0:
-            # in case of zero-size dimension return the boundaries between which simulation falls
-            ind = np.searchsorted(bound_coords, center, side="right")
-
-            # in case when the center coincides with the right most boundary
-            if ind >= len(bound_coords):
-                ind = len(bound_coords) - 1
-
-            return bound_coords[ind - 1 : ind + 1]
-
-        else:
-            bound_coords = bound_coords[bound_coords <= bound_max]
-            bound_coords = bound_coords[bound_coords >= bound_min]
-
-            # if not extending to simulation bounds, repeat beginning and end
-            dl_min = dl[0]
-            dl_max = dl[-1]
-            while bound_coords[0] - dl_min >= bound_min:
-                bound_coords = np.insert(bound_coords, 0, bound_coords[0] - dl_min)
-            while bound_coords[-1] + dl_max <= bound_max:
-                bound_coords = np.append(bound_coords, bound_coords[-1] + dl_max)
-
-            # in case a `custom_offset` is provided, it's possible the bounds were numerically within
-            # the simulation bounds but were still chopped off, which is fixed here
-            if self.custom_offset is not None:
-                if np.isclose(bound_coords[0] - dl_min, bound_min):
-                    bound_coords = np.insert(bound_coords, 0, bound_coords[0] - dl_min)
-                if np.isclose(bound_coords[-1] + dl_max, bound_max):
-                    bound_coords = np.append(bound_coords, bound_coords[-1] + dl_max)
-
-            return bound_coords
+        return self._postprocess_unaligned_grid(
+            axis=axis,
+            simulation_box=structures[0].geometry,
+            machine_error_relaxation=self.custom_offset is not None,
+            bound_coords=bound_coords,
+        )
 
 
 class AutoGrid(GridSpec1d):
@@ -454,7 +529,7 @@ def _make_coords_initial(
         return np.array(bound_coords)
 
 
-GridType = Union[UniformGrid, CustomGrid, AutoGrid]
+GridType = Union[UniformGrid, CustomGrid, AutoGrid, CustomGridBoundaries]
 
 
 class GridSpec(Tidy3dBaseModel):
@@ -657,6 +732,14 @@ def make_grid(
         coords = Coords(**coords_dict)
         return Grid(boundaries=coords)
 
+    @classmethod
+    def from_grid(cls, grid: Grid) -> GridSpec:
+        """Import grid directly from another simulation, e.g. ``grid_spec = GridSpec.from_grid(sim.grid)``."""
+        grid_dict = {}
+        for dim in "xyz":
+            grid_dict["grid_" + dim] = CustomGridBoundaries(coords=grid.boundaries.to_dict[dim])
+        return cls(**grid_dict)
+
     @classmethod
     def auto(
         cls,

tidy3d/components/simulation.py

@@ -37,7 +37,7 @@
 from .geometry.utils import flatten_groups, traverse_geometries
 from .geometry.utils_2d import get_bounds, get_thickened_geom, snap_coordinate_to_grid, subdivide
 from .grid.grid import Coords, Coords1D, Grid
-from .grid.grid_spec import AutoGrid, CustomGrid, GridSpec, UniformGrid
+from .grid.grid_spec import AutoGrid, GridSpec, UniformGrid
 from .lumped_element import LumpedElementType
 from .medium import (
     AbstractCustomMedium,
@@ -4297,6 +4297,7 @@ def subsection(
         monitors: Tuple[MonitorType, ...] = None,
         remove_outside_structures: bool = True,
         remove_outside_custom_mediums: bool = False,
+        include_pml_cells: bool = False,
         **kwargs,
     ) -> Simulation:
         """Generate a simulation instance containing only the ``region``.
@@ -4327,6 +4328,9 @@ def subsection(
             Remove structures outside of the new simulation domain.
         remove_outside_custom_mediums : bool = True
             Remove custom medium data outside of the new simulation domain.
+        include_pml_cells : bool = False
+            Keep PML cells in simulation boundaries. Note that retained PML cells will be converted
+            to regular cells, and the simulation domain boundary will be moved accordingly.
         **kwargs
             Other arguments passed to new simulation instance.
         """
@@ -4336,7 +4340,10 @@ def subsection(
             raise SetupError("Requested region does not intersect simulation domain")
 
         # restrict to the original simulation domain
-        new_bounds = Box.bounds_intersection(self.bounds, region.bounds)
+        if include_pml_cells:
+            new_bounds = Box.bounds_intersection(self.simulation_bounds, region.bounds)
+        else:
+            new_bounds = Box.bounds_intersection(self.bounds, region.bounds)
         new_bounds = [list(new_bounds[0]), list(new_bounds[1])]
 
         # grid spec inheritace
@@ -4345,11 +4352,7 @@ def subsection(
         elif isinstance(grid_spec, str) and grid_spec == "identical":
             # create a custom grid from existing one
             grids_1d = self.grid.boundaries.to_list
-            new_grids = [
-                CustomGrid(dl=tuple(np.diff(grids_1d[dim])), custom_offset=grids_1d[dim][0])
-                for dim in range(3)
-            ]
-            grid_spec = GridSpec(grid_x=new_grids[0], grid_y=new_grids[1], grid_z=new_grids[2])
+            grid_spec = GridSpec.from_grid(self.grid)
 
             # adjust region bounds to perfectly coincide with the grid
             # note, sometimes (when a box already seems to perfrecty align with the grid)

tidy3d/plugins/mode/mode_solver.py

@@ -1536,6 +1536,7 @@ def reduced_simulation_copy(self):
             boundary_spec=new_bspec,
             remove_outside_custom_mediums=True,
             remove_outside_structures=True,
+            include_pml_cells=True,
         )
 
         return self.updated_copy(simulation=new_sim)