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.

In [1]:
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()
CONFIGURATION FILE {'root': 'C:\\Program Files\\Lumerical\\FDTD\\api\\python', 'lumapi': 'C:\\Program Files\\Lumerical\\FDTD\\api\\python'}

User defined parameters

In [2]:
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:

In [3]:
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)
In [4]:
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:

In [5]:
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)
In [6]:
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.

In [7]:
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.

In [8]:
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=' ')
Accurate interface detection enabled
c:\users\jpond\appdata\local\programs\python\python37\lib\site-packages\matplotlib\figure.py:445: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.
  % get_backend())
Accurate interface detection enabled
Initializing super optimization
Wavelength range of source object will be superseded by the global settings.
Wavelength range of source object will be superseded by the global settings.
Running scipy optimizer
bounds = [[0.1 1. ]
 [0.1 1. ]
 [0.1 1. ]
 [0.1 1. ]
 [0.1 1. ]
 [0.1 1. ]
 [0.1 1. ]
 [0.1 1. ]
 [0.1 1. ]
 [0.1 1. ]]
start = [0.225      0.26388889 0.30277778 0.34166667 0.38055556 0.41944444
 0.45833333 0.49722222 0.53611111 0.575     ]
Running forward solves
FOM = 0.6929737002498921
Running forward solves
FOM = 0.7344183730030833
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 0 results
Plots updated with optimization 1 iteration 0 results
Saved frame
Running forward solves
FOM = 0.366577161281213
Running forward solves
FOM = 0.5057628350665708
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running forward solves
FOM = 0.8976685961997988
Running forward solves
FOM = 0.8626135277485157
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 1 results
Plots updated with optimization 1 iteration 1 results
Saved frame
Running forward solves
FOM = 0.9184145034186335
Running forward solves
FOM = 0.9017391545560013
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 2 results
Plots updated with optimization 1 iteration 2 results
Saved frame
Running forward solves
FOM = 0.914431843465571
Running forward solves
FOM = 0.9339108595256067
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 3 results
Plots updated with optimization 1 iteration 3 results
Saved frame
Running forward solves
FOM = 0.9316427286475351
Running forward solves
FOM = 0.9323592720756777
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 4 results
Plots updated with optimization 1 iteration 4 results
Saved frame
Running forward solves
FOM = 0.935984677816464
Running forward solves
FOM = 0.9345160072625376
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 5 results
Plots updated with optimization 1 iteration 5 results
Saved frame
Running forward solves
FOM = 0.9391325459433196
Running forward solves
FOM = 0.9362859098259368
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 6 results
Plots updated with optimization 1 iteration 6 results
Saved frame
Running forward solves
FOM = 0.9404556732509005
Running forward solves
FOM = 0.9393897702698981
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 7 results
Plots updated with optimization 1 iteration 7 results
Saved frame
Running forward solves
FOM = 0.9438766404967568
Running forward solves
FOM = 0.9434738981333566
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 8 results
Plots updated with optimization 1 iteration 8 results
Saved frame
Running forward solves
FOM = 0.9551783613854415
Running forward solves
FOM = 0.9353408407015839
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running forward solves
FOM = 0.947490806662743
Running forward solves
FOM = 0.9446003136587704
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running forward solves
FOM = 0.9509543512886547
Running forward solves
FOM = 0.9426161123304521
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 9 results
Plots updated with optimization 1 iteration 9 results
Saved frame
Running forward solves
FOM = 0.9559896077037915
Running forward solves
FOM = 0.9462176564638285
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 10 results
Plots updated with optimization 1 iteration 10 results
Saved frame
Running forward solves
FOM = 0.9537573604865681
Running forward solves
FOM = 0.958394795083825
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 11 results
Plots updated with optimization 1 iteration 11 results
Saved frame
Running forward solves
FOM = 0.9585648379261811
Running forward solves
FOM = 0.9673268215724918
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 12 results
Plots updated with optimization 1 iteration 12 results
Saved frame
Running forward solves
FOM = 0.9713385989725506
Running forward solves
FOM = 0.9554944557122483
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running forward solves
FOM = 0.9692680982545432
Running forward solves
FOM = 0.9657038218052519
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 13 results
Plots updated with optimization 1 iteration 13 results
Saved frame
Running forward solves
FOM = 0.9678336757403456
Running forward solves
FOM = 0.9722061950252814
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 14 results
Plots updated with optimization 1 iteration 14 results
Saved frame
Running forward solves
FOM = 0.9712096195964588
Running forward solves
FOM = 0.9748690771072419
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 15 results
Plots updated with optimization 1 iteration 15 results
Saved frame
Running forward solves
FOM = 0.9744634125495509
Running forward solves
FOM = 0.9751776264596861
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 16 results
Plots updated with optimization 1 iteration 16 results
Saved frame
Running forward solves
FOM = 0.9754635621842564
Running forward solves
FOM = 0.9748405269747723
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 17 results
Plots updated with optimization 1 iteration 17 results
Saved frame
Running forward solves
FOM = 0.9761712725996641
Running forward solves
FOM = 0.9746389666379829
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 18 results
Plots updated with optimization 1 iteration 18 results
Saved frame
Running forward solves
FOM = 0.9758137631503605
Running forward solves
FOM = 0.975377401244028
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 19 results
Plots updated with optimization 1 iteration 19 results
Saved frame
Running forward solves
FOM = 0.9751974558998519
Running forward solves
FOM = 0.9760240291535777
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 20 results
Plots updated with optimization 1 iteration 20 results
Saved frame
Running forward solves
FOM = 0.9753458033372067
Running forward solves
FOM = 0.9759353418288329
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 21 results
Plots updated with optimization 1 iteration 21 results
Saved frame
Running forward solves
FOM = 0.9755861857325558
Running forward solves
FOM = 0.9756815072420051
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running forward solves
FOM = 0.9754318751282209
Running forward solves
FOM = 0.9758559465512243
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 22 results
Plots updated with optimization 1 iteration 22 results
Saved frame
Running forward solves
FOM = 0.975729838504347
Running forward solves
FOM = 0.9755125823739006
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running forward solves
FOM = 0.975483920828763
Running forward solves
FOM = 0.9758087591950921
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 23 results
Plots updated with optimization 1 iteration 23 results
Saved frame
Running forward solves
FOM = 0.9755518463938421
Running forward solves
FOM = 0.9757365083763704
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running forward solves
FOM = 0.9755064094798527
Running forward solves
FOM = 0.9757892446232429
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 24 results
Plots updated with optimization 1 iteration 24 results
Saved frame
Running forward solves
FOM = 0.9755253662569988
Running forward solves
FOM = 0.975772983026233
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Running adjoint solves
Calculating gradients
Getting d eps: dx = 1e-10
..........
Plots updated with optimization 0 iteration 25 results
Plots updated with optimization 1 iteration 25 results
Saved frame
FINAL FOM = 1.9512983492832319
FINAL PARAMETERS = [0.23927992 0.26714977 0.29743278 0.50978788 0.65446577 0.635906
 0.68331354 0.6381805  0.65502561 0.62360218]

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

In [12]:
#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.

In [13]:
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.

In [14]:
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.

In [17]:
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()
Minimum transmission over bandwidth is 0.938234165446647
Minimum transmission over bandwidth factor is 0.8916490595367699
Bandwidth factor is 1.199999999999999
final design score is 1.069978871444123

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!