This notebook will introduce inverse design using an example optimization task of a waveguide mode converter.
First, we will import the required python modules. Note that here we are importing several in addition to the ones we used in Notebook 01.
import numpy as np
import autograd.numpy as npa
import copy
import matplotlib as mpl
mpl.rcParams['figure.dpi']=100
import matplotlib.pylab as plt
from autograd.scipy.signal import convolve as conv
from skimage.draw import circle
import ceviche
from ceviche import fdfd_ez, jacobian
from ceviche.optimizers import adam_optimize
from ceviche.modes import insert_mode
import collections
# Create a container for our slice coords to be used for sources and probes
Slice = collections.namedtuple('Slice', 'x y')
As we saw in Notebook 01, the ceviche
package has a built-in method insert_mode()
that allows different modes to be inserted as sources.
Below, we demonstrate how this functionality can be used to excite the first and second order modes of a straight waveguide:
# Define simulation parameters (see above)
omega = 2*np.pi*200e12
dl = 25e-9
Nx = 200
Ny = 120
Npml = 20
# Define permittivity for a straight waveguide
epsr = np.ones((Nx, Ny))
epsr[:,50:67] = 12
# Source position
src_y = np.arange(20,100)
src_x = 30 * np.ones(src_y.shape, dtype=int)
# Source for mode 1
source1 = insert_mode(omega, dl, src_x, src_y, epsr, m=1)
# Source for mode 2
source2 = insert_mode(omega, dl, src_x, src_y, epsr, m=2)
# Run the simulation exciting mode 1
simulation = ceviche.fdfd_ez(omega, dl, epsr, [Npml, Npml])
Hx, Hy, Ez = simulation.solve(source1)
# Visualize the electric field
ax = ceviche.viz.real(Ez, outline=epsr, cmap='RdBu_r')
ax.plot(src_x,src_y,'k')
# Run the simulation exciting mode 2
simulation = ceviche.fdfd_ez(omega, dl, epsr, [Npml, Npml])
Hx, Hy, Ez = simulation.solve(source2)
# Visualize the electric field
ax = ceviche.viz.real(Ez, outline=epsr, cmap='RdBu_r')
ax.plot(src_x,src_y,'k')
plt.show()
Our toy optimization problem will be to design a device that converts an input in the first-order mode into an output as the second-order mode. First, we define the parameters of our device and optimization:
# Angular frequency of the source in Hz
omega=2*np.pi*200e12
# Spatial resolution in meters
dl=40e-9
# Number of pixels in x-direction
Nx=100
# Number of pixels in y-direction
Ny=100
# Number of pixels in the PMLs in each direction
Npml=20
# Initial value of the structure's relative permittivity
epsr_init=12.0
# Space between the PMLs and the design region (in pixels)
space=10
# Width of the waveguide (in pixels)
wg_width=12
# Length in pixels of the source/probe slices on each side of the center point
space_slice=8
# Number of epochs in the optimization
Nsteps=100
# Step size for the Adam optimizer
step_size=1e-2
We now define some utility functions for initialization and optimization:
def init_domain(Nx, Ny, Npml, space=10, wg_width=10, space_slice=5):
"""Initializes the domain and design region
space : The space between the PML and the structure
wg_width : The feed and probe waveguide width
space_slice : The added space for the probe and source slices
"""
# Parametrization of the permittivity of the structure
bg_epsr = np.zeros((Nx, Ny))
epsr = np.zeros((Nx, Ny))
# Region within which the permittivity is allowed to change
design_region = np.zeros((Nx, Ny))
# Input waveguide
bg_epsr[0:int(Npml+space),int(Ny/2-wg_width/2):int(Ny/2+wg_width/2)] = epsr_init
# Input probe slice
input_slice = Slice(x=np.array(Npml+1),
y=np.arange(int(Ny/2-wg_width/2-space_slice), int(Ny/2+wg_width/2+space_slice)))
# Output waveguide
bg_epsr[int(Nx-Npml-space)::,int(Ny/2-wg_width/2):int(Ny/2+wg_width/2)] = epsr_init
# Output probe slice
output_slice = (Slice(x=np.array(Nx-Npml-1),
y=np.arange(int(Ny/2-wg_width/2-space_slice), int(Ny/2+wg_width/2+space_slice))))
design_region[Npml+space:Nx-Npml-space, Npml+space:Ny-Npml-space] = 1
epsr[Npml+space:Nx-Npml-space, Npml+space:Ny-Npml-space] = epsr_init
return epsr, bg_epsr, design_region, input_slice, output_slice
def mask_combine_epsr(epsr, bg_epsr, design_region):
"""Utility function for combining the design region epsr and the background epsr
"""
return epsr*design_region + bg_epsr*(design_region==0).astype(np.float)
def viz_sim(epsr, source, slices=[]):
"""Solve and visualize a simulation with permittivity 'epsr'
"""
simulation = fdfd_ez(omega, dl, epsr, [Npml, Npml])
Hx, Hy, Ez = simulation.solve(source)
fig, ax = plt.subplots(1, 2, constrained_layout=True, figsize=(6,3))
ceviche.viz.real(Ez, outline=epsr, ax=ax[0], cbar=False)
for sl in slices:
ax[0].plot(sl.x*np.ones(len(sl.y)), sl.y, 'b-')
ceviche.viz.abs(epsr, ax=ax[1], cmap='Greys');
plt.show()
return (simulation, ax)
def mode_overlap(E1, E2):
"""Defines an overlap integral between the simulated field and desired field
"""
return npa.abs(npa.sum(npa.conj(E1)*E2))
We can visualize what our starting device looks like and how it behaves. Our device is initialized by the init_domain()
function which was defined several cells above.
# Initialize the parametrization rho and the design region
epsr, bg_epsr, design_region, input_slice, output_slice = \
init_domain(Nx, Ny, Npml, space=space, wg_width=wg_width, space_slice=space_slice)
epsr_total = mask_combine_epsr(epsr, bg_epsr, design_region)
# Setup source
source = insert_mode(omega, dl, input_slice.x, input_slice.y, epsr_total, m=1)
# Setup probe
probe = insert_mode(omega, dl, output_slice.x, output_slice.y, epsr_total, m=2)
# Simulate initial device
simulation, ax = viz_sim(epsr_total, source, slices=[input_slice,output_slice])
# get normalization factor (field overlap before optimizing)
_, _, Ez = simulation.solve(source)
E0 = mode_overlap(Ez, probe)
/Users/twhughes/anaconda3/lib/python3.7/site-packages/matplotlib/contour.py:1243: UserWarning: No contour levels were found within the data range. warnings.warn("No contour levels were found"
We will now define our objective function. This is a scalar-valued function which our optimizer uses to improve the device's performance.
Our objective function will consist of maximizing an overlap integral of the field in the output waveguide of the simulated device and the field of the waveguide's second order mode. The function takes in a single argument, epsr
and returns the value of the overlap integral. The details of setting the permittivity and solving for the fields happens inside the objective function.
def objective(epsr):
"""Objective function called by optimizer
1) Takes the epsr distribution as input
2) Runs the simulation
3) Returns the overlap integral between the output wg field
and the desired mode field
"""
epsr = epsr.reshape((Nx, Ny))
simulation.eps_r = mask_combine_epsr(epsr, bg_epsr, design_region)
_, _, Ez = simulation.solve(source)
return mode_overlap(Ez, probe) / E0
Now we can run our optimization. Note that we need to define our simulation object outside of the objective function and before the optimization is started.
# Simulate initial device
simulation, ax = viz_sim(epsr_total, source, slices=[input_slice, output_slice])
# Compute the gradient of the objective function using revere-mode differentiation
objective_jac = jacobian(objective, mode='reverse')
# Maximize the objective function using an ADAM optimizer
(epsr_optimum, loss) = adam_optimize(objective, epsr.flatten(), objective_jac,
Nsteps=Nsteps, direction='max', step_size=step_size)
# Simulate and show the optimal device
epsr_optimum = epsr_optimum.reshape((Nx, Ny))
epsr_optimum_total = mask_combine_epsr(epsr_optimum, bg_epsr, design_region)
simulation, ax = viz_sim(epsr_optimum_total, source, slices=[input_slice, output_slice])
Epoch: 1/100 | Duration: 0.09 secs | Value: 1.000000e+00 Epoch: 2/100 | Duration: 0.08 secs | Value: 3.984201e+03 Epoch: 3/100 | Duration: 0.09 secs | Value: 7.957030e+03 Epoch: 4/100 | Duration: 0.09 secs | Value: 1.191068e+04 Epoch: 5/100 | Duration: 0.08 secs | Value: 1.583605e+04 Epoch: 6/100 | Duration: 0.08 secs | Value: 1.972413e+04 Epoch: 7/100 | Duration: 0.09 secs | Value: 2.356618e+04 Epoch: 8/100 | Duration: 0.09 secs | Value: 2.735380e+04 Epoch: 9/100 | Duration: 0.08 secs | Value: 3.107895e+04 Epoch: 10/100 | Duration: 0.09 secs | Value: 3.473394e+04 Epoch: 11/100 | Duration: 0.09 secs | Value: 3.831152e+04 Epoch: 12/100 | Duration: 0.09 secs | Value: 4.180494e+04 Epoch: 13/100 | Duration: 0.08 secs | Value: 4.520808e+04 Epoch: 14/100 | Duration: 0.08 secs | Value: 4.851556e+04 Epoch: 15/100 | Duration: 0.08 secs | Value: 5.172283e+04 Epoch: 16/100 | Duration: 0.09 secs | Value: 5.482626e+04 Epoch: 17/100 | Duration: 0.08 secs | Value: 5.782317e+04 Epoch: 18/100 | Duration: 0.08 secs | Value: 6.071186e+04 Epoch: 19/100 | Duration: 0.08 secs | Value: 6.349160e+04 Epoch: 20/100 | Duration: 0.08 secs | Value: 6.616260e+04 Epoch: 21/100 | Duration: 0.08 secs | Value: 6.872598e+04 Epoch: 22/100 | Duration: 0.08 secs | Value: 7.118368e+04 Epoch: 23/100 | Duration: 0.08 secs | Value: 7.353849e+04 Epoch: 24/100 | Duration: 0.08 secs | Value: 7.579393e+04 Epoch: 25/100 | Duration: 0.08 secs | Value: 7.795418e+04 Epoch: 26/100 | Duration: 0.08 secs | Value: 8.002397e+04 Epoch: 27/100 | Duration: 0.09 secs | Value: 8.200852e+04 Epoch: 28/100 | Duration: 0.08 secs | Value: 8.391338e+04 Epoch: 29/100 | Duration: 0.08 secs | Value: 8.574435e+04 Epoch: 30/100 | Duration: 0.08 secs | Value: 8.750735e+04 Epoch: 31/100 | Duration: 0.09 secs | Value: 8.920834e+04 Epoch: 32/100 | Duration: 0.08 secs | Value: 9.085323e+04 Epoch: 33/100 | Duration: 0.08 secs | Value: 9.244773e+04 Epoch: 34/100 | Duration: 0.08 secs | Value: 9.399738e+04 Epoch: 35/100 | Duration: 0.08 secs | Value: 9.550740e+04 Epoch: 36/100 | Duration: 0.08 secs | Value: 9.698266e+04 Epoch: 37/100 | Duration: 0.09 secs | Value: 9.842768e+04 Epoch: 38/100 | Duration: 0.08 secs | Value: 9.984654e+04 Epoch: 39/100 | Duration: 0.08 secs | Value: 1.012429e+05 Epoch: 40/100 | Duration: 0.08 secs | Value: 1.026201e+05 Epoch: 41/100 | Duration: 0.08 secs | Value: 1.039809e+05 Epoch: 42/100 | Duration: 0.08 secs | Value: 1.053277e+05 Epoch: 43/100 | Duration: 0.08 secs | Value: 1.066625e+05 Epoch: 44/100 | Duration: 0.08 secs | Value: 1.079871e+05 Epoch: 45/100 | Duration: 0.08 secs | Value: 1.093027e+05 Epoch: 46/100 | Duration: 0.09 secs | Value: 1.106104e+05 Epoch: 47/100 | Duration: 0.08 secs | Value: 1.119109e+05 Epoch: 48/100 | Duration: 0.08 secs | Value: 1.132046e+05 Epoch: 49/100 | Duration: 0.08 secs | Value: 1.144919e+05 Epoch: 50/100 | Duration: 0.08 secs | Value: 1.157728e+05 Epoch: 51/100 | Duration: 0.08 secs | Value: 1.170471e+05 Epoch: 52/100 | Duration: 0.08 secs | Value: 1.183146e+05 Epoch: 53/100 | Duration: 0.08 secs | Value: 1.195750e+05 Epoch: 54/100 | Duration: 0.08 secs | Value: 1.208279e+05 Epoch: 55/100 | Duration: 0.08 secs | Value: 1.220726e+05 Epoch: 56/100 | Duration: 0.08 secs | Value: 1.233086e+05 Epoch: 57/100 | Duration: 0.08 secs | Value: 1.245354e+05 Epoch: 58/100 | Duration: 0.08 secs | Value: 1.257522e+05 Epoch: 59/100 | Duration: 0.08 secs | Value: 1.269585e+05 Epoch: 60/100 | Duration: 0.08 secs | Value: 1.281537e+05 Epoch: 61/100 | Duration: 0.08 secs | Value: 1.293372e+05 Epoch: 62/100 | Duration: 0.08 secs | Value: 1.305083e+05 Epoch: 63/100 | Duration: 0.08 secs | Value: 1.316667e+05 Epoch: 64/100 | Duration: 0.08 secs | Value: 1.328117e+05 Epoch: 65/100 | Duration: 0.08 secs | Value: 1.339429e+05 Epoch: 66/100 | Duration: 0.08 secs | Value: 1.350599e+05 Epoch: 67/100 | Duration: 0.08 secs | Value: 1.361624e+05 Epoch: 68/100 | Duration: 0.08 secs | Value: 1.372500e+05 Epoch: 69/100 | Duration: 0.08 secs | Value: 1.383225e+05 Epoch: 70/100 | Duration: 0.08 secs | Value: 1.393796e+05 Epoch: 71/100 | Duration: 0.08 secs | Value: 1.404211e+05 Epoch: 72/100 | Duration: 0.08 secs | Value: 1.414469e+05 Epoch: 73/100 | Duration: 0.08 secs | Value: 1.424569e+05 Epoch: 74/100 | Duration: 0.08 secs | Value: 1.434509e+05 Epoch: 75/100 | Duration: 0.08 secs | Value: 1.444290e+05 Epoch: 76/100 | Duration: 0.08 secs | Value: 1.453911e+05 Epoch: 77/100 | Duration: 0.08 secs | Value: 1.463372e+05 Epoch: 78/100 | Duration: 0.08 secs | Value: 1.472673e+05 Epoch: 79/100 | Duration: 0.09 secs | Value: 1.481814e+05 Epoch: 80/100 | Duration: 0.09 secs | Value: 1.490797e+05 Epoch: 81/100 | Duration: 0.09 secs | Value: 1.499621e+05 Epoch: 82/100 | Duration: 0.09 secs | Value: 1.508287e+05 Epoch: 83/100 | Duration: 0.09 secs | Value: 1.516797e+05 Epoch: 84/100 | Duration: 0.10 secs | Value: 1.525150e+05 Epoch: 85/100 | Duration: 0.09 secs | Value: 1.533349e+05 Epoch: 86/100 | Duration: 0.09 secs | Value: 1.541395e+05 Epoch: 87/100 | Duration: 0.09 secs | Value: 1.549287e+05 Epoch: 88/100 | Duration: 0.09 secs | Value: 1.557028e+05 Epoch: 89/100 | Duration: 0.09 secs | Value: 1.564619e+05 Epoch: 90/100 | Duration: 0.10 secs | Value: 1.572060e+05 Epoch: 91/100 | Duration: 0.10 secs | Value: 1.579354e+05 Epoch: 92/100 | Duration: 0.10 secs | Value: 1.586501e+05 Epoch: 93/100 | Duration: 0.10 secs | Value: 1.593504e+05 Epoch: 94/100 | Duration: 0.10 secs | Value: 1.600362e+05 Epoch: 95/100 | Duration: 0.09 secs | Value: 1.607078e+05 Epoch: 96/100 | Duration: 0.10 secs | Value: 1.613652e+05 Epoch: 97/100 | Duration: 0.10 secs | Value: 1.620088e+05 Epoch: 98/100 | Duration: 0.11 secs | Value: 1.626385e+05 Epoch: 99/100 | Duration: 0.11 secs | Value: 1.632545e+05 Epoch: 100/100 | Duration: 0.10 secs | Value: 1.638571e+05
At the end of the optimization we can see our final device. From the field pattern, we can easily observe that the device is doing what we intend: the even mode enters from the left and exits as the odd mode on the right.
However, an additional observation is that our device's permittivity changes continuously. This is not ideal if we wanted to fabricated our device. We're also not constraining the minimum and maximum values of $\epsilon_r$. Thus, we need to consider alternative ways of parameterizing our device.