This section covers the details of shape and topology based inverse design using lumopt. As prerequisite please familiarize yourself with photonic inverse design and understand the basics of the lumopt workflow.
The following examples can help provide context to the information presented here.
Shape
Shape based optimization allows users to take strong existing designs and improve their performance. Constraints can be easily included in the parameterization such as design rules for manufacturing. To do this the user needs to parameterize their geometry, through a shape function \( s(p) \) and provide bounds on the \( p \) values. This is supported with a handful of geometry classes FunctionDefinedPolygon which accepts a python function of parameters, and ParamterizedGeometry which accepts a python API function. FDTD is used to calculate the \( FOM( s(p) ) \).
Figure 1: Shape parameters and the shape function
The Adjoint method motivation section described how the parameters enter the equations through the permittivity. When users define a parameterized geometry an additional step, we refer to as the d_eps_calculation, is required to determine how varying each parameter changes the permittivity in Maxwell’s equations and the subsequent shape derivative. This is accomplished by meshing the simulation and returning the \( \epsilon (x,y,z) \), but without running the solver. Since this must be done for each parameter, including more parameters comes at some extra cost; although, this should be small compared to running the simulation; however, in large 3D simulations this meshing step can be intensive even if it is brief. If you need assistance improving the performance, and parallelization of the d_eps calculation please reach out to Lumerical support. For more information on the shape derivative approximation method see Chpt 5 of Owen Miller's thesis.
It is important to consider that the choice of parameters may be too restrictive, and that gradient based techniques are inherently dependent on initial conditions. It not always obvious what a better choice of parameters and initial conditions may be so some experimentation may be required to achieve the greatest performance.
Shape Geometry Classes
For each geometry class the user must pass the shape function, parameter bounds, geometry specifications and material properties. It is important to verify that your geometry varies smoothly and does not self-intersect for any combination of acceptable parameters as this can produce discontinuities in the gradient calculation. Spline points are not bounded in the same way the parameters are which can be problematic. When end points are clamped this may lead to discontinuities, employing "not-a-knot spline" usually avoids these errors. By visualizing different shapes realized from possible sets of bounded parameters most of these errors can be recognized and avoided.
Functiondefinedpolygon
This python class takes a user defined function of the optimization parameters. The function defines a polygon by returning a set of vertices, as a numpy array of coordinate pairs np.array([(x0,y0),…,(xn,yn)]), where the vertices must be ordered in a counter-clockwise direction.
from lumopt.geometries.polygon import FunctionDefinedPolygon
class FunctionDefinedPolygon(func,
initial_params,
bounds,
z,
depth,
eps_out,
eps_in,
edge_precision,
dx,
deps_num_threads)
: func : function
function(parameters) - Returns \( s(p) \) an np.array([(x0,y0),…,(xn,yn)]) of coordinate tuples ordered CCW.
Typically, a subset of polygon points \( p \) will be used as parameters with a spline technique to up-sample these values and create a smoother curve. This user defined function is a common source of error.
: initial_params: np.array
Parameters in a form accepted by the func, that define the initial geometry. Size = [1, n_param]
: bounds: np.array
A set of tuples that define the min and max of each parameter. Size = [2, n_param]
: z : float
Center of polygon along the z-axis. Default = 0.0
: depth: float
Span of polygon along the z-axis.
: eps_out: float
permittivity of the material around the polygon.
: eps_in: float
Permittivity of the polygon material.
: edge_precision: int
Number of quadrature points along each edge for computing the FOM gradient using the shape derivative approximation method.
: dx: float > 0.0
Step size for computing the FOM gradient using permittivity perturbations. Default = 1e-5
: deps_num_threads: int, optional
Number of threads to be uses in the d_eps calculation. This utilizes the job manager, and changes your resource configuration. By default this will be performed serially, but it is possible to utilize more cores simultaneously by setting this greater than 1. This is an advanced, experimental feature with some some additional overhead associated; thus, testing and benchmarking is required. In the event of an unexpected termination one will have to manually fix the update. Resource configuration elements and controls. Default = 1
ParameterizedGeometry
from lumopt.geometries.parameterized_geometry import ParameterizedGeometry
class ParameterizedGeometry(func,
initial_params,
bounds,
dx,
deps_num_threads,
obj_bounds)
: func : python API function
function(parameters, fdtd, only_update, (optional arguments))
This function is less concise, but more intuitive and flexible than FunctionDefinePolygon. It allows users to call the API methods, in the same way you would use lsf script commands to set-up and update the geometry rather than passing the (x,y) vertices. The flag 'only_update' is used to avoid frequent recreations of the parameterized geometry which reduces the performance. When the flag is true, it is assumed that the geometry was already added at least once to the CAD. For an example of it's usage see Inverse design of grating coupler.
: initial_params: np.array
Parameters in a form accepted by the func, that define the initial geometry. Size = [1, n_param]
: bounds: np.array
A set of tuples that define the min and max of each parameter. Size = [2, n_param]
: dx: float > 0.0
Step size for computing the FOM gradient using permittivity perturbations. Default = 1e-5
: deps_num_threads: int, optional
Number of threads to be uses in the d_eps calculation. This utilizes the job manager, and changes your resource configuration. By default this will be performed serially, but it is possible to utilize more cores simultaneously by setting this greater than 1. This is an advanced, experimental feature with some some additional overhead associated; thus, testing and benchmarking is required. In the event of an unexpected termination one will have to manually fix the update. Resource configuration elements and controls. Default = 1
: obj_bounds: python API function, optional
When calculating the gradient due to the geometry change, the software would re-mesh the structure to identify the geometry variation in one iteration. In some cases, only a small region would be changed. This function could be applied to customize a specific region so that re-meshing is only focusing on this small region, speeding up the re-mesh time consequently. The result of this function should be a np.array like [[xmin,xmax],[ymin,ymax]], which defines the bounding box of the re-mesh region. Default = none
Shape Considerations
From 2D to 3D
Shape optimization in 2D is supported through MODE varFDTD and FDTD. Using varFDTD will automatically collapse the material and geometric contributions in the z-direction to more accurately reflect the effective indices. FDTD simply uses the material index, ignoring the 3rd dimension. For in-plane integrated photonic components varFDTD provides the most straightforward route to accurate permittivity values. One can use FDTD in this situation as well, but determining the effective indices is up to the user. In situations where the propagation is not entirely in plane, such as with the grating coupler example, FDTD is required. To use varFDTD set the use_varFDTD = True in the optimization class.
In shape optimization it is straightforward to move from 2D to 3D since all the supported geometry classes are transferable. The only update required is to change the base file monitors and simulation region from 2D to 3D FDTD. Bulk material indices should be specified for eps_in and eps_out as the geometric contributions are automatically incorporated into the update equations. We have found strong correlation between high performing devices in 2D and device performance in 3D; however, expect the FOM to be reduced when moving to 3D as new loss mechanisms are considered. Due to the complexity and time required for 3D optimizations we recommend starting in 2D and then importing those optimized structures into 3D for a final abbreviated optimization run.
GDS Export
When exporting the optimized geometry we recommend using the encrypted GDS export script, outlined and available for download here:
Topology
In topology optimization the user does not need to parameterize the geometry, which simplifies the workflow; however, the simulation needs to perform a few more complicated steps to produce physically realizable devices. To set up topology optimization simply define the optimizable footprint and materials in the relevant 2D or 3D geometry class, outline below. Each FDTD mesh cell in the discretized optimization volume becomes a parameter. This technique often scales to thousands of parameters, but with the adjoint method we still only need 2 simulations per iteration to calculate the gradient. Furthermore, since the permittivity values are the optimization parameters, there is a direct link to Maxwell’s equations so there is no need to perform a d_eps as calculation.
Figure 2: Topology optimization schematic
In the initialization section we describe important aspects of the setup, and how to provide initial conditions. Once the optimization is running it proceeds in two or three steps. First the greyscale phase where the parameters vary continuously between the core and cladding index. Next, a binarization phase is run to force the mesh cells to take either the core or cladding index. Finally, an optional design for manufacturing step is implemented. In this step a penalty function is added to the FOM. The magnitude of the penalty function is determined by the degree to which the minimum feature size constraint is violated. Finally, we will discuss how to move from 2D to 3D, and export your designs to GDS .
Topology Geometry Classes
TopologyOptimization2D
This class is used in 2D simulations.
from lumopt.geometries.topology import TopologyOptimization2D
class TopologyOptimization2D(params,
eps_min,
eps_max,
x, y, z,
filter_R,
eta,
beta,
min_feature_size)
: params : np.array
Initial parameters that define the initial geometry. Size = [n_x, n_y]
: eps_min : float
Permittivity of the core material
: eps_max: float
Permittivity of the cladding
: x: np.array
The x position of the parameters and FDTD unit cells in [m].
: y : np.array
The y positions of the parameters and FDTD unit cells in [m].
: z : float
Z position of the topology object in [m]. Default 0.0
: filter_R: float
The radius of the Heaviside filter in [m]. Outlined in greyscale. Default = 200e-9
: eta : float
Eta parameter of the Heaviside filter. Provides a shift of the binarization threshold towards eps_min if eta <0.5, and towards eps_max if eta >0.5. Default = 0.5
: beta: float
Beta parameter of the Heaviside filter. This is the starting value used, but the beta parameter is ramped up to ~1000 during the binarization phase. Default = 1
: min_feature_size: float, optional
Minimum feature size in the DFM phase. Features smaller than this value will be penalized in the final DFM stage. Must be filter_R < min_feature_size < 2.0 * filter_R. A value of 0.0 will disable the DFM phase. Default = 0.0
TopologyOptimization3DLayered
This class is used in 3D simulations, and operates in the same way as the 2D geometry with the exception that z is now a required nump array of points, and the optimization region is a volume. To reduce the field results to the XY plane there is an integration along z. The pattern of the layer is optimized for, but not each mesh cell individually: ie it is not a full 3D voxelated topological geometry, that feature is currently not supported.
from lumopt.geometries.topology import TopologyOptimization3DLayered
class TopologyOptimization3DLayered(params,
eps_min,
eps_max,
x, y, z,
filter_R,
eta,
beta,
min_feature_size)
: params : np.array
Initial parameters that define the initial geometry. Size = [n_x, n_y, n_z]
: eps_min : float
Permittivity of the cladding
: eps_max: float
Permittivity of the core material
: x: np.array
The x position of the parameters and FDTD unit cells in [m].
: y : np.array
The y positions of the parameters and FDTD unit cells in [m].
: z : np.array
The z position of the parameters and FDTD unit cells in [m].
: filter_R: float
The radius of the Heaviside filter in [m]. Outlined in greyscale. Default = 200e-9
: eta : float
Eta parameter of the Heaviside filter. Provides a shift of the binarization threshold towards eps_min if eta <0.5, and towards eps_max if eta >0.5. Default = 0.5
: beta: float
Beta parameter of the Heaviside filter. This is the starting value used, but the beta parameter is ramped up to ~1000 during the binarization phase. Default = 1
: min_feature_size: float, optional
Minimum feature size in the DFM phase. Features smaller than this value will be penalized in the final DFM stage. Must be filter_R < min_feature_size < 2.0 * filter_R. A value of 0.0 will disable the DFM phase. Default = 0.0
Running
This is section describes the various stages of topology optimization.
Figure 3: Topology optimization steps
Initialization
It is important the mesh defined in the FDTD simulation precisely matches the optimization parameters defined in Python. To do this make sure a uniform mesh is defined over the optimizable footprint. The python parameters may be defined as follows.
#Footprint and pixel size
size_x = 3000 #< Length of the device (in nm). Longer devices typically lead to better performance
delta_x = 25 #< Size of a pixel along x-axis (in nm)
size_y = 3000 #<Extent along the y-axis (in nm). When using symmetry this should half the size.
delta_y = 25 #< Size of a pixel along y-axis (in nm)
size_z = 240 #< Size of device in z direction
delta_z = 40 #< Size of a pixel along z-axis (in nm)
x_points=int(size_x/delta_x)+1
y_points=int(size_y/delta_y)+1
z_points=int(size_z/delta_z)+1
#Convert to SI position vectors
x_pos = np.linspace(-size_x/2,size_x/2,x_points)*1e-9
y_pos = np.linspace(-size_y/2,size_y/2,y_points)*1e-9
z_pos = np.linspace(-size_z/2,size_z/2,z_points)*1e-9
In the lumerical script that generates the base file ensure that the grid matches. Alternative methods such as a python function or base file are also allowed, but the grid and parameter values will also need to match.
#Define footprint
opt_size_x=3e-6;
opt_size_y=3e-6;
opt_size_z=240e-9;
#Define pixel size
dx = 25e-9;
dy = 25e-9;
dz = 40e-9;
addmesh;
set('name','opt_fields_mesh');
set('override x mesh',true); set('dx',dx);
set('override y mesh',true); set('dy',dy);
set('override z mesh',true); set('dy',dz);
set('x',0); set('x span',opt_size_x);
set('y',0); set('y span',opt_size_y);
set('z',0); set('z span',opt_size_z);
To set the initial conditions you can pass an np.array the same size as the parameters. If you want start with a uniform permittivity value across the optimizable domain this is straightforward.
##Initial conditions
#initial_cond = np.ones((x_points,y_points)) #< Start with the domain filled with eps_max
#initial_cond = 0.5*np.ones((x_points,y_points)) #< Start with the domain filled with (eps_max+eps_min)/2
#initial_cond = np.zeros((x_points,y_points)) #< Start with the domain filled with eps_min
Alternatively, define an initial structure in the base simulation, and pass None. In this case the optimizer meshes the simulation to determine the initial parameter set. This should be a structure group containing the relevant geometry called initial_guess.
Figure 4: Topology optimization set-up
Greyscale
In the python optimization each cell corresponds to a parameter \( \rho \in [0,1] \) which is linearly mapped to the \( \epsilon \) values in FDTD. A circular top hat convolution with a user specified radius is used to remove sharp corner, small holes or islands which can’t be manufactured through convolution. Depending on you manufacturing process, typical values for the filter radius will be between 50 and 500nm.
Figure 5: Top-hat convolution diagram
In this stage of the optimization the FOM error tends to very quickly reduce but may take 100’s of iterations to reach a minima. The values obtained during greyscale are typically very good, but the filtering process is imperfect, and there will be intermediate eps values and unproducible features via lithographic means. The following steps are used to produce realizable devices.
Binarization
During binarization the Heaviside filter changes the mapping of the python elements to FDTD permittivity. In the greyscale phase the beta parameter is 1, producing a linear mapping. During binarization the beta parameter is incrementally ramped up. The Heaviside filter is defined as follows, and ramping up the beta parameter converts it from a linear mapping to a step function.
$$\epsilon_{i}=\epsilon_{l o w}+\frac{\tanh (\beta \eta)+\tanh \left(\beta\left(\rho_{i}-\eta\right)\right)}{\tanh (\beta \eta)+\tanh (\beta(1-\eta))} \Delta \epsilon$$
Figure 6: Heaviside filter
Ramping \( \beta \) up will perturbs the system producing sharp peaks in the FOM. It is followed by a set of optimization iterations to adjust to the perturbation, and reduce the FOM error. It is possible to change the number of iterations between increments of Beta, by changing the opt.continuation_max parameter.
Figure 7: FOM during optimization steps
Design For Manufacturing (optional)
Even when using the Heaviside filter for spatial filtering, the topology optimization designs tend to create features that break process design rules. To ensure manufacturability we explored ways of explicitly enforcing specific minimum feature size constraints. Our implementation uses a technique published by Zhou et al. in 2015[1]. This paper was focused structural mechanics, so we had to make some small adjustments to make it work reliably for photonic devices.
The main idea is to come up with an indicator function which is minimal if the minimum feature size constraints are fulfilled everywhere. This term is then added as a penalty term to the FOM during optimization. To implement this phase pass a non-zero min_feature_size constraint to the geometry class.
Topology Considerations
From 2D to 3D
All dimensional considerations from shape apply to topology. We recommend using effective indices as calculated in varFDTD as the permittivity values in 2D opt, but one should work with FDTD for all topology optimizations.
When going from 2D to 3D in topology the performance does not correlate as well, as in shape based designs. That being said good 2D results demonstrate the possibility of producing strong 3D designs which is important. Performance improvements in 3D tend to proceed more quickly if one imports a design that is not fully binarized, which relaxes the abillity of the 3D optimizer to vary parameters. Importing partially binarized 2D designs into 3D is outlined below.
## First, we specify in which path to find the 2d optimization
path2d = os.path.join(cur_path,"[Name of 2D opt folder]")
## Next, check the log-file to figure out which iteration we should load from.
# It is usually best to load something mid-way through the binarization
convergence_data = np.genfromtxt(os.path.join(path2d, 'convergence_report.txt'), delimiter=',')
## Find the first row where the beta value (3rd column) is larger than beta_threshold (here 10)
beta_threshold = 10
beta_filter = convergence_data[:,2]>beta_threshold
convergence_data = convergence_data[beta_filter,:] #< Trim data down
iteration_number = int(convergence_data[0,0]-2)
## Load the 2d parameter file for the specified iteration
geom2d = TopologyOptimization2D.from_file(os.path.join(path2d, 'parameters_{}.npz'.format(iteration_number) ), filter_R=filter_R, eta=0.5)
startingParams = geom2d.last_params #< Use the loaded parameters as starting parameters for 3d
startingBeta = geom2d.beta #< Also start with the same beta value. One could start with 1 as well.
Export to GDS
Since the topology designed devices do not have a well defined shape we developed an alternate method to extract the GDSII designs using contours. This functionality is outlined in the following KX page.
Further Resources
Webinars
- Feb 2019 - Broadband shape
- Aug 2019 - Updates to Shape, Grating Couplers, and Intro to Topology Optimization
- April 2020 - Update to Topology Optimization
Examples
- GDS Pattern Extraction Using getcontours
- Inverse Design of an Efficient Narrowband Grating Coupler
- Broadband Grating Coupler Design Bandwidth vs Peak Efficiency
- Inverse Design of a Y-Splitter using Topology Optimization
- Topology Optimization of an O,C Band Splitter with min_feature_size Constraint