Simulation setup and ports¶

A bifrost simulation is an ordinary fdtdx simulation in which the device, the sources and the monitors come from the StackDevice: it knows where its ports are, so you ask for a source or monitor on a port instead of computing coordinates. The example below sets one up and places it without running the solve; it finishes in seconds on a CPU.

Build, place and inspect¶

The device and the numerical settings:

import fdtdx
import jax
import gdsfactory as gf
import bifrost as bf
from gdsfactory.technology import LayerLevel, LayerStack

gf.gpdk.PDK.activate()

component = gf.components.straight(length=5)  # a 5 um long, 500 nm wide waveguide with ports o1 and o2

layer_stack = LayerStack(
    layers={"core": LayerLevel(layer=gf.gpdk.LAYER.WG, zmin=0.0, thickness=0.22, material="core", mesh_order=1)}
)
materials = {
    "core": fdtdx.Material(permittivity=3.48**2),
    "clad": fdtdx.Material(permittivity=1.444**2),
}
stack = bf.StackDevice.from_layer_stack(
    component, layer_stack, materials, background="clad", zrange_layers=["core"],
    extend_ports=(("o1", 2.0), ("o2", 2.0)),          # the leads end at the box faces ...
    # ... which are 1.6 um from the guide in y and 1.59 um from the core in z. The box is then
    # 9 x 3.7 x 3.4 um: a whole number of cells at the 50 nm grid below.
    boundary_offset=(0, 0, 1.6, -1.6, 1.59, -1.59),
    name="stack",
)

RESOLUTION = 0.05  # um, coarse enough for a quick check
config = stack.get_simulation_config(                                # fdtdx.SimulationConfig, in the device's precision
    time=100e-15,                                                    # seconds of simulated time
    grid=fdtdx.core.grid.UniformGrid(spacing=RESOLUTION * 1e-6),     # the cell edge, in meters
)
wave = fdtdx.WaveCharacter(wavelength=1.55e-6)                       # fdtdx works in meters; bifrost's lengths are um
pulse = fdtdx.GaussianPulseProfile(
    center_wave=wave, spectral_width=fdtdx.WaveCharacter(frequency=0.05 * wave.get_frequency())
)

• Units. bifrost's lengths (boundary_offset, extend_ports, port_offset) are in µm. fdtdx's are in meters: SimulationConfig(time=...) is in seconds, UniformGrid(spacing=...) and WaveCharacter(wavelength=...) are in meters. The * 1e-6 conversions are the usual source of a factor-of-a-million mistake.
• The grid. UniformGrid(spacing=...) is the cell edge. fdtdx rounds the box to a whole number of cells in every direction, so keep its sides on the grid (see Layer stacks and materials).
• The pulse. A Gaussian pulse with a spectral width of 5 % of the center frequency ramps up smoothly and lets one run cover a band; the monitors below record the center wavelength.
• The precision. stack.get_simulation_config(...) is fdtdx.SimulationConfig with the device's dtype filled in; see Precision.

The box, the source and the monitors, each placed on a port:

volume, constraints = stack.get_simulation_volume()   # the box, from the device's bounding box
objects = [volume, stack]
te0 = dict(mode_index=0, filter_pol="te")             # the fundamental TE mode
for obj, placing in (
    # a mode source 0.5 um outside port o1 (negative port_offset moves it out of the component)
    stack.create_mode_plane_source_at(
        "o1", (0.5, 0.3), wave_character=wave, temporal_profile=pulse, port_offset=-0.5, name="source", **te0,
    ),
    # mode monitors on both ports; "in" counts the mode travelling into the component
    stack.create_mode_overlap_detector_at(
        "o1", (0.5, 0.3), direction="in", wave_characters=(wave,), scaling_mode="pulse", name="in", **te0,
    ),
    stack.create_mode_overlap_detector_at(
        "o2", (0.5, 0.3), wave_characters=(wave,), scaling_mode="pulse", name="out", **te0,
    ),
):
    objects.append(obj)
    constraints += placing

The absorbing boundary and the placement:

boundaries, placing = fdtdx.boundary_objects_from_config(
    fdtdx.BoundaryConfig.from_uniform_bound(thickness=round(1.0 / RESOLUTION)), volume   # 1 um of PML
)
objects += boundaries.values()
constraints += placing

key = jax.random.PRNGKey(0)
objects, arrays, params, config, _ = fdtdx.place_objects(objects, config, constraints, key)
print("grid cells:", objects.volume.grid_shape)
for monitor in (*objects.sources, *objects.detectors):
    print(f"{monitor.name:8} {monitor.grid_shape} cells at {[(s.start, s.stop) for s in monitor.grid_slice]}")

Output:

grid cells: (180, 74, 68)
source   (1, 41, 23) cells at [(30, 31), (16, 57), (22, 45)]
in       (1, 41, 23) cells at [(40, 41), (16, 57), (22, 45)]
out      (1, 41, 23) cells at [(140, 141), (16, 57), (22, 45)]

The box is 180 × 74 × 68 cells of 50 nm. Each slice is (start, stop) cell indices along x, y and z: each plane is one cell thick along its normal and 41 × 23 cells (2.05 × 1.15 µm) across, the 0.5 µm guide plus 0.5 wavelengths per side, and the 0.22 µm core plus 0.3 wavelengths per side. The source is 10 cells (0.5 µm) before the input monitor, which sits on the port at x = 0.

Sources and monitors on ports¶

stack.get_simulation_volume() returns the simulation box and its placement constraints. Each stack.create_<object>_at(port, port_padding, ...) method returns an fdtdx object together with the constraints that place it on that port, centered and normal to the port's direction. Add both to your lists:

for obj, placing in (stack.create_mode_plane_source_at("o1", (0.5, 0.3), ...),):
    objects.append(obj)
    constraints += placing

The arguments every method shares:

The methods and the fdtdx objects they make:

Other fdtdx fields go in as keyword arguments. Boundaries, regular fdtdx devices and other objects are added as usual.

Place, apply and run¶

After fdtdx.place_objects, bf.apply_stack writes the stack into the permittivity arrays, in place of fdtdx.apply_params (the waveguide bend runs this end to end):

objects, arrays, params, config, _ = fdtdx.place_objects(objects, config, constraints, key)
arrays, objects, _ = bf.apply_stack(arrays, objects, {**params, **bf.init_stack_params(objects)}, key)
_, arrays = fdtdx.run_fdtd(arrays=arrays, objects=objects, config=config, key=key)

fdtdx's own Device objects only mix two materials, so fdtdx.apply_params doesn't know about a StackDevice. apply_stack writes every stack into the arrays first and then calls fdtdx.apply_params, so regular devices and sources are still updated, and devices win where they overlap a stack. bf.init_stack_params(objects) supplies the stack's own parameters (vertices, heights, thicknesses), which is what you differentiate; see Gradients.

Graded grids¶

The cells can have different sizes along every axis: fine cells where the layers and the edges are, coarse ones in the cladding around them. bf.graded_edges places the edges of one axis and StackDevice.rectilinear_grid turns them into fdtdx's grid over the box, here graded along z:

half = 0.5 * stack.partial_real_shape[2] * 1e6  # the box's half height (um)
z = bf.graded_edges(
    stack.z_center - half, stack.z_center + half, 0.06,  # at most 60 nm cells
    fine=[(-0.05, 0.27, 0.02)],    # 20 nm around the 220 nm core
    snap=[0.0, 0.22],              # the core's faces on cell edges
    uniform_ends=(10, 10),         # 10 cells of 60 nm under each PML, whose thickness fdtdx counts in cells
)
config = stack.get_simulation_config(time=100e-15, grid=stack.rectilinear_grid(0.06, z_um=z))
source, placing = stack.create_mode_plane_source_at("o1", (0.5, 0.5), wave, grid=config)

• x and y may be graded too, with any averaging: pass x_um / y_um edges to rectilinear_grid (from graded_edges, e.g. fine across a gap and coarse along the waveguides, or uniform ones of a size other than spacing_um). Every column of voxels is painted on its own lateral cell, a rectangle of the cell's aspect in a frame scaled the same way along x and y, so normals, sidewall slopes and round corners stay Euclidean; the corner model and xy_subsamples work per column, and with "subpixel_full" a junction voxel's closure solves its cell problem on the voxel's own box. Rough surfaces on graded x or y are drawn on an even lattice at the finest cell and interpolated to the columns, so they come out slightly smoothed in the coarse ones. Keep the cells within a factor 16 of square (a placement warning says when they are not).
• The smallest cell sets the time step of the whole simulation. Grading saves cells and memory (also the gradient's checkpoints), not time steps. graded_edges keeps neighbouring cells within a factor growth (1.25) and does not make slivers: a snap coordinate within half a cell of another edge is dropped, with a warning.
• Pass grid= to the create_*_at methods (the config or the grid). fdtdx converts a plane's physical size to cells from the bottom of the domain, which is wrong where the cells differ; with grid= the plane's faces snap to the nearest edges at its own height. The device must fill the simulation volume (get_simulation_volume, any boundary_offset); a grid that spans more than its box is refused.
• Build the config with get_simulation_config. fdtdx scales the fields by the cell widths in the dtype of the grid's edges, so a float32 simulation needs float32 edges; rectilinear_grid makes them in the device's dtype and get_simulation_config(grid=...) converts those of a RectilinearGrid built otherwise. Placing a device with a config whose graded edges are in another dtype raises a ValueError.
• Fill fractions stay exact. A level's part of a cell is its overlap with the cell, and each field component averages over the cells between neighbouring cell centers (the dual cells) along the axes fdtdx staggers it on: Ex along y and z, Ey along x and z, Ez along x and y.
• Optimizations compile once on a fixed graded grid, like on a uniform one; a graded grid's step costs what its cells cost. In benchmarks/rectilinear a sloped silicon strip graded in y and z gives the uniform 20 nm grid's effective index within 7e-5 with 27 % of its cells, and the inverse-design bend graded to 80 nm away from the guide gives the uniform 40 nm grid's transmission within 4.4e-4 with 63 %. Keep the fine cells over the mode's evanescent tail, not only over the core: cells that grow too early (120 nm in the bend's silica) cost accuracy.

Precision¶

The whole simulation runs in the StackDevice's dtype, jnp.float32 (the default) or jnp.float64, and the device is the only place to set it:

stack = bf.StackDevice.from_layer_stack(
    component, layer_stack, materials, background="clad", zrange_layers=["core"], dtype=jnp.float64
)
config = stack.get_simulation_config(time=100e-15, grid=fdtdx.core.grid.UniformGrid(spacing=50e-9))

• get_simulation_config builds the fdtdx.SimulationConfig with the device's dtype and passes every other argument (time, grid, backend, gradient_config, ...) on. Passing dtype= to it is an error. A config you build yourself must use the device's dtype as well, or placing the device raises a ValueError.
• The create_*_at monitors follow the device: phasor, mode-overlap and diffractive monitors are complex64 or complex128, field, energy and Poynting-flux monitors float32 or float64, unless you pass dtype= yourself. fdtdx's sources already follow the config.
• Nothing needs switching around the run. Importing gdsfactory turns JAX's 64-bit mode on, and bf.apply_stack keeps fdtdx's absorbing-boundary (PML) coefficients in the device's dtype, so fdtdx.run_fdtd works in either precision. Scripts written before this wrapped the run in jax.enable_x64(False); that still works for float32, but a float64 device refuses to be placed while the 64-bit mode is off.
• float64 doubles the memory of every array and is slow on GPUs built for single precision. fdtdx 0.6.3 still computes a few auxiliary quantities in float32 (the time window of phasor monitors and the face weights of mode-overlap monitors), so a float64 run has float64 fields and permittivity, not every intermediate in double precision.

Troubleshooting¶

• Configured grid shape ... does not match simulation volume shape: a side of the box is exactly halfway between two cell counts. Move boundary_offset by a fraction of a cell, or change the grid spacing.
• A source or monitor is in the wrong place or clipped. Print monitor.grid_slice as above and compare with the box's cell count. Check the sign of port_offset (negative moves the plane out of the component) and direction.
• The result looks a million times off. A length went into fdtdx in µm. fdtdx wants meters.