diff --git a/CHANGELOG.md b/CHANGELOG.md index 5c3bf9e6f1..921f3da623 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/docs/api/discretization.rst b/docs/api/discretization.rst index b9a3be7d96..a15172083d 100644 --- a/docs/api/discretization.rst +++ b/docs/api/discretization.rst @@ -11,6 +11,7 @@ Discretization tidy3d.AutoGrid tidy3d.UniformGrid tidy3d.CustomGrid + tidy3d.CustomGridBoundaries tidy3d.Coords tidy3d.FieldGrid tidy3d.YeeGrid diff --git a/tests/test_components/test_grid_spec.py b/tests/test_components/test_grid_spec.py index 49d36c2387..2b4d8c4718 100644 --- a/tests/test_components/test_grid_spec.py +++ b/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)) diff --git a/tidy3d/__init__.py b/tidy3d/__init__.py index a527535740..85825ab606 100644 --- a/tidy3d/__init__.py +++ b/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", diff --git a/tidy3d/components/grid/grid_spec.py b/tidy3d/components/grid/grid_spec.py index e315d7734d..0fa2dac1f5 100644 --- a/tidy3d/components/grid/grid_spec.py +++ b/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,8 +355,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 ------- @@ -250,7 +362,7 @@ def _make_coords_initial( 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, diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py index f93d08e40c..4afce4f42e 100644 --- a/tidy3d/components/simulation.py +++ b/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) diff --git a/tidy3d/plugins/mode/mode_solver.py b/tidy3d/plugins/mode/mode_solver.py index b5399a0ca7..6bbc585545 100644 --- a/tidy3d/plugins/mode/mode_solver.py +++ b/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)