pyoptics.fdtd

Finite-difference time-domain (FDTD) solver for Maxwell’s equations in 1D and 2D.

The FDTD method (Yee 1966) samples E and H on staggered grids in space and time and leapfrogs them through the curl equations. It is rigorous for any structure that can be drawn on the grid: waveguides, resonators, gratings, scatterers of any size, metals. A single run with a short pulse gives the response over a whole band.

Building blocks

• Grid: domain size, cell size and time step (normalized units: c = 1, lengths and wavelengths in the same unit, times in that unit over c).

• Shapes (Slab, Rect, Circle, Annulus, Polygon, unions and differences) and Medium (permittivity, conductivity), drawn on the grid with sub-cell averaging of the permittivity.

• Boundaries: PML (convolutional PML), “pec” walls, “periodic”.

• Sources: PointSource, LineSource (soft currents), TFSF (plane wave in a total-field / scattered-field box), with the time profiles GaussianPulse and ContinuousWave.

• Monitors: Probe, Snapshots, DFTField, Flux.

• Simulation puts them together; pyoptics.fdtd.plot draws structures, fields and animations.

Polarization: Ez and Hz

2D simulations are invariant along z and split into two independent polarizations, named here after the field component along z:

• mode="Ez": E along z (out of the plane), fields Ez, Hx, Hy,

• mode="Hz": H along z, fields Ex, Ey, Hz.

The literature names them inconsistently. In the FDTD literature (Taflove and Hagness) Ez is “TM” (H transverse to z) and Hz “TE”. For waveguides, gratings and scattering, the names refer to the plane of incidence instead: E perpendicular to it (here: E along z) is “TE” or s, so there Ez is “TE” and Hz “TM”. pyoptics.litho.rcwa and pyoptics.analytical.cylinder_field use the latter: their TE case is the Ez mode here.

Example

>>> from pyoptics import fdtd
>>> grid = fdtd.Grid((8.0, 8.0), 1/30)
>>> sim = fdtd.Simulation(grid, structure=[(fdtd.Circle((0, 0), 1.0), fdtd.Medium.from_index(1.5))],
...                       sources=[fdtd.TFSF((-2, 2, -2, 2), fdtd.GaussianPulse(1.0))])
>>> dft = fdtd.DFTField(1.0, components=["Ez"])
>>> sim.monitors.append(dft)
>>> sim.run(time=40)

This is the first milestone of the design in proposals/fdtd.md: NumPy kernels and non-dispersive media. Planned: numba kernels, dispersive media, mode sources and S-parameters, Bloch boundaries, 3D.

class pyoptics.fdtd.Grid(size, spacing, center=None)

Bases: object

Regular Yee grid in 1D or 2D.

Parameters:

sizenumber or tuple

Extent of the computational domain along each axis, including the absorbing layers. A number gives a 1D grid.

spacingnumber or tuple

Cell size, the same for all axes if a number. The number of cells is round(size/spacing).

centernumber or tuple, optional

Centre of the domain, 0 by default.

axis(a, offset=0.0)

Coordinates of the nodes along axis a, shifted by offset cells.

coords(component, mode='Ez')

1D coordinate arrays (one per axis) of a field component.

index(position, component, mode='Ez')

Index tuple of the node of a component closest to a position.

time_step(courant=0.99)

Time step as a fraction courant of the stability limit 1/sqrt(sum(1/spacing**2)).

property extent

(xmin, xmax[, ymin, ymax]) of the domain, for imshow().

class pyoptics.fdtd.Shape

Bases: object

Base class: a region of space, defined by contains().

class pyoptics.fdtd.Union(*shapes)

Bases: Shape

Points inside any of the shapes.

class pyoptics.fdtd.Difference(a, b)

Bases: Shape

Points inside a but not inside b.

class pyoptics.fdtd.Everywhere

Bases: Shape

The whole domain.

class pyoptics.fdtd.Slab(lo, hi, axis=0)

Bases: Shape

Layer lo <= coordinate < hi along an axis (0: x, 1: y), unbounded in the other directions. Works in 1D and 2D.

class pyoptics.fdtd.Rect(center, size)

Bases: Shape

Axis-parallel rectangle (2D) given by its centre and size.

class pyoptics.fdtd.Circle(center, radius)

Bases: Shape

Disc (2D), the cross section of a cylinder along z.

class pyoptics.fdtd.Annulus(center, r_inner, r_outer)

Bases: Shape

Ring between two radii (2D), e.g. a ring resonator.

class pyoptics.fdtd.Polygon(vertices)

Bases: Shape

Simple polygon (2D, convex or concave) given by its vertices.

class pyoptics.fdtd.Medium(eps=1.0, sigma=0.0)

Bases: object

Non-dispersive, possibly conductive medium.

Parameters:

epsnumber

Relative permittivity (>= 1 for stability with the default time step).

sigmanumber, optional

Conductivity in normalized units (absorption).

classmethod from_index(n, k=0.0, wavelength=None)

Medium with the complex index n + ik at one vacuum wavelength (exact there; the conductivity model makes k vary with wavelength).

complex_eps(wavelength)

Complex relative permittivity at a vacuum wavelength.

index(wavelength=None)

Complex refractive index n + ik (at a wavelength if conductive).

class pyoptics.fdtd.PML(cells=12, order=3, reflection=1e-08, kappa_max=1.0, alpha_max=0.0)

Bases: object

Absorbing layer of the given number of cells on both ends of an axis.

Parameters:

cellsint

Thickness in cells.

ordernumber

Polynomial grading order m of sigma and kappa.

reflectionnumber

Theoretical reflection R0 at normal incidence; sets sigma_max = -(m + 1)*ln(R0)/(2*thickness).

kappa_maxnumber

Maximum real coordinate stretch (>= 1); > 1 helps with evanescent and grazing waves.

alpha_maxnumber

Complex frequency shift at the inner edge (in units of 1/time, decreasing linearly to 0 at the wall); absorbs low frequencies less and so avoids late-time drift for long runs.

coefficients(n_cells, spacing, dt, offset)

b, c and 1/kappa for the left and the right layer along an axis of n_cells cells, at nodes shifted by offset cells. Each is an array of length cells; the left layer covers indices 0..cells - 1, the right one n_cells - cells..n_cells - 1.

class pyoptics.fdtd.GaussianPulse(wavelength, bandwidth=0.2, delay=None, amplitude=1.0)

Bases: object

Sine carrier under a Gaussian envelope (no DC component).

s(t) = amplitude*exp(-(t - delay)**2/(2 tau**2))*sin(omega_0 (t - delay))

Parameters:

wavelengthnumber

Vacuum wavelength of the carrier.

bandwidthnumber

Spectral standard deviation relative to the carrier frequency; tau = 1/(2 pi bandwidth f_0).

delaynumber, optional

Centre of the pulse, 5 tau by default (the pulse starts smoothly).

property end

Time after which the pulse is negligible (< 4e-6 of its peak).

class pyoptics.fdtd.ContinuousWave(wavelength, ramp=None, amplitude=1.0)

Bases: object

sin(omega t), switched on smoothly over a ramp time (raised cosine).

Parameters:

wavelengthnumber

rampnumber, optional

Rise time, 3 periods by default.

class pyoptics.fdtd.PointSource(position, profile, component=None, amplitude=1.0)

Bases: _CurrentSource

Point dipole (a line current along z in 2D, a sheet in 1D).

Parameters:

positionnumber or tuple

profilecallable

Time profile, e.g. GaussianPulse or ContinuousWave.

componentstr, optional

Field component the current drives: Ez (default in 1D and the Ez mode), Ex or Ey (electric), Hz (magnetic current, default in the Hz mode).

amplitudenumber

class pyoptics.fdtd.LineSource(position, profile, axis=0, component=None, extent=None, transverse=None, amplitude=1.0)

Bases: _CurrentSource

Current sheet across the domain at a position along an axis (2D), or a sheet in 1D. With periodic boundaries along the line it launches plane waves to both sides.

Parameters:

positionnumber

Coordinate of the line along axis.

profilecallable

Time profile.

axisint

Normal of the line (0: a line x = position, 1: y = position).

componentstr, optional

As for PointSource.

extent(lo, hi), optional

Part of the line that radiates (default: all of it).

transversecallable, optional

Amplitude across the line as a function of the coordinate along it, e.g. a Gaussian beam profile.

amplitudenumber

class pyoptics.fdtd.TFSF(box, profile, n_background=1.0)

Bases: Source

Plane wave travelling towards +x, injected on the boundary of a box (total-field / scattered-field formulation, 2D, both modes).

Inside the box the total field is computed, outside only the scattered field: without scatterers the field outside the box is zero up to round-off. The incident wave comes from an auxiliary 1D grid with the same cell size and time step, whose numerical dispersion along x is identical to that of the 2D grid.

Parameters:

box(xmin, xmax, ymin, ymax)

Total field region; it should lie in the background medium near its edges and stay clear of the PML.

profilecallable

Time profile of the incident field (its E at the source of the auxiliary grid; the wave reaches x = xmin a bit later).

n_backgroundnumber

Refractive index of the background medium.

incident(x, t=None)

Incident Ez (Ez mode) or Ey (Hz mode) of the auxiliary grid at x, at the current time of the simulation (nearest node).

class pyoptics.fdtd.Probe(position, component=None)

Bases: Monitor

Time series of one field component at the node nearest to a point.

Attributes after a run: times, values (arrays).

spectrum(wavelengths)

DFT of the recorded series at vacuum wavelengths (convention as DFTField).

class pyoptics.fdtd.Snapshots(component=None, every=10, region=None, stride=1, dtype=<class 'numpy.float16'>, keep_last=None, start=0.0, path=None)

Bases: Monitor

Decimated frames of one field component.

Parameters:

componentstr, optional

Default: Ez in the Ez mode, Hz in the Hz mode.

everyint

Record every every steps.

regiontuple, optional

Crop, ((xmin, xmax), (ymin, ymax)).

strideint

Spatial subsampling.

dtypenumpy dtype

Storage type of the frames; float16 halves the memory of float32 and is plenty for plots.

keep_lastint, optional

Keep only the last frames (ring buffer).

startnumber

Time at which recording starts.

pathstr, optional

Stream the frames to this HDF5 file (needs h5py) instead of keeping them in memory.

close()

Close the HDF5 file of a streaming recorder.

stack()

The recorded frames as a SnapshotStack (read back from the file when streaming).

class pyoptics.fdtd.SnapshotStack(times, coords, frames, component, mode='Ez')

Bases: object

Frames of one component: times (n,), coords (one array per axis), frames (n, …) and the component name.

save(path)

Save to .npz or (with h5py) .h5.

class pyoptics.fdtd.DFTField(wavelengths, components=None, region=None, start=0.0, stride=1)

Bases: Monitor

Running DFT of field components at vacuum wavelengths.

Parameters:

wavelengthsnumber or array

componentssequence of str, optional

Default: all components of the mode.

regiontuple, optional

((xmin, xmax), (ymin, ymax)); default: the whole grid.

startnumber

Time from which on to accumulate (e.g. after a CW source has reached the steady state).

strideint

Spatial subsampling.

field(component)

Complex amplitudes (n_wavelengths, …) of a component.

data()

The result as DFTFieldData (e.g. to save it).

steady_state(component)

Complex amplitudes normalized to a continuous wave: for a CW run accumulated over whole periods, this is the phasor of the field (F*2/duration).

class pyoptics.fdtd.DFTFieldData(wavelengths, fields, coords, mode='Ez')

Bases: object

Complex field amplitudes at some wavelengths: fields[component] has the shape (n_wavelengths, …), coords[component] one array per axis.

class pyoptics.fdtd.Flux(position, wavelengths, axis=0, extent=None)

Bases: Monitor

Spectral power flux through a line (2D) or a point (1D).

The line is normal to axis at position, optionally limited to extent along the other axis; the flux counts positive along +axis. H is interpolated to the E nodes. flux() returns 1/2 Re(sum E H*) dA per wavelength (arbitrary units: divide by the flux of a reference run).

flux()

Power flux per wavelength (n_wavelengths,).

fields()

DFT amplitudes of E and H (interpolated) along the line.

class pyoptics.fdtd.Simulation(grid, structure=(), sources=(), monitors=(), mode='Ez', boundaries='pml', background=Medium(eps=1, sigma=0), courant=0.99, dt=None, subpixel=None, dtype=<class 'numpy.float64'>)

Bases: object

FDTD simulation in 1D or 2D.

Parameters:

gridGrid

structuresequence of (Shape, Medium)

Drawn in order on top of the background, later ones on top.

sourcessequence of Source

monitorssequence of Monitor

mode“Ez” or “Hz”

2D polarization: Ez has the electric field along z (out of the plane, fields Ez, Hx, Hy), Hz the magnetic field along z (Ex, Ey, Hz). Ignored in 1D, which has Ez and Hy. See the module docstring of pyoptics.fdtd for the TE/TM names.

boundariesstr, PML or sequence of them

Per axis “pml” (or a PML instance), “pec” or “periodic”; a single entry applies to all axes. Default: PML() everywhere.

backgroundMedium

courantnumber

Time step as a fraction of the stability limit.

dtnumber, optional

Explicit time step (overrides courant).

subpixelint, optional

Sub-samples per cell and axis for averaging the permittivity across material boundaries (1: staircase). Default 64 in 1D, 8 in 2D. The average is arithmetic, which is the right one for E parallel to the boundaries (1D, and the Ez mode in 2D).

dtypenumpy dtype

float64 (default) or float32.

property time

Time of the electric field (H is half a step behind).

rasterize(component)

Permittivity and conductivity at the nodes of a component, averaged over subpixel**ndim samples per cell.

setup()

Rasterize the structure, allocate the fields and set up sources and monitors (done by run() if needed).

step()

One time step: H from n - 1/2 to n + 1/2, E from n to n + 1.

run(time=None, steps=None, callback=None, every=None, progress=False)

Advance the simulation by a time span or a number of steps.

Parameters:

timenumber, optional

Simulated time span (rounded up to whole steps).

stepsint, optional

callbackcallable, optional

callback(sim), called every every steps (e.g. live plots).

everyint, optional

progressbool

Print the number of steps and the run time at the end.

energy()

Electromagnetic energy 1/2 sum(eps E**2 + H**2) dV.

E and H are half a step apart, so this oscillates slightly around the true energy. With track_energy = True (set before stepping), the E term is eps E(n)*E(n + 1) instead, which makes it the quantity that the Yee scheme conserves exactly in lossless media inside PEC walls.

permittivity(component=None)

Rasterized permittivity at the nodes of a component (default: the first electric component).

pml_cells()

Number of PML cells per axis (0 where there is none).