Inverse design of y-branch

Use the built-in parametric optimization (PO) method for inverse design of a silicon-on-insulator Y-branch to achieve the maximum broadband transmission. The varFDTD 2.5D solver now supports this technique allowing for realistic material models and monolithic planar devices. The adjoint method is the basis for the efficient optimization algorithm which significantly reduces the number of iterations required for convergence in comparison with other optimization methods such as particle swarm. The optimal design is then exported into a GDS file for further simulation and/or fabrication.

download example

Overview

Understand the simulation workflow and key results

Lumerical's inverse design capability provides unparalleled optimization performance by combining the power of gradient-based optimization routine with efficiencies found in fundamental properties of Maxwell equations. This example will demonstrate how to use the inverse design method to generate a Y branch shape that delivers near-perfect transmission efficiency. Moreover, we will demonstrate how to modify the example with your parameters so you can reuse this approach for your own design.

This example draws extensively from the LumOpt framework:

• Git
• Documentation

Step 1: Define base simulation parameters

The goal for this initial step is to define the static inputs of the varFDTD simulation that will be used for the adjoint optimization. Key simulation settings for the input/output waveguides, monitors for field gradient and figure of merit (FOM) calculations, source and mesh are parametrized for ease of running the optimization routine. It is possible to provide these necessary parameters with either

• MODE project file
• lumerical script (.lsf format)
• Python function handle.

In this example we will demonstrate how to use the python interoperability API to define a function that works with the LumOpt optimization routines, allowing users to work exclusively within python.

Step 2: Define the optimizable geometry

The second step is to define a parametrized polygon geometry that will represent the body of the Y branch connecting the input and output waveguide (the optimized geometry). This is done with a python function that accepts a set of optimization parameters and generates a corresponding polygon geometry in the base simulation. Initial optimization parameter values and their limits will also get defined in this step.

Step 3: Run the parametric optimization in 2.5D

Run the optimization python script with the base simulation, and parametrized polygon object as an input. This will be a 2.5D simulation for the fast delivery of a near optimal design which can be used as the initial guess for the 3D optimization in the next step. To find the best set of parameters that define the optimal Y branch shape, the optimization routine will use the specified field monitors to calculate:

• the field gradients: the gradients of the field due to perturbation of permittivity resulting from slight changes in geometry of the optimized region
• the gradient of the figure of merit (FOM): the gradient of the mode overlap to the fundamental TE mode of the waveguides as a result of the shape changes

After optimization, the optimized Y branch component shape will be exported into GDSII format which can be used for further simulation and/or fabrication (mask design).

Step 4: Run the parametric optimization in 3D [optional]

As an optional step, the best solution found by the 2.5D optimization can be used as an initial guess for a 3D optimization to verify the 2.5D simulation results and to further improve the performance of the design if possible. If one does not have a varFDTD license it is possible to start with the 3D optimization; however, without an initial solution, the optimization may progress slowly. We have included symmetric boundary conditions to reduce the sim time by a factor of four.

Run and results

Instructions for running the model and discussion of key results

Step 1: Define base simulation parameters

The resulting base simulation which includes the varFDTD simulation region, (non-dispersive) materials, source, monitors, mesh override for optimizable geometry, and input and output waveguides is shown in the figure below.

Step 2: Define the optimizable geometry

This step is not necessary if you want to run the example as is, however it is required for modification of the optimizable geometry based on your needs.

Below is a screenshot of the optimization python script showing the lines where you can define the initial parameters and their bound values. The second figure shows the optimizable geometry within the based simulation with the default initial parameters.

Step 3: Run the parametric optimization in 2.5D

A command prompt (in windows operating system) or terminal (in macOS or Linux) window will open showing the progress of the optimization. A screenshot of the command prompt window at the end of the optimization is shown below which contains information about the optimization progress.

During the optimization, a few new sessions of varFDTD will get opened by the optimization routine to perform simulations based on the guesses made by the optimization algorithm. These CAD windows can be hidden by specifying 'hide_fdtd_cad=true' in the optimizer. Another window will also report the optimization progress providing various plots for evolution of the optimization geometry, figure of merit, parameters and other related information. The figure below shows the optimization report window at the end of the optimization process.  The geometry of the final solution can be seen in the plot and it is shown that the FOM will have little improvement after about 25 iterations which suggests the number of iterations can be further reduced for future runs. The field gradients (due to the changes in geometry) and forward direction field distributions within the device’s structure can also be seen in the plot. In addition, it is clear that the changes in the optimized parameter values become less intense after a few iterations since the FOM gets closer to the optimum value. The changes in FOM gradient (based on changes in the parameter values) also decrease on average as seen in gradient evolution plot as the iteration continues.

At the end of optimization, the final parameter values will be saved in a text file (2D_parameters.txt) in the same folder as the example files which will be used in the 3D optimization for the next step. In addition, copies of the simulation files for each iteration, a text file containing the history of parameters and FOM values for each iteration and still images of the optimization report window for each iteration will be saved in a folder (with a name starting with “opts”) in the same folder as example files for future reference. Each optimization run creates a new folder with the same name but a numbered suffix to avoid overwriting previous runs data. Finally, we obtain a 2.5D design with FOM of about 0.9811.

Here is the geometry optimized in 2.5D:

The optimization python script file will also export the optimal design to a GDS file at the end of optimization. The GDS file will be saved in the same folder as the example files. The example files come with an encrypted script files with .lsfx extension that is needed for the GDS export operation. They should be saved in the same location as the python script file.

Step 4: Run the parametric optimization in 3D [optional]

Looking at the results of the 3D optimization will reveal that the Y-branch’s transmission has been predicted to be about 94% by 3D simulation (which is more realistic than the 2D case with the predicted value of 98%) and has been improved to over 95% at the end of optimization.

The following plot shows the final broadband transmission obtained from the optimized 3D structure between 1300 and 1800nm:

Important model settings

Description of important objects and settings used in this model

Definition of optimizable geometry : The optimizable geometry is defined as a polygon that has fixed values for the x coordinates of its boundary points and the points’ y coordinates can be modified to obtain the optimal geometry. These y coordinates are defined as the nodes of a cubic spline (as shown in the figure in the appendix section), which makes it possible to perform cubic spline interpolation on the points in order to obtain a smooth polygon geometry.

Optimization fields monitor position : The optimization fields DFT monitor (opt_field) is used to collect the field data within the optimizable geometry which is then used to calculate field gradients used in optimization algorithm. As such the location of this monitor is of great importance and should cover the entire optimizable geometry.

Figure of merit fields monitor position : Since this monitor is used to calculate the figure of merit required for optimization (mode overlap to the fundamental TE mode of the waveguides), it should be located within the output waveguides of the Y-branch.

Dimensions of optimizable geometry vs base simulation : It is important to make sure that the dimensions of the optimizable geometry and the input output waveguides setup in the base simulation match so that the Y-branch geometry is properly formed.

Height of the geometry objects : The height (depth) of the geometry objects should match in both base simulation (input, output waveguides) and the optimizable geometry (polygon). This is especially important for the 3D simulation to ensure the correct geometry setup. For the waveguides, this is set through the initialization method (by setting their z span) and for the polygon, it is adjusted as the input argument “depth” when calling the splitter function to form the optimizable geometry.

Optimization parameters bounds and initial values : The range over which the optimization algorithm is allowed to vary the parameters and also the initial values for those parameters as a starting point are defined as shown in the lines below:

initial_points_x = np.linspace(-1.0e-6, 1.0e-6, 10)
initial_points_y = np.linspace(0.25e-6, 0.6e-6, initial_points_x.size)
bounds = [(0.2e-6, 0.8e-6)] * initial_points_y.size

Simulation bandwidth : The inverse design PO method supports broadband simulations. The bandwidth is defined using start and stop wavelengths complemented by the number of frequency points in the given interval. The wavelength specification is passed to the optimization script when defining the FOM. Therefore, the wavelength settings defined in the base simulation script will be overridden.

Additionally, each frequency point can have its target FOM specified independently. This means that the optimization routine can aim to match the custom FOM profile as a function of wavelength rather than simply maximizing the FOM across the full bandwidth.

Here is a code section from this Y-branch example that specifies the bandwidth from 1300 to 1800nm with 21 frequency points:

wavelengths = Wavelengths(start = 1300e-9, stop = 1800e-9, points = 21)

The bandwidth settings is then passed to the optimization routine during FOM definition. Notice that we specified FOM target of 1 for all frequency points to maximize the transmission across the desired bandwidth:

fom = ModeMatch(monitor_name = 'fom',
                mode_number = 'fundamental mode',
                direction = 'Forward',
                target_T_fwd = lambda wl: np.ones(wl.size),
                norm_p = 1)

Maximum number of iterations : While the algorithm has the capability to stop the optimization once the gradient of the figure of merit falls below a certain threshold, the “max_iter” variable used for defining the optimization algorithm can be used to limit the number of iterations the algorithm is able to perform.

Symmetry : As of 2019a r4 the inverse design PO method supports symmetric boundaries. This can be particularly useful for 3D simulations where it can decrease the simulation time by a factor of 4.

Updating the model with your parameters

Instructions for updating the model based on your device parameters

Geometry : If you need to have your own geometry defined for the Y-branch, including SOI device layer thickness (waveguide height), output waveguides spacing, device length (spacing between input and output waveguides) and waveguide width, corresponding changes should be made in base simulation setup script and/or optimization setup python script. This might need changes in objects span and location in the base script, the polygon’s “depth” parameter in Python script which its value is also passed to the GDS export portion of the script. The source, simulation region, mesh override, and field gradient and FOM monitors dimensions should also be adjusted accordingly to make sure they cover the entire structure properly. The coordinates of the initial values and bounds of the optimizable geometry at points where it is connected to the output waveguide might need to be adjusted to ensure a proper connection between the optimizable region and the waveguides. In addition, the range of x coordinates of the optimizable geometry might need adjustment in the optimization setup python script where the splitter function is defined as shown below (points_x variable) to ensure its proper connection to the waveguides:

def splitter(params):
    points_x = np.concatenate(([initial_points_x.min() - 0.01e-6], initial_points_x, [initial_points_x.max() + 0.01e-6]))
    points_y = np.concatenate(([initial_points_y.min()], params, [initial_points_y.max()]))
    n_interpolation_points = 100
    polygon_points_x = np.linspace(min(points_x), max(points_x), n_interpolation_points)
    interpolator = sp.interpolate.interp1d(points_x, points_y, kind = 'cubic')
    polygon_points_y = interpolator(polygon_points_x)
    polygon_points_up = [(x, y) for x, y in zip(polygon_points_x, polygon_points_y)]
    polygon_points_down = [(x, -y) for x, y in zip(polygon_points_x, polygon_points_y)]
    polygon_points = np.array(polygon_points_up[::-1] + polygon_points_down)
    return polygon_points

It is recommended to set the coordinates such that there is a slight overlap between the optimizable region and the edges of the waveguides to ensure they are connected in the simulation.

Materials : The materials included in the simulation (the material making the optimizable geometry and the material surrounding it) are defined in Step 1 via lsf, project file or python function. In this example, we have taken the index of silicon and silicon dioxide at the desired simulation wavelength of 1550nm. One should make sure that the material defined (or given index) in the initialization is the same as that passed to the optimizer when defining the geometry as shown below;

eps_in = Material(name = 'Si: non-dispersive', mesh_order = 2)
eps_out = Material(name = 'SiO2: non-dispersive', mesh_order = 3)
depth = 220.0e-9
polygon = FunctionDefinedPolygon(func = splitter,
    initial_params = initial_points_y,
    bounds = bounds,
    z = 0.0,
      depth = depth,
    eps_out = eps_out,
     eps_in = eps_in,
    edge_precision = 5,
     dx = 1.0e-9)

With varFDTD it is now possible to use the material library to define the material properties; however, the optimization technique has not been tested with dispersive materials; thus, the effective index is solved at the center frequency and passed to the optimizer.

Taking the model further

Information and tips for users that want to further customize the model

This example provides predefined polygon geometry that is already parametrized and later optimized to create the final Y branch shape. While the polygon geometry allows some basic customization such as the boundaries in the y-direction, an advanced user can modify the python function “splitter” to achieve customize the properties of the parametrized geometry (sometimes called pCell).

Number of spline nodes : The number of optimized parameters corresponds to the number of spline nodes that control the geometry shape. A higher number of nodes will allow more nuanced control of the overall shape with a smaller minimum curve radius. On the other hand, the manufacturing process could limit the minimum feature size and a large number of parameters can also translate into longer optimization time. Changing the number of spline nodes requires modifications in the splitter function as well as in the bounds array size and initial parameters passed to the optimization routine.

Custom pCell : There is a virtually unlimited number of possibilities for how to translate an array of scalar parameters into a polygon shape. Therefore, advanced Python users can utilize the robustness of the language and write the splitter function from scratch. The only limitation is that the input is a scalar array and the output is a set of polygon points. The example below shows an example function that parametrizes not only the body of the Y branch, but also the output waveguides:

def splitter(params = np.linspace(0.25e-6, 2e-6, 20)):
    ''' Defines a taper where the paramaters are the y coordinates of the nodes of a cubic spline. '''
    points_x = np.concatenate(([-2.51e-6], np.linspace(-2.5e-6,2.5e-6,20), [2.51e-6]))
    points_y = np.concatenate(([0.25e-6], params, [2e-6]))
    n_interpolation_points = 100
    px = np.linspace(min(points_x), max(points_x), n_interpolation_points)
    interpolator = sp.interpolate.interp1d(points_x, points_y, kind = 'cubic')
    py = interpolator(px)
    py = np.minimum(2.5e-6, py)
    py = np.maximum(np.concatenate((np.ones(50)*0.2e-6, np.ones(50)*0.53e-6)), py)
    px = np.concatenate((px, px[40::][::-1]))
    py = np.concatenate((py, py[40::][::-1]-0.5e-6))
    polygon_points_up = [(x, y) for x, y in zip(px, py)]
    polygon_points_down = [(x, -y) for x, y in zip(px, py)]
    polygon_points = np.array(polygon_points_up[::-1] + polygon_points_down)
    return polygon_points

Video demonstrating the optimization progress of the modified pCell:

Additional resources

Additional documentation, examples and training material

See also

• Python API
• Y Branch PSO
• Y Branch S parameter extraction
• Symmetric boundary conditions
• Inverse design of waveguide crossing
• Inverse design of grating coupler

Related Ansys Innovation Courses

Appendix: Parametric Optimization (Adjoint method)

Additional background information and theory

In a typical device optimization, the goal is to minimize or maximize a figure or merit \( F(p) \) that depends on a set of \( N \) design parameters \( p=(p_1,...,p_n) \). The figure of merit can be any quantity used to benchmark the device’s performance and the design parameters can be any quantity that alters the device’s response. To evaluate \( F(p) \) for a given set of values of \( p \), a physics based modeling tool - such as FDTD - is used to map values of \( p \) to values of \( F \). The model usually contains a description of the device’s geometry, materials and all relevant environmental conditions. In a parameterized shape optimization, the goal is to allow some parts of the device’s geometry to vary to minimize or maximize the figure of merit. For this purpose, a shape function \( f(p) \) is typically defined to specify part of the device’s geometry; an example is shown in the following figure:

Y-splitter with a parameterized junction boundary $f(p)$ defined using connected splines.

The shape function is typically user-defined and must be fed to the modeling tool to construct the device’s geometry. Once the figure of merit \( F(f(p)) \) can be evaluated for a given shape function, a minimizer can make multiple evaluations of the figure of merit to find the optimal shape. There is a very large body of minimization tools and techniques that can be used to try to find the optimal shape parameters. When selecting a minimization technique, it is important to consider that running a device simulation to evaluate the figure of merit can be a slow operation. So, the best-suited minimization techniques are those that limit the number of evaluations of the figure of merit. With a good starting point, gradient-based methods can quickly find the optimal parameters with a relatively small number of figure of merit evaluations provided that it is also possible to evaluate the gradient of the figure of merit:

$$∇_pF=(\frac{\partial F}{\partial p_1},⋯,\frac{\partial F}{\partial p_n})$$

To evaluate the gradient without having to perform one simulation for each partial derivative, a technique called the adjoint method can be employed. This technique exploits the properties of linear partial differential equations to evaluate the gradient using only two simulations independently of the number of optimization parameters. The Python-based shape optimization module in FDTD (lumopt) allows users to define a shape function and run the corresponding device simulation to evaluate the figure of merit. The module also implements the adjoint method to provide the gradient. In this way, any of the gradient-based minimizers in SciPy’s optimize module can be used to perform a fully automated optimization.