This tutorial illustates the adjoint solver for DFT fields. Specifically, we use the FourierFields
adjoint solver to design a metalens that maximizes the average intensity over three frequencies at a spot near the metalens. In a similar example based on the Near2FarFields
adjoint solver, the focal spot is far from the metalens. Compared with Near2FarFields
, the FourierFields
adjoint solver has some advantages, including higher efficiency for fields at nearby points and greater flexibility for surrounding media (not limited to homogeneous media surrounding far-field regions when Near2FarFields
is applied). On the other hand, the disadvantage is that FourierFields
is less efficient for faraway points, which require a large computational cell.
At first, let us import meep
, its adjoint module meep.adjoint
, and other packages for automatic differentiation, nonlinear optimization, and making figures.
import meep as mp
import meep.adjoint as mpa
import numpy as np
from autograd import numpy as npa
from autograd import tensor_jacobian_product, grad
import nlopt
from matplotlib import pyplot as plt
from matplotlib.patches import Circle
mp.verbosity(0)
Si = mp.Medium(index=3.4)
Air = mp.Medium(index=1.0)
Basic setup
resolution = 30
design_region_width = 4
design_region_height = 2
pml_size = 1.0
Sx = 2 * pml_size + design_region_width
Sy = 2 * pml_size + design_region_height + 5
cell_size = mp.Vector3(Sx, Sy)
wavelengths = np.array([1.5, 1.55, 1.6])
frequencies = 1/wavelengths
nf = len(frequencies) # number of frequencies
minimum_length = 0.09 # minimum length scale (microns)
eta_i = 0.5 # blueprint (or intermediate) design field thresholding point (between 0 and 1)
eta_e = 0.55 # erosion design field thresholding point (between 0 and 1)
eta_d = 1 - eta_e # dilation design field thresholding point (between 0 and 1)
filter_radius = mpa.get_conic_radius_from_eta_e(minimum_length, eta_e)
design_region_resolution = int(resolution)
pml_layers = [mp.PML(pml_size)]
fcen = 1 / 1.55
width = 0.2
fwidth = width * fcen
source_center = [0, -(design_region_height / 2 + 1.5), 0]
source_size = mp.Vector3(Sx, 0, 0)
src = mp.GaussianSource(frequency=fcen, fwidth=fwidth, is_integrated=True)
source = [mp.Source(src, component=mp.Ez, size=source_size, center=source_center)]
Nx = int(round(design_region_resolution * design_region_width)) + 1
Ny = int(round(design_region_resolution * design_region_height)) + 1
design_variables = mp.MaterialGrid(mp.Vector3(Nx, Ny), Air, Si, grid_type="U_MEAN")
design_region = mpa.DesignRegion(
design_variables,
volume=mp.Volume(
center=mp.Vector3(),
size=mp.Vector3(design_region_width, design_region_height, 0),
),
)
def mapping(x, eta, beta):
# filter
filtered_field = mpa.conic_filter(
x,
filter_radius,
design_region_width,
design_region_height,
design_region_resolution,
)
# projection
projected_field = mpa.tanh_projection(filtered_field, beta, eta)
# interpolate to actual materials
return projected_field.flatten()
geometry = [
mp.Block(
center=design_region.center, size=design_region.size, material=design_variables
)
]
sim = mp.Simulation(
cell_size=cell_size,
boundary_layers=pml_layers,
geometry=geometry,
sources=source,
default_material=Air,
resolution=resolution,
)
To use the mpa.FourierFields
adjoint solver, we need to specify the monitor where DFT fields are evaluated. In this example, we try to focus the light at three different frequencies. The objective function is simply the intensity at the focal spot averaged over those frequencies.
monitor_position, monitor_size = mp.Vector3(0,design_region_height/2+1.5), mp.Vector3(0.01,0)
FourierFields = mpa.FourierFields(sim,mp.Volume(center=monitor_position,size=monitor_size),mp.Ez,yee_grid=True)
ob_list = [FourierFields]
def J(fields):
return npa.mean(npa.abs(fields[:,1]) ** 2) # The index 1 corresponds to the point at the center of our monitor.
opt = mpa.OptimizationProblem(
simulation=sim,
objective_functions=[J],
objective_arguments=ob_list,
design_regions=[design_region],
frequencies=frequencies,
maximum_run_time=2000,
)
opt.plot2D(True)
To apply nlopt
for nonlinear optimization, the following function is defined. This function also illustrates the structure in the design region during optimization and records the evaluation history of the objective function.
evaluation_history = []
cur_iter = [0]
def f(v, gradient, beta):
print("Current iteration: {}".format(cur_iter[0] + 1))
f0, dJ_du = opt([mapping(v, eta_i, beta)]) # compute objective and gradient
if gradient.size > 0:
gradient[:] = tensor_jacobian_product(mapping, 0)(
v, eta_i, beta, np.sum(dJ_du, axis=1)
) # backprop
evaluation_history.append(np.real(f0))
plt.figure()
ax = plt.gca()
opt.plot2D(
False,
ax=ax,
plot_sources_flag=False,
plot_monitors_flag=False,
plot_boundaries_flag=False,
)
circ = Circle((2, 2), minimum_length / 2)
ax.add_patch(circ)
ax.axis("off")
plt.show()
cur_iter[0] = cur_iter[0] + 1
return np.real(f0)
The optimizer runs in multiple loops with increasing beta, which gradually turns on the thresholding. When beta changes, the cost function also changes, which requires the optimizer to be reset.
algorithm = nlopt.LD_MMA
n = Nx * Ny # number of parameters
# Initial guess
x = np.ones((n,)) * 0.5
# lower and upper bounds
lb = np.zeros((Nx * Ny,))
ub = np.ones((Nx * Ny,))
cur_beta = 4
beta_scale = 1.5
num_betas = 9
update_factor = 12
ftol = 1e-5
for iters in range(num_betas):
solver = nlopt.opt(algorithm, n)
solver.set_lower_bounds(lb)
solver.set_upper_bounds(ub)
solver.set_max_objective(lambda a, g: f(a, g, cur_beta))
solver.set_maxeval(update_factor)
solver.set_ftol_rel(ftol)
x[:] = solver.optimize(x)
cur_beta = cur_beta * beta_scale
Current iteration: 1 Starting forward run... Starting adjoint run...
/home/mawc/install/envs/mp/lib/python3.10/site-packages/meep/adjoint/filter_source.py:152: RuntimeWarning: invalid value encountered in true_divide l2_err = np.sum(np.abs(H - H_hat.T)**2 / np.abs(H)**2)
Calculating gradient...
Current iteration: 2 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 3 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 4 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 5 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 6 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 7 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 8 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 9 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 10 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 11 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 12 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 13 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 14 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 15 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 16 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 17 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 18 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 19 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 20 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 21 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 22 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 23 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 24 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 25 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 26 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 27 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 28 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 29 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 30 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 31 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 32 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 33 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 34 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 35 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 36 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 37 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 38 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 39 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 40 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 41 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 42 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 43 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 44 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 45 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 46 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 47 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 48 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 49 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 50 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 51 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 52 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 53 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 54 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 55 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 56 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 57 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 58 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 59 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 60 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 61 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 62 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 63 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 65 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 66 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 67 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 68 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 69 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 70 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 71 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 72 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 73 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 74 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 75 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 76 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 77 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 78 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 79 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 80 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 81 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 82 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 83 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 84 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 85 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 86 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 87 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 88 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 89 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 90 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 91 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 92 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 93 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 94 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 95 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 96 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 97 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 98 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 99 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 100 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 101 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 102 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 103 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 104 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 105 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 106 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 107 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 108 Starting forward run... Starting adjoint run... Calculating gradient...
plt.figure()
plt.plot(evaluation_history, "o-")
plt.grid(True)
plt.xlabel("Iteration")
plt.ylabel("FOM")
plt.show()
We can plot the optimized geometry.
opt.update_design([mapping(x, eta_i, cur_beta/beta_scale)]) # cur_beta/beta_scale is the final beta in the optimization.
plt.figure()
ax = plt.gca()
opt.plot2D(
False,
ax=ax,
plot_sources_flag=False,
plot_monitors_flag=False,
plot_boundaries_flag=False,
)
circ = Circle((2, 2), minimum_length / 2)
ax.add_patch(circ)
ax.axis("off")
plt.show()
To check the performance of our final structure, we plot $|E_z|^2$ at those three frequencies. Because our objective function only concerns the average of intensities, the optimizer does not necessarily maximize the performance at every frequency.
f0, dJ_du = opt([mapping(x, eta_i, cur_beta / beta_scale)], need_gradient=False)
intensities = np.abs(opt.get_objective_arguments()[0][:,1]) ** 2
plt.figure()
plt.plot(wavelengths, intensities, "-o")
plt.grid(True)
plt.xlabel("Wavelength (μm)")
plt.ylabel("|Ez|^2 intensity (a.u.)")
plt.show()
Starting forward run...
Let us launch continuous-wave sources at individual frequencies and illustrate the fields.
opt.sim = mp.Simulation(
cell_size=mp.Vector3(Sx, Sy),
boundary_layers=pml_layers,
geometry=geometry,
sources=source,
default_material=Air,
resolution=resolution,
)
src = mp.ContinuousSource(frequency=frequencies[0], fwidth=fwidth, is_integrated=True)
source = [mp.Source(src, component=mp.Ez, size=source_size, center=source_center)]
opt.sim.change_sources(source)
plt.figure(figsize=(16, 10))
plt.subplot(1,3,1)
plt.title('wavelength = '+str(wavelengths[0]))
opt.sim.run(until=300)
opt.sim.plot2D(fields=mp.Ez)
opt.sim.reset_meep()
opt.sim = mp.Simulation(
cell_size=mp.Vector3(Sx, Sy),
boundary_layers=pml_layers,
geometry=geometry,
sources=source,
default_material=Air,
resolution=resolution,
)
src = mp.ContinuousSource(frequency=frequencies[1], fwidth=fwidth, is_integrated=True)
source = [mp.Source(src, component=mp.Ez, size=source_size, center=source_center)]
opt.sim.change_sources(source)
plt.subplot(1,3,2)
plt.title('wavelength = '+str(wavelengths[1]))
opt.sim.run(until=300)
opt.sim.plot2D(fields=mp.Ez)
opt.sim.reset_meep()
opt.sim = mp.Simulation(
cell_size=mp.Vector3(Sx, Sy),
boundary_layers=pml_layers,
geometry=geometry,
sources=source,
default_material=Air,
resolution=resolution,
)
src = mp.ContinuousSource(frequency=frequencies[2], fwidth=fwidth, is_integrated=True)
source = [mp.Source(src, component=mp.Ez, size=source_size, center=source_center)]
opt.sim.change_sources(source)
plt.subplot(1,3,3)
plt.title('wavelength = '+str(wavelengths[2]))
opt.sim.run(until=300)
opt.sim.plot2D(fields=mp.Ez)
opt.sim.reset_meep()