Tutorial example of inverse design of a taper¶
This example demonstrates inverse design using lumopt and FDTD Solutions. The example will design the taper/split section of a Y splitter for TE polarization modes. The shape to be optimized meets an input waveguide on the left side and two output waveguides on the right side. The sides of the taper is drawn with N points (currently N=10) evenly spaced in the X axis and the Y coordinates are determined by the optimization. A cubic spline is drawn through the N points to make a smooth shape. The device is made symmetric about the x axis. The design we will optimize is shown in the figure below, and the optimization region is shown by the dashed blue box.
The output waveguides are close together and strongly coupled at the output of the splitter contained in the dashed blue box. The figure of merit that we choose is the symmetric coupled mode of the 2 individual waveguide modes. We know that if we slowly separate the waveguides the power will transfer to the 2 TE modes with minimal loss or reflection. This example includes code to add output waveguides following a Bezier curve that separates the modes. This results in a complete splitter design.
The total size constraints for our splitter are shown in the figure. These contraints are fairly agressive to make it a challenging design problem!
We first optimize the device in 2D FDTD using an effective index approach. (We have precalculated that the slab effective index for 220nm Silicon encapsulated by n=1.44 glass is 2.85 using MODE Solutions.) This approximate approach is very accurate in this case. To confirm we run the optimization in 3D, but seed that optimization with the 2D result so it terminates very quickly.
The complete splitter design is saved into a gdsii mask file so it could be opened in a layout tool and fabricated. We also show how to do a full S-parameter extraction of the final device to obtain a compact model that could be used in circuit simulation.
Import required modules¶
There are quite a few, but you can just start your project with these.
You may have to change the path to the Lumerical FDTD folder if working on Linux/Mac. You don't need this line if you have configured Python to find the lumopt and lumapi modules.
import sys
sys.path.insert(0, "C:\\Program Files\\Lumerical\\FDTD\\api\\python\lumopt")
sys.path.insert(0, "C:\\Program Files\\Lumerical\\FDTD\\api\\python")
import os
import numpy as np
import scipy as sp
from lumopt import CONFIG
from lumopt.utilities.load_lumerical_scripts import load_from_lsf
from lumopt.utilities.wavelengths import Wavelengths
from lumopt.utilities.materials import Material
from lumopt.geometries.polygon import FunctionDefinedPolygon
from lumopt.figures_of_merit.modematch import ModeMatch
from lumopt.optimizers.generic_optimizers import ScipyOptimizers
from lumopt.optimization import Optimization
import lumapi
#os.chdir(os.path.dirname(__file__))
projdir = os.getcwd()
User defined parameters¶
Np = 10 # the number of control points
bc_type = 'clamped' # choose 'clamped' or 'not-a-knot' or 'natural' - see https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.interpolate.CubicSpline.html
wavelength1 = 1530e-9 # C+L band
wavelength2 = 1650e-9
#wavelength1 = 1230e-9 # O+C band
#wavelength2 = 1565e-9
wavelengths = Wavelengths(start = wavelength1, stop = wavelength2, points = 11)
bounds = [(0.1e-6, 1e-6)]*Np
rerun_optimization_2d = True
run_optimization_3d = False
delta = 20e-9 # determines the overetch/underetch values. Set to 0 for optimizing only the nominal design
inspect_base_project = True
display_example_polygon = True
Load the base project¶
The base project is defined with Lumerical script. This defines your input/output waveguides, input source, figure of merit measurement location, materials and simulation parameters. Basically everything but the shape to be optimized.
We will open an FDTD window and show the base project, but this is not a necessary step.
When building your own example it be easiest to start with one of the examples and change it. The Lumerical script follows a very simple pattern of creating simulation objects and setting properties. It is fairly straightfforward to change the vertical layer structure, change the device footprint or reposition the input/output waveguides.
The example will open an FDTD Solutions window and show you the base simulation so you can inspect it before continuing. It looks something like this:
script = load_from_lsf('splitter_base_TE.lsf')
if inspect_base_project:
#OPTIONAL: open FDTD session and inspect the base project
fdtd = lumapi.FDTD()
fdtd.eval(script)
if inspect_base_project:
#wait = input("PRESS ENTER TO CONTINUE.")
fdtd.close()
Define parametric geometry¶
This is the main part of the user input. The input to this funtion is a list of the parameter values. You can make as many parameters as you like, adjoint optimization is very efficient! The output of your function is a list of vertices of a polygon that will be drawn in FDTD Solutions. The parameterized geometry defines your design space. Choosing the parameterization carefully can limit the design space to devices that will manufacture well. You can also use the bounds to limit it from going into undesirable design space, for example limiting to a minimum feature size. It is best if your shape can't intersect itself.
For this design the edge of the taper is defined with a cubic spline with N nodes (N is the length of the input arguemnt params). The X values are just evenly spaced between right and left of taper. The Y values are the optimization parameters. Upper and lower edges are identical to get an even splitting ratio. The taper with is 500nm on the input side and 1.2um on the output, matching the input/output waveguides.
This example will open an FDTD Solutions window and show you the shape for a specific choice of parameters. It should look like this:
def taper_splitter(params, delta, bc_type):
Nx = len(params)
dx = 2e-6/Nx
points_x = np.concatenate(([-1e-6], np.linspace(-1e-6+0.5*dx,1e-6-0.5*dx,Nx), [1e-6]))
points_y = np.concatenate(([0.225e-6], params, [0.575e-6]))
px = np.linspace(min(points_x), max(points_x), 100)
interpolator = sp.interpolate.CubicSpline(points_x, points_y, bc_type = bc_type)
interpolator_prime = interpolator.derivative(nu=1)
py = interpolator(px)
pyp = interpolator_prime(px)
theta = np.arctan(pyp)
theta[0] = 0.
theta[-1] = 0.
px2 = px-delta*np.sin(theta)
py2 = py+delta*np.cos(theta)
px2[px2<px[0]] = px[0]
px2[px2>px[-1]] = px[-1]
polygon_points_up = [(x, y) for x, y in zip(px2, py2)]
polygon_points_down = [(x, -y) for x, y in zip(px2, py2)]
polygon_points = np.array(polygon_points_up[::-1] + polygon_points_down)
return polygon_points
if display_example_polygon:
test_params = np.fromstring('2.49111789e-07 2.69971496e-07 3.13218940e-07 5.58600573e-07 6.60844448e-07 6.41800230e-07 6.94548758e-07 6.36098641e-07 6.66422617e-07 6.32802984e-07',sep=' ')
test_vertices = taper_splitter(test_params,0,bc_type)
fdtd = lumapi.FDTD()
fdtd.addpoly(vertices=test_vertices)
if display_example_polygon:
#wait = input("PRESS ENTER TO CONTINUE.")
fdtd.close()
The 2D optimization definition¶
Below we have a function to assemble a 2d optimization. For first testing you may want to reduce the max iterations.
We created this as a functionb because we want to combine 2 optimizations to create our robust design.
def make_opt(initial_params, bounds, offset, bc_type, wavelengths):
#pass parameters to base script
base_script = ('offset=%.14e;\n' % offset) + script
func = lambda p:taper_splitter(p, offset, bc_type)
geometry = FunctionDefinedPolygon(func = func, initial_params = initial_params, bounds = bounds, z = 0.0, depth = 220e-9, eps_out = 1.44 ** 2, eps_in = 2.85 ** 2, edge_precision = 5, dx = 0.1e-9)
######## DEFINE FIGURE OF MERIT ########
fom = ModeMatch(monitor_name = 'fom', mode_number = 2, direction = 'Forward')
######## DEFINE OPTIMIZATION ALGORITHM ########
optimizer = ScipyOptimizers(max_iter = 25, method = 'L-BFGS-B', scaling_factor = 1e6)
######## PUT EVERYTHING TOGETHER ########
return Optimization(base_script = base_script, wavelengths = wavelengths, fom = fom, geometry = geometry, optimizer = optimizer, hide_fdtd_cad = True, use_deps = True)
Run the 2D optimization¶
Below we assemble and run the 2D optimization. If delta is 0, we use only one, otherwise, we add two optimizations to create a single fom from the overetched and underetched designs.
Optionally, if you set rerun_optimization_2d, you can simply use previously optimized parameters. These can be cut and pasted from the text output file easily.
if rerun_optimization_2d:
initial_params = np.linspace(0.225e-6, 0.575e-6, Np)
if delta==0:
opt = make_opt(initial_params,bounds,0,bc_type,wavelengths)
else:
opt = make_opt(initial_params,bounds,-delta,bc_type,wavelengths) + make_opt(initial_params,bounds,delta,bc_type,wavelengths)
######## RUN THE OPTIMIZER ########
result2d = opt.run()
result2d = np.array(result2d[1])*1e-6
os.chdir(projdir)
else:
#result_string = '2.49111789e-07 2.69971496e-07 3.13218940e-07 5.58600573e-07 6.60844448e-07 6.41800230e-07 6.94548758e-07 6.36098641e-07 6.66422617e-07 6.32802984e-07'
result_string = '2.39112467e-07 2.67451880e-07 2.97245655e-07 5.11033222e-07 6.55044127e-07 6.35547994e-07 6.83456742e-07 6.37961641e-07 6.54524012e-07 6.24185751e-07'
result2d = np.fromstring(result_string,sep=' ')
2D Expected results¶
Save the splitter as a GDSII file¶
To make a useful splitter we need to separate the output waveguides so the modes aren't coupled. A function is provided to add the output arms using a Bezier curve. We choose the lenght of the curve to limit loss and reflections. This part requires some trial and error or intuition. We could do this much easier with inverse design and that is an exercise left for the student!
Port markers are also added
#Cubic Bezier curve, with constant width
def fork_arm(length, offset, width, dx, dy):
t=np.linspace(0, 1, 100)
x = (3*(1-t)**2*t)*length/2 + (3*(1-t)*t**2)*length/2 + (t**3)*length
y = (3*(1-t)*t**2)*offset + (t**3)*offset
xp = (3*(1-t)**2)*length/2 + (3*t**2)*length/2
yp = (6*(1-t)*t)*offset
m = yp/xp
xt = np.cos(np.arctan(m) + np.pi/2)*width/2 + x
yt = np.sin(np.arctan(m) + np.pi/2)*width/2 + y
xb = np.cos(np.arctan(m) - np.pi/2)*width/2 + x
yb = np.sin(np.arctan(m) - np.pi/2)*width/2 + y
x = np.concatenate((xt, np.flip(xb,0))) + dx
y = np.concatenate((yt, np.flip(yb,0))) + dy
return np.column_stack((x,y))
def create_gds(fname,params,delta):
with lumapi.FDTD(hide=True) as f:
gds = f.gdsopen(fname)
f.gdsbegincell(gds, 'splitter')
v = taper_splitter(params,delta,bc_type)
v[:,0] += 1.1e-6
f.gdsaddpoly(gds, 1, v)
v = fork_arm(3.75e-6, 650e-9, 450e-9+2.*delta, 2.1e-6, 350e-9)
f.gdsaddpoly(gds, 1, v)
v = fork_arm(3.75e-6, -650e-9, 450e-9+2.*delta, 2.1e-6, -350e-9)
f.gdsaddpoly(gds, 1, v)
f.gdsaddrect(gds, 1, 0.05e-6, 0, 0.1e-6, 450e-9+2.*delta)
f.gdsaddrect(gds, "1:10", 5.85e-6, 1e-6, 200e-9, 450e-9+2.*delta)
f.gdsaddrect(gds, "1:10", 5.85e-6, -1e-6, 200e-9, 450e-9+2.*delta)
f.gdsaddrect(gds, "1:10", 0, 0, 200e-9, 450e-9+2.*delta)
f.gdsendcell(gds)
f.gdsclose(gds)
create_gds('splitter_2d.gds',result2d,0.)
if delta!=0.:
create_gds('splitter_2d_plus_delta.gds',result2d,delta)
create_gds('splitter_2d_minus_delta.gds',result2d,-delta)
GDSII output¶
An nice looking taper, shown here from KLayout!
It will also create gds files for the over and under etched designs (if you chose delta != 0). Here are both shown in KLayout.
Run a 3D optimization to refine resign¶
The 3D optimization is very simiular to 2D. The same parametric shape is used, the base project is the same but we pass in an additional variable to tell it that the 'sim_dimension' is '3D'.
This step is optional since it can take some time to run. Currently the maximum iterations are limited to 5. It is seeded with the best 2D results it shouldn't take many iterations if an improvement is possible.
def make_opt_3d(initial_params, bounds, offset, bc_type, wavelengths):
#pass parameters to base script
base_script = ('offset=%.14e;\n' % offset) + script
silicon = Material(name = 'Si (Silicon) - Palik', mesh_order = 2)
func = lambda p:taper_splitter(p, offset, bc_type)
geometry = FunctionDefinedPolygon(func = func, initial_params = initial_params, bounds = bounds, z = 0.0, depth = 220e-9, eps_out = 1.44 ** 2, eps_in = silicon, edge_precision = 5, dx = 0.1e-9)
######## DEFINE FIGURE OF MERIT ########
fom = ModeMatch(monitor_name = 'fom', mode_number = 1, direction = 'Forward')
######## DEFINE OPTIMIZATION ALGORITHM ########
optimizer = ScipyOptimizers(max_iter = 5, method = 'L-BFGS-B', scaling_factor = 1e6)
######## PUT EVERYTHING TOGETHER ########
return Optimization(base_script = base_script, wavelengths = wavelengths, fom = fom, geometry = geometry, optimizer = optimizer, hide_fdtd_cad = True, use_deps = True)
if run_optimization_3d:
script = "sim_dimension='3D';\n" + load_from_lsf('splitter_base_TE.lsf')
initial_params = result2d
if delta==0:
opt = make_opt_3d(initial_params,bounds,0,bc_type,wavelengths)
else:
opt = make_opt_3d(initial_params,bounds,-delta,bc_type,wavelengths) + make_opt_3d(initial_params,bounds,delta,bc_type,wavelengths)
######## RUN THE OPTIMIZER ########
result3d = opt.run()
result3d = np.array(result3d[1])*1e-6
os.chdir(projdir)
else:
result3d = result2d;
Create GDS for final simulation¶
Here we create the GDS from the 3D simulation. If you skipped the 3D simulaton then the 2D results are used.
create_gds('splitter_3d.gds',result3d,0.)
if delta!=0.:
create_gds('splitter_3d_plus_delta.gds',result3d,delta)
create_gds('splitter_3d_minus_delta.gds',result3d,-delta)
Extract S parameters and figure of merit¶
Here we setup an S parameter simulation of the full structure (including output arms) in 3D using the gds files we created.
def extract_parameters_3d(fname,gdsfile,delta=0.,calculate_S_parameters=True):
extract_script = load_from_lsf('splitter_base_model_extract.lsf')
S = 0
with lumapi.FDTD() as f:
f.putv("delta",delta)
f.eval(extract_script)
f.setglobalsource('wavelength start',wavelength1);
f.setglobalsource('wavelength stop',wavelength2);
f.setnamed("FDTD::ports","monitor frequency points",100);
f.gdsimport(gdsfile, "splitter", "1:0",'Si (Silicon) - Palik', -110e-9, 110e-9);
f.save(fname)
if calculate_S_parameters:
f.runsweep()
S = f.getsweepresult("s-parameter sweep","S matrix");
f.exportsweep("s-parameter sweep","splitter.dat")
return S
S = extract_parameters_3d("splitter_model_extract.fsp","splitter_3d.gds",0.,True)
if delta!=0.:
extract_parameters_3d("splitter_model_extract_plus_delta.fsp","splitter_3d_plus_delta.gds",delta,False)
extract_parameters_3d("splitter_model_extract_minus_delta.fsp","splitter_3d_minus_delta.gds",-delta,False)
# display and plot the final results!
Slambda = S['lambda']
Sfull = S['S']
S21 = Sfull[1,0,:]
S11 = Sfull[0,0,:]
S22 = Sfull[1,1,:]
Tmin = min(2.*abs(S21)**2)
pi = np.arccos(-1.)
Tmin_factor = np.arctan( 100*(Tmin-0.91))/pi+0.5
bw_factor = (max(Slambda)-min(Slambda))/100e-9
design_score = Tmin_factor*bw_factor;
print('Minimum transmission over bandwidth is',Tmin)
print('Minimum transmission over bandwidth factor is',Tmin_factor)
print('Bandwidth factor is',bw_factor[0])
print('final design score is',design_score[0])
import matplotlib.pyplot as plt
plt.figure(1)
plt.subplot(211)
plt.plot(Slambda*1e9,2.*abs(S21)**2)
plt.ylabel("Transmission")
plt.subplot(212)
plt.plot(Slambda*1e9,20.*np.log10(abs(S11)),label='S11')
plt.plot(Slambda*1e9,20.*np.log10(abs(S22)),label='S22')
plt.legend()
plt.xlabel("lambda (nm)")
plt.ylabel("Reflections (dB)")
plt.show()
View some results in FDTD Solutions¶
We did create fsp files for the nominal device (and overetch, underetch if relevant). You can open and run them and plot some results. Note that the comparison of over etch, under etch and nominal shows we could do well to increase the width of our waveguides a little at these longer wavelengths!
