GDS Import#
This tutorial demonstrates how to set up a simulation based on importing a GDS file. There are two examples: (1) computing the S-parameters of a four-port network using a silicon directional coupler and (2) finding the modes of a ring resonator. These two component devices are used in photonic integrated circuits to split/combine and filter an input signal. For more information on directional couplers and ring resonators, see Section 4.1 of Silicon Photonics Design by Chrostowski and Hochberg.
The GDS layout files are read using the gdstk Python package, which must be installed separately (e.g., pip install gdsktk). gdstk is a library for creating and reading GDSII/OASIS files. The polygons it extracts from a given layer are converted into Meep Prism objects (the device geometry) and Volume objects (for source and flux/mode-monitor regions). The two examples below define a few small helper functions for this conversion which can be reused in your own scripts:
import gdstk
import meep as mp
def get_gds_cell(fname):
"""Returns the (single) top-level cell of the GDS file `fname`."""
return gdstk.read_gds(fname).top_level()[0]
def get_gds_prisms(material, cell, layer, datatype=0, zmin=0.0, zmax=0.0):
"""Returns a list of `mp.Prism`s, one for each polygon on (`layer`, `datatype`)."""
prisms = []
for poly in cell.get_polygons(layer=layer, datatype=datatype):
vertices = [mp.Vector3(x, y, zmin) for x, y in poly.points]
prisms.append(
mp.Prism(
vertices,
height=zmax - zmin,
axis=mp.Vector3(0, 0, 1),
material=material,
)
)
return prisms
def get_gds_vol(cell, layer, datatype=0, zmin=0.0, zmax=0.0):
"""Returns an `mp.Volume` spanning the bounding box of (`layer`, `datatype`)."""
polygons = cell.get_polygons(layer=layer, datatype=datatype)
xs = [x for poly in polygons for x, y in poly.points]
ys = [y for poly in polygons for x, y in poly.points]
xmin, xmax = min(xs), max(xs)
ymin, ymax = min(ys), max(ys)
center = mp.Vector3(0.5 * (xmin + xmax), 0.5 * (ymin + ymax), 0.5 * (zmin + zmax))
size = mp.Vector3(xmax - xmin, ymax - ymin, zmax - zmin)
dims = 2 if (zmin == 0 and zmax == 0) else 3
return mp.Volume(center=center, size=size, dims=dims)
S-Parameters of a Directional Coupler#
The directional coupler as well as the source and mode monitor geometries are described by the GDS file examples/coupler.gds. A snapshot of this file viewed using KLayout is shown below. The figure labels have been added in post processing. The design consists of two identical strip waveguides which are positioned close together via an adiabatic taper such that their modes couple evanescently. There is a source (labelled "Source") and four mode monitors (labelled "Port 1", "Port 2", etc.). The input pulse from Port 1 is split in two and exits through Ports 3 and 4. The design objective is to find the separation distance which maximizes the outgoing power in Port 4 at a wavelength of 1.55 μm. More generally, though not included in this example, it is possible to have two additional degrees of freedom: (1) the length of the straight waveguide section where the two waveguides are coupled and (2) the length of the tapered section (the taper profile is described by a hyperbolic tangent (tanh) function).
The GDS file is adapted from the SiEPIC EBeam PDK with four major modifications:
-
the computational cell is centered at the origin of the \(xy\) plane and defined on layer 0
-
the source and four mode monitors are defined on layers 1-5
-
the lower and upper branches of the coupler are defined on layers 31 and 32
-
the straight waveguide sections are perfectly linear
Note that rather than being specified as part of the GDS file, the volume regions of the source and flux monitors could have been specified in the simulation script.
The simulation script is in examples/coupler.py.
import argparse
import gdstk
import meep as mp
gds_file = 'coupler.gds'
CELL_LAYER = 0
PORT1_LAYER = 1
PORT2_LAYER = 2
PORT3_LAYER = 3
PORT4_LAYER = 4
SOURCE_LAYER = 5
UPPER_BRANCH_LAYER = 31
LOWER_BRANCH_LAYER = 32
default_d = 0.3
t_oxide = 1.0
t_Si = 0.22
t_air = 0.78
dpml = 1
cell_thickness = dpml + t_oxide + t_Si + t_air + dpml
oxide = mp.Medium(epsilon=2.25)
silicon = mp.Medium(epsilon=12)
fcen = 1 / 1.55
df = 0.2 * fcen
def get_gds_cell(fname):
"""Returns the (single) top-level cell of the GDS file `fname`."""
return gdstk.read_gds(fname).top_level()[0]
def get_gds_prisms(material, cell, layer, datatype=0, zmin=0.0, zmax=0.0):
"""Returns a list of `mp.Prism`s, one for each polygon on (`layer`, `datatype`)."""
prisms = []
for poly in cell.get_polygons(layer=layer, datatype=datatype):
vertices = [mp.Vector3(x, y, zmin) for x, y in poly.points]
prisms.append(
mp.Prism(
vertices,
height=zmax - zmin,
axis=mp.Vector3(0, 0, 1),
material=material,
)
)
return prisms
def get_gds_vol(cell, layer, datatype=0, zmin=0.0, zmax=0.0):
"""Returns an `mp.Volume` spanning the bounding box of (`layer`, `datatype`)."""
polygons = cell.get_polygons(layer=layer, datatype=datatype)
xs = [x for poly in polygons for x, y in poly.points]
ys = [y for poly in polygons for x, y in poly.points]
xmin, xmax = min(xs), max(xs)
ymin, ymax = min(ys), max(ys)
center = mp.Vector3(0.5 * (xmin + xmax), 0.5 * (ymin + ymax), 0.5 * (zmin + zmax))
size = mp.Vector3(xmax - xmin, ymax - ymin, zmax - zmin)
dims = 2 if (zmin == 0 and zmax == 0) else 3
return mp.Volume(center=center, size=size, dims=dims)
def main(args):
cell_zmax = 0.5 * cell_thickness if args.three_d else 0
cell_zmin = -0.5 * cell_thickness if args.three_d else 0
si_zmax = 0.5 * t_Si if args.three_d else 10
si_zmin = -0.5 * t_Si if args.three_d else -10
# read cell size, volumes for source region and flux monitors,
# and coupler geometry from the GDS file using gdstk
gds_cell = get_gds_cell(gds_file)
upper_branch = get_gds_prisms(silicon, gds_cell, UPPER_BRANCH_LAYER, zmin=si_zmin, zmax=si_zmax)
lower_branch = get_gds_prisms(silicon, gds_cell, LOWER_BRANCH_LAYER, zmin=si_zmin, zmax=si_zmax)
cell = get_gds_vol(gds_cell, CELL_LAYER, zmin=cell_zmin, zmax=cell_zmax)
p1 = get_gds_vol(gds_cell, PORT1_LAYER, zmin=si_zmin, zmax=si_zmax)
p2 = get_gds_vol(gds_cell, PORT2_LAYER, zmin=si_zmin, zmax=si_zmax)
p3 = get_gds_vol(gds_cell, PORT3_LAYER, zmin=si_zmin, zmax=si_zmax)
p4 = get_gds_vol(gds_cell, PORT4_LAYER, zmin=si_zmin, zmax=si_zmax)
src_vol = get_gds_vol(gds_cell, SOURCE_LAYER, zmin=si_zmin, zmax=si_zmax)
# displace upper and lower branches of coupler (as well as source and flux regions)
if args.d != default_d:
delta_y = 0.5*(args.d - default_d)
delta = mp.Vector3(y=delta_y)
p1.center += delta
p2.center -= delta
p3.center += delta
p4.center -= delta
src_vol.center += delta
cell.size += 2*delta
for np in range(len(lower_branch)):
lower_branch[np].center -= delta
for nv in range(len(lower_branch[np].vertices)):
lower_branch[np].vertices[nv] -= delta
for np in range(len(upper_branch)):
upper_branch[np].center += delta
for nv in range(len(upper_branch[np].vertices)):
upper_branch[np].vertices[nv] += delta
geometry = upper_branch+lower_branch
if args.three_d:
oxide_center = mp.Vector3(z=-0.5 * t_oxide)
oxide_size = mp.Vector3(cell.size.x, cell.size.y, t_oxide)
oxide_layer = [mp.Block(material=oxide, center=oxide_center, size=oxide_size)]
geometry = geometry + oxide_layer
sources = [
mp.EigenModeSource(
src=mp.GaussianSource(fcen,fwidth=df),
volume=src_vol,
eig_parity=mp.NO_PARITY if args.three_d else mp.EVEN_Y+mp.ODD_Z,
)
]
sim = mp.Simulation(
resolution=args.res,
cell_size=cell.size,
boundary_layers=[mp.PML(dpml)],
sources=sources,
geometry=geometry
)
mode1 = sim.add_mode_monitor(fcen, 0, 1, mp.ModeRegion(volume=p1))
mode2 = sim.add_mode_monitor(fcen, 0, 1, mp.ModeRegion(volume=p2))
mode3 = sim.add_mode_monitor(fcen, 0, 1, mp.ModeRegion(volume=p3))
mode4 = sim.add_mode_monitor(fcen, 0, 1, mp.ModeRegion(volume=p4))
sim.run(until_after_sources=100)
# S parameters
p1_coeff = sim.get_eigenmode_coefficients(
mode1, [1], eig_parity=mp.NO_PARITY if args.three_d else mp.EVEN_Y+mp.ODD_Z
).alpha[0,0,0]
p2_coeff = sim.get_eigenmode_coefficients(
mode2, [1], eig_parity=mp.NO_PARITY if args.three_d else mp.EVEN_Y+mp.ODD_Z
).alpha[0,0,1]
p3_coeff = sim.get_eigenmode_coefficients(
mode3, [1], eig_parity=mp.NO_PARITY if args.three_d else mp.EVEN_Y+mp.ODD_Z
).alpha[0,0,0]
p4_coeff = sim.get_eigenmode_coefficients(
mode4, [1], eig_parity=mp.NO_PARITY if args.three_d else mp.EVEN_Y+mp.ODD_Z
).alpha[0,0,0]
# transmittance
p2_trans = abs(p2_coeff) ** 2 / abs(p1_coeff) ** 2
p3_trans = abs(p3_coeff) ** 2 / abs(p1_coeff) ** 2
p4_trans = abs(p4_coeff) ** 2 / abs(p1_coeff) ** 2
print(
f"trans:, {args.d:.2f}, {p2_trans:.6f}, {p3_trans:.6f}, {p4_trans:.6f}"
)
if __name__ == '__main__':
parser = argparse.ArgumentParser()
parser.add_argument(
"-res", type=int, default=50, help="resolution (default: 50 pixels/um)"
)
parser.add_argument(
"-d", type=float, default=0.1, help="branch separation (default: 0.1 um)"
)
parser.add_argument(
"--three_d",
action="store_true",
default=False,
help="d calculation? (default: False)"
)
args = parser.parse_args()
main(args)
For a given waveguide separation distance (\(d\)), the simulation computes the transmittance of Ports 2, 3, and 4. The transmittance is the square of the S-parameter which itself is equivalent to the mode coefficient. There is an additional mode monitor at Port 1 to compute the input power from the adjacent eigenmode source; this is used for normalization when computing the transmittance. The eight layers of the GDS file are each converted to a Simulation object: the upper and lower branches of the coupler are defined as a collection of Prisms, the rectilinear regions of the source and flux monitor as a Volume and FluxRegion. The size of the cell in the \(y\) direction is dependent on \(d\). The default dimensionality is 2d. (Note that for a 2d cell the Prism objects returned by get_gds_prisms must have a finite height. The finite height of Volume objects returned by get_gds_vol are ignored in 2d.) An optional input parameter (three_d) converts the geometry to 3d by extruding the coupler geometry in the \(z\) direction and adding an oxide layer beneath similar to a silicon on insulator (SOI) substrate. A schematic of the coupler design in 3d generated using MayaVi is shown below.
The coupler properties are computed for a range of separation distances from 0.02 to 0.30 μm with increments of 0.02 μm from the shell command line:
for d in $(seq 0.02 0.02 0.30); do
mpirun -np 2 python coupler.py -d ${d} |tee -a directional_coupler.out;
done
grep trans: directional_coupler.out |cut -d , -f2- > directional_coupler.dat;
The transmittance results converted into insertion loss for Ports 3 and 4 are shown in the figure below. (There is essentially no flux into Port 2 and thus \(|S_{21}|^2\) is not shown in the figure.) When the two waveguide branches are sufficiently separated (\(d\) > 0.2 μm), practically all of the input power remains in the top branch and is transferred to Port 3. A small amount of the input power is lost due to scattering into radiative modes within the light cone in the tapered sections where the translational symmetry of the waveguide is broken. This is why the power in Port 3 never reaches exactly 100%. For separation distances of less than approximately 0.2 μm, evanescent coupling of the modes from the top to the lower branch begins to transfer some of the input power to Port 4. For \(d\) of 0.13 μm, the input signal is split evenly into Ports 3 and 4. For \(d\) of 0.06 μm, the input power is transferred completely to Port 4. Finally, for \(d\) of less than 0.06 μm, the evanescent coupling becomes rapidly ineffective and the signal again remains mostly in Port 3.
These quantitative results can also be verified qualitatively using the field profiles shown below for \(d\) of 0.06, 0.13, and 0.30 μm. To generate these images, the pulse source is replaced with a continuous wave (CW) and the fields are time stepped for a sufficiently long run time until they have reached steady state. The array slicing routines get_epsilon and get_efield_z are then used to obtain the dielectric and field data over the entire cell.
sources = [mp.EigenModeSource(src=mp.ContinuousSource(fcen,fwidth=df),
volume=src_vol,
eig_parity=mp.EVEN_Y+mp.ODD_Z)]
sim = mp.Simulation(resolution=res,
cell_size=cell.size,
boundary_layers=[mp.PML(dpml)],
sources=sources,
geometry=geometry)
sim.run(until=400) # arbitrary long run time to ensure that fields have reached steady state
eps_data = sim.get_epsilon()
ez_data = np.real(sim.get_efield_z())
import matplotlib.pyplot as plt
plt.figure()
plt.plot2D(fields=mp.Ez,
plot_sources_flag=False,
plot_monitors_flag=False,
plot_boundaries_flag=False)
plt.axis('off')
plt.show()
The field profiles confirm that for \(d\) of 0.06 μm (Figure 1), the input signal in Port 1 of the top branch is almost completely transferred to Port 4 of the bottom branch. For \(d\) of 0.13 μm (Figure 2), the input signal is split evenly between the two branches. Finally, for \(d\) of 0.30 μm (Figure 3), there is no longer any evanescent coupling and the signal remains completely in the top branch. The waveguide regions with no fields in Ports 3 and 4 are PML.
When computing the reflection coefficient |S11|2, is it necessary to perform a separate normalization run to obtain the incident fields?#
No (generally). In the single-run calculation of the reflection coefficent \(|S_{11}|^2\) which is based on the back-scattered fields in Port 1 (due to the finite taper/bend length which breaks translational symmetry) given the forward-propagating fields of an eigenmode source also in Port 1, slight discretization errors in the eigenmode-coefficient extraction (as described in paragraph 3 of Section 4.2.2 of this book chapter) will result in a "noise floor" below which the reflection cannot be measured in this way. This "noise floor" only applies at a fixed resolution — as the resolution is increased, the discretization error in the mode-coefficient calculation goes away, and \(|S_{11}|^2\) should approach the "true" reflection from the taper/bend. This is demonstrated in the figure below which shows a plot of the \(S_{11}\) and \(S_{21}\) reflectance versus resolution. (In these types of calculations, it is important that the source and mode monitor in the same port be separated by at least several pixels in order to avoid any overlap due to the grid discretization.)
In the limit of infinite resolution, the discretization error is removed and the reflectance for \(S_{11}\) and \(S_{21}\) converge to their "true" values of ~10-6 and ~10-8, respectively. (Note that the back-scattered fields in Port 2 are two orders of magnitude smaller than those in Port 1 because the input fields in the upper branch of the directional coupler must cross into the lower branch to reach Port 2.) In this example, \(|S_{21}|^2\) requires a resolution of at least ~150 to minimize discretization errors. The discretization errors due to the eigenmode-coefficient extraction can be greatly reduced by using a separate normalization run to compute the incident fields for just a straight waveguide (i.e., no taper/bend) which are then subtracted from the Fourier-transformed fields in Port 1 and 2 of the directional coupler. This procedure is similar to those involving flux calculations. (Alternatively, for single-mode waveguides, the mode-coefficient calculation can be replaced entirely with just computing the Poynting flux in the ports. This approach is more accurate at lower resolutions.) For practical applications, however, reflectance values less than 40 dB (e.g., for telecom multi-path interference tolerances) are often considered negligible. On the other hand, there may be theoretical investigations where trying to resolve such small reflections could be important. (As reflections approach 10-15, the limits of floating-point precision will eventually limit accuracy even for the normalization approach.)
Importing a GDS Layer using a Tuple#
In the directional coupler example above, individual layers of the GDS file were referenced by a single layer number (i.e., 1, 2, 31, 32, etc.) with an implicit data type of 0. However, there are certain GDS files in which a layer is referenced using a 2-tuple of a layer number and a data type (e.g., (37,4)). gdstk supports this natively: Cell.get_polygons accepts both a layer and a datatype argument. The example below uses the same helper functions defined at the top of this tutorial, passing the data type explicitly.
import gdstk
import meep as mp
# load the top-level cell of the GDS file
gds_cell = gdstk.read_gds(gds_file).top_level()[0]
# define cell size and center from the bounding box of the entire layout
(xmin, ymin), (xmax, ymax) = gds_cell.bounding_box()
cell_center = mp.Vector3(0.5 * (xmin + xmax), 0.5 * (ymin + ymax))
design_geometry = []
for poly in gds_cell.get_polygons(layer=37, datatype=4):
# define vertices relative to the center of the cell
vertices = [mp.Vector3(x, y) - cell_center for x, y in poly.points]
design_geometry.append(mp.Prism(vertices=vertices,
height=0.5,
axis=mp.Vector3(0, 0, +1),
material=mp.Medium(index=3.5)))
design_geometry.append(mp.Prism(vertices=vertices,
height=0.5,
axis=mp.Vector3(0, 0, -1),
material=mp.Medium(index=3.5)))
Note that for each polygon in the GDS layer, there are two Prism objects: one extending in the \(+z\) direction and the other in \(-z\) with a combined height of 1.0 μm.
Modes of a Ring Resonator#
The next example is similar to Tutorial/Basics/Modes of a Ring Resonator and consists of two parts: (1) creating the ring resonator geometry using gdstk and (2) finding its modes using Harminv. The cell, geometry, source, and monitor are defined on separate layers within the same GDS file. This example demonstates that gdstk can be used to create a GDS file and to read it back in for the Meep simulation.
The simulation script is in examples/ring_gds.py.
import gdstk
from matplotlib import pyplot as plt
import meep as mp
import numpy as np
# core and cladding materials
Si = mp.Medium(index=3.4)
SiO2 = mp.Medium(index=1.4)
# layer numbers for GDS file
RING_LAYER = 0
SOURCE0_LAYER = 1
SOURCE1_LAYER = 2
MONITOR_LAYER = 3
SIMULATION_LAYER = 4
resolution = 50 # pixels/μm
dpml = 1 # thickness of PML
zmin = 0 # minimum z value of simulation domain (0 for 2D)
zmax = 0 # maximum z value of simulation domain (0 for 2D)
def get_gds_cell(fname):
"""Returns the (single) top-level cell of the GDS file `fname`."""
return gdstk.read_gds(fname).top_level()[0]
def get_gds_prisms(material, cell, layer, datatype=0, zmin=0.0, zmax=0.0):
"""Returns a list of `mp.Prism`s, one for each polygon on (`layer`, `datatype`)."""
prisms = []
for poly in cell.get_polygons(layer=layer, datatype=datatype):
vertices = [mp.Vector3(x, y, zmin) for x, y in poly.points]
prisms.append(
mp.Prism(
vertices,
height=zmax - zmin,
axis=mp.Vector3(0, 0, 1),
material=material,
)
)
return prisms
def get_gds_vol(cell, layer, datatype=0, zmin=0.0, zmax=0.0):
"""Returns an `mp.Volume` spanning the bounding box of (`layer`, `datatype`)."""
polygons = cell.get_polygons(layer=layer, datatype=datatype)
xs = [x for poly in polygons for x, y in poly.points]
ys = [y for poly in polygons for x, y in poly.points]
xmin, xmax = min(xs), max(xs)
ymin, ymax = min(ys), max(ys)
center = mp.Vector3(0.5 * (xmin + xmax), 0.5 * (ymin + ymax), 0.5 * (zmin + zmax))
size = mp.Vector3(xmax - xmin, ymax - ymin, zmax - zmin)
dims = 2 if (zmin == 0 and zmax == 0) else 3
return mp.Volume(center=center, size=size, dims=dims)
def create_ring_gds(radius, width):
lib = gdstk.Library()
ring_cell = lib.new_cell(f"ring_resonator_r{radius}_w{width}")
# Draw the ring
ring_cell.add(
gdstk.ellipse(
(0, 0),
radius + width / 2,
inner_radius=radius - width/2,
layer=RING_LAYER,
)
)
# Draw the first source
ring_cell.add(
gdstk.rectangle((radius - width, 0), (radius + width, 0), layer=SOURCE0_LAYER)
)
# Draw the second source
ring_cell.add(
gdstk.rectangle((-radius - width, 0), (-radius + width, 0), layer=SOURCE1_LAYER)
)
# Draw the monitor location
ring_cell.add(
gdstk.rectangle(
(radius - width / 2, 0), (radius + width / 2, 0), layer=MONITOR_LAYER
)
)
# Draw the simulation domain
pad = 2 # padding between waveguide and edge of PML
ring_cell.add(
gdstk.rectangle(
(-radius - width / 2 - pad, -radius-width/2-pad),
(radius + width / 2 + pad, radius + width / 2 + pad),
layer=SIMULATION_LAYER,
)
)
filename = f"ring_r{radius}_w{width}.gds"
lib.write_gds(filename)
return filename
def find_modes(filename,wvl=1.55,bw=0.05):
# Read in the ring structure using gdstk
gds_cell = get_gds_cell(filename)
geometry = get_gds_prisms(Si, gds_cell, RING_LAYER, zmin=-100, zmax=100)
cell = get_gds_vol(gds_cell, SIMULATION_LAYER, zmin=zmin, zmax=zmax)
src_vol0 = get_gds_vol(gds_cell, SOURCE0_LAYER, zmin=zmin, zmax=zmax)
src_vol1 = get_gds_vol(gds_cell, SOURCE1_LAYER, zmin=zmin, zmax=zmax)
mon_vol = get_gds_vol(gds_cell, MONITOR_LAYER, zmin=zmin, zmax=zmax)
fcen = 1 / wvl
df = bw * fcen
src = [
mp.Source(
mp.GaussianSource(fcen, fwidth=df),
component=mp.Hz,
volume=src_vol0,
),
mp.Source(
mp.GaussianSource(fcen, fwidth=df),
component=mp.Hz,
volume=src_vol1,
amplitude=-1,
)
]
sim = mp.Simulation(
cell_size=cell.size,
geometry=geometry,
sources=src,
resolution=resolution,
boundary_layers=[mp.PML(dpml)],
default_material=SiO2
)
h = mp.Harminv(mp.Hz, mon_vol.center, fcen, df)
sim.run(mp.after_sources(h), until_after_sources=100)
fig, ax = plt.subplots()
sim.plot2D(
ax=ax,
fields=mp.Hz,
eps_parameters={'contour':True}
)
fig.savefig("ring_fields.png", bbox_inches="tight", dpi=150)
wvl = np.array([1 / m.freq for m in h.modes])
Q = np.array([m.Q for m in h.modes])
sim.reset_meep()
return wvl, Q
if __name__ == "__main__":
filename = create_ring_gds(2.0, 0.5)
wvls, Qs = find_modes(filename, 1.55, 0.05)
for w, Q in zip(wvls, Qs):
print(f"mode: {w}, {Q}")
Note the absence of symmetries even though, in principle, the ring geometry and the two line sources satisfy two mirror symmetry planes through the \(x\) (even) and \(y\) (odd) axes. This omission is due to the fact that the ring geometry created using gdstk and imported from the GDS file is actually a Prism consisting of a discrete number of vertices (rather than two overlapping Cylinders as in Tutorial/Basics/Modes of a Ring Resonator). Discretization artifacts of the ring geometry slightly break its mirror symmetry. (Attempting to use symmetries in this case yields unpredictable results.)
For this ring geometry, Harminv finds a mode with wavelength 1.5490604 μm and \(Q\) of 124691.308. The \(H_z\) field profile is shown below. As expected, due to the large \(Q\) the mode is tightly confined to the ring and exhibits little radiative loss.
