The adjoint solver in meep now supports the adjoint simulation for near-to-far fields transformation. We present a simple optimization of metalens using the epigraph formulation.
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.quiet(quietval=True)
Si = mp.Medium(index=3.4)
SiO2 = mp.Medium(index=1.44)
Basic setup
design_region_width = 32
design_region_height = 2
pml_size = 1.0
resolution = 30
Sx = design_region_width
Sy = 2*pml_size + design_region_height + 5
cell_size = mp.Vector3(Sx,Sy)
nf = 2
frequencies = np.array([1/1.52, 1/1.58])
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(1*resolution)
pml_layers = [mp.PML(pml_size, direction=mp.Y)]
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)
source = [mp.EigenModeSource(src,
eig_band = 1,
size = source_size,
center=source_center)]
Nx = int(design_region_resolution*design_region_width)
Ny = int(design_region_resolution*design_region_height)
design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,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),
mp.Block(center=design_region.center, size=design_region.size, material=design_variables, e1=mp.Vector3(x=-1))] # design region
kpoint = mp.Vector3()
sim = mp.Simulation(cell_size=cell_size,
boundary_layers=pml_layers,
k_point=kpoint,
geometry=geometry,
sources=source,
default_material=SiO2,
symmetries=[mp.Mirror(direction=mp.X)],
resolution=resolution)
We will have two objective functions, one for each far point. However, we only need one mpa.Near2FarFields
objective. We also need to provide a near-field monitor, from which the field at far point will be calculated. Only single monitor is supported right now, so the monitor needs to extend to the entire cell to capture all outgoing fields.
When evaluated, mpa.Near2FarFields will return a numpy array with shape npoints by nfreq by 6 numpy array, where the third axis corresponds to the field components $E_x, E_y, E_z, H_x, H_y, H_z$, in that order. We will specify a objective as a function 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 different points.
far_x = [mp.Vector3(0,40,0)]
NearRegions = [mp.Near2FarRegion(center=mp.Vector3(0,design_region_height/2+1.5), size=mp.Vector3(design_region_width+1.5,0), weight=+1)]
FarFields = mpa.Near2FarFields(sim, NearRegions ,far_x)
ob_list = [FarFields]
def J1(alpha):
return -npa.abs(alpha[0,:,2])**2
opt = mpa.OptimizationProblem(
simulation = sim,
objective_functions = [J1],
objective_arguments = ob_list,
design_regions = [design_region],
frequencies=frequencies,
decay_by = 1e-4,
decay_fields=[mp.Ez],
maximum_run_time = 2000
)
Our objective function that we pass to nlopt is rather simple. We'll introduce a dummy parameter t
. The goal of the optimization problem will be to simply minimize the value of t
. The gradient of this functional is rather straightforward.
evaluation_history = []
cur_iter = [0]
def f(x, grad):
t = x[0] # "dummy" parameter
v = x[1:] # design parameters
if grad.size > 0:
grad[0] = 1
grad[1:] = 0
return t
The key to the epigraph formulation (and most nonlinear optimization problems) lies in the constraints. We'll define one constraint for every frequency point of interest ($f_i$). We'll define our constraint as
$$f_i < t$$where $t$ is the previosuly defined dummy parameter. Each constraint function is then defined by
$$ c_i = f_i-t $$within some tolerance.
More detail about this formulation can be found in the nlopt documentation.
def c(result,x,gradient,eta,beta):
print("Current iteration: {}; current eta: {}, current beta: {}".format(cur_iter[0],eta,beta))
t = x[0] # dummy parameter
v = x[1:] # design parameters
f0, dJ_du = opt([mapping(v,eta,beta)])
# Backprop the gradients through our mapping function
my_grad = np.zeros(dJ_du.shape)
for k in range(opt.nf):
my_grad[:,k] = tensor_jacobian_product(mapping,0)(v,eta,beta,dJ_du[:,k])
# Assign gradients
if gradient.size > 0:
gradient[:,0] = -1 # gradient w.r.t. "t"
gradient[:,1:] = my_grad.T # gradient w.r.t. each frequency objective
result[:] = np.real(f0) - t
# store results
evaluation_history.append(np.real(f0))
# visualize
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
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,))
# insert dummy parameter bounds and variable
x = np.insert(x,0,0) # our initial guess for the worst error
lb = np.insert(lb,0,-np.inf)
ub = np.insert(ub,0,0)
cur_beta = 64
beta_scale = 2
num_betas = 4
update_factor = 15
ftol = 1e-5
for iters in range(num_betas):
solver = nlopt.opt(algorithm, n+1)
solver.set_lower_bounds(lb)
solver.set_upper_bounds(ub)
solver.set_min_objective(f)
solver.set_maxeval(update_factor)
solver.set_ftol_rel(ftol)
solver.add_inequality_mconstraint(lambda r,x,g: c(r,x,g,eta_i,cur_beta), np.array([1e-3]*nf))
x[:] = solver.optimize(x)
cur_beta = cur_beta*beta_scale
Current iteration: 0; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run...
/usr/local/lib/python3.8/site-packages/meep/adjoint/filter_source.py:91: 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: 1; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 2; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 3; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 4; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 5; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 6; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 7; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 8; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 9; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 10; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 11; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 12; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 13; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 14; current eta: 0.5, current beta: 64 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 15; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 16; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 17; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 18; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 19; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 20; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 21; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 22; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 23; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 24; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 25; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 26; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 27; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 28; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 29; current eta: 0.5, current beta: 128 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 30; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 31; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 32; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 33; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 34; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 35; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 36; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 37; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 38; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 39; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 40; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 41; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 42; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 43; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 44; current eta: 0.5, current beta: 256 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 45; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 46; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 47; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 48; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 49; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 50; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 51; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 52; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 53; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 54; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 55; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 56; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 57; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 58; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
Current iteration: 59; current eta: 0.5, current beta: 512 Starting forward run... Starting adjoint run... Calculating gradient...
lb = -np.min(evaluation_history,axis=1)
ub = -np.max(evaluation_history,axis=1)
mean = -np.mean(evaluation_history,axis=1)
num_iters = lb.size
plt.figure()
plt.fill_between(np.arange(num_iters),ub,lb,alpha=0.3)
plt.plot(mean,'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[1:],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 use a CW source at the desired frequency and plot the fields after the struture.
opt.sim = mp.Simulation(cell_size=mp.Vector3(Sx,90),
boundary_layers=pml_layers,
k_point=kpoint,
geometry=geometry,
sources=source,
default_material=SiO2,
resolution=resolution)
src = mp.ContinuousSource(frequency=1/1.52,fwidth=fwidth)
source = [mp.EigenModeSource(src,
eig_band = 1,
size = source_size,
center=source_center)]
opt.sim.change_sources(source)
opt.sim.run(until=400)
plt.figure(figsize=(10,20))
opt.sim.plot2D(fields=mp.Ez)
<AxesSubplot:xlabel='X', ylabel='Y'>
opt.sim = mp.Simulation(cell_size=mp.Vector3(Sx,90),
boundary_layers=pml_layers,
k_point=kpoint,
geometry=geometry,
sources=source,
default_material=SiO2,
resolution=resolution)
src = mp.ContinuousSource(frequency=1/1.58,fwidth=fwidth)
source = [mp.EigenModeSource(src,
eig_band = 1,
size = source_size,
center=source_center)]
opt.sim.change_sources(source)
opt.sim.run(until=200)
plt.figure(figsize=(10,20))
opt.sim.plot2D(fields=mp.Ez)
<AxesSubplot:xlabel='X', ylabel='Y'>