diff --git a/tests/test_grid.py b/tests/test_grid.py index 99a61dd36f..978045675c 100644 --- a/tests/test_grid.py +++ b/tests/test_grid.py @@ -54,6 +54,20 @@ def test_sim_grid(): assert np.all(b == np.array([-2, -1, 0, 1, 2])) +def test_sim_pml_grid(): + + sim = td.Simulation( + size=(4, 4, 4), + grid_size=(1, 1, 1), + pml_layers=(td.PML(num_layers=2), td.Absorber(num_layers=2), td.StablePML(num_layers=2)), + ) + + for c in sim.grid.centers.dict().values(): + assert np.all(c == np.array([-3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5])) + for b in sim.grid.boundaries.dict().values(): + assert np.all(b == np.array([-4, -3, -2, -1, 0, 1, 2, 3, 4])) + + def test_sim_discretize_vol(): sim = td.Simulation( diff --git a/tidy3d/components/monitor.py b/tidy3d/components/monitor.py index 89b46417df..5d7d6286b2 100644 --- a/tidy3d/components/monitor.py +++ b/tidy3d/components/monitor.py @@ -57,8 +57,8 @@ class TimeMonitor(Monitor, ABC): @pydantic.validator("stop", always=True) def stop_greater_than_start(cls, val, values): - """ make sure stop is greater than or equal to start""" - start = values.get('start') + """make sure stop is greater than or equal to start""" + start = values.get("start") if val and val < start: raise SetupError("Monitor start time is greater than stop time.") return val diff --git a/tidy3d/components/pml.py b/tidy3d/components/pml.py index 16e32474cd..345291c395 100644 --- a/tidy3d/components/pml.py +++ b/tidy3d/components/pml.py @@ -1,11 +1,59 @@ """ Defines profile of Perfectly-matched layers (absorber) """ -from typing import Union +from typing import Union, Literal import pydantic from .base import Tidy3dBaseModel +"""TODO: better docstrings.""" + + +class AbsorberParams(Tidy3dBaseModel): + """Specifies parameters for an Absorber or PML. Sigma is in units of 2*EPSILON_0/dt.""" + + sigma_order: pydantic.NonNegativeInt + sigma_min: pydantic.NonNegativeFloat + sigma_max: pydantic.NonNegativeFloat + + +class PMLParams(AbsorberParams): + """Extra parameters needed for complex frequency-shifted PML. Kappa is dimensionless, alpha + is in the same units as sigma.""" + + kappa_order: pydantic.NonNegativeInt + kappa_min: pydantic.NonNegativeFloat + kappa_max: pydantic.NonNegativeFloat + alpha_order: pydantic.NonNegativeInt + alpha_min: pydantic.NonNegativeFloat + alpha_max: pydantic.NonNegativeFloat + + +AbsorberPs = AbsorberParams(sigma_order=3, sigma_min=0.0, sigma_max=6.4) +StandardPs = PMLParams( + sigma_order=3, + sigma_min=0.0, + sigma_max=1.5, + kappa_order=3, + kappa_min=1.0, + kappa_max=3.0, + alpha_order=1, + alpha_min=0.0, + alpha_max=0.0, +) +StablePs = PMLParams( + sigma_order=3, + sigma_min=0.0, + sigma_max=1.0, + kappa_order=3, + kappa_min=1.0, + kappa_max=5.0, + alpha_order=1, + alpha_min=0.0, + alpha_max=0.9, +) + + class AbsorberSpec(Tidy3dBaseModel): """Specifies the absorber along a single dimension.""" @@ -19,6 +67,8 @@ class PML(AbsorberSpec): ---------- num_layers : ``int``, optional Number of layers of PML to add to + and - boundaries, default = 12. + pml_params : :class:PMLParams + Parameters of the complex frequency-shifted absorption poles. Example ------- @@ -26,6 +76,7 @@ class PML(AbsorberSpec): """ num_layers: pydantic.NonNegativeInt = 12 + parameters: PMLParams = StandardPs class StablePML(AbsorberSpec): @@ -42,6 +93,7 @@ class StablePML(AbsorberSpec): """ num_layers: pydantic.NonNegativeInt = 40 + parameters: Literal[StablePs] = StablePs class Absorber(AbsorberSpec): @@ -58,6 +110,7 @@ class Absorber(AbsorberSpec): """ num_layers: pydantic.NonNegativeInt = 40 + parameters: AbsorberParams = AbsorberPs PMLTypes = Union[PML, StablePML, Absorber, None] diff --git a/tidy3d/components/simulation.py b/tidy3d/components/simulation.py index 2f065a5c2f..b1efd5d843 100644 --- a/tidy3d/components/simulation.py +++ b/tidy3d/components/simulation.py @@ -9,7 +9,7 @@ from mpl_toolkits.axes_grid1 import make_axes_locatable from descartes import PolygonPatch -from .types import Symmetry, Ax, Shapely, FreqBound +from .types import Symmetry, Ax, Shapely, FreqBound, Numpy from .geometry import Box from .grid import Coords1D, Grid, Coords from .medium import Medium, MediumType, eps_complex_to_nk @@ -271,16 +271,25 @@ def plot_symmetries( return ax @property - def pml_thicknesses( - self, - ) -> List[Tuple[float, float]]: - """thicknesses (um) of PML in all three axis and directions (+, -)""" - grid_sizes = self.grid.cell_sizes - pml_thicknesses = [] + def num_pml_layers(self) -> List[Tuple[float, float]]: + """Number of PML layers in all three axes and directions (-, +).""" + num_layers = [] for pml_axis, pml_layer in enumerate(self.pml_layers): - sizes_axis = grid_sizes.dict()["xyz"[pml_axis]] - num_layers = pml_layer.num_layers - pml_thicknesses.append((num_layers * sizes_axis[0], num_layers * sizes_axis[-1])) + 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 + + @property + def pml_thicknesses(self) -> List[Tuple[float, float]]: + """Thicknesses (um) of PML in all three axes and directions (-, +)""" + num_layers = self.num_pml_layers + pml_thicknesses = [] + for boundaries in self.grid.boundaries.dict().values(): + thick_l = boundaries[num_layers[0]] - boundaries[0] + thick_r = boundaries[-1] - boundaries[-1 - num_layers[1]] + pml_thicknesses.append((thick_l, thick_r)) return pml_thicknesses @add_ax_if_none @@ -400,16 +409,34 @@ def tmesh(self) -> Coords1D: def grid(self) -> Grid: """:class:`Grid` interface to the spatial locations in Simulation""" cell_boundary_dict = {} - for key, dl, center, size in zip("xyz", self.grid_size, self.center, self.size): + 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: num_cells = int(np.floor(size / dl)) + # Make sure there's at least one cell + num_cells = max(num_cells, 1) size_snapped = dl * num_cells if size_snapped != size: log.warning(f"dl = {dl} not commensurate with simulation size = {size}") 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) cell_boundary_dict[key] = bound_coords boundaries = Coords(**cell_boundary_dict) return Grid(boundaries=boundaries) + @staticmethod + def add_pml_to_bounds(num_layers: Tuple[int, int], bounds: Numpy): + """Append PML pixels at the beginning and end of bounds.""" + if bounds.size < 2: + return bounds + + first_step = bounds[1] - bounds[0] + last_step = bounds[-1] - bounds[-2] + add_left = bounds[0] - first_step * np.arange(num_layers[0], 0, -1) + add_right = bounds[-1] + last_step * np.arange(1, num_layers[1] + 1) + new_bounds = np.concatenate((add_left, bounds, add_right)) + + return new_bounds + @property def wvl_mat_min(self) -> float: """minimum wavelength in the material""" diff --git a/tidy3d/log.py b/tidy3d/log.py index 7968681ea7..0bf7df6c8e 100644 --- a/tidy3d/log.py +++ b/tidy3d/log.py @@ -9,7 +9,7 @@ log = logging.getLogger("rich") -class Tidy3dError(Exception): +class Tidy3DError(Exception): """any error in tidy3d""" def __init__(self, message: str = None): @@ -18,25 +18,25 @@ def __init__(self, message: str = None): super().__init__(self, message) -class ValidationError(Tidy3dError): +class ValidationError(Tidy3DError): """error when constructing tidy3d components""" -class SetupError(Tidy3dError): +class SetupError(Tidy3DError): """error regarding the setup of the components (outside of domains, etc)""" -class FileError(Tidy3dError): +class FileError(Tidy3DError): """error reading or writing to file""" -class WebError(Tidy3dError): +class WebError(Tidy3DError): """error with the webAPI""" -class AuthenticationError(Tidy3dError): +class AuthenticationError(Tidy3DError): """error authenticating a user through webapi webAPI""" -class DataError(Tidy3dError): +class DataError(Tidy3DError): """error accessing data"""