photonfdtd 0.12.0

Local-only FDTD + waveguide mode solver for integrated photonics (NumPy core; optional JAX and compiled Rust/CUDA backends).

photonfdtd

A small, fully-local open-source FDTD + waveguide mode solver written in Python/NumPy. It implements a Yee-grid FDTD time-stepper with CPML absorbing boundaries, dispersive materials, anisotropic subpixel smoothing, and a full-vectorial waveguide mode solver, with an API intentionally similar in spirit to Tidy3D.

This is alpha-stage software (v0.1). The pieces that exist are tested and correct; the pieces that don't, don't. See Status below.

Install

pip install photonfdtd

or from source, either directly from GitHub:

pip install "git+https://github.com/ngpaladi/photonfdtd"

or from a checkout:

git clone https://github.com/ngpaladi/photonfdtd
cd photonfdtd
pip install -e .

The PyPI wheel is pure Python (NumPy core; photonfdtd[jax] for the differentiable XLA backend). The optional compiled Rust / CUDA backends (backend="rust", backend="rust-cuda") are built separately from the Rust sources shipped in the sdist — see benchmarks/README.md.

Releasing (maintainers)

Releases publish to PyPI automatically via .github/workflows/publish.yml using PyPI Trusted Publishing (OIDC — no API token stored). To cut a release:

1. Bump the version in pyproject.toml and src/photonfdtd/__init__.py (kept in sync by tests/test_version.py).
2. Publish a GitHub Release tagged vX.Y.Z matching that version. The workflow runs the tests, guards tag-vs-version, builds the sdist + wheel, and publishes.

One-time setup before the first release: on PyPI, add a trusted publisher for the project (Publishing → add a GitHub Actions publisher) with owner ngpaladi, repo photonfdtd, workflow publish.yml, environment pypi. Until the project exists, use PyPI's pending publisher form with the same values.

Quick start

A point dipole radiating in 2D vacuum

import photonfdtd as pf

lam0 = 1.0e-6
freq0 = pf.C_0 / lam0
dx = lam0 / 20

grid = pf.Grid(size=(4e-6, 4e-6), cell_size=dx, pml_layers=(12, 12, 0))
src = pf.PointDipole(
    position=(0.0, 0.0),
    component="Ez",
    waveform=pf.GaussianPulse(freq0=freq0, fwhm=10e-15),
)
mon = pf.FieldMonitor(name="snap", components=("Ez",), interval=20)

sim = pf.Simulation(grid, sources=[src], monitors=[mon], run_time=200e-15)
result = sim.run()

result.fields["snap"]["Ez"] is a (n_frames, ny, nz) array of field snapshots.

Solving a slab waveguide mode

import photonfdtd as pf

clad = pf.Medium.from_index(1.0)
core = pf.Medium.from_index(2.0)
slab = pf.Box(center=(0.0, 0.0), size=(0.5e-6, 0.3e-6), medium=core)

ms = pf.ModeSolver(
    size=(4e-6, 3e-6),
    cell_size=20e-9,
    structures=[slab],
    wavelength=1.55e-6,
    num_modes=2,
)
result = ms.solve()
print(result.n_eff)            # array of effective indices

What v0.2 actually does

• Yee-grid FDTD in 1D / 2D / 3D, with the Courant-stable time step selected automatically.
• CPML absorbing boundaries (Roden & Gedney 2000, kappa = 1) on any number of axes.
• Isotropic dielectric media, non-dispersive or dispersive. Non-dispersive media stamp a single eps_r per cell; dispersive media (DispersiveMedium, with .lorentz / .drude / .sellmeier constructors) carry Lorentz, Drude and Sellmeier poles advanced by the auxiliary-differential-equation (ADE) method during stepping. A cited material library (pf.silica(), pf.silicon(), pf.silicon_nitride(), pf.lithium_niobate(), pf.gold(), pf.silver()) reproduces each source's published index — see docs/material_data.md.
• Anisotropic subpixel smoothing (Simulation(subpixel=True)): partially filled interface cells get an effective diagonal permittivity tensor (harmonic mean normal to the interface, arithmetic mean tangential), removing the staircase error and restoring second-order boundary accuracy (Farjadpour et al. 2006). See docs/accuracy.rst.
• Geometry primitives: axis-aligned Box and arbitrary-polygon PolySlab (a polygon in xy extruded between two z bounds).
• Sources: soft point-dipole (PointDipole), distributed line/area mode injection (ModeSource), an energy-normalised SinglePhotonSource whose amplitude is set so the launched wavepacket carries exactly $h,\nu$ of total electromagnetic energy, and a moving-charge ChargedParticle current source that emits Cherenkov radiation when it outruns the local phase velocity (see examples/04_cherenkov.py).
• Monitors: time-domain field snapshots (FieldMonitor), Poynting flux (FluxMonitor), and a frequency-domain DFTMonitor that accumulates a running Fourier transform at chosen frequencies so storage scales with the number of frequencies rather than the number of timesteps (routinely 50-1000x smaller than an equivalent time-domain monitor, exact at those frequencies).
• Memory efficiency for large volumes: per-array float32/float64 precision, the permittivity grid released during stepping, a shared curl scratch buffer that avoids per-step full-domain temporaries, and FieldMonitor(compression=...), which streams snapshots to disk as per-frame-scaled, quantised, compressed blocks - keeping RAM flat regardless of recording length and shrinking stored data ~10-30x (8-bit) or ~5-7x (16-bit) versus an uncompressed float64 monitor.
• Out-of-core stepping (sim.run(out_of_core=True, tile_cells=...)): for a volume whose fields do not fit in RAM, the six field arrays, ce_field and the CPML state are memory-mapped to disk and each timestep sweeps the domain in slabs along one axis, so peak RAM is bounded by the tile size, not the grid. Reproduces the in-core result to machine precision. NumPy backend, point sources and FieldMonitor (incl. compression=) supported; see photonfdtd.outofcore.
• Differentiable PIC inverse design. A substrate/etched-core/cladding stack driven by a 2-D etch density ρ(x,y) (EtchedCore) maps — conic filter → tanh projection → anisotropic subpixel sidewalls → 3-D ε — entirely in JAX, so value_and_grad_density returns d(loss)/d(ρ) through the full time evolution (topology optimization; gradient matches finite differences to ~1e-8). Ports are launched one-way with UniModeSource (equivalence-principle current sheets). See docs/accuracy.rst.
• Mode-decomposition S-parameters (s_parameters, mode_amplitudes): project the frequency-domain fields on a memory-light port plane (any axis via DFTMonitor(plane_axis=..., plane_position=...)) onto a solved mode to get complex forward/backward modal amplitudes and hence S-parameters. A mode down a straight lossless guide gives |S21|≈1; the overlap is backend-agnostic (NumPy/CuPy/JAX) so |S_ij|² is differentiable through the adjoint — the PIC inverse-design objective. See docs/accuracy.rst.
• Full-vectorial waveguide mode solver (ModeSolver): a finite-difference eigenproblem for the transverse fields on the Yee grid that distinguishes TE from TM and captures high-index-contrast boundary effects. Validated against the exact symmetric-slab dispersion to better than 0.01% on both TE0 and TM0 (Zhu & Brown 2002; Fallahkhair et al. 2008).
• gdsfactory adapter (from_gdsfactory) that reads a layout Component, maps its layers onto user-supplied materials, and returns a pre-built Simulation.
• Optional acceleration backends for time-stepping: a Numba JIT path and a CuPy GPU path (use_gpu=True). The GPU path uses only generic CuPy array ops, so it runs on either an NVIDIA GPU (CuPy's CUDA build, cupy-cuda12x) or an AMD GPU (CuPy's ROCm build, cupy-rocm-5-0). Both are optional; without them the plain NumPy core runs.
• Differentiable JAX backend (use_jax=True): a pure-functional reimplementation of the Yee + CPML step under jax.lax.scan, JIT-compiled through XLA (CPU/GPU/TPU from one path) and matching the in-core solver to machine precision. Because the stepper is a pure function of the permittivity, pf.jax_value_and_grad_eps(sim, loss) returns the gradient of any monitor-based scalar w.r.t. eps_r — a time-domain adjoint for gradient-based inverse design / topology optimization. See photonfdtd.jaxbackend. (pip install "photonfdtd[jax]".)

Running on GPU

JAX is the GPU backend. With use_jax=True the forward run and the differentiable / reversible adjoints run on the GPU automatically once a CUDA jaxlib is present — no code change: pip install "jax[cuda12]". Validated on an NVIDIA RTX 4080 (CUDA 12): forward matches NumPy to ~1e-16, and both adjoints run on-device. If JAX falls back to CPU with a "cuSPARSE not found"-style error, the nvidia CUDA libraries just aren't on the loader path — add them, e.g. export LD_LIBRARY_PATH=$(python -c "import nvidia,glob,os;print(':'.join(glob.glob(os.path.dirname(nvidia.__file__)+'/*/lib')))"). Large differentiable runs fit in GPU memory thanks to gradient checkpointing and the O(1)-memory reversible adjoint (photonfdtd.reversible).

The CuPy backend (use_gpu=True) is a deprecated legacy path — superseded by JAX for GPU work, but retained (optional; pip install cupy-cuda12x) for AMD/ROCm GPUs and for GPU/host/disk out-of-core (use_gpu=True + run(out_of_core=True)), where fields live on disk and each tile is processed on the GPU so peak device memory is one tile — a volume larger than GPU memory still runs (JAX cannot stream to disk).

What v0.2 does not do (yet)

Items below are intentional out-of-scope and tracked in the issue tracker:

• Dispersive media (Lorentz / Drude / Debye) - planned next, via the auxiliary-differential-equation method.
• Full-vectorial 2D mode solving with E and H couplings.
• Anisotropic media.
• Sub-cell averaging for polygon edges.
• Total-field / scattered-field plane wave injection.

Roadmap

• v0.3: dispersive media (Lorentz pole-residue), full-vectorial mode solver, anisotropic materials.
• v0.4: Numba/JAX backend for ~10-100x speed-up.

License

MIT.