"""Automatic simulation setup: grid, PML, boundaries, and budget checks.
The goal is one function that turns ``(wavelength, domain, physics hints)``
into a trustworthy :class:`~fdtdx.config.SimulationConfig` with the same
"behind the scenes" defaults Tidy3D provides:
* resolution from points-per-wavelength (or an explicit target ``dl``)
* per-axis cell counts adjusted to integers — Tidy3D adjusts ``dl`` per axis
so the domain is covered by an integer number of cells; FDTDX mirrors that
with a :class:`~fdtdx.core.grid.QuasiUniformGrid`
* PML thickness in cells, uniform by default, per-face overridable
* a voxel-budget check with a clear error instead of an OOM crash
The trust argument: if one automatic function reproduces every golden in the
benchmark suite, its defaults are trustworthy. See ``benchmarks/`` and the
parity runners, which use these helpers.
"""
from __future__ import annotations
import math
import warnings
from collections.abc import Sequence
from itertools import pairwise
from typing import TYPE_CHECKING, Literal
import jax.numpy as jnp
import numpy as np
from fdtdx.config import SimulationConfig
from fdtdx.constants import c, eps0
from fdtdx.core.grid import QuasiUniformGrid, RectilinearGrid
if TYPE_CHECKING:
from fdtdx.objects.boundaries.initialization import BoundaryConfig
#: Default voxel budget (matches the benchmark hardware: RTX 3080 12GB).
DEFAULT_MAX_VOXELS = 40_000_000
#: Default points per wavelength for the auto resolution.
DEFAULT_PPW = 20.0
#: Default PML thickness in cells.
DEFAULT_PML_LAYERS = 10
BoundaryType = Literal["pml", "periodic", "pec", "pmc", "bloch"]
def snap_half_integer(scaled: float) -> int:
"""Ceil ``scaled`` — the cell-count rule Tidy3D uses.
Tidy3D covers the domain with ``ceil(size / dl)`` cells and adjusts the
per-axis spacing to ``size / n`` (verified empirically: 0.22um at 50nm ->
5 cells of 44nm, 2.5um at 40nm -> 63 cells of 39.68nm, 1.5um at 40nm ->
38 cells of 39.47nm). Ceil is monotonic, so float representation error
(e.g. ``2.5um / 40nm = 62.50000000000001``) cannot flip the result. The
relative epsilon (1e-5) absorbs the float32 roundoff of a resolved grid's
spacing — the spacing is derived from float32 edge differences, whose
catastrophic cancellation at large edge magnitudes reaches ~1e-6..1e-5
relative error (e.g. a 100-cell 50nm grid reports 4.999993e-08) — while
leaving genuine non-integer ratios (4.4, 62.5, ...) untouched.
Args:
scaled: The ``length / spacing`` ratio.
Returns:
The integer cell count covering the length.
"""
return math.ceil(scaled * (1.0 - 1e-5))
[docs]
def auto_pml_layers(
physical_thickness: float,
spacing: float | tuple[float, float, float],
) -> int:
"""Choose a uniform PML cell count that preserves physical thickness.
Refining a mesh must not silently make its absorber thinner. For an
anisotropic grid, the smallest spacing controls the uniform layer count,
so every face is at least ``physical_thickness`` thick (up to the same
floating-point guard used by :func:`snap_half_integer`).
"""
if physical_thickness <= 0:
raise ValueError(
f"PML physical thickness must be positive, got {physical_thickness}."
)
spacings = spacing if isinstance(spacing, tuple) else (spacing,)
if any(value <= 0 for value in spacings):
raise ValueError(f"Grid spacings must be positive, got {spacings}.")
return max(1, snap_half_integer(physical_thickness / min(spacings)))
[docs]
def auto_grid(
wavelength: float,
domain_size: tuple[float, float, float],
*,
ppw: float = DEFAULT_PPW,
dl: float | None = None,
max_voxels: int = DEFAULT_MAX_VOXELS,
max_refractive_index: float = 1.0,
min_feature_size: float | None = None,
min_cells_per_feature: float = 3.0,
budget_behavior: Literal["raise", "coarsen"] = "raise",
) -> QuasiUniformGrid:
"""Build a grid for ``domain_size`` at ``wavelength / ppw`` resolution.
Each axis gets an integer cell count (``ceil(length / dl)`` with a small
roundoff guard) and the per-axis spacing is adjusted to
``length / n`` — the same adjustment Tidy3D's uniform grid performs, so a
domain of e.g. 0.22um at 40nm becomes 6 cells of 0.0367um instead of a
fractional 5.5 cells.
Args:
wavelength: Free-space wavelength in metres. Used only when ``dl`` is
not given.
domain_size: Physical ``(Lx, Ly, Lz)`` extent in metres.
ppw: Points per wavelength *inside the highest-index material* when
``max_refractive_index`` is supplied. Defaults to 20.
dl: Explicit target cell width in metres. Overrides ``ppw``.
max_voxels: Voxel budget. Raises when exceeded (instead of an OOM
crash later) unless ``budget_behavior="coarsen"``.
max_refractive_index: Largest refractive index relevant to the scene.
The wavelength-derived target spacing is
``wavelength / (ppw * max_refractive_index)``. Defaults to 1 so
existing free-space behavior is preserved.
min_feature_size: Optional smallest geometrical feature that must be
represented. At least ``min_cells_per_feature`` cells are targeted.
min_cells_per_feature: Target cells across ``min_feature_size``.
Defaults to 3.
budget_behavior: ``"raise"`` preserves the requested accuracy and
fails before allocation. ``"coarsen"`` chooses the finest global
spacing that fits ``max_voxels`` and emits a warning. This explicit
warning is important because global coarsening may under-resolve a
requested material wavelength or feature.
Returns:
A :class:`~fdtdx.core.grid.QuasiUniformGrid` with per-axis spacings
adjusted so every axis has an integer cell count.
Raises:
ValueError: If the voxel count exceeds ``max_voxels``.
"""
if max_refractive_index < 1.0:
raise ValueError(f"max_refractive_index must be >= 1, got {max_refractive_index}.")
if ppw <= 0:
raise ValueError(f"Points per wavelength must be positive, got {ppw}.")
if min_cells_per_feature <= 0:
raise ValueError(f"min_cells_per_feature must be positive, got {min_cells_per_feature}.")
if min_feature_size is not None and min_feature_size <= 0:
raise ValueError(f"min_feature_size must be positive, got {min_feature_size}.")
if budget_behavior not in {"raise", "coarsen"}:
raise ValueError(f"budget_behavior must be 'raise' or 'coarsen', got {budget_behavior!r}.")
if dl is None:
material_dl = wavelength / (ppw * max_refractive_index)
feature_dl = math.inf if min_feature_size is None else min_feature_size / min_cells_per_feature
dl = min(material_dl, feature_dl)
if dl <= 0:
raise ValueError(f"Cell width must be positive, got {dl}.")
counts = [snap_half_integer(length / dl) for length in domain_size]
counts = [max(n, 1) for n in counts]
n_voxels = math.prod(counts)
if n_voxels > max_voxels:
if budget_behavior == "raise":
raise ValueError(
f"Auto grid for domain {domain_size} at dl={dl:.3g} needs {n_voxels:,} voxels "
f"(shape {tuple(counts)}), exceeding the budget of {max_voxels:,}. "
"Decrease ppw (coarser), reduce the domain, use symmetry/local refinement, "
"or raise max_voxels."
)
requested_dl = dl
# Solve the continuous cubic estimate first, then move upward until
# the integer ceil rule also fits the exact budget.
dl = max(dl, (math.prod(domain_size) / max_voxels) ** (1.0 / 3.0))
counts = [max(snap_half_integer(length / dl), 1) for length in domain_size]
while math.prod(counts) > max_voxels:
dl *= 1.001
counts = [max(snap_half_integer(length / dl), 1) for length in domain_size]
warnings.warn(
f"Auto grid coarsened from requested dl={requested_dl:.3g} m to about {dl:.3g} m "
f"to fit {math.prod(counts):,}/{max_voxels:,} voxels. Check wavelength and feature convergence.",
UserWarning,
stacklevel=2,
)
spacings = tuple(length / n for length, n in zip(domain_size, counts, strict=True))
return QuasiUniformGrid(dx=spacings[0], dy=spacings[1], dz=spacings[2])
[docs]
def auto_interface_aligned_grid(
domain_size: tuple[float, float, float],
target_spacing: float | tuple[float, float, float],
*,
interfaces: tuple[Sequence[float], Sequence[float], Sequence[float]] = ((), (), ()),
center: tuple[float, float, float] = (0.0, 0.0, 0.0),
max_voxels: int = DEFAULT_MAX_VOXELS,
) -> RectilinearGrid:
"""Build a rectilinear mesh whose cells end exactly at material interfaces.
Each interval between the domain boundary and supplied interface planes is
divided into the fewest equal cells whose width does not exceed
``target_spacing``. This removes mesh-origin changes in thin films while
keeping the mesh quasi-uniform and the cell count predictable. Interface
coordinates are absolute physical coordinates in the same frame as
``center``.
This policy is particularly useful for layered photonics: pass the slab,
etch, and substrate planes on their normal axis and let subpixel smoothing
handle curved or oblique in-plane boundaries.
"""
if isinstance(target_spacing, tuple):
spacings = target_spacing
else:
spacings = (target_spacing, target_spacing, target_spacing)
if any(length <= 0 for length in domain_size):
raise ValueError(f"Domain sizes must be positive, got {domain_size}.")
if any(spacing <= 0 for spacing in spacings):
raise ValueError(f"Target spacings must be positive, got {spacings}.")
edge_arrays: list[np.ndarray] = []
for axis in range(3):
lower = center[axis] - 0.5 * domain_size[axis]
upper = center[axis] + 0.5 * domain_size[axis]
tolerance = 1e-10 * domain_size[axis]
interior = sorted(float(value) for value in interfaces[axis])
for value in interior:
if not lower + tolerance < value < upper - tolerance:
raise ValueError(
f"Interface {value:.6g} on axis {axis} must lie strictly inside "
f"the domain [{lower:.6g}, {upper:.6g}]."
)
anchors = [lower]
for value in interior:
if value - anchors[-1] > tolerance:
anchors.append(value)
anchors.append(upper)
pieces: list[np.ndarray] = []
for segment_index, (start, stop) in enumerate(pairwise(anchors)):
count = max(1, snap_half_integer((stop - start) / spacings[axis]))
segment = np.linspace(start, stop, count + 1, dtype=np.float64)
if segment_index:
segment = segment[1:]
pieces.append(segment)
edge_arrays.append(np.concatenate(pieces))
shape = tuple(len(edges) - 1 for edges in edge_arrays)
voxels = math.prod(shape)
if voxels > max_voxels:
raise ValueError(
f"Interface-aligned grid needs {voxels:,} voxels (shape {shape}), "
f"exceeding the budget of {max_voxels:,}. Increase target_spacing, "
"reduce the domain, or raise max_voxels."
)
return RectilinearGrid.custom(
x_edges=jnp.asarray(edge_arrays[0]),
y_edges=jnp.asarray(edge_arrays[1]),
z_edges=jnp.asarray(edge_arrays[2]),
)
[docs]
def auto_boundary_config(
pml_layers: int = DEFAULT_PML_LAYERS,
boundary_types: dict[str, BoundaryType] | None = None,
structures: Sequence[object] | None = None,
domain_size: tuple[float, float, float] | None = None,
wavelength: float | None = None,
stabilize_evanescent: bool = False,
pml_alpha_fraction: float = 0.20,
) -> BoundaryConfig:
"""Uniform PML boundary config with per-face overrides and inference.
Defaults to PML on every face — the Tidy3D default. Pass
``boundary_types`` to override individual faces, e.g.
``{"min_x": "periodic", "max_x": "periodic"}``.
When ``structures`` and ``domain_size`` are given, any axis whose extent is
completely filled by a structure is inferred as **PEC** on both faces —
the convention the Tidy3D goldens use (z-extent = slab thickness with the
slab filling it). Explicit ``boundary_types`` always win over inference.
Args:
pml_layers: PML thickness in cells on every PML face.
boundary_types: Per-face boundary type overrides keyed by
``"min_x"/"max_x"/"min_y"/"max_y"/"min_z"/"max_z"``.
structures: Scene objects (e.g. waveguides, slabs) used to infer PEC
faces. Any object whose ``partial_real_shape`` fills an axis of
``domain_size`` marks that axis PEC.
domain_size: Physical ``(Lx, Ly, Lz)`` extent in metres. Required when
``structures`` is given.
wavelength: Reference free-space wavelength. Required when
``stabilize_evanescent`` is enabled.
stabilize_evanescent: Use a complex-frequency-shifted PML suitable for
resonators, waveguides, and photonic-crystal terminations with
strong evanescent or grazing fields. Material crossing a PML must
also be continued with :func:`fdtdx.extend_material_to_pml` after
geometry parameters are applied.
pml_alpha_fraction: CFS strength as a fraction of
``2 pi f epsilon_0`` at the PML interface. The 0.20 default was
selected by a late-time stability sweep of an Ez photonic-crystal
surface cavity; smaller values left a growing PML mode.
Returns:
A :class:`~fdtdx.objects.boundaries.initialization.BoundaryConfig`.
"""
# Imported lazily: the boundaries package imports from the fdtdx root, so a
# module-level import here would create a circular import.
from fdtdx.objects.boundaries.initialization import BoundaryConfig
inferred: dict[str, BoundaryType] = {}
if structures is not None:
if domain_size is None:
raise ValueError("domain_size is required when structures is given.")
# The SimulationVolume is the domain itself — it always fills every
# axis and must not trigger PEC inference.
from fdtdx.objects.static_material.static import SimulationVolume
for obj in structures:
if isinstance(obj, SimulationVolume):
continue
shape = getattr(obj, "partial_real_shape", None)
if shape is None:
continue
for axis, length in enumerate(domain_size):
obj_len = shape[axis]
if obj_len is not None and abs(obj_len - length) <= 1e-9 * max(1.0, length):
face_min, face_max = ("min_x", "max_x", "min_y", "max_y", "min_z", "max_z")[2 * axis : 2 * axis + 2]
inferred.setdefault(face_min, "pec")
inferred.setdefault(face_max, "pec")
merged = dict(inferred)
if boundary_types:
merged.update(boundary_types)
if pml_alpha_fraction <= 0:
raise ValueError("pml_alpha_fraction must be positive")
if stabilize_evanescent and (wavelength is None or wavelength <= 0):
raise ValueError("a positive wavelength is required to stabilize evanescent PML fields")
alpha_start = None
alpha_end = None
if stabilize_evanescent:
assert wavelength is not None
alpha_start = pml_alpha_fraction * 2.0 * np.pi * c / wavelength * eps0
alpha_end = 0.0
return BoundaryConfig.from_uniform_bound(
thickness=pml_layers,
boundary_type="pml",
override_types=merged,
alpha_start=alpha_start,
alpha_end=alpha_end,
)
[docs]
def auto_config(
wavelength: float,
domain_size: tuple[float, float, float],
time: float,
*,
ppw: float = DEFAULT_PPW,
dl: float | None = None,
max_voxels: int = DEFAULT_MAX_VOXELS,
max_refractive_index: float = 1.0,
min_feature_size: float | None = None,
min_cells_per_feature: float = 3.0,
budget_behavior: Literal["raise", "coarsen"] = "raise",
dtype: jnp.dtype = jnp.float32,
courant_factor: float = 0.99,
) -> SimulationConfig:
"""One-call simulation config: auto grid + defaults.
Args:
wavelength: Free-space wavelength in metres (resolution basis).
domain_size: Physical ``(Lx, Ly, Lz)`` extent in metres.
time: Total simulation time in seconds.
ppw: Points per wavelength. Defaults to 20.
dl: Explicit target cell width in metres. Overrides ``ppw``.
max_voxels: Voxel budget (see :func:`auto_grid`).
max_refractive_index: Highest scene refractive index used to resolve
wavelength in material.
min_feature_size: Optional smallest geometry feature to resolve.
min_cells_per_feature: Target cells across the smallest feature.
budget_behavior: Whether an over-budget request raises or explicitly
coarsens with a warning.
dtype: Field dtype. Defaults to float32.
courant_factor: Courant safety factor. Defaults to 0.99.
Returns:
A :class:`~fdtdx.config.SimulationConfig` ready for
:func:`fdtdx.place_objects`. Pair with :func:`auto_boundary_config`
for the boundary objects.
"""
grid = auto_grid(
wavelength=wavelength,
domain_size=domain_size,
ppw=ppw,
dl=dl,
max_voxels=max_voxels,
max_refractive_index=max_refractive_index,
min_feature_size=min_feature_size,
min_cells_per_feature=min_cells_per_feature,
budget_behavior=budget_behavior,
)
return SimulationConfig(
grid=grid,
time=time,
dtype=dtype,
courant_factor=courant_factor,
)
