Meep FDTD#
Meep is a free and open source finite-difference time-domain (FDTD) software package for electromagnetic simulations spanning a broad range of applications.
You can install Meep, NLOPT and the mode solver MPB (which Meep uses to compute S-parameters and launch mode sources) using Conda
conda install -c conda-forge pymeep=*=mpi_mpich_* nlopt -y
To update the installed package, you can do:
conda update pymeep=*=mpi_mpich_* -y
The Conda packages are available for Linux, macOS, or Windows WSL. For more details on installing Meep, see the user manual.
The gdsfactory gmeep plugin computes the transmission spectrum for planar photonic components.
One advantage of using the gmeep plugin is that you only need to define your component geometry once using gdsfactory. The geometry is automatically imported into Meep, so there is no need to define the geometry separately for Meep.
For extracting S-parameters, gmeep automatically swaps the source between ports to compute the full S-matrix. This process involves:
adding monitors to each port of the device
extending the ports into the PML absorbing boundary layers
running the simulation and computing elements of the S-matrix in post-processing using the correct ratio of S-parameters. The port monitors are used to compute the discrete-time Fourier-transformed (DTFT) fields which are then decomposed into complex mode coefficients known as S-parameters. The S-parameters can be computed over a range of frequencies.
The resolution is in units of pixels/μm. As a general rule, you should run with at least resolution=30 for 1/30 μm/pixel (33 nm/pixel)
Note that most of the examples use resolution=20 in order to run fast.
Here are some examples of how to extract S-parameters in Meep in planar devices.
For Sparameters we follow the syntax o1@0,o2@0 where o1 is the input port @0 mode 0 (usually fundamental TE mode) and o2@0 refers to output port o2 mode 0.
To speed up the simulations below we define the file to save filepath and automatically loads the simulation results when the file already exists.
top view
________________________________
| |
| xmargin_left | port_extension
|<---------> port_margin ||<-->
o2_|___________ _________||_o3
| \ / |
| \ / |
| ====== |
| / \ |
o1_|___________/ \__________|_o4
| | <-------->|
| |ymargin_bot xmargin_right|
| | |
|___|___________________________|
side view
________________________________
| | |
| | |
| zmargin_top |
|xmargin_left | |
|<---> _____ _|___ |
| | | | | |
| | | | | |
| |_____| |_____| |
| | |
| | |
| |zmargin_bot |
| | |
|_______|_______________________|
Serial Simulation (single core)#
Running Meep using a single CPU core can be slow as the single core needs to update all the simulation grid points sequentially at each time step.
import gdsfactory as gf
import matplotlib.pyplot as plt
import numpy as np
from autograd import tensor_jacobian_product
from gdsfactory.generic_tech import get_generic_pdk
from meep import MaterialGrid, Medium, Vector3, Volume
from meep.adjoint import (
DesignRegion,
get_conic_radius_from_eta_e,
)
import gplugins
import gplugins.gmeep as gm
gf.config.rich_output()
PDK = get_generic_pdk()
PDK.activate()
Using MPI version 4.1, 1 processes
2024-05-10 10:17:41.883 | INFO | gplugins.gmeep:<module>:39 - Meep '1.28.0' installed at ['/home/runner/micromamba/lib/python3.11/site-packages/meep']
c = gf.components.straight(length=2)
c.plot()
run=False only plots the simulations for you to review that is set up correctly
sp = gm.write_sparameters_meep(c, run=False, ymargin_top=3, ymargin_bot=3, is_3d=False)
2024-05-10 10:17:42.347 | INFO | gdsfactory.component:_write_library:2003 - Wrote to '/tmp/gdsfactory/straight_length2.gds'
2024-05-10 10:17:42.363 | WARNING | gdsfactory.klive:show:49 - UserWarning: Could not connect to klive server. Is klayout open and klive plugin installed?
2024-05-10 10:17:42.389 | WARNING | meep:_get_epsilon_grid:4442 - ComplexWarning: Casting complex values to real discards the imaginary part
help(gm.write_sparameters_meep)
Help on function write_sparameters_meep in module gplugins.gmeep.write_sparameters_meep:
write_sparameters_meep(component: 'ComponentSpec', port_source_names: 'list[str] | None' = None, port_source_modes: 'dict[str, list]' = None, port_modes: 'list[int]' = None, port_symmetries: 'PortSymmetries | None' = None, resolution: 'int' = 30, wavelength_start: 'float' = 1.5, wavelength_stop: 'float' = 1.6, wavelength_points: 'int' = 50, dirpath: 'PathType | None' = None, layer_stack: 'LayerStack | None' = None, port_margin: 'float' = 2, port_monitor_offset: 'float' = -0.1, port_source_offset: 'float' = -0.1, filepath: 'Path | None' = None, overwrite: 'bool' = False, animate: 'bool' = False, animate_center: 'tuple' = None, animate_size: 'tuple' = None, lazy_parallelism: 'bool' = False, run: 'bool' = True, dispersive: 'bool' = False, xmargin: 'float' = 0, ymargin: 'float' = 0, zmargin: 'float' = 0, xmargin_left: 'float' = 0, xmargin_right: 'float' = 0, ymargin_top: 'float' = 0, ymargin_bot: 'float' = 0, zmargin_top: 'float' = 0, zmargin_bot: 'float' = 0, decay_by: 'float' = 0.001, is_3d: 'bool' = False, z: 'float' = 0, plot_args: 'dict | None' = None, only_return_filepath_sim_settings=False, verbosity: 'int' = 0, **settings) -> 'dict[str, np.ndarray]'
Returns Sparameters and writes them to npz filepath.
Simulates each time using a different input port (by default, all of them)
unless you specify port_symmetries:
.. code::
port_symmetries_crossing = {
"o1@0,o1@0": ["o2@0,o2@0", "o3@0,o3@0", "o4@0,o4@0"],
"o2@0,o1@0": ["o1@0,o2@0", "o3@0,o4@0", "o4@0,o3@0"],
"o3@0,o1@0": ["o1@0,o3@0", "o2@0,o4@0", "o4@0,o2@0"],
"o4@0,o1@0": ["o1@0,o4@0", "o2@0,o3@0", "o3@0,o2@0"],
}
- Only simulations using the outer key port names will be run
- The associated value is another dict whose keys are the S-parameters computed
when this source is active
- The values of this inner Dict are lists of s-parameters whose values are copied
.. code::
top view
________________________________
| |
| xmargin_left | port_extension
|<---------> port_margin ||<-->
o2_|___________ _________||_o3
| \ / |
| \ / |
| ====== |
| / \ |
o1_|___________/ \__________|_o4
| | <-------->|
| |ymargin_bot xmargin_right|
| | |
|___|___________________________|
side view
________________________________
| | |
| | |
| zmargin_top |
|xmargin_left | |
|<---> _____ _|___ |
| | | | | |
| | | | | |
| |_____| |_____| |
| | |
| | |
| |zmargin_bot |
| | |
|_______|_______________________|
Args:
component: to simulate.
resolution: in pixels/um (30: for coarse, 100: for fine).
port_source_names: list of ports to excite. Defaults to all.
port_symmetries: Dict to specify port symmetries, to save number of simulations.
dirpath: directory to store Sparameters.
layer_stack: contains layer to thickness, zmin and material.
Defaults to active pdk.layer_stack.
port_margin: margin on each side of the port.
port_monitor_offset: offset between Component and monitor port in um.
port_source_offset: offset between Component and source port in um.
filepath: to store pandas Dataframe with Sparameters in npz format.
Defaults to dirpath/component_.npz.
overwrite: overwrites stored Sparameter npz results.
animate: saves a MP4 images of the simulation for inspection, and also
outputs during computation. The name of the file is the source index.
lazy_parallelism: toggles the flag "meep.divide_parallel_processes" to
perform the simulations with different sources in parallel.
By default MPI just runs the same copy of the Python script everywhere,
with the C++ under MEEP actually being parallelized.
divide_parallel_processes allows us to logically split this one calculation
into (in this case "cores") subdivisions.
The only difference in the scripts is that a different integer n
is returned depending on the subdivision it is running in.
So we use that n to select different sources, and each subdivision calculates
its own Sparams independently. Afterwards, we collect all
results in one of the subdivisions (if rank == 0).
run: runs simulation, if False, only plots simulation.
dispersive: use dispersive models for materials (requires higher resolution).
xmargin: left and right distance from component to PML.
xmargin_left: west distance from component to PML.
xmargin_right: east distance from component to PML.
ymargin: top and bottom distance from component to PML.
ymargin_top: north distance from component to PML.
ymargin_bot: south distance from component to PML.
zmargin_top: +z distance from component to PML.
zmargin_bot: -z distance from component to PML.
is_3d: if True runs in 3D (much slower).
z: for 2D plot.
plot_args: if animate or not run, customization keyword arguments passed to
`plot2D()` (i.e. `labels`, `eps_parameters`, `boundary_parameters`, `field_parameters`, etc.)
keyword Args:
extend_ports_length: to extend ports beyond the PML (um).
zmargin_top: thickness for cladding above core (um).
zmargin_bot: thickness for cladding below core (um).
tpml: PML thickness (um).
clad_material: material for cladding.
wavelength_start: wavelength min (um).
wavelength_stop: wavelength max (um).
wavelength_points: wavelength steps.
dfcen: delta frequency.
port_source_name: input port name.
port_margin: margin on each side of the port (um).
distance_source_to_monitors: in (um).
port_source_offset: offset between source Component port and source MEEP port.
port_monitor_offset: offset between Component and MEEP port monitor.
material_name_to_meep: map layer_stack names with meep material database name
or refractive index. dispersive materials have a wavelength dependent index.
Returns:
sparameters in a Dict (wavelengths, s11a, o1@0,o2@0, ...)
where `a` is the angle in radians and `m` the module.
As you’ve noticed we added ymargin_top and ymargin_bot to ensure we have enough distance to the PML
You can also do this directly with another version of the function that adds ymargin_top and ymargin_bot
c = gf.components.straight(length=2)
sp = gm.write_sparameters_meep(c, run=False, is_3d=False)
2024-05-10 10:17:42.543 | INFO | gdsfactory.component:_write_library:2003 - Wrote to '/tmp/gdsfactory/straight_length2.gds'
Because components with left-right ports are very common write_sparameters_meep y_margin = 3um
c = gf.components.taper(length=2.0, width1=0.5, width2=1)
c.plot()
sp = gm.write_sparameters_meep(c, run=False)
2024-05-10 10:17:42.881 | INFO | gdsfactory.component:_write_library:2003 - Wrote to '/tmp/gdsfactory/taper_length2p0_width21.gds'
2.5D Simulation#
For faster simulations you can do an effective mode approximation, to compute the mode of the slab and run a 2D simulation to speed your simulations
core_material = gplugins.get_effective_indices(
core_material=3.4777,
clad_materialding=1.444,
nsubstrate=1.444,
thickness=0.22,
wavelength=1.55,
polarization="te",
)[0]
core_material
2.8494636999424405
sp = gm.write_sparameters_meep(
c,
resolution=20,
is_3d=False,
material_name_to_meep=dict(si=core_material),
filepath="data/taper.npz",
)
2024-05-10 10:17:43.155 | INFO | gplugins.gmeep.write_sparameters_meep:write_sparameters_meep:396 - Simulation loaded from PosixPath('data/taper.npz')
gplugins.plot.plot_sparameters(sp)
gplugins.plot.plot_sparameters(sp, keys=("o2@0,o1@0",))
For a small taper length, the matrix element S\(_{21}\) (transmission in Port 2 given a source in Port 1) is around 0 dB which is equivalent to ~100% transmission.
Port Symmetries#
You can skip redundant simulations in reciprocal devices. If the device looks the same going from in -> out as out -> in, we just need to run one simulation.
c = gf.components.bend_euler(radius=3)
c.plot()
2024-05-10 10:17:43.392 | WARNING | gdsfactory.cross_section:validate_radius:223 - UserWarning: min_bend_radius 3 < CrossSection.radius_min 5.0.
sp = gm.write_sparameters_meep_1x1_bend90(c, run=False)
2024-05-10 10:17:43.565 | INFO | gdsfactory.component:_write_library:2003 - Wrote to '/tmp/gdsfactory/bend_euler_radius3.gds'
Warning: grid volume is not an integer number of pixels; cell size will be rounded to nearest pixel.
sp = gm.write_sparameters_meep_1x1_bend90(c, run=True, filepath="data/bend90_meep.npz")
2024-05-10 10:17:43.981 | INFO | gplugins.gmeep.write_sparameters_meep:write_sparameters_meep:396 - Simulation loaded from PosixPath('data/bend90_meep.npz')
list(sp.keys())
['o1@0,o1@0', 'o2@0,o1@0', 'o1@0,o2@0', 'o2@0,o2@0', 'wavelengths']
gplugins.plot.plot_sparameters(sp)
gplugins.plot.plot_sparameters(sp, keys=("o2@0,o1@0",), logscale=False)
gplugins.plot.plot_sparameters(sp, keys=("o2@0,o1@0",))
c = gf.components.crossing()
c.plot()
Here are the port symmetries for a crossing
port_symmetries_crossing = {
"o1": {
"o1@0,o1@0": ["o2@0,o2@0", "o3@0,o3@0", "o4@0,o4@0"],
"o2@0,o1@0": ["o1@0,o2@0", "o3@0,o4@0", "o4@0,o3@0"],
"o3@0,o1@0": ["o1@0,o3@0", "o2@0,o4@0", "o4@0,o2@0"],
"o4@0,o1@0": ["o1@0,o4@0", "o2@0,o3@0", "o3@0,o2@0"],
}
}
sp = gm.write_sparameters_meep(
c,
resolution=20,
ymargin=0,
port_symmetries=gm.port_symmetries.port_symmetries_crossing,
run=False,
)
2024-05-10 10:17:44.577 | INFO | gdsfactory.component:_write_library:2003 - Wrote to '/tmp/gdsfactory/crossing.gds'
sp = gm.write_sparameters_meep(
c,
resolution=20,
ymargin=0,
port_symmetries=gm.port_symmetries.port_symmetries_crossing,
run=True,
filepath="data/meep_crossing.npz",
)
2024-05-10 10:17:45.151 | INFO | gplugins.gmeep.write_sparameters_meep:write_sparameters_meep:396 - Simulation loaded from PosixPath('data/meep_crossing.npz')
gm.plot.plot_sparameters(sp)
gm.plot.plot_sparameters(sp, keys=("o3@0,o1@0",))
As you can see this crossing looks beautiful but is quite lossy (9 dB @ 1550 nm)
Multimode Simulations#
You can also simulate structures that source and measure multiple modes. We define the port_source_modes as a dictionary with the keys as the source keys and the values as list of mode numbers to output (starting at 0 for TE00). Similarly on the output ports, port_modes is a list of modes to measure at every output port.
c = gf.components.straight(length=5, width=2)
sp = gm.write_sparameters_meep(
c,
run=False,
ymargin_top=3,
ymargin_bot=3,
is_3d=False,
resolution=20,
)
2024-05-10 10:17:45.492 | INFO | gdsfactory.component:_write_library:2003 - Wrote to '/tmp/gdsfactory/straight_length5_width2.gds'
sp = gm.write_sparameters_meep(
c,
run=True,
ymargin_top=3,
ymargin_bot=3,
is_3d=False,
port_source_names=["o1"],
port_source_modes={"o1": [0, 1]},
port_modes=[0, 1],
resolution=20,
overwrite=False,
)
gm.plot.plot_sparameters(sp, with_simpler_labels=False)
2024-05-10 10:17:45.684 | INFO | gdsfactory.component:_write_library:2003 - Wrote to '/tmp/gdsfactory/straight_length5_width2.gds'
2024-05-10 10:17:45.710 | WARNING | meep:create_structure:4436 - ComplexWarning: Casting complex values to real discards the imaginary part
2024-05-10 10:17:45.711 | WARNING | meep:_set_materials:4439 - ComplexWarning: Casting complex values to real discards the imaginary part
2024-05-10 10:17:50.391 | INFO | gplugins.gmeep.write_sparameters_meep:write_sparameters_meep:583 - Write simulation results to PosixPath('/tmp/gdsfactory/sparameters/straight/straight_length5_width2_89835eccf5dbba18d0f4a75d146af9e6.npz')
2024-05-10 10:17:50.407 | INFO | gplugins.gmeep.write_sparameters_meep:write_sparameters_meep:585 - Write simulation settings to PosixPath('/tmp/gdsfactory/sparameters/straight/straight_length5_width2_89835eccf5dbba18d0f4a75d146af9e6.yml')
Parallel Simulation (multicore/MPI)#
Meep supports distributed-memory parallelization via MPI which can be used to provide a significant speedup compared to serial calculations.
By default MPI just runs the same copy of the Python script everywhere, with the C++ under MEEP actually being parallelized.
divide_parallel_processes allows us to logically split this one calculation into (in this case “cores”) subdivisions.
The only difference in the scripts is that a different integer n is returned depending on the subdivision it is running in.
So we use that n to select different sources, and each subdivision calculates its own Sparams independently. Afterwards, we collect all the results in one of the subdivisions (if rank == 0)
As a demonstration, lets try to reproduce the results of the directional coupler results from the Meep manual which indicates that to obtain a 3 dB (50%/50%) splitter you need a separation distance of 130 nm over a coupler length of 8 μm.
help(gf.components.coupler)
Help on function coupler in module gdsfactory.components.coupler:
coupler(gap: 'float' = 0.236, length: 'float' = 20.0, coupler_symmetric: 'ComponentFactory' = <function coupler_symmetric at 0x7fb5e852f7e0>, coupler_straight: 'ComponentFactory' = <function coupler_straight at 0x7fb5e852f600>, dy: 'float' = 4.0, dx: 'float' = 10.0, cross_section: 'CrossSectionSpec' = 'xs_sc') -> gdsfactory.component.Component
Symmetric coupler.
Args:
gap: between straights in um.
length: of coupling region in um.
coupler_symmetric: spec for bend coupler.
coupler_straight: spec for straight coupler.
dy: port to port vertical spacing in um.
dx: length of bend in x direction in um.
cross_section: spec (CrossSection, string or dict).
.. code::
dx dx
|------| |------|
o2 ________ ______o3
\ / |
\ length / |
======================= gap | dy
/ \ |
________/ \_______ |
o1 o4
coupler_straight coupler_symmetric
c = gf.components.coupler(length=8, gap=0.13)
c.plot()
gm.write_sparameters_meep(component=c, run=False)
2024-05-10 10:17:50.735 | INFO | gdsfactory.component:_write_library:2003 - Wrote to '/tmp/gdsfactory/coupler_gap0p13_length8.gds'
Warning: grid volume is not an integer number of pixels; cell size will be rounded to nearest pixel.
<meep.simulation.Simulation object at 0x7fb59ff7ca90>
filepath = gm.write_sparameters_meep_mpi(
component=c,
cores=4,
resolution=30,
filepath='data/meep_coupler.npz'
)
sp = np.load(filepath)
gplugins.plot.plot_sparameters(sp, keys=["o1@0,o3@0", "o1@0,o4@0"])
Batch Simulations#
You can also run a batch of multicore simulations.
Given a list of write_sparameters_meep keyword arguments jobs launches them in different cores using MPI where each simulation runs with “cores_per_run” cores.
If there are more simulations than cores each batch runs sequentially.
help(gm.write_sparameters_meep_batch)
Help on function write_sparameters_meep_batch in module gplugins.gmeep.write_sparameters_meep_batch:
write_sparameters_meep_batch(jobs: 'list[dict]', cores_per_run: 'int' = 2, total_cores: 'int' = 4, temp_dir: 'Path' = PosixPath('/home/runner/.gdsfactory/sp/temp'), delete_temp_files: 'bool' = True, dirpath: 'Path | None' = None, layer_stack: 'LayerStack | None' = None, **kwargs) -> 'list[Path]'
Write Sparameters for a batch of jobs using MPI and returns results filepaths.
Given a list of write_sparameters_meep keyword arguments `jobs` launches them in
different cores using MPI where each simulation runs with `cores_per_run` cores.
If there are more simulations than cores each batch runs sequentially.
Args
jobs: list of Dicts containing the simulation settings for each job.
for write_sparameters_meep.
cores_per_run: number of processors to assign to each component simulation.
total_cores: total number of cores to use.
temp_dir: temporary directory to hold simulation files.
delete_temp_files: deletes temp_dir when done.
dirpath: directory to store Sparameters.
layer_stack: contains layer to thickness, zmin and material.
Defaults to active pdk.layer_stack.
keyword Args:
resolution: in pixels/um (30: for coarse, 100: for fine).
port_symmetries: Dict to specify port symmetries, to save number of simulations.
dirpath: directory to store Sparameters.
port_margin: margin on each side of the port.
port_monitor_offset: offset between monitor GDS port and monitor MEEP port.
port_source_offset: offset between source GDS port and source MEEP port.
filepath: to store pandas Dataframe with Sparameters in CSV format..
animate: saves a MP4 images of the simulation for inspection, and also
outputs during computation. The name of the file is the source index.
lazy_parallelism: toggles the flag "meep.divide_parallel_processes" to
perform the simulations with different sources in parallel.
dispersive: use dispersive models for materials (requires higher resolution).
xmargin: left and right distance from component to PML.
xmargin_left: west distance from component to PML.
xmargin_right: east distance from component to PML.
ymargin: top and bottom distance from component to PML.
ymargin_top: north distance from component to PML.
ymargin_bot: south distance from component to PML.
extend_ports_length: to extend ports beyond the PML
layer_stack: Dict of layer number (int, int) to thickness (um).
zmargin_top: thickness for cladding above core.
zmargin_bot: thickness for cladding below core.
tpml: PML thickness (um).
clad_material: material for cladding.
is_3d: if True runs in 3D.
wavelength_start: wavelength min (um).
wavelength_stop: wavelength max (um).
wavelength_points: wavelength steps.
dfcen: delta frequency.
port_source_name: input port name.
port_margin: margin on each side of the port.
distance_source_to_monitors: in (um) source goes before.
port_source_offset: offset between source GDS port and source MEEP port.
port_monitor_offset: offset between monitor GDS port and monitor MEEP port.
Returns:
filepath list for sparameters numpy saved files (wavelengths, o1@0,o2@0, ...).
c = gf.components.straight(length=3.1)
gm.write_sparameters_meep(c, ymargin=3, run=False)
2024-05-10 10:17:57.165 | INFO | gdsfactory.component:_write_library:2003 - Wrote to '/tmp/gdsfactory/straight_length3p1.gds'
<meep.simulation.Simulation object at 0x7fb59e9e1890>
c1_dict = {"component": c, "ymargin": 3}
jobs = [
c1_dict,
]
filepaths = gm.write_sparameters_meep_batch_1x1(
jobs=jobs,
cores_per_run=4,
total_cores=8,
lazy_parallelism=True,
is_3d=False,
filepath='data/meep_straight3.npz'
)
sp = np.load(filepaths[0])
gplugins.plot.plot_sparameters(sp)
Adjoint Optimization#
gdsfactory extends Meep’s Adjoint Optimization features to optimize and generate primitive photonic components.
This example is based on this Meep Adjoint Optimization tutorial
We define some useful variables that we will need later. We can leave out many of the small design parameters by defining a minimum length and applying that to a filter and using that as a constraint in our optimization.
design_region_width = 2.5
design_region_height = 3.5
eta_e = 0.55
minimum_length = 0.1
filter_radius = get_conic_radius_from_eta_e(minimum_length, eta_e)
eta_i = 0.5
eta_d = 1 - eta_e
resolution = 20
design_region_resolution = int(5 * resolution)
Nx = int(design_region_resolution * design_region_width)
Ny = int(design_region_resolution * design_region_height)
pml_size = 1.0
waveguide_length = 1.5
waveguide_width = 0.5
Sx = 2 * pml_size + 2 * waveguide_length + design_region_width
Sy = 2 * pml_size + design_region_height + 0.5
cell_size = (Sx, Sy)
SiO2 = Medium(index=1.44)
Si = Medium(index=3.4)
We define the design region, design variables, and the component to optimize.
design_variables = MaterialGrid(Vector3(Nx, Ny), SiO2, Si, grid_type="U_MEAN")
design_region = DesignRegion(
design_variables,
volume=Volume(
center=Vector3(),
size=Vector3(design_region_width, design_region_height, 0),
),
)
c = gf.Component("mmi1x2")
arm_separation = 1.0
straight1 = c << gf.components.taper(6.5, width2=1)
straight1.move(straight1.ports["o2"], (-design_region_width / 2.0, 0))
straight2 = c << gf.components.taper(6.5, width1=1, width2=0.5)
straight2.move(straight2.ports["o1"], (design_region_width / 2.0, 1))
straight3 = c << gf.components.taper(6.5, width1=1, width2=0.5)
straight3.move(straight3.ports["o1"], (design_region_width / 2.0, -1))
c.add_port("o1", port=straight1.ports["o1"])
c.add_port("o2", port=straight2.ports["o2"])
c.add_port("o3", port=straight3.ports["o2"])
c.plot()
We define the objective function, and obtain the optimization object.
def mapping(x, eta, beta):
# filter
filtered_field = conic_filter(
x,
filter_radius,
design_region_width,
design_region_height,
design_region_resolution,
)
# projection
projected_field = tanh_projection(filtered_field, beta, eta)
projected_field = (
npa.fliplr(projected_field) + projected_field
) / 2 # up-down symmetry
# interpolate to actual materials
return projected_field.flatten()
seed = 240
np.random.seed(seed)
x0 = mapping(
np.random.rand(Nx * Ny),
eta_i,
128,
)
def J(source, top, bottom):
power = npa.abs(top / source) ** 2 + npa.abs(bottom / source) ** 2
return npa.mean(power)
opt = gm.get_meep_adjoint_optimizer(
c,
J,
[design_region],
[design_variables],
x0,
resolution=resolution,
cell_size=(
Sx + 2 + design_region_width + 2 * pml_size,
design_region_height + 2 * pml_size + 1.5,
),
tpml=1.0,
extend_ports_length=0,
port_margin=1,
port_source_offset=-5.5,
port_monitor_offset=-5.5,
symmetries=[mp.Mirror(direction=mp.Y)],
wavelength_points=10,
)
We’ll define a simple objective function that returns the gradient, and records the figure of merit. We’ll plot the new geometry after each iteration.
evaluation_history = []
cur_iter = [0]
def f(v, gradient, cur_beta):
print(f"Current iteration: {cur_iter[0] + 1}")
f0, dJ_du = opt([mapping(v, eta_i, cur_beta)])
plt.figure()
ax = plt.gca()
opt.plot2D(
False,
ax=ax,
plot_sources_flag=False,
plot_monitors_flag=False,
plot_boundaries_flag=False,
)
plt.show()
if gradient.size > 0:
gradient[:] = tensor_jacobian_product(mapping, 0)(
v, eta_i, cur_beta, np.sum(dJ_du, axis=1)
)
evaluation_history.append(np.max(np.real(f0)))
cur_iter[0] = cur_iter[0] + 1
return np.real(f0)
We can define bitmasks to describe the boundary conditions.
# Define spatial arrays used to generate bit masks
x_g = np.linspace(-design_region_width / 2, design_region_width / 2, Nx)
y_g = np.linspace(-design_region_height / 2, design_region_height / 2, Ny)
X_g, Y_g = np.meshgrid(x_g, y_g, sparse=True, indexing="ij")
# Define the core mask
left_wg_mask = (X_g == -design_region_width / 2) & (np.abs(Y_g) <= waveguide_width)
top_right_wg_mask = (X_g == design_region_width / 2) & (
np.abs(Y_g + arm_separation) <= waveguide_width
)
bottom_right_wg_mask = (X_g == design_region_width / 2) & (
np.abs(Y_g - arm_separation) <= waveguide_width
)
Si_mask = left_wg_mask | top_right_wg_mask | bottom_right_wg_mask
# Define the cladding mask
border_mask = (
(X_g == -design_region_width / 2)
| (X_g == design_region_width / 2)
| (Y_g == -design_region_height / 2)
| (Y_g == design_region_height / 2)
)
SiO2_mask = border_mask.copy()
SiO2_mask[Si_mask] = False
We can then finally run the optimizer, and visualize the optimized component.
n = Nx * Ny # number of parameters
# Initial guess
x = np.ones((n,)) * 0.5
x[Si_mask.flatten()] = 1 # set the edges of waveguides to silicon
x[SiO2_mask.flatten()] = 0 # set the other edges to SiO2
# lower and upper bounds
lb = np.zeros((Nx * Ny,))
lb[Si_mask.flatten()] = 1
ub = np.ones((Nx * Ny,))
ub[SiO2_mask.flatten()] = 0
cur_beta = 4
beta_scale = 2
num_betas = 7
update_factor = 12
run_optimization = False
if run_optimization:
for iters in range(num_betas):
print("current beta: ", cur_beta)
if iters != num_betas - 1:
x[:] = gm.run_meep_adjoint_optimizer(
n,
lambda a, g: f(a, g, cur_beta),
x,
lower_bound=lb,
upper_bound=ub,
maxeval=update_factor,
)
else:
optimized_component = gm.run_meep_adjoint_optimizer(
n,
lambda a, g: f(a, g, cur_beta),
x,
lower_bound=lb,
upper_bound=ub,
maxeval=update_factor,
get_optimized_component=True,
opt=opt,
threshold_offset_from_max=0.09,
)
cur_beta = cur_beta * beta_scale
optimized_component.plot()
final_figure_of_merit = 10 * np.log10(
0.5 * np.array(evaluation_history[-1])
) # around -3.7 dB
The final optimized structure should look like this:
