The adjoint solver in meep now supports the adjoint simulation for near-to-far fields transformation. We present a simple optimization to illustrate this feature.
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
from scipy import special, signal
mp.verbosity(0)
Si = mp.Medium(index=3.4)
Air = mp.Medium(index=1.0)
Using MPI version 3.1, 1 processes
Basic setup
design_region_width = 15
design_region_height = 2
pml_size = 1.0
resolution = 30
Sx = 2 * pml_size + design_region_width
Sy = 2 * pml_size + design_region_height + 5
cell_size = mp.Vector3(Sx, Sy)
nf = 3
frequencies = np.array([1 / 1.5, 1 / 1.55, 1 / 1.6])
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(design_region_width, 0, 0)
src = mp.GaussianSource(frequency=fcen, fwidth=fwidth)
source = [mp.Source(src, component=mp.Ez, size=source_size, center=source_center)]
Nx = int(design_region_resolution * design_region_width) + 1
Ny = int(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)
projected_field = (
npa.flipud(projected_field) + projected_field
) / 2 # left-right symmetry
# interpolate to actual materials
return projected_field.flatten()
geometry = [
mp.Block(
center=design_region.center, size=design_region.size, material=design_variables
),
# mp.Block(center=design_region.center, size=design_region.size, material=design_variables, e1=mp.Vector3(x=-1))
#
# The commented lines above impose symmetry by overlapping design region with the same design variable. However,
# currently there is an issue of doing that; instead, we use an alternative approach to impose symmetry.
# See https://github.com/NanoComp/meep/issues/1984 and https://github.com/NanoComp/meep/issues/2093
]
kpoint = mp.Vector3()
sim = mp.Simulation(
cell_size=cell_size,
boundary_layers=pml_layers,
geometry=geometry,
sources=source,
default_material=Air,
symmetries=[mp.Mirror(direction=mp.X)],
resolution=resolution,
)
To use the mpa.Near2FarFields
objective, we need to specify the far point(s) of interest; and provide near2far field monitor(s) mp.Near2FarRegion
, from which the field at far point will be calculated. The monitor(s) has to capture all outgoing fields.
When evaluated, mpa.Near2FarFields
will return a numpy array with shape (num_of_points, nfreq, 6), where the third axis corresponds to the field components $E_x, E_y, E_z, H_x, H_y, H_z$, in that order. In general, we specify the mpa.Near2FarFields
objectives functions of the field components at frequencies of interest at points of interest. In this case, we would like to optimize $|E_z|^2$, and focus the fields of different frequency at the same point (0,15,0).
In this example, we try to focus the light at three different frequencies. Our objective is simply the mean of the fields at those frequencies.
far_x = [mp.Vector3(0, 15, 0)]
NearRegions = [
mp.Near2FarRegion(
center=mp.Vector3(0, design_region_height / 2 + 1.5),
size=mp.Vector3(design_region_width, 0),
weight=+1,
)
]
FarFields = mpa.Near2FarFields(sim, NearRegions, far_x)
ob_list = [FarFields]
def J1(FF):
return npa.mean(npa.abs(FF[0, :, 2]) ** 2)
opt = mpa.OptimizationProblem(
simulation=sim,
objective_functions=[J1],
objective_arguments=ob_list,
design_regions=[design_region],
frequencies=frequencies,
maximum_run_time=2000,
)
opt.plot2D(True)
evaluation_history = []
cur_iter = [0]
def f(v, gradient, cur_beta):
print("Current iteration: {}".format(cur_iter[0] + 1))
f0, dJ_du = opt([mapping(v, eta_i, cur_beta)]) # compute objective and gradient
if gradient.size > 0:
gradient[:] = tensor_jacobian_product(mapping, 0)(
v, eta_i, cur_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)
We'll now run our optimizer in loop. The loop will increase beta and reset the optimizer, which is important since the cost function changes.
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 = 2
num_betas = 6
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... 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...
plt.figure()
plt.plot(evaluation_history, "o-")
plt.grid(True)
plt.xlabel("Iteration")
plt.ylabel("FOM")
plt.show()
We can plot our results and see the resulting geometry.
opt.update_design([mapping(x, eta_i, cur_beta)])
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, it turns out that the optimizer focused on maximizing the performance at one frequency, and the structure behaved poorly at other two frequencies. The better way is to use the epigraph formulation, which is illustrated in the next tutorial.
f0, dJ_du = opt([mapping(x, eta_i, cur_beta // 2)], need_gradient=False)
frequencies = opt.frequencies
intensities = np.abs(opt.get_objective_arguments()[0][0, :, 2]) ** 2
Starting forward run...
plt.figure()
plt.plot(1 / frequencies, intensities, "-o")
plt.grid(True)
plt.xlabel("Wavelength (microns)")
plt.ylabel("|Ez|^2 Intensities")
plt.show()