# Antenna Inverse Design#

In this notebook we perform inverse design of an antenna using a medium ranging between air and an experimentally measured gold from tidy3dโs material library.

We use topology optimization to design the patterning on this sheet of gold to maximizes the intensity enhancement in a central location, while respecting fabrication penalties.

This uses the native automatic differentiation support in tidy3d >= 2.7.0 through autograd.

## Set up#

We start by importing the packages we need, including autograd and the autograd wrapper for numpy.

[1]:

import matplotlib.pylab as plt

import tidy3d as td
import tidy3d.web as web


Weโll set a few helper constants for the notebook.

[3]:

# whether to run checks for setting up sim and normalization enhancement
run_pre_sims = True

um = 1e0
nm = 1e-3


Then we define the parameters setting our source spectrum and FDTD spatial resolution.

[4]:

wavelength = 0.910

freq = td.C_0 / wavelength
fwidth = freq / 10
run_time = 20 / fwidth

dl = 3 * nm

# min_steps_per_wvl = 30


Next we define our material parameters.

We import the gold model from our material library, which is a td.PoleResidue medium from a fit to measured data.

[5]:

# select one of the material library variants
print(tuple(td.material_library["Au"].variants.keys()))

('Olmon2012crystal', 'Olmon2012stripped', 'Olmon2012evaporated', 'Olmon2012Drude', 'JohnsonChristy1972', 'RakicLorentzDrude1998')

[6]:

medium_SiO2 = td.Medium(permittivity=1.44**2)
medium_gold = td.material_library["Au"]["JohnsonChristy1972"]
eps_background = 1.0

[7]:

wvls = np.linspace(0.5, 1.0, 101)
freqs = td.C_0 / wvls
eps_complex = medium_gold.eps_model(freqs)
eps_real, sigma = medium_gold.eps_complex_to_eps_sigma(eps_complex, freq=freq)
plt.plot(wvls, eps_real, label="permittivity (real)")
plt.plot(wvls, np.imag(eps_complex), label="permittivity (imag)")
plt.plot(wvls, sigma, label="conductivity")
plt.plot(wvls, np.ones_like(eps_real), label="1", linestyle="--")
plt.plot(wvls, np.zeros_like(eps_real), label="0", linestyle="--")
plt.xlabel("wavelength (um)")
plt.ylabel("relative permittivity")
plt.legend()
plt.show()

[8]:

density_values = np.linspace(0, 1, 101)
eps_complex = 0j * np.zeros_like(density_values)
for i, d in enumerate(density_values):
new_eps_inf = 1.0 + d * (medium_gold.eps_inf - 1)
new_poles = [(a, c * d) for (a, c) in medium_gold.poles]
new_medium = medium_gold.updated_copy(eps_inf=new_eps_inf, poles=new_poles)
eps_complex[i] = new_medium.eps_model(freq)

eps_real, sigma = medium_gold.eps_complex_to_eps_sigma(eps_complex, freq=freq)
plt.plot(density_values, eps_real, label="permittivity (real)")
plt.plot(wvls, eps_complex.imag, label="permittivity (imag)")
plt.plot(density_values, sigma, label="conductivity")
plt.plot(density_values, np.ones_like(eps_real), label="1", linestyle="--")
plt.plot(density_values, np.zeros_like(eps_real), label="0", linestyle="--")
plt.xlabel("density value")
plt.ylabel("relative permittivity")
plt.legend()
plt.show()


Next we define the geometric parmaeters describing our device.

We have a metal rectanglular slab sitting on top of a SiO2 substrate with air on top.

The rectangular slab has a hole in the center, where weโll be measuring the field intensity enhancement.

[9]:

size_design_x = 500 * nm
size_design_y = 500 * nm
thick_metal = 50 * nm
thick_sub = 0.8 * wavelength
thick_air = 0.8 * wavelength
buffer_xy = 0.0 * wavelength



We derive the total sizes needed for our simulation in all dimensions.

[10]:

size_x = size_design_x + 2 * buffer_xy
size_y = size_design_y + 2 * buffer_xy
size_z = thick_sub + thick_metal + thick_air


As well as set some variables to help us keep track of location within the simulation.

[11]:

top_z = size_z / 2.0

center_air_z = top_z - thick_air / 2.0
center_metal_z = top_z - thick_air - thick_metal / 2.0
center_sub_z = top_z - thick_air - thick_metal - thick_sub / 2.0


Next, we define our static components, such as the substrate, source, and monitors.

We will include the pattern gold film later on, but we set the geometry as a td.Box here so we can use it later.

All of these components will be combined into our base td.Simulation, which defines our setup without any gold film.

[12]:

design_region_geometry = td.Box(
center=(0, 0, center_metal_z),
size=(td.inf, td.inf, thick_metal),
)

substrate = td.Structure(
geometry=td.Box(
center=(0, 0, -size_z + thick_sub),
size=(td.inf, td.inf, size_z),
),
medium=medium_SiO2,
)

plane_wave = td.PlaneWave(
center=(0, 0, center_air_z),
size=(td.inf, td.inf, 0),
source_time=td.GaussianPulse(freq0=freq, fwidth=fwidth),
pol_angle=0,
direction="-",
)

point_monitor = td.FieldMonitor(
size=(0, 0, 0),
center=(0, 0, center_metal_z),
freqs=[freq],
name="point",
)

field_monitor_xy = td.FieldMonitor(
center=(0, 0, center_metal_z),
size=(td.inf, td.inf, 0),
freqs=[freq],
name="field_xy",
)

eps_monitor_xy = td.PermittivityMonitor(
center=(0, 0, center_metal_z),
size=(td.inf, td.inf, 0),
freqs=[freq],
name="eps_xy",
)

field_monitor_xz = td.FieldMonitor(
center=(0, 0, 0),
size=(td.inf, 0, td.inf),
freqs=[freq],
name="field_xz",
)

sim_no_antenna = td.Simulation(
size=(size_x, size_y, size_z),
run_time=run_time,
structures=[substrate],
sources=[plane_wave],
monitors=[point_monitor, field_monitor_xy, field_monitor_xz, eps_monitor_xy],
grid_spec=td.GridSpec.uniform(dl=dl),
boundary_spec=td.BoundarySpec.pml(x=False, y=False, z=True),
# grid_spec=td.GridSpec.auto(wavelength=wavelength, min_steps_per_wvl=min_steps_per_wvl),
)

[13]:

# shift the plot positions a little so the field monitors don't overlap plots
shift_plot = 0.01

f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 4), tight_layout=True)
sim_no_antenna.plot(x=0.0 + shift_plot, ax=ax1)
sim_no_antenna.plot(y=0.0 + shift_plot, ax=ax2)
sim_no_antenna.plot(z=center_metal_z + shift_plot, ax=ax3)
plt.show()


## Antenna Parameterization#

Next, weโll define our antenna as a td.Structure, which we add to the simulation.

We will be optimizing this pattern with respect to an array of parameters. So it makes sense to write our components from here on as functions so we can update them as they change.

We start by defining some variables to define our design parameterization, such as the pixel size, feature size, and projection (binarization) strength.

[14]:

# resolution of the design region pixels
# pixel_size = 10 * nm
pixel_size = dl

# radius of the circular filter (um) (higher = larger features)

# projection strengths (higher = more binarized)
beta_project = 10
beta_penalty = 10


Next, we write a function to get the density (between 0 and 1) of gold on our film as a function of our optimization parameters.

Note: there are many helper functions in tidy3d.plugins.autograd that are used for inverse design. These have very close analogues in the adjoint plugin, but are compatible with autograd and have extra features added.

[15]:

from tidy3d.plugins.autograd.invdes import make_filter_and_project, get_kernel_size_px

def get_density(params: np.ndarray, beta: float = beta_project) -> np.ndarray:
"""Generic function to get the etch density as a function of design parameters, using filter and projection."""
filter_project = make_filter_and_project(filter_size, beta=beta)
return filter_project(params)


Next, we need to generate the td.Structure containing the gold film by modifying the material properties of the gold model from our material library.

We use a td.CustomPoleResidue to define a spatially-varying dispersive medium. The relative permittivity at infinite frequency (eps_inf) is set between that of air (1) and that of the gold eps_inf linearly, depending on the density value at each point. The numerator of the poles are also multiplied by the density so that they interpolate between values of 0 for density of 0 and the values for gold when density is 1.

We apply the mask over the central region to set it to air and combine everything together into a td.Structure.

[16]:

def make_antenna_from_density(density: np.ndarray) -> td.Structure:
"""Make a td.Structure containing a td.CustomPoleResidue corresponding to a density array."""

rmin, rmax = sim_no_antenna.bounds
coords = {}
for key, pt_min, pt_max, num_pts in zip("xyz", rmin, rmax, density.shape):
coord_edges = np.linspace(pt_min, pt_max, num_pts + 1)
coord_centers = (coord_edges[1:] + coord_edges[:-1]) / 2.0
coords[key] = coord_centers

np.ones_like(density),
coords=coords,
)

eps_inf_scaled_array = eps_background + density_masked * (medium_gold.eps_inf - eps_background)

poles_arrays_scaled = [
]

eps_inf_scaled = td.ScalarFieldDataArray(
eps_inf_scaled_array,
coords=coords,
)

poles_scaled = []
for a_array, c_array in poles_arrays_scaled:
a_scaled = td.ScalarFieldDataArray(a_array, coords=coords)
c_scaled = td.ScalarFieldDataArray(c_array, coords=coords)
poles_scaled.append((a_scaled, c_scaled))

medium_gold_scaled = td.CustomPoleResidue(
eps_inf=eps_inf_scaled, poles=poles_scaled, interp_method="linear"
)

return td.Structure(geometry=design_region_geometry, medium=medium_gold_scaled)


Next we write a quick helper function to create the td.Structure as a function of the optimization parameters.

[17]:

def make_antenna(params: np.ndarray, beta: float = beta_project) -> td.Structure:
"""Make the antenna structure."""
density = get_density(params, beta=beta)
return make_antenna_from_density(density)


We also write a function that generates a td.Simulation with the structure added to our original simulation, along with a mesh override structure to ensure even meshing in the design region.

[18]:

def make_sim_with_antenna(
params: np.ndarray, include_field_mnts: bool = True, beta: float = beta_project
) -> td.Simulation:
"""Make a simulation as a function of the density parameters for the antenna regions."""

antenna = make_antenna(params, beta=beta)

# add uniform mesh override structures to simulation (if desired)
design_region_mesh = td.MeshOverrideStructure(
geometry=antenna.geometry.updated_copy(size=(td.inf, td.inf, thick_metal)),
dl=[dl] * 3,
enforce=True,
)

sim = sim_no_antenna.updated_copy(
structures=list(sim_no_antenna.structures) + [antenna],
)

if not include_field_mnts:
sim = sim.updated_copy(monitors=[point_monitor])

return sim


Letโs test this out by generating a set of initial parameters for the optimization and calling our code.

[19]:

# some variables we might need later

nx = int(size_design_x // pixel_size)
ny = int(size_design_y // pixel_size)

# some intial parameters to test with
def make_symmetric_x(arr: np.ndarray) -> np.ndarray:
"""make an array symmetric in x."""
return (arr + np.flipud(arr)) / 2.0

params0 = make_symmetric_x(0.5 * np.ones((nx, ny, 1)))

[20]:

sim_antenna_params0 = make_sim_with_antenna(params0)

[21]:

f, axes = f, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 4), tight_layout=True)
sim_antenna_params0.plot(x=0.0 + shift_plot, monitor_alpha=0.0, ax=ax1)
sim_antenna_params0.plot(y=0.0 + shift_plot, monitor_alpha=0.0, ax=ax2)
sim_antenna_params0.plot(z=center_metal_z + shift_plot, monitor_alpha=0.0, ax=ax3)

for ax in axes:
ax.set_aspect("equal")

plt.show()


If we want, we can also visualize the initial fields.

[22]:

if run_pre_sims:
f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)
sim_data.plot_field("field_xy", field_name="E", val="abs^2", ax=ax1)
sim_data.plot_field("field_xz", field_name="E", val="abs^2", ax=ax2)
plt.show()

18:39:03 EDT Created task 'check fields' with task_id

             View task using web UI at
4-432e-94b9-2b691b85c131'.

18:39:05 EDT status = success

18:39:08 EDT loading simulation from simulation_data.hdf5


Next, we are ready to put everything together into an objective function to optimize.

We start by writing functions to compute the intensity at our measurement point through a tidy3d simulation.

[23]:

def extract_intensity(sim_data: td.SimulationData) -> float:
"""Grab the intensity from a simulation data."""
return np.sum(sim_data.get_intensity("point").values)

def measure_intensity(
params: np.ndarray,
verbose=False,
beta: float = beta_project,
) -> float:
"""Measure intensity as a function of design parameters."""
sim = make_sim_with_antenna(params, beta=beta, include_field_mnts=False)
return extract_intensity(sim_data)

if run_pre_sims:
intensity0 = measure_intensity(
)
else:
intensity0 = 2083.6917

print(f"Intensity without structure = {intensity0:.4f} (au)")

18:39:09 EDT Created task 'antenna_normalize' with task_id

             View task using web UI at
7-4b2a-848e-27e595572629'.

18:39:10 EDT status = success

18:39:11 EDT loading simulation from simulation_data.hdf5

Intensity without structure = 2083.6917 (au)


We then define a penalty function to discourage density patterns that have small feature sizes below the radius we defined.

[24]:

from tidy3d.plugins.autograd.invdes import make_erosion_dilation_penalty

penalty_fn = make_erosion_dilation_penalty(filter_size, beta=beta_penalty)

def penalty(params: np.ndarray, beta: float = None) -> float:
"""Define the erosion dilation invariance penalty for Si density parameters."""
beta_kwargs = {}
if beta is not None:
beta_kwargs["beta"] = beta
density = get_density(params, **beta_kwargs)
pen_val = penalty_fn(density)
return pen_val


## Objective Function#

Finally, we define our objective function as the ratio of our measured intensity compared to the intensity with no gold film minus the fabrication penalty, which is normalized between 0 and 1.

[25]:

# let's us grab and print autograd values while they're in the objective function

def intensity_enhancement(
params: np.ndarray, task_name: str = "antenna_intensity", verbose=False
) -> float:
return intensity_with_params / intensity0

def objective(
params, beta: float = beta_project, verbose: bool = False, penalty_only: bool = False
) -> float:
if penalty_only:
enhancement_factor = 0.0
else:
print(f"\tenhancement = {getval(enhancement_factor):.2e}")

# penalty_value = 0
penalty_value = penalty(params, beta=beta)
print(f"\tpenalty val = {getval(penalty_value):.2e}")
objective_value = enhancement_factor * (1.0 - penalty_value)
print(f"\tobjective = {getval(objective_value):.2e}")
return objective_value


We can use autograd to simply compute a function that returns the value of our objective along with the gradient, to feed to our optimizer.

[26]:

val_grad_fn = autograd.value_and_grad(objective)


We can take this opportunity to run through the function to verify everything, see our initial objective function value, and visualize the gradients of our objective w.r.t. the initial parameters.

[27]:

if run_pre_sims:
print(f"starting objective function value = {val}")
im1 = plt.imshow(np.flipud(np.squeeze(grad)).T, cmap="PiYG", vmax=vmag1, vmin=-vmag1)
plt.colorbar(im1)

             Created task 'antenna' with task_id

             View task using web UI at
4-44bd-afdd-bdb3ef6b35d1'.

18:39:12 EDT status = success

18:39:17 EDT loading simulation from simulation_data.hdf5

        enhancement = 1.20e+00
penalty val = 9.98e-01
objective = 2.18e-03

             Created task 'antenna_adjoint' with task_id

             View task using web UI at
9-48ce-a49e-3237e5d44e70'.

18:39:18 EDT status = success

18:39:23 EDT loading simulation from simulation_data.hdf5

starting objective function value = 0.0021820654597771005


## Optimization#

Next, we will use optax (an open source optimization package) to optimize this objective function with respect to our parameter array.

Note that the parameters are clipped between 0 and 1.

We will also visualize the density of the gold regions as we go.

[28]:

import optax

# hyperparameters
num_steps = 20
learning_rate = 0.01

beta_min = 1.0
beta_max = 55.0

def plot_density(density: np.ndarray, ax=None) -> None:
"""Plot the density of the device."""
if ax is None:
_, ax = plt.subplots(figsize=(2, 2))
arr = np.flipud(1 - density.squeeze().T)
plt.imshow(arr, cmap="gray", vmin=0, vmax=1, interpolation="none")
plt.axis("off")
plt.show()

# initialize adam optimizer with starting parameters (all combined)
params = params0
opt_state = optimizer.init(params)

# store history
objective_values_history = []
params_history = [params]
beta_history = []

# optimization loop
for i in range(num_steps):
print(f"step = {i + 1}")

beta = beta_min + i / num_steps * (beta_max - beta_min)
beta_history.append(beta)

# plot the densities, to monitor
plot_density(get_density(params, beta=beta))

# outputs
# print(f"\tobjective = {val:.4e}")
print(f"\tbeta = {beta:.4e}")

# compute and apply updates to the optimizer based on gradient (-1 sign to maximize obj_fn)
params = np.array(params)

# clip the parameters between 0 and 1
params = np.clip(params, 0.0, 1.0)

# save history
objective_values_history.append(val)
params_history.append(params)

step = 1

        enhancement = 1.20e+00
penalty val = 9.98e-01
objective = 2.18e-03
beta = 1.0000e+00
step = 2

        enhancement = 2.24e+00
penalty val = 9.98e-01
objective = 4.34e-03
beta = 3.7000e+00
step = 3

        enhancement = 4.15e+00
penalty val = 9.97e-01
objective = 1.33e-02
beta = 6.4000e+00
step = 4

        enhancement = 7.73e+00
penalty val = 9.90e-01
objective = 7.96e-02
beta = 9.1000e+00
step = 5

        enhancement = 1.22e+01
penalty val = 9.55e-01
objective = 5.56e-01
beta = 1.1800e+01
step = 6

        enhancement = 2.06e+01
penalty val = 8.20e-01
objective = 3.72e+00
beta = 1.4500e+01
step = 7

        enhancement = 3.75e+01
penalty val = 5.83e-01
objective = 1.56e+01
beta = 1.7200e+01
step = 8

        enhancement = 6.72e+01
penalty val = 4.88e-01
objective = 3.44e+01
beta = 1.9900e+01
step = 9

        enhancement = 1.30e+02
penalty val = 4.65e-01
objective = 6.94e+01
beta = 2.2600e+01
step = 10

        enhancement = 2.76e+02
penalty val = 4.66e-01
objective = 1.47e+02
beta = 2.5300e+01
step = 11

        enhancement = 3.23e+02
penalty val = 4.79e-01
objective = 1.68e+02
beta = 2.8000e+01
step = 12

        enhancement = 4.25e+02
penalty val = 4.27e-01
objective = 2.43e+02
beta = 3.0700e+01
step = 13

        enhancement = 4.78e+02
penalty val = 3.73e-01
objective = 3.00e+02
beta = 3.3400e+01
step = 14

        enhancement = 6.31e+02
penalty val = 3.26e-01
objective = 4.25e+02
beta = 3.6100e+01
step = 15

        enhancement = 6.11e+02
penalty val = 2.90e-01
objective = 4.34e+02
beta = 3.8800e+01
step = 16

        enhancement = 6.07e+02
penalty val = 2.63e-01
objective = 4.48e+02
beta = 4.1500e+01
step = 17

        enhancement = 6.11e+02
penalty val = 2.42e-01
objective = 4.63e+02
beta = 4.4200e+01
step = 18

        enhancement = 6.26e+02
penalty val = 2.24e-01
objective = 4.86e+02
beta = 4.6900e+01
step = 19

        enhancement = 6.48e+02
penalty val = 2.07e-01
objective = 5.14e+02
beta = 4.9600e+01
step = 20

        enhancement = 6.86e+02
penalty val = 2.04e-01
objective = 5.46e+02
beta = 5.2300e+01


## Analyze Results#

Letโs grab and evaluate the final parameters, since we exited the loop without doing this.

[29]:

params_final = params_history[-1]
beta_final = beta_history[-1]
objective_value_final = objective(params_final, beta=beta_final)
objective_values_history.append(objective_value_final)

        enhancement = 6.97e+02
penalty val = 2.23e-01
objective = 5.42e+02

[30]:

plt.plot(objective_values_history)
plt.xlabel("iterations")
plt.ylabel("objective function")
plt.show()

[ ]:




We can also generate a figure looking at the final fields and density pattern.

[31]:

sim_final = make_sim_with_antenna(params_final)

20:12:34 EDT Created task 'antenna final' with task_id

             View task using web UI at
4-40c5-91e9-46e7e5a5b6ea'.

20:12:44 EDT status = queued

             To cancel the simulation, use 'web.abort(task_id)' or
Terminating the Python script will not stop the job running on the
cloud.

20:12:51 EDT status = preprocess

20:12:54 EDT Maximum FlexCredit cost: 0.563. Use 'web.real_cost(task_id)' to get
the billed FlexCredit cost after a simulation run.

             starting up solver

             running solver

20:14:41 EDT early shutoff detected at 20%, exiting.

             status = postprocess

20:14:46 EDT status = success

             View simulation result at
4-40c5-91e9-46e7e5a5b6ea'.

20:14:47 EDT loading simulation from simulation_data.hdf5

[32]:

f, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(12, 8), tight_layout=False)

ax1.plot(objective_values_history)
ax1.set_xlabel("iterations")
ax1.set_ylabel("objective function")

density = get_density(params_final, beta=beta)
ax2.imshow(np.flipud(1 - density.squeeze().T), cmap="grey")
ax2.set_title("antenna density pattern")

sim_data_final.plot_field("field_xy", field_name="E", val="abs^2", ax=ax3)
sim_data_final.plot_field("field_xz", field_name="E", val="abs^2", ax=ax4)
ax3.set_title("intensity antenna (xy) plane")
ax4.set_title("intensity side (xz) plane")

plt.show()