PhotonHub quickstart¶

Before you start

This guide requires photonhub 0.1.4 or later. Earlier releases lack ph.Domain, ph.Mesh, ph.Port, ph.GaussianBeam, the *_rad angle names, cloud wait_timeout_s and cloud submit(max_usd=...). With the environment below activated, check yours with python -m pip show photonhub and upgrade with python -m pip install -U photonhub. Local runs also need the separate phsolver binary; see Install.

PhotonHub is 3D FDTD for photonics, made by LeapField. Its Python client lets you describe a validated simulation, estimate cost, run on a local CPU or a configured GPU, and read results as labelled xarray arrays.

import photonhub as ph

This guide starts with a small CPU run. The examples check the workflow. Their short durations and coarse meshes do not establish converged physics results. For worked studies with stated validation status, see the notebook gallery.

Install¶

Python 3.11 or newer is required. Check the version before creating .venv. If Python is missing or the version is older than 3.11, install Python 3.11 or newer, open a new terminal, and check again. If .venv was created with an older Python, delete it before trying again; creating it over the old directory keeps the old interpreter. On macOS or Linux, run rm -rf .venv. On Windows, run rmdir /s /q .venv in Command Prompt or Remove-Item -Recurse -Force .venv in PowerShell. The .venv folder is hidden in Finder because its name starts with a dot.

The client is a normal Python package.

On macOS or Linux, check:

python3 --version

If the newer Python has a versioned name (for example python3.12), use that name in the version check and environment command instead of python3.

python3 -m venv .venv
source .venv/bin/activate
python -m pip install photonhub

On Windows, in Command Prompt or PowerShell, check:

py -3 --version

If py is not found, install Python from python.org with its launcher. py -3 --version must report 3.11 or newer.

py -3 -m venv .venv

Run Python as .venv\Scripts\python.exe, or activate the environment with .venv\Scripts\activate.bat in Command Prompt or .venv\Scripts\Activate.ps1 in PowerShell (which the default script policy may refuse); the path form always works.

.venv\Scripts\python.exe -m pip install photonhub

For this published package, pip on an older Python reports No matching distribution found for photonhub.

On Windows, run saved examples with .venv\Scripts\python.exe your_script.py.

The base install covers the core workflow and 2D plots. GDS import needs photonhub[gds]. Interactive 3D plots need photonhub[viz]; install that extra before the 3D section below.

To run locally you also need the solver binary phsolver. Installing or importing the SDK does not fetch an engine. An explicit install-solver command, in photonhub 0.1.4 or later, can download one. The headless guide shows it and a manual route.

ph.find_solver() checks, in order: an explicit solver_path=, PHOTONHUB_SOLVER, a recorded install or link, PATH, then an eligible build/phsolver in a source checkout. Invalid explicit paths and environment overrides raise SolverRunError instead of falling through.

print(ph.find_solver())  # the selected path, or None if none is available

Build a `Simulation`¶

A Simulation is the complete scene: domain size, grid, run controls, sources, monitors, and boundaries. Construction performs structural and cross-field validation, so malformed geometry or an out-of-domain monitor fails before a solver runs. The SDK also resolves the mesh, cell counts and timestep. It checks grid and stability constraints. sim.check_runnable() checks run readiness, including sources and boundary-layer overlap. The selected engine remains authoritative through phsolver validate; device memory and backend-specific limits also need engine validation or cloud preflight.

import numpy as np

f0 = 1.934e14                              # center frequency (Hz)

sim = ph.Simulation(
    size_um=(3.0, 3.0, 3.0),               # domain extents, microns
    grid=ph.UniformMesh(dl_um=0.06),   # uniform grid, 60 nm cells
    run=ph.RunSpec(n_steps=1500),          # or run_time_s=...
    pml_num_layers=10,                     # PML on every face is the default
    sources=[
        ph.PointDipole(
            center_um=(1.5, 1.5, 1.5),
            polarization="Ez",
            source_time=ph.GaussianPulse(freq0_hz=f0, fwidth_hz=0.5 * f0),
        ),
    ],
    monitors=[
        ph.PowerMonitor(name="flux", axis="z", position_um=2.1,
                       freqs_hz=tuple(np.linspace(1.2e14, 2.7e14, 31))),
        ph.SnapshotMonitor(name="snapshot", fields=("Ez",)),
    ],
)

Building blocks you will reach for most:

Geometry and materials: Box, Sphere, Cylinder and Polygon, each wrapped in a Structure with a Medium. Cylinder supports annular sectors with inner_radius_um, angle_start_rad and angle_stop_rad. Angles are in radians. Polygon describes an extruded polygon. Structures apply in list order, with the last one winning overlaps. The domain has a Background medium.
Sources: PointDipole, PlaneWave and ModeSource. A plane wave requires periodic or Bloch transverse axes; an oblique plane wave needs Bloch ones. The analysis functions build profiled mode and Gaussian-beam excitations as source planes.
Monitors: PowerMonitor measures power through a plane over frequency. ProfileMonitor records complex fields in a region. TimeMonitor records a time series, and SnapshotMonitor records full-domain fields. Mode monitors use profile fields to measure modal power.
Meshes: UniformMesh uses constant spacing. GradedMesh and auto_mesh refine selected regions. CPU and GPU paths support graded meshes, but tensor_full and contour_full subpixel methods require uniform meshes. See Known limitations.
Declarative setup: Domain, Mesh and RunSpec(transits=...) fit the domain, mesh and duration to the device and band. The ports section below shows this form. Explicit size_um, grid and duration also work.

Estimate the cost before you run¶

sim.cost_estimate() resolves cell counts, timestep and step count in Python. It estimates memory, output size, cost and wall time using planning inputs. It does not reserve a device or promise a runtime or price.

estimate = sim.cost_estimate()
print(estimate.summary())
print(estimate.num_cells, estimate.num_steps)

The work unit is Tcell-step: num_cells * num_steps / 1e12. The planning rate and throughput are overridable. Use a server quote before cloud submission. See Cloud GPU runs.

Auto-shutoff watches field-energy decay. On a capable local CPU solver with frequency monitors, it also estimates the remaining DFT or flux tail and waits until that estimate falls below dft_shutoff (its limits are under Silicon PIC devices); GPU, cloud and older solvers use the energy check alone. A broadband dipole pulse can leave a static residue, so a run may use its whole step cap. Narrow the pulse or set the cap deliberately. Reaching a cap is not evidence of convergence.

Run locally¶

ph.run_local writes the spec to a directory, runs phsolver, and loads the outputs into a RunResult:

dipole_sim = sim
data = ph.run_local(dipole_sim, output_dir="dipole-run", device="cpu")
dipole_data = data
print(data.monitor_names)                  # e.g. ['snapshot', 'flux']
print(data.steps_run, data.run["mcells_per_s"])

While it runs, a live status line is drawn to stderr with progress, the field decay against the shutoff threshold (so you can see how close the run is to finishing early), whether the sources are still injecting or ringing down, stability, throughput, and an ETA. Pass quiet=True to silence it, or a progress= callback to consume the raw event dicts yourself.

The default device is CPU. A GPU run requires a GPU-enabled solver or the cloud service. The packaged headless solver is CPU-only. CPU and GPU comparisons use numerical tolerances; exact equality is not a general promise.

Construction and preflight can raise ValueError before the solver starts. SolverRunError covers solver discovery, failed execution, divergence, a local timeout= (which stops the subprocess), and unreadable solver outputs. Filesystem and user callback errors can still propagate separately.

Local, cloud and desktop runs share the RunResult structure. To compare numerical results, match the engine revision, mesh, materials, subpixel settings, boundaries, precision, sources, monitors and run controls. data.provenance records solver_version and git_sha. A schema version alone does not identify the engine. Cloud API and solver images are versioned separately, and an older engine may return plausible, different numbers.

Read results¶

Each monitor is an xarray.DataArray with named dimensions and physical coordinates:

flux = data["flux"]                        # dims ('f',), f in Hz, plus a wlen_um coordinate (microns)
snap = data["snapshot"]                    # dims ('t','component','z','y','x')

ez = snap.isel(t=-1).sel(component="Ez")   # final-step Ez field
plane = ez.sel(z=1.5, method="nearest")    # a slice, in microns

Frequency-domain monitors (flux, DFT) are the continuous-wave response to the sources as declared: a PowerMonitor reads watts (a beam or mode launched with power_watts=1 reads 1 W on a full plane below it), a ProfileMonitor reads V/m and A/m. With symmetry= planes the simulation models part of the device, and these are still the whole device's watts, what the simulation without the planes reads: a launch centred on a symmetry plane puts half its power_watts into the modeled part per plane, and a PowerMonitor counts the mirrored parts back in (.attrs["symmetry_factor"]; data.wire(name) is the modeled part alone). The convention is recorded in .attrs["normalization"]. Ratios are amplitude-free either way; the cross-validation examples (Fresnel, Mie, waveguide) show how to turn paired runs into R, T, C_sca, and n_eff.

Sweeps and async¶

For parameter sweeps, ph.Batch runs a named collection together, prices the whole batch up front, and returns a dict-like BatchResults:

sim_a = sim                                          # the scene built above
sim_b = sim.with_changes(run=ph.RunSpec(n_steps=800))    # an edit is validated like a new build
batch = ph.Batch({"thin": sim_a, "thick": sim_b})
per_sim, total_usd = batch.quote()
results = batch.run(max_workers=2)         # 2 concurrent local subprocesses
for name, sim_data in results.items():     # successful runs only
    print(name, sim_data.steps_run)
print(results.errors)                       # named failures

A single run can also go async with ph.submit(sim).result(). The same sweep contract is available through ph.cloud.Batch when a compatible service is configured; placement and concurrency are service-defined. See cloud.md.

Silicon PIC devices: shapes in, mode-resolved transmission out¶

The device workflow combines geometry, a waveguide mode source and mode monitors. It returns mode-resolved transmission. This is one-mode-at-a-time power transmission (not an S-matrix), and it is broadband: a single Gaussian-pulse run gives the whole T(λ) spectrum via the running DFT.

A parametric component library¶

ph.library builds common PIC devices in a short call: straight, bend, taper, crossing, coupler, ring. Each returns a Component bundling its structures (drop into a Simulation) with its ports (the planes where a mode source / monitor attaches):

wg = ph.library.straight(length_um=10.0)   # one Box + 'in'/'out' ports
res = ph.library.ring(radius_um=5.0, gap_um=0.2)   # ring Cylinder + bus Box
for port in wg.ports:
    print(port.name, port.center_um, port.axis, port.width_um)

Builders place their geometry around center_um, the origin by default; with domain= (below) the box is fitted around the device wherever it sits.

Ports: mode launch and mode-resolved readout¶

A port is a plane on a waveguide where a mode is launched or read out. Give the Simulation its ports, the wavelengths to read out and the port to drive, and it solves each port's mode on its own grid at construction, launches the input mode a third of a wavelength behind the input plane (source_offset_um; a length, not a cell count, so the readout does not drift as you refine the mesh), and records the modal power at every port. data.transmission(name) is then the modal power into that port over the launched power, an xarray over f with a wlen_um coordinate:

# A 0.5 x 0.22 um silicon strip along z in oxide cladding, drawn around the origin.
si = ph.Medium(permittivity=3.48**2)
strip = ph.Structure(geometry=ph.Box(center_um=(0.0, 0.0, 0.0), size_um=(0.5, 0.22, 4.0)),
                     medium=si)
ports = [ph.Port("in", (0.0, 0.0, -0.6), "z", 0.5, thickness_um=0.22, thickness_axis="y"),
         ph.Port("out", (0.0, 0.0, 0.6), "z", 0.5, thickness_um=0.22, thickness_axis="y")]

sim = ph.Simulation(
    domain=ph.Domain(clearance_um=0.8),  # fit around the strip, with PML outside
    mesh=ph.Mesh(dl_um=0.1, uniform=True),  # coarse uniform mesh for this workflow check
    run=ph.RunSpec(n_steps=200),       # short workflow check, not a converged spectrum
    background=ph.Background(permittivity=1.44**2),
    structures=[strip],
    wlens_um=(1.50, 1.55, 1.60),       # the readout wavelengths; the pulse follows
    ports=ports,
    source="in",                       # TE0 into the input port
)
port_sim = sim
data = ph.run_local(port_sim, output_dir="port-run", device="cpu")
t = data.transmission("out")           # modal power ratio, dims ('f',) with wlen_um

domain= fits the box around what you drew: two thirds of a wavelength of cladding beyond the geometry on the closed sides (ph.Domain(clearance_um=...) sets it), ten cells of the guide's own cell beyond the launch and readout planes on the sides the guides leave through (port_margin_um), the PML outside; a guide that stops short of the wall is continued through it in its own width, thickness and medium. The coordinates stay yours: sim.origin_um records where the wire's corner sits in your frame, and every plot and every result coordinate is in your frame. mesh= builds the graded mesh at cells_per_wlen in the highest index (uniform=True for a uniform one, refine= for local overrides, dl_um= for a spacing given directly), and transits= gives the run length as passes over the domain; with no run at all, or RunSpec(shutoff=...) alone, the run is capped at 40 transits. It may end at auto-shutoff or reach the cap. With a frequency monitor, a capable local CPU solver also estimates how much its recorded values may still change. It uses dft_shutoff=1e-4 by default and never stops before the energy rule would. The estimate assumes one geometrically decaying mode per frequency. Single-mode cavity tests met the requested tolerance, but beating modes and late arrivals can stop early. A coupled-cavity measurement missed by 242 times its 1e-4 tolerance. A later stop is not always more accurate: at an unlucky beat phase the error can be about twice the energy-only error.

Measured single-mode stops took roughly 12 to 24 amplitude decay times, so the default 40-transit cap can still end a slow resonance. A tail that falls onto a static residue after block 16,384 can require up to 16,384 more progress blocks. A slowly decaying component can hold a single-bin monitor to its cap: the measured two-source run took 8,000,000 steps, about 14.7 times its energy-only stop. The measured 4|4 cavity at 1e-4 rose from 2,256,600 to 4,096,700 steps with no accuracy gain. The cap warning identifies the frequency that held the run and compares its estimated tail with the threshold. Increase transits or run_time_s for a longer run, or set dft_shutoff=1e-3 or 0 if a looser estimate or the energy-only rule fits your task. Inspect the decay and demonstrate duration convergence before interpreting spectra. size_um, grid and run_time_s remain the hand-built forms; a scene built with them is in the wire's corner frame. A structure's medium and the background may be materials-library entries (ph.materials.Si, ph.materials.SiO2): the simulation takes their index at the band centre and warns when it moves by more than 0.5 % across wlens_um.

A mirror-symmetric device is described whole. With domain=, a symmetry=(0, -1, 0) (an electric mirror normal to y; +1 is a magnetic one) names the plane through the device's centre, and the simulation runs the half on the far side: ports in the mirrored half are read through their images, a monitor crossing the plane is recorded on the kept half and comes back expanded to the full device, and every plot shows the whole device. On a hand-built size_um the symmetry field keeps its wire meaning, the half domain with the mirror on the low face.

A port names its plane (the centre's component along axis), the guide's width and thickness, and the mode as a family and an index (mode="TE0", "TM1"). The mode window defaults to the core plus a pad sized so the mode's intensity is 30 dB below its peak at the window's edge, with a 20 % margin, at the longest wavelength the port reads; the fitted domain leaves room for it. window_um=(half_w_um, half_v_um) sets it explicitly. ph.mode_window_um(width, thickness, n_core, n_clad, wlen_max) gives the rule's window from the guide alone, before the port's own solve refines it. Which side of the plane faces the domain wall is inferred from where the port sits relative to the device; pass out_direction="+" or "-" when the device is symmetric about the port. Declarative ports resolve to explicit solver sources and monitors at construction. The manual pipeline below exposes those pieces. For reflection, read the driven port backward: sim.port_monitors["in"].mode_power(data, direction="-"). That sim is the one that ran. For a result opened again later, read ports through the result: data.transmission("out"), or data.simulation.port_monitors. A simulation rebuilt from your script can have another grid, and its port_monitors are not checked against the saved fields.

The pieces underneath¶

You can also build sources and monitors with photonhub.analysis. Solve a grid-consistent mode on a geometry-only simulation, create its source plane, and project each monitor's fields onto the mode. This is useful when a declarative port cannot express the excitation you need.

from photonhub.analysis import (
    mode_launch, mode_monitor, solve_yee_mode, solve_yee_mode_bank,
    transmission_spectrum,
)

# Shared scene: a 0.5 x 0.22 um silicon strip along z in oxide cladding.
si = ph.Medium(permittivity=3.48**2)
base = dict(
    size_um=(3.0, 3.0, 4.0),
    grid=ph.UniformMesh(dl_um=0.1),
    run=ph.RunSpec(n_steps=200),  # workflow check only
    pml_num_layers=6,
    background=ph.Background(permittivity=1.44**2),
    structures=[ph.Structure(
        geometry=ph.Box(center_um=(1.5, 1.5, 2.0), size_um=(0.5, 0.22, 4.0)),
        medium=si)],
)

pulse = ph.GaussianPulse(freq0_hz=1.934e14, fwidth_hz=4e13)
freqs = (0.95 * 1.934e14, 1.934e14, 1.05 * 1.934e14)

# 1. Solve TE0 on the simulation grid
geometry_sim = ph.Simulation(**base)
window = dict(h_center_um=1.5, v_center_um=1.5,
              half_w_um=0.8, half_v_um=0.6, dl_um=0.1)
mode = solve_yee_mode(geometry_sim, "z", 0.8, 1.55, "TE", 0, **window)
modes = solve_yee_mode_bank(geometry_sim, "z", 0.8, freqs, "TE", 0, **window)

# 2. mode_launch returns the dipoles that form a source plane
sources = mode_launch(geometry_sim, mode, axis="z", position_um=0.8,
                      source_time=pulse, modes_by_freq=modes)
mm_in = mode_monitor(geometry_sim, mode, axis="z", position_um=1.2,
                     freqs_hz=freqs, name="in", modes_by_freq=modes)
mm_out = mode_monitor(geometry_sim, mode, axis="z", position_um=3.2,
                      freqs_hz=freqs, name="out", modes_by_freq=modes)

# 3. assemble with the launch sources + each monitor's field monitor, then run
sim  = ph.Simulation(**base, sources=sources,
                     monitors=[mm_in.field_monitor, mm_out.field_monitor])
data = ph.run_local(sim, output_dir="mode-run", device="cpu")

t = transmission_spectrum(mm_out, mm_in, data)   # T over f, with a wlen_um coordinate

The older ModeSolver → mode_source path remains available for scalar-mode compatibility, including the launches in the inverse-design examples. New full-vector simulations should use solve_yee_mode / mode_launch.

Gaussian beams (fibre / free-space excitation)¶

For an edge coupler, a grating coupler, or anything else fed from off-chip, the excitation is a Gaussian beam, not a guided mode. Declare it as the simulation's source and the source plane is built at construction, at the wavelength wlens_um names:

sim = ph.Simulation(**{**base, "size_um": (4.0, 6.0, 6.0)}, wlens_um=1.55,
                    source=ph.GaussianBeam("x", 1.0, mfd_um=2.0, polarization="Ez", n=1.45))

gaussian_beam_source is the function underneath; it returns the same kind of list of dipoles that form a source plane mode_launch does:

from photonhub.analysis import gaussian_beam_source

# A small spot keeps this construction example small. The source plane must
# fit in the interior, outside the absorbing boundary layers.
beam_geometry_sim = ph.Simulation(**{**base, "size_um": (4.0, 6.0, 6.0)})

sources = gaussian_beam_source(
    beam_geometry_sim, axis="x", position_um=1.0, source_time=pulse,
    mfd_um=2.0,               # 1/e^2 intensity diameter in microns
    polarization="Ez",         # the E component on the x-normal plane
    n=1.45,                    # the medium the beam lives in (oxide cladding)
    power_watts=1.0)           # launched beam power -> T reads an absolute fraction

Beyond the at-waist spot, the beam carries its own phase, so it also expresses what a real optical setup does:

waist_um=(wh, wv) / mfd_um=(...): an elliptical (lensed-fibre) spot.
waist_distance_um: the waist need not be on the injection plane. Negative focuses the beam inside the domain; positive launches an already-diverging beam.
angle_theta_rad / angle_phi_rad: off-normal injection, with both angles in radians.
pol_angle_rad: linear polarization angle in radians in the plane.
freqs_hz: a broadband launch (one source plane per frequency), worth it when the profile actually moves across the band (offset waist, tilt, dispersive n).

gaussian_beam(...) alone returns the beam as a VectorMode, which is what you want as the reference of a mode_monitor for a chip-to-fibre coupling readout, or as an operand of mode_overlap. Symmetry planes are honoured automatically.

Older spellings angle_start, angle_stop, angle_theta, angle_phi and pol_angle are compatibility aliases. Use the _rad names for new code.

transmission_spectrum(mm_out, mm_in, data) ratios modal powers at the two planes and cancels the common source normalization. For reflection, mm_in.mode_power(data, direction="-") reads backward modal power. Refine the mesh and extend the run before using either as a physics result.

Visualize¶

Simulation carries sim.plot (the scene), sim.plot_index (the meshed permittivity on a cut plane; it draws Box/Sphere/Cylinder/Polygon geometry), and sim.plot_3d (an interactive view, via the optional photonhub[viz] extra). RunResult.plot_field renders a recorded field slice.

ax = dipole_sim.plot(z=1.5)             # scene: structures + sources + monitors + PML
dipole_sim.plot_index(z=1.5, grid=True)   # meshed index with the grid overlay

For the optional 3D view, install the visualization extra first:

python -m pip install "photonhub[viz]"

On Windows, in Command Prompt or PowerShell:

.venv\Scripts\python.exe -m pip install "photonhub[viz]"

fig = dipole_sim.plot_3d()

Two module-level helpers cover the mode pipeline:

from photonhub.viz import plot_mode, plot_spectrum

plot_mode(mode)                 # FDE mode-profile heatmap (Ex/Ey, n_eff, λ)
plot_spectrum(t)                # transmission T vs wavelength (nm)
# a {label: {freq_hz: T}} mapping draws several labelled traces, e.g. for a
# coupler's two output ports:  plot_spectrum({"through": t_through, "cross": t_cross})

All 2D plotters return a matplotlib Axes, accept ax=, and never call plt.show().

Where to go next¶

The Examples publish selected studies with their settings.
The analysis guide covers mode solving, S-parameters, resonances, far field and EME.
Troubleshooting explains recovery and support. Known limitations lists beta restrictions and engine-version skew.

For beta help, contact beta@leapfield.ai.