Inverse design¶
Discrete adjoint topology and shape optimisation. A DesignRegion holds the optimisable permittivity, a WidthProfile renders a few widths into it, a PolygonProfile renders a polygon's vertices as exact area fractions, and an objective such as ModeTransmission defines what to maximise.
Design region¶
dataclass
¶
DesignRegion(
dl_um: float,
cells: Tuple[
Tuple[int, int], Tuple[int, int], Tuple[int, int]
],
shape: Tuple[int, int, int],
eps_min: float = 1.0,
eps_max: float = 12.25,
parametrization: str = "boxes",
box_face_offset_cells: float = 0.25,
merge_boxes: bool = False,
)
A rectangular topology-optimization region tiled into shape pixels.
The region occupies integer Yee-cell ranges [i0,i1) x [j0,j1) x [k0,k1)
on a uniform grid of pitch dl_um; shape = (npx, npy, npz) pixels must
divide each range evenly. Each pixel is a density rho in [0, 1] mapped to
a relative permittivity eps = eps_min + rho*(eps_max - eps_min). The
controls become quarter-cell-offset pixel boxes by default, so no Yee E
sample lies on any pixel face. box_face_offset_cells=0 retains integer
faces and uses exact engine-style face ownership. parametrization="trilinear" selects one
trilinear custom medium instead and is CPU only. eps_min >= 1.
On a one-cell periodic axis (a quasi-2D scene) the quarter-cell faces of
cells=(0, 1) hold the half-cell sample, which the E component along the
axis reads, and miss the sample at 0, which the two across it read. A
gradient refuses such a region when a source or the objective names a
component that misses, and warns when it cannot tell, as with a mode
source. cells=(-1, 1) holds both samples in every parametrization.
box_face_offset_cells=0 holds them too, but it moves the faces on
every axis, which changes the design.
merge_boxes=True emits the same pixel model with far fewer structures:
one box at eps(0) over the whole region, then one box per run of
pixels with bitwise equal density along the region's longest axis
(adjacent equal runs merged), and none for a density-0 pixel. With
quarter-cell faces no Yee sample lies on any face, so the engine samples
exactly the permittivity of the one-box-per-pixel layout at every point,
and the gradient is unchanged. The engine scans every structure at every
sample, so its setup time grows with the number of structures: on a
quasi-2D CPU test with 259200 pixel boxes (720 x 360, 0.025 um cells) the
one-box-per-pixel layout took 483 s to set up. A shape rendering that
is 0 or 1 over most of the region needs only a few thousand merged
boxes. Requires parametrization="boxes" and
box_face_offset_cells=0.25.
Build with on_grid from physical microns; the raw constructor takes
cell indices.
property
¶
Number of density parameters: the product of the pixel-grid dimensions.
property
¶
Number of solver cells along each axis of one design pixel.
property
¶
Realized physical extent along each axis, in microns.
property
¶
Center of the realized design region, in microns.
classmethod
¶
on_grid(
*,
center_um: Tuple[float, float, float],
size_um: Tuple[float, float, float],
dl_um: float,
shape: Tuple[int, int, int],
eps_min: float = 1.0,
eps_max: float = 12.25,
parametrization: str = "boxes",
box_face_offset_cells: float = 0.25,
merge_boxes: bool = False,
) -> "DesignRegion"
Snap a physical region (centre/size in microns) to integer cells.
The low corner snaps to round((center-size/2)/dl) and the span to
round(size/dl) cells, then rounded UP so each axis divides into its
pixel count. size_um and center_um describe the snapped
integer-cell extent; box faces shift by box_face_offset_cells.
structures(rho: ndarray) -> Tuple[Structure, ...]
Material structures for this region's parametrization.
trilinear_structure(rho: ndarray) -> Tuple[Structure, ...]
One CPU-only trilinear custom-medium box for the density grid.
The control at node (i,j,k) is density
rho[min(i,nx-1),min(j,ny-1),min(k,nz-1)]. The box faces are a
quarter cell past integer planes. No Yee E node lies on a face.
box_structures(rho: ndarray) -> Tuple[Structure, ...]
Pixel boxes with quarter-cell faces unless integer faces are selected.
With merge_boxes the merged layout of merged_box_structures.
merged_box_structures(
rho: ndarray,
) -> Tuple[Structure, ...]
The merge_boxes layout: a box at eps(0) over the region, then
one box per run of bitwise equal densities along the longest axis,
equal runs on adjacent lines of the next longest axis merged, and
nothing for a pixel whose permittivity equals eps(0).
Every pixel's permittivity is the same float as in the one-box-per-pixel layout, and with quarter-cell faces every Yee sample lies inside exactly one of the merged boxes or only in the first, so the engine samples the two layouts identically.
monitor(
freq_hz,
name: str = "design_region",
*,
simulation: Optional[Simulation] = None,
) -> ProfileMonitor
The DFT monitor recording all three E-components over the region (quarter-cell faces so the components co-snap).
freq_hz is one frequency, or a sequence or array of frequencies
for a BroadbandModeTransmission objective, which records the
region at every one of its freqs_hz in the same run.
Pass simulation (any simulation with the run's size, grid and
boundaries) for a quasi-2D run: on an axis that the simulation
realizes as one plain periodic cell the monitor becomes the plane at
a quarter cell, which the engine snaps to that cell for every
component, instead of a box reaching outside the one-cell domain.
Other axes are unchanged. value_and_gradient passes the
forward simulation for the adjoint run.
dataclass
¶
WidthProfile(
region: DesignRegion,
axis: str,
width_axis: str,
positions_um: Tuple[float, ...],
guide_center_um: float,
edge_um: Optional[float] = None,
)
A guide whose width varies along one axis, rendered with a soft edge.
density and jacobian take widths_um, where
widths_um[j] is the full width at positions_um[j] along axis.
Between controls the width is interpolated linearly; before the first
and after the last control it is held. The guide is centered at
guide_center_um along width_axis and fills the region along the
third axis. A pixel whose box center lies at s along axis and
t along width_axis takes density S(u) with
u = (w(s)/2 - |t - guide_center_um|)/edge_um + 1/2 and S the
quintic step 6u^5 - 15u^4 + 10u^3 on [0, 1]. The edge is
edge_um wide and centered on the nominal boundary; outside it the
density is exactly 0 or 1. u is linear in the widths and S has
continuous first and second derivatives, so the density is twice
continuously differentiable in the widths (its third derivative jumps
where u crosses 0 or 1), and jacobian returns its exact
derivative for parametric_value_and_gradient and
optimize_parametric.
The edge sets how evenly a width acts on the density. Only pixel centers
inside the edge band move, so d(sum rho)/dw follows S' of those
centers. With an edge of one pixel pitch along width_axis exactly
one center lies in the band: as the half-width moves across one pitch
that sensitivity swings from 0 to 1.92 times its mean, and it is exactly
zero whenever the nominal edge lies on a pixel face (the edge pixels sit
at u = 0 and 1). The width is then dead, and an optimizer meets
plateaus and strong curvature between them. On a 3D CPU test (0.044 um
uniform mesh, one cell per pixel, subpixel off, non-dispersive Si in
SiO2, PML 6 on every axis, 3000 fixed steps, shutoff 0) a one-cell edge
made dT/dt of a uniform guide thickness swing from -0.099 per um to
exactly 0 over one cell of half-thickness; a two-cell edge kept it
between -0.021 and -0.062. An edge of two pitches holds the sensitivity
within 0.94 to 1.05 times its mean, and three within 0.99 to 1.01. An
edge below the pitch can hold no pixel center at all: on a region of two
cells per pixel a one-cell edge left the density exactly 0 or 1, and the
gradient exactly zero, over 49% of a two-cell sweep of the half-width.
The default edge is the larger of two cells and the pitch along
width_axis, which is two pitches on one cell per pixel. On coarser
pixels it is one pitch, so pass edge_um of two pitches for an
optimization there. A smaller edge_um than the pitch warns.
A width whose edge reaches the region's first or last pixel along
width_axis is clipped there, and its gradient is lost; density
warns when either boundary pixel has nonzero density.
Build the region with on_grid for one cell per pixel along
axis and width_axis, or pass any pixel-box region. Use
density as expand and jacobian as jacobian.
Symmetric or composite shapes are built from these two methods with the
chain rule.
classmethod
¶
on_grid(
*,
center_um: Tuple[float, float, float],
size_um: Tuple[float, float, float],
dl_um: float,
axis: str,
width_axis: str,
positions_um: Tuple[float, ...],
guide_center_um: float,
eps_min: float = 1.0,
eps_max: float = 12.25,
edge_um: Optional[float] = None,
) -> "WidthProfile"
Profile on a region of one cell per pixel along axis and
width_axis and one pixel across the third axis.
The region snaps as DesignRegion.on_grid does.
Density grid of shape region.shape for the given widths.
Warns when the first or last pixel along width_axis has nonzero
density: the guide's edge is then cut by the region boundary.
Exact d density / d widths, shape (region.n_params, num_widths)
in flat density(...).ravel() order.
dataclass
¶
PolygonProfile(
region: "DesignRegion",
vertices_um: Tuple[Tuple[float, float], ...],
plane: Tuple[str, str] = ("x", "y"),
movable: Optional[Tuple[int, ...]] = None,
move: Optional[str] = None,
mirror: Optional[Tuple[str, float]] = None,
)
A polygon rendered as exact area fractions on a pixel-box region.
The polygon lies in the plane of plane = (a, b), two of x, y and z,
and is extruded through the region along the third axis. Each pixel of
region takes as its density the exact fraction of its box (face
offset included) that the polygon covers, the continuous material
averaging of Michaels and Yablonovitch, Opt. Express 26, 31717 (2018),
Eq. 14, for two materials. Moving a vertex changes the density only in
the pixels its two edges cross, continuously everywhere (the render has
no staircase jumps) and with a continuous derivative, except for a move
across a pixel face of an edge lying on that face (see below); the second
derivative jumps where the vertex crosses a pixel face or one of its
edges sweeps over a pixel corner, so a central difference across such a
point converges at first order only.
vertices_um lists the base vertices as (a, b) pairs in microns.
Without mirror they are the whole polygon, closed from the last
vertex back to the first. With mirror = (axis, value), axis one
of the plane's axes, they are an open chain and the polygon is the
chain followed by its mirror image about axis = value in reverse
order, so the shape is symmetric by construction; a chain vertex and
its image move together, mirrored. The parameters are vertex
displacements in microns, zero at the base polygon: for each vertex
index in movable (default every vertex, or every chain vertex), one
displacement along each axis in move (default both plane axes, in
plane order), so p.reshape(len(movable), len(move)) with the
default order. Vertices not listed stay fixed, for example the ends of
a taper where it meets fixed guides.
density renders, jacobian returns the exact
d density / d p as a scipy.sparse matrix of shape
(region.n_params, num_params) in flat density(...).ravel() order;
pass the two as expand and jacobian to
parametric_value_and_gradient and optimize_parametric,
which chain the sparse matrix without densifying it.
difference_jacobian is the forward difference of the exact
render at a tiny step, for checking. The area of a pixel is not
differentiable when an edge lies exactly on a pixel face and moves
across it; jacobian then returns the mean of the two one-sided
derivatives. A base vertex within 1e-9 pixel of a face is placed on
it.
The region must use pixel boxes with quarter-cell faces (the
DesignRegion default); build it with on_grid, which
merges equal pixels into a few boxes (merge_boxes) so a large region
stays cheap to set up. A polygon that leaves the region is cut at its
boundary; structures outside the region are the caller's. The polygon
must be simple: planar_density, density and
jacobian test every pair of non-adjacent edges and raise,
naming the first pair, when two touch or cross (for 400 vertices a
small fraction of the render's cost); pass check_simple=False to
skip the test when the polygon is known to stay simple. They also raise
when a pixel's covered fraction leaves [0, 1]. With mirror the
test covers the chain and its image, so a chain crossing the axis
raises too.
classmethod
¶
on_grid(
*,
center_um: Tuple[float, float, float],
size_um: Tuple[float, float, float],
dl_um: float,
vertices_um,
plane: Tuple[str, str] = ("x", "y"),
movable: Optional[Tuple[int, ...]] = None,
move: Optional[str] = None,
mirror: Optional[Tuple[str, float]] = None,
eps_min: float = 1.0,
eps_max: float = 12.25,
) -> "PolygonProfile"
Profile on a region of one cell per pixel in the plane and one pixel across the third axis, with equal pixels merged into boxes.
The region snaps as DesignRegion.on_grid does. For a
quasi-2D run (one periodic cell along the third axis) center the
region on 0 along that axis with size 2 * dl_um: its
quarter-cell faces at -0.75 and 1.25 cells then hold every Yee
sample of the single plane, and region.monitor(freq_hz,
simulation=sim) collapses its monitor onto that plane.
The same profile with count fixed vertices added on each end
edge of the chain, spacing_um apart, so the curvature penalty
sees the corners where the chain meets its fixed ends.
For a mirrored chain whose first and last vertices are fixed, for
example far inside the guides a taper joins: the new vertices lie on
the edge from vertex 1 back toward vertex 0, at k * spacing_um
from vertex 1 for k = 1 .. count, and likewise from the
second-to-last vertex toward the last. The polygon at p = 0 and
the order of p are unchanged. A move of vertex 1 or of the
second-to-last vertex now tilts only the short edge to the nearest
new vertex instead of the whole end edge, so the render away from
p = 0 and those vertices' Jacobian columns change; every other
column is the same. With spacing_um about the movable vertices'
own spacing, the radius at each junction no longer averages over a
neighbour far away (see curvature_radii).
The full polygon at p, shape (num_vertices, 2) in plane order.
The vertices listed in vertices_um at p (the chain, without
its mirror image), shape (len(vertices_um), 2).
Covered fraction of each pixel in the plane, shape (n_a, n_b).
Raises if two non-adjacent edges touch or cross (skipped with
check_simple=False) or a pixel's fraction leaves [0, 1].
Density grid of shape region.shape: the exact covered fraction
of each pixel box, the same along the third axis. Raises as
planar_density does.
Exact d density / d p, a scipy.sparse CSR matrix of shape
(region.n_params, num_params) in flat density(...).ravel()
order, per micron of displacement. Raises if two non-adjacent edges
touch or cross, unless check_simple=False.
Forward difference of the exact render, (density(p + h e_j) -
density(p)) / h, as a sparse matrix like jacobian.
step_um defaults to 1e-8 of a cell. Entries below 1e-6 of a
pixel of area per pixel of displacement are rounding of the column
sums and are dropped. For checking jacobian; it renders
the polygon num_params + 1 times and does not test that it is
simple. A step that reverses the polygon's orientation, as one
larger than the thickness of a sliver can, gives a wrong column: on
a sliver 2^-30 pixel thick the default step is off by 160%, where
jacobian stays exact.
curvature_penalty of the vertex chain at p and its
gradient with respect to p.
With mirror the chain is vertices_um (the image has the same
radii and is not counted again); without it, the closed polygon,
every vertex with its two neighbours. Fixed vertices count as
neighbours and as points. Use it as penalty of
optimize_parametric::
penalty=lambda p: profile.curvature_penalty(
p, r0_um=0.15, alpha=..., kappa_per_um=...)
curvature_penalty(
points_um,
*,
r0_um: float,
alpha: float,
kappa_per_um: float,
) -> Tuple[float, ndarray]
Logistic penalty on small radii of curvature along a vertex chain.
P = alpha * sum_i 1 / (1 + exp(-kappa (r0 - R_i))) over the
interior points of the chain, with R_i from
curvature_radii: a radius well below r0_um costs about
alpha, one well above costs nearly nothing, and kappa_per_um
sets how sharp the step is (1/kappa is its width in um). This is
the penalty of Michaels and Yablonovitch, Opt. Express 26, 31717 (2018),
Appendix B, with their fit read as parameter values -1, 0 and 1 at the
three points. Returns (P, dP/dpoints), the gradient exact with
shape points_um.shape; a straight run of three points has an
infinite radius and contributes neither value nor gradient.
Each R_i assumes the point's two neighbours are roughly equally far
from it (see curvature_radii): on a circle, neighbours at -0.1
and +0.5 rad read 1.77 times the true radius. Where a movable chain
meets a fixed vertex far away, as a taper meets a fixed guide, the
corner there is therefore barely penalized; PolygonProfile.pad_ends
adds fixed vertices along the guides at the chain's own spacing.
Objectives¶
dataclass
¶
TransmissionPort(
mode: Union[Mode, VectorMode],
axis: str,
position_um: float,
out_direction: str,
name: str,
center_um: Optional[Tuple[float, float]] = None,
thickness_axis: Optional[str] = None,
)
One modal port plane for ModeTransmission.
out_direction points away from the design. The input port's incident
wave travels in the opposite direction. The plane is snapped by
photonhub.analysis.mode_monitor, and its four-field monitor is
cropped clear of transverse PML, as for ModePower: on each
transverse PML axis the aperture starts 2.25 cells after the low PML slab
and ends 2.75 cells before the high one. Both faces must sit at the same
quarter-cell phase so that all four tangential components snap to the same
cells, so the two margins differ by half a cell and the aperture center
lies a quarter cell below the domain center. The overlap therefore differs
from a full-plane port when the mode's tail reaches the crop; the
objective warns when more than 1e-3 of the reference mode's sampled power
lies outside the aperture.
sport(simulation: Simulation, freq_hz: float) -> SPort
Build the S-matrix port and its PML-clear four-field monitor.
dataclass
¶
ModeTransmission(
input_port: TransmissionPort,
output_port: TransmissionPort,
freq_hz: float,
incident_amplitude: Optional[complex] = None,
)
Single-frequency modal transmission between two finite-aperture ports.
By default T = |b_out|**2 / |a_in|**2 reads the incident wave
a_in and the outgoing wave b_out from the same forward run. This
is the S-matrix element |S_out,in|**2 of this objective's own ports,
which are cropped clear of transverse PML (see TransmissionPort);
a full-plane port on the same run reads a slightly different value. The
adjoint differentiates both waves exactly. Include monitors in
build_forward.
How well T measures power transmission depends on the reference
modes. The incident projection also picks up part of the design's
reflection through any mismatch between the reference mode and the
solver's own discrete mode, and that leak changes with the plane position
rather than fading with distance from the design or from the source. Two
reflective 3D CPU tests measured it: the adjoint validation record (0.05
um mesh, eight input planes 0.75 to 7.75 cells before the design) and an
independent review scene (0.04 um mesh, nine planes 0.75 to 15.75 cells
before the design), both on a uniform Yee mesh with subpixel off,
non-dispersive Si in SiO2, six PML layers on every axis, 3000 fixed
steps and shutoff 0. Across the planes the incident power differed from
a design-free run by up to 5.7% with scalar modes in both tests, 0.42%
and 1.0% with Yee-staggered modes solved on the smoothed cross-section
(solve_mode_on_cross_section with its default supersample=8),
and 0.16% and 0.24% with Yee-staggered modes solved on the staircase
permittivity of the subpixel-off run (supersample=1); the same-run
T spread 6.4% and 12.4%, 1.7% and 1.4%, and 0.11% and 0.16%. The
adjoint runs subpixel off, so solve port modes with supersample=1.
Constructing the objective warns when a port's mode is not
Yee-staggered, or records a smoothed solve (supersample other than
1 in its solve_params).
The geometry guard requires the input readout rows to lie strictly
between the source and the design and the output rows to lie beyond the
design. With staircase Yee modes nothing more was needed in either test:
the design-free guide read within 0.11% of 1 at every plane from 1.25
cells past the source. Scalar and smoothed readings carry a
plane-dependent error, several percent with scalar modes and up to about
1% with smoothed ones, that distance from the source does not remove: in
the review scene the scalar design-free guide read 0.985 to 1.036 over
planes 1.5 to 16.5 cells past the source, beating without settling. The
remedy is staircase Yee modes (supersample=1), not a larger gap.
Pass incident_amplitude to hold the denominator fixed:
T = |b_out|**2 / |a_ref|**2, where a_ref is
reference_incident_amplitude of a design-free run with the same
source, mesh and input port. The adjoint then keeps only the output-port
term. Prefer it for strongly reflective designs, so that an optimizer
cannot raise T by steering the reflection into the incident reading:
at input planes 0.75 cells and 6.75 or 7.75 cells before the design, the
same-run incident term was 6.8-22% of the gradient norm with scalar
modes, 0.7-3.8% with smoothed and 0.37-1.6% with staircase Yee modes
over the two tests. The fixed amplitude removes that reflection term from
the gradient, not the plane-dependent error of a scalar or smoothed mode,
which also enters a_ref. It assumes the source launches the same wave
with and without the design, which a mode source does up to the
reflection that the boundary behind the source returns: 0.17% of the
incident power with six PML layers in the review scene and 0.003% with
twelve. The same-run a includes that return and a_ref cannot, so
fixed-incident T carries it (0.27% higher with twelve layers there),
and it is part of the staircase residual above.
class-attribute
instance-attribute
¶
ports(simulation: Simulation) -> Tuple[SPort, SPort]
Return the exact finite-aperture SPorts used by this objective.
The ports of the two most recent simulations are cached by identity, so repeated readouts of one run do not rebuild the monitors or re-solve the automatic per-frequency set of modes. The crop warning fires once per port name and aperture (the cropped monitor's center and size, and the grid step), so a later simulation with a smaller cross-section or more PML layers is checked again.
monitors(
simulation: Simulation,
) -> Tuple[ProfileMonitor, ProfileMonitor]
The input and output four-field monitors for build_forward.
waves(
data: RunResult, simulation: Optional[Simulation] = None
) -> Tuple[complex, complex]
Return same-run incident input a and outgoing output b.
reference_incident_amplitude(
data: RunResult, simulation: Optional[Simulation] = None
) -> complex
Read a_ref from a design-free run with the same input port.
Run a normalization simulation with the same source, mesh, run
window, input monitor and input guide, without the design structures,
and pass the returned complex amplitude as incident_amplitude to a
new objective with the same ports.
amplitude(
data: RunResult, simulation: Optional[Simulation] = None
) -> complex
Complex b/a, using the fixed reference when supplied.
value(
data: RunResult, simulation: Optional[Simulation] = None
) -> float
Power ratio of the output wave to the chosen incident reading.
transpose_sources(simulation: Simulation, data: RunResult)
One point-current cloud for the exact differential of transmission.
dataclass
¶
BroadbandModeTransmission(
input_port: TransmissionPort,
output_port: TransmissionPort,
freqs_hz: Tuple[float, ...],
incident_amplitudes: Optional[
Tuple[complex, ...]
] = None,
reduction: Union[
str, Callable[[ndarray], Tuple[float, ndarray]]
] = "mean",
weights: Optional[Tuple[float, ...]] = None,
softmin_beta: float = 20.0,
)
Modal transmission at several frequencies, reduced to one figure of merit.
T_k = |b_k|**2 / |a_k|**2 is read at every frequency of freqs_hz
from ONE forward run: b_k is the outgoing wave of the output port and
a_k the fixed incident_amplitudes[k], or the same-run incident wave
at the input port when incident_amplitudes is None. The ports are the
finite-aperture TransmissionPort planes of
ModeTransmission, recording all frequencies. A port mode solved
with solve_mode_on_cross_section carries its solve settings, and
each frequency is then read against the mode re-solved at that
frequency; the re-solved modes are kept for later simulations of the
same objective. Read every a_k from one design-free broadband run
with reference_incident_amplitudes.
The figure of merit is F = reduction(T):
"mean":sum_k w_k T_k / sum_k w_k, withw=weights(default all 1);"sum":sum_k w_k T_k;"softmin":-(1/beta) log(sum_k w_k exp(-beta T_k) / sum_k w_k)withbeta=softmin_beta, a smooth minimum that lies between the smallestT_kand the weighted mean and approaches the minimum asbetagrows;- a callable taking the array
Tand returning(F, dF/dT).
value_and_gradient returns F and dF/drho from one forward
and one adjoint solve whatever the number of frequencies. Its adjoint run
drives every port sample with several delayed copies of the forward
pulse whose amplitudes are solved so the summed current has, at each
f_k, exactly the single-frequency adjoint current scaled by
dF/dT_k; the adjoint window is longer than the forward one by the
span of those delays, about 1/df for the smallest frequency spacing
df. Include monitors and region.monitor(freqs_hz) in
build_forward, and give the forward simulation one pulse that covers
every frequency (each within 3.7 fwidth_hz of its freq0_hz, where
its spectrum is at least 1e-3 of the peak).
Launch each frequency's own mode when T_k should match one run per
wavelength: mode_launch(..., modes_by_freq=...) with a mode solved at
every frequency. A launch of one profile injects the centre frequency's
mode at every frequency, while each T_k is read against the port
mode of its own frequency. On a quasi-2D CPU taper (0.025 um uniform
mesh, subpixel off, non-dispersive, 12 PML layers on x and y, 8000 fixed
steps, shutoff 0, 21 wavelengths from 1.35 to 1.85 um, one pulse centred
at 1.55 um) the one-run T_k differed from one run per wavelength by
up to 1.3% at the ends of the band with the centre mode launched, and by
at most 5.2e-4 with each frequency's own mode.
class-attribute
instance-attribute
¶
class-attribute
instance-attribute
¶
ports(simulation: Simulation) -> Tuple[SPort, SPort]
The two finite-aperture ports, recording every frequency.
The ports of the two most recent simulations are cached by identity, and each port's per-frequency modes by port name, grid and domain size, so readouts of later simulations on the same mesh re-solve no mode.
monitors(
simulation: Simulation,
) -> Tuple[ProfileMonitor, ProfileMonitor]
The input and output four-field monitors for build_forward.
waves(
data: RunResult, simulation: Optional[Simulation] = None
) -> Tuple[ndarray, ndarray]
Same-run incident waves a and outgoing waves b, one per frequency.
reference_incident_amplitudes(
data: RunResult, simulation: Optional[Simulation] = None
) -> Tuple[complex, ...]
Read every a_ref(f_k) from one design-free broadband run.
Run a normalization simulation with the same source, mesh, run
window, input monitor and input guide, without the design structures,
and pass the returned amplitudes as incident_amplitudes to a new
objective with the same ports and frequencies.
amplitudes(
data: RunResult, simulation: Optional[Simulation] = None
) -> ndarray
Complex b_k/a_k per frequency, using the fixed references when supplied.
transmissions(
data: RunResult, simulation: Optional[Simulation] = None
) -> ndarray
T_k per frequency, in freqs_hz order.
(F, dF/dT) of the reduction at the given T_k.
value(
data: RunResult, simulation: Optional[Simulation] = None
) -> float
The figure of merit F.
dataclass
¶
ModePower(
mode: Union[Mode, VectorMode],
axis: str,
position_um: float,
freq_hz: float,
direction: str = "+",
name: str = "mode_out",
center_um: Optional[Tuple[float, float]] = None,
thickness_axis: Optional[str] = None,
beta: complex = BETA_MODE,
)
Maximize the power coupled into a guided mode at an output port.
J = |c|^2, where c is the P_mode-normalized complex modal amplitude
of the recorded plane (mode_amplitude). For port-to-port
transmission, use ModeTransmission with suitable port modes.
Warning
J is a RELATIVE objective, not a calibrated power
transmission. |c|^2 == T would require the launch to read
c_in == 1 on an input plane, but mode_source(power_watts=1)
normalizes to SI watts, so J carries a scene/grid-dependent
positive scale (measured ~4/dl_um on SOI strip scenes: J = 23.4
at dl = 0.05 um where T ~ 0.29). Maximizing J still maximizes T;
optimization and relative comparisons are unaffected, but do NOT
report J as transmission; use the S-matrix path
(photonhub.analysis.smatrix), which normalizes by the driven
port's incident amplitude. The accuracy of that reading depends on
the reference mode and the placement of its input plane.
The adjoint uses electric and magnetic point currents built from the transpose of the actual four-field overlap, including transverse Yee co-location for scalar modes. The recording monitor excludes transverse PML rows.
monitor(simulation: Simulation) -> ProfileMonitor
Four-field output plane with transverse PML rows excluded.
The finite-aperture overlap differs from the former full-plane objective; it keeps the adjoint current cloud in the physical interior where the reciprocal Yee operator applies.
amplitude(data: RunResult) -> complex
The complex normalized modal amplitude c on the output plane.
NOTE: this reads WITHOUT the longitudinal Yee de-stagger that
ModeMonitor.mode_power applies by default (mode_amplitude's
destagger_dl is left off). The objective and its adjoint
coefficient conj(c) use the SAME reading, so the gradient stays
self-consistent; only the absolute |c|² differs from mode_power's
default reading by the O(beta*dl/2) stagger term.
value(data: RunResult) -> float
Figure of merit |c|^2, relative modal power (not a calibrated
transmission; see the class warning).
adjoint_source(forward_sim: Simulation)
Legacy backward mode launch for direct callers.
Gradient evaluation uses transpose_sources instead.
adjoint_coeff(data: RunResult) -> complex
The complex excitation coefficient conj(c).
transpose_sources(simulation: Simulation, data: RunResult)
Point-current cloud equal to the transpose of the modal readout.
Returns (sources, coefficient). The first electric dipole has
positive unit amplitude, fixing the engine's DFT normalization.
dataclass
¶
PointIntensity(
probe_um: Tuple[float, float, float],
freq_hz: float,
component: str = "Ez",
name: str = "probe",
beta: complex = BETA,
)
Maximize |E_comp(probe)|^2 at a single point and frequency, the
simplest adjoint objective: its adjoint source is a single point dipole at
the probe, polarized along component, with post-multiplied coefficient
conj(u) (u = the forward phasor at the probe).
A focusing/concentrator design (push energy to a focal point) is exactly this objective.
monitor() -> ProfileMonitor
Build the single-point DFT monitor for this objective component and frequency.
amplitude(data: RunResult) -> complex
The complex forward objective phasor u = E_comp(probe).
adjoint_source(forward_sim: Simulation) -> PointDipole
Unit point dipole at the probe, reusing the forward pulse so the normalization matches (the objective coefficient is applied in post-processing, not in the source amplitude).
Gradients¶
value_and_gradient(
build_forward: BuildForward,
region: DesignRegion,
objective: Objective,
rho: ndarray,
*,
device: Optional[str] = None,
run_kwargs: Optional[dict] = None,
beta: Optional[complex] = None,
monitor_name: str = "design_region",
_warn_beta: bool = True,
) -> GradientResult
Return the objective and physical dJ/drho from two solves.
The returned value and gradient use the caller's simulation with subpixel
smoothing off, with a warning when the caller's simulation enables it.
The adjoint run injects the transpose of the objective readout and records
the E field over the design. beta is a compatibility override; the
default is derived from the leapfrog frequency and vacuum permittivity.
Supplying beta is deprecated. Both solves use the caller's run settings,
including the default stop rule. The gradient is exact for the recorded
window; its error relative to a converged gradient follows the truncation
of the forward results. With the energy-only stop (shutoff 1e-5) that
was about 0.2% on a measured point-source scene, and more on resonant
designs. On a local CPU solver the default stop also waits for an estimate
of the remaining change in the recorded frequency-domain results
(RunSpec.dft_shutoff) when the solver supports it; the estimate is not
a guarantee, and GPU runs and older solvers stop on the field energy
alone. For resonant designs and validation, set shutoff=0 with an
explicit n_steps or run_time_s after checking convergence. A
simulation fitted with domain= is refused; give size_um and
grid.
Quarter-cell pixel boxes support CPU and GPU; the optional trilinear
custom medium supports CPU only.
With a BroadbandModeTransmission the forward run reads every
T_k and the design fields at every frequency, and ONE adjoint run
carries the adjoint currents of all frequencies at once (see that class);
the result is a BroadbandGradientResult whose value is F
and grad is dF/drho. The adjoint run is longer than the forward
run by the span of its delayed pulses, adjoint_extra_steps in the
result; the forward run's other settings, its stop rule included, carry
over. beta is refused for it. The broadband path is validated on the
CPU.
build_adjoint_simulation(
forward_simulation: Simulation,
forward: RunResult,
region: DesignRegion,
objective: Objective,
*,
monitor_name: str = "design_region",
) -> Tuple[Simulation, complex]
Return the runnable adjoint simulation and its gradient coefficient.
Run the returned simulation with run_local, then pass both results
and the returned complete coefficient to assemble_gradient
with its default beta. The coefficient already includes the full
discrete physical factor.
The forward simulation must use subpixel off.
value_and_gradient uses these same pieces automatically.
assemble_gradient(
region: DesignRegion,
forward: RunResult,
adjoint: RunResult,
coeff: complex,
freq_hz: float,
*,
monitor_name: str = "design_region",
beta: complex = BETA,
) -> ndarray
Contract the field product with exact trilinear control weights.
beta is the derived physical coefficient supplied by the driver.
The default is retained solely for old direct assembly callers.
finite_difference_check(
build_forward: BuildForward,
region: DesignRegion,
objective: Objective,
rho: ndarray,
h: float,
*,
indices: Optional[Tuple[int, ...]] = None,
device: Optional[str] = None,
run_kwargs: Optional[dict] = None,
) -> FiniteDifferenceResult
Compare the discrete adjoint with central differences of real solves.
Each checked control takes two additional forward simulations, at
rho[i] + h and rho[i] - h. The objective value, every plus and
minus value and every solve length are recorded. Keep the density
perturbations inside the physical [0, 1] range.
dataclass
¶
GradientResult(
value: float,
grad: ndarray,
forward: RunResult,
adjoint: RunResult,
amplitude: complex,
forward_steps_run: Optional[int],
adjoint_steps_run: Optional[int],
shutoff: float,
forward_simulation: Optional[Simulation] = None,
adjoint_simulation: Optional[Simulation] = None,
)
Objective and gradient for the subpixel-off model.
shutoff is shared by the forward and adjoint solves. Their actual
lengths are forward_steps_run and adjoint_steps_run; they can differ
when the caller's shutoff ends either solve early. forward_simulation
and adjoint_simulation are the exact specifications sent to the solver.
class-attribute
instance-attribute
¶
forward_simulation: Optional[Simulation] = None
class-attribute
instance-attribute
¶
adjoint_simulation: Optional[Simulation] = None
dataclass
¶
BroadbandGradientResult(
value: float,
grad: ndarray,
forward: RunResult,
adjoint: RunResult,
amplitude: complex,
forward_steps_run: Optional[int],
adjoint_steps_run: Optional[int],
shutoff: float,
forward_simulation: Optional[Simulation] = None,
adjoint_simulation: Optional[Simulation] = None,
freqs_hz: Tuple[float, ...] = (),
transmissions: Optional[ndarray] = None,
amplitudes: Optional[ndarray] = None,
dF_dT: Optional[ndarray] = None,
gradients_by_freq: Optional[ndarray] = None,
taps: int = 0,
tap_span_s: float = 0.0,
tap_condition: float = float("nan"),
adjoint_extra_steps: int = 0,
adjoint_dipoles: int = 0,
)
Bases: GradientResult
Objective and gradient of a BroadbandModeTransmission.
value is F and grad is dF/drho, from one forward and one
adjoint solve. amplitude holds b_k/a_k per frequency, as does
amplitudes. transmissions are the T_k and dF_dT the
reduction's slopes. gradients_by_freq[k] is frequency k's part of
grad, already weighted by dF_dT[k]; the parts sum to grad.
The adjoint run drives each port sample with taps delayed pulses
spread over tap_span_s seconds, adjoint_dipoles dipoles in all,
and runs adjoint_extra_steps steps longer than the forward run.
tap_condition is the condition number of the tap solve.
dataclass
¶
FiniteDifferenceResult(
indices: Tuple[int, ...],
h: float,
adjoint_gradient: ndarray,
finite_difference_gradient: ndarray,
cosine: float,
norm_ratio: float,
relative_l2: float,
forward_steps_run: Optional[int],
adjoint_steps_run: Optional[int],
finite_difference_steps_run: Tuple[
Tuple[Optional[int], Optional[int]], ...
],
value: Optional[float] = None,
finite_difference_values: Tuple[
Tuple[float, float], ...
] = (),
jacobian: Optional[ndarray] = None,
)
Adjoint and central finite-difference gradients from real solver runs.
indices selects flat controls: densities for
finite_difference_check, shape parameters for
parametric_finite_difference_check. value is the objective at
the checked point. finite_difference_values and
finite_difference_steps_run hold the plus and minus objective values
and run lengths in index order. Checking all controls costs two extra
forward solves per control, in addition to the adjoint pair. For a
parameter check, jacobian is the drho/dp the adjoint gradient was
chained through.
Shape-parameter gradients¶
parametric_value_and_gradient(
build_forward: BuildForward,
region: DesignRegion,
objective: Objective,
p: ndarray,
expand: Callable[[ndarray], ndarray],
*,
jacobian: Optional[Callable[[ndarray], ndarray]] = None,
fd_step: float = 0.001,
device: Optional[str] = None,
run_kwargs: Optional[dict] = None,
monitor_name: str = "design_region",
check_jacobian: bool = True,
) -> ParametricGradientResult
Return the objective and dJ/dp for shape parameters p.
expand(p) renders the density grid (region.shape, values in
[0, 1]). One forward and one adjoint solve give dJ/drho as in
value_and_gradient, and the chain rule gives
dJ/dp = (drho/dp)^T dJ/drho without further solves. The cost is two
solves whatever the number of parameters.
jacobian(p) supplies drho/dp exactly, as
WidthProfile.jacobian does, with shape
(region.n_params, len(p)) in flat density order or
region.shape + (len(p),). It may be a scipy.sparse matrix of
the first shape, as PolygonProfile.jacobian returns; it is then
chained and returned sparse, never densified. Without it, column j of drho/dp is
(expand(p + fd_step e_j) - expand(p - fd_step e_j)) / (2 fd_step),
with fd_step in the parameter's own unit. Two cases:
- A smooth
expand: the difference is accurate to second order infd_step, and the error scales with(fd_step / edge)^2for an edge of widthedge. On the two 3D shape-validation scenes, rendered withWidthProfileon a 0.044 um mesh with its default two-cell edge, the default 1e-3 um step moveddJ/dpby 1.5e-5 and 2.3e-4 relative to the exact Jacobian, within the resolution of central differences of solves there (1.2e-4 and 6.4e-4). - A piecewise-linear
expand(np.clip,np.interpin the parameters,np.maximumof two renders): away from a kink the difference is exact up to rounding, but for every density whose kink lies withinfd_stepofpit returns a weighted mean of the two one-sided slopes (their plain mean at the kink), which is the slope on neither side. At the kink the objective has no gradient, and an optimizer stepping either way sees a different slope from the one reported. Prefer a smooth rendering, keeppaway from kinks, or supplyjacobian.
Two checks render expand again before any solve, at p and along
every parameter at multiples of fd_step/8 out to 1.25 fd_step,
and warn without changing the result. Without jacobian, a kink
within fd_step of p makes the forward and backward differences
disagree by more than the O(fd_step) a smooth rendering allows; it
shows as a slope change concentrated at one sample, which a twice
continuously differentiable expand cannot produce. The warning names
the parameters. With jacobian, the supplied matrix is compared with
the central difference at fd_step/8, away from such kinks and
allowing for its truncation error, which catches a wrong unit, sign or
flattening order before the 2N + 2 solves of a finite-difference
check. Both warn above 1% of the largest one-sided slope of expand.
Neither warns on WidthProfile renderings, including widths
whose edge pixels sit at the ends of their soft edge. The checks render
one parameter at a time, so their memory stays within a few times the
Jacobian's; for a sparse jacobian the comparison also runs one
column at a time and keeps about 30 renders of expand whatever the
number of parameters (about 67 MB and 21 s for the 400 vertex
displacements of a 736 x 384 pixel PolygonProfile on a laptop
CPU). check_jacobian=False skips them, and a render that fails
beyond p +/- fd_step skips them too.
Check the result against real solves with
parametric_finite_difference_check. Run settings, the subpixel
model, the stop rule and the refusals are those of
value_and_gradient.
parametric_finite_difference_check(
build_forward: BuildForward,
region: DesignRegion,
objective: Objective,
p: ndarray,
expand: Callable[[ndarray], ndarray],
h: float,
*,
jacobian: Optional[Callable[[ndarray], ndarray]] = None,
fd_step: float = 0.001,
indices: Optional[Tuple[int, ...]] = None,
device: Optional[str] = None,
run_kwargs: Optional[dict] = None,
check_jacobian: bool = True,
) -> FiniteDifferenceResult
Compare the chained parameter gradient with central differences of real solves.
Each checked parameter takes two forward simulations, at p[i] + h
and p[i] - h, each re-rendered with expand; h is in the
parameter's own unit (microns for widths). The adjoint side is
parametric_value_and_gradient with the same jacobian,
fd_step and check_jacobian. The objective value, every plus and minus value, every
solve length and the Jacobian used are recorded. Keep every rendered
density in [0, 1] and keep p +/- h clear of kinks in expand,
where central differences of the objective are not a derivative.
dataclass
¶
ParametricGradientResult(
value: float,
grad: ndarray,
params: ndarray,
rho: ndarray,
jacobian: ndarray,
jacobian_method: str,
density_gradient: GradientResult,
)
Objective and its gradient with respect to shape parameters.
grad is dJ/dp = jacobian.T @ density_gradient.grad. jacobian
is drho/dp with shape (region.n_params, len(params)) in flat
density order, a scipy.sparse matrix when the caller's
jacobian(p) returned one (as PolygonProfile.jacobian does),
and jacobian_method says whether the caller supplied
it or it is a central difference of expand. rho is
expand(params). density_gradient is the underlying
GradientResult, with the forward and adjoint simulations and
their run lengths.
Optimisation¶
optimize(
build_forward: BuildForward,
region: DesignRegion,
objective: Objective,
rho0: ndarray,
*,
n_iters: int = 30,
method: str = "lbfgs",
step: float = 0.05,
bounds: Tuple[float, float] = (0.0, 1.0),
maximize: bool = True,
device: Optional[str] = None,
run_kwargs: Optional[dict] = None,
beta: Optional[complex] = None,
param_map: Optional[
Callable[[ndarray], ndarray]
] = None,
callback: Optional[
Callable[[int, GradientResult, ndarray], None]
] = None,
) -> OptimizeResult
Topology optimization driven by the adjoint gradient.
A BroadbandModeTransmission objective optimizes its figure of
merit over all of its frequencies at the same cost per evaluation, one
forward and one adjoint solve; best is then a
BroadbandGradientResult.
method is "lbfgs" (default, SciPy L-BFGS-B) or "adam"
(normalized-gradient Adam); see _descend. L-BFGS-B is the default
here as in optimize_parametric: in an optimizer comparison with the
exact gradient (one 3D CPU run per optimizer on one mesh) it reached the
higher transmission on a 64-pixel topology crossing and led Adam at equal
evaluations from the 6th evaluation on. Adam was the default before that
comparison, for its scale-free update; pass method="adam" to keep it.
Every L-BFGS-B evaluation is one SciPy asks for, the first at rho0;
each iteration takes one or more (its line search), so n_iters
iterations can take more than n_iters evaluations. For a FEW smooth
shape parameters, prefer optimize_parametric. n_iters bounds
the optimizer iterations. Densities are box-constrained to bounds.
Maximizes objective by default; returns the BEST design seen. step
is used only by Adam.
device ("cpu" / "gpu" / "gpu:N") runs every iteration's two
solves on that backend (forwarded to run_local).
param_map (a (region.shape) -> (region.shape) array map) constrains the
design: the structure is built from param_map(rho) and the gradient is
mapped through it, a proper projected gradient. For a LINEAR, self-adjoint
projection (e.g. a symmetry-averaging map, which is its own adjoint) this is
exact; best and rho are reported in the mapped (constrained) space.
Use it to enforce device symmetry or a density filter. Returned values and
gradients use the subpixel-off model. beta= is deprecated. Both solves
retain the caller's run settings and are exact for their recorded window;
the error against a converged gradient follows the truncation of the
forward results (see value_and_gradient). For resonant designs and
validation, set shutoff=0 with an explicit n_steps or
run_time_s after checking convergence.
optimize_parametric(
build_forward: BuildForward,
region: DesignRegion,
objective: Objective,
p0: ndarray,
expand: Callable[[ndarray], ndarray],
*,
n_iters: int = 30,
method: str = "lbfgs",
step: float = 0.02,
bounds: Optional[Tuple[float, float]] = None,
maximize: bool = True,
device: Optional[str] = None,
run_kwargs: Optional[dict] = None,
beta: Optional[complex] = None,
fd_step: float = 0.001,
jacobian: Optional[Callable[[ndarray], ndarray]] = None,
callback: Optional[
Callable[[int, GradientResult, ndarray], None]
] = None,
check_jacobian: bool = True,
penalty: Optional[
Callable[[ndarray], Tuple[float, ndarray]]
] = None,
) -> ParametricResult
Adjoint shape optimization of a few geometric parameters p
instead of a free density per pixel.
expand(p) renders the density grid (region.shape, values in
[0, 1]), for example a guide's widths rendered by
WidthProfile.density. Each evaluation is
parametric_value_and_gradient: one forward and one adjoint solve
give dJ/drho, and the chain rule through drho/dp gives
dJ/dp without further solves, so the cost per evaluation does not
grow with the number of pixels or parameters, nor with the number of
frequencies of a BroadbandModeTransmission. jacobian(p) supplies
drho/dp exactly (WidthProfile.jacobian); without it a
central difference of expand with step fd_step is used, which is
wrong for densities within fd_step of a kink of expand. Both
paths run the no-solve checks of parametric_value_and_gradient
on the first evaluation, at p0, and warn on a kink or on a supplied
Jacobian that disagrees with expand; check_jacobian=False skips
them. Check the gradient at a
representative p with parametric_finite_difference_check.
method is "lbfgs" (default, SciPy L-BFGS-B) or "adam"; see
_descend. Every L-BFGS-B evaluation is one SciPy asks for, the
first at p0; each iteration takes one or more (its line search).
bounds (lo, hi) box-constrains every parameter. Returns the best
parameters seen. Maximizes objective by default. Values and
gradients use the subpixel-off model. beta= is deprecated. Both
solves keep the caller's run settings and are exact for their recorded
window; the error against a converged gradient follows the truncation
of the forward results (see value_and_gradient). For resonant
designs and validation, set shutoff=0 with an explicit n_steps
or run_time_s after checking convergence.
penalty(p) returns a value and its gradient with respect to p,
for example PolygonProfile.curvature_penalty; it needs no
solve. The optimizer then works on the merit J - penalty when
maximizing and J + penalty when minimizing, so the penalty always
counts against the design. history keeps J and
penalty_history the penalty, each per evaluation, and the best
evaluation is the one with the best merit.
dataclass
¶
OptimizeResult(
rho: ndarray,
history: List[float],
grads: List[ndarray],
best: GradientResult,
)
Optimization history and best value for the subpixel-off model.
history, grads, and best all refer to that same model even if
the caller's constructed simulation enabled subpixel smoothing.
dataclass
¶
ParametricResult(
params: ndarray,
rho: ndarray,
history: List[float],
params_history: List[ndarray],
best: GradientResult,
grads: List[ndarray] = list(),
penalty_history: List[float] = list(),
)
Shape-optimization history and the best parameters seen.
history, params_history and grads hold every evaluation in
order, the first at the starting point; each evaluation is two solves.
With L-BFGS-B they include its line-search trials, so the last
evaluation need not be the best. grads are the
chained dJ/dp of each evaluation, before any scaling inside the
optimizer.
best is the density-level GradientResult of the best
evaluation. With a penalty, history and grads stay the
objective's own, penalty_history holds the penalty of each
evaluation, and the best evaluation is the one with the best penalized
merit (objective minus penalty when maximizing, plus when minimizing).