Standard grating coupler: parameter optimization

import os
import numpy as np
import matplotlib.pyplot as plt

import gdsfactory as gf
from gdsfactory.technology import LogicalLayer, LayerLevel, LayerStack
import luminescent as lumi

# SOI grating-coupler dimensions [um]
w_wg = 0.5
taper_length = 15.0
n_periods = 25
fill_factor = 0.5
fiber_angle = 6.9 # in clad or 10deg in air
taper_angle = 40.0
period=.61

# SOI dimensions [um]: 220 nm silicon with a 130 nm residual slab
h_core = 0.22
h_slab = 0.13
h_teeth = h_core - h_slab

# wavelengths [um]
wavelength = 1.55
bandwidth = 0.1
nwls = 11
wavelengths = np.linspace(wavelength - bandwidth / 2, wavelength + bandwidth / 2, nwls)

# Gaussian fiber mode dimensions [um]
D_beam = 10.4  # mode field diameter, SMF-28
Dlim = 1.5 * D_beam  # beam frame width
fiber_offset = 0.6


def setup(path, wavelength, wavelengths,n_periods, period,fill_factor, taper_length, taper_angle):
    boundaries = ["PML", ["PEC", "PML"], "PML"]
    nres = 6
    stop_decay = 1e-3
    PML_decay = [0.01, 0.1, 0.01]
    shot_every = 10

    l_wg = 2
    margin = 1.

    source_offset = fiber_offset/2
    zmargin = .6
    zmin =  - zmargin
    zmax = h_core + fiber_offset
    z_port = zmax - source_offset

    nclad = 1.444
    materials = {
        "Si": lumi.material(3.478),
        "SiO2": lumi.material(nclad),
    }
    background = "SiO2"

    WG = (1, 0)
    SLAB = (2, 0)
    BBOX = (3, 0)
    layer_stack = LayerStack(
        layers={
            "core": LayerLevel(
                layer=LogicalLayer(layer=WG),
                thickness=h_core,
                zmin=0,
                material="Si",
                mesh_order=200,
            ),
            "slab": LayerLevel(
                layer=LogicalLayer(layer=SLAB),
                thickness=h_slab,
                zmin=0,
                material="Si",
                mesh_order=100,
            ),
        }
    )

    c = gf.Component()
    gc = c << gf.components.grating_coupler_elliptical_uniform(
        n_periods=n_periods,
        period=period,
        fill_factor=fill_factor,
        taper_length=taper_length,
        taper_angle=taper_angle,
        wavelength=wavelength,
        fiber_angle=fiber_angle,
        nclad=nclad,
        layer_slab=SLAB,
        cross_section=gf.cross_section.cross_section(width=w_wg, layer=WG),
    )
    wg = c << gf.components.straight(
        length=l_wg,
        cross_section=gf.cross_section.cross_section(width=w_wg, layer=WG),
    )
    wg.connect("o2", gc.ports["o1"])
    c << gf.components.bbox(component=c, layer=BBOX, top=margin, bottom=margin)
    full_bbox = c.bbox_np().tolist()

    clipper = gf.Component()
    clipper_ref = clipper << gf.components.rectangle(
        size=(full_bbox[1][0] - full_bbox[0][0], full_bbox[1][1]),
        layer=WG,
    )

    half_component = gf.Component()
    clipper_ref.movex(full_bbox[0][0])
    for layer in (WG, SLAB, BBOX):
        clipped = gf.boolean(
            c,
            clipper,
            operation="and",
            layer=layer,
            layer1=layer,
            layer2=WG,
        )
        half_component << clipped
    half_component.add_port("wg", port=gc.ports["o1"])
    half_component.flatten()
    c = half_component
    c.plot()

    bbox = c.bbox_np().tolist()
    bbox[0][1] = 0
    bbox[0].append(zmin)
    bbox[1].append(zmax)

# sources, monitors
    frames = [
        lumi.plane_frame(
            "wg",
            start=[-w_wg / 2 - margin, zmin],
            stop=[0, zmax],
            boundaries=[("PML", "PEC"), "PML"],
            wls_mode_solver=[wavelength],
        ),
        lumi.plane_frame(
            "fiber",
            modes=[
                [
                    lumi.gaussian_mode(
                        D_beam,
                        nclad,
                        dx=wavelength / nclad / nres,
                        theta=np.deg2rad(fiber_angle),
                        wl=wavelength,
                        half="y+",
                        Dlim=Dlim,
                    )
                ]
            ],
            start=[-Dlim / 2, 0],
            stop=[Dlim / 2, Dlim / 2],
        ),
    ]

    fiber_origin = [
        gc.ports["o2"].center[0] + np.sin(np.deg2rad(fiber_angle)) * fiber_offset,
        0,
        z_port,
    ]
    fiber_basis = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
    sources = [
        lumi.source(
            "fiber",
            fiber_origin,
            fiber_basis,
            frame="fiber",
            wl=wavelength,
            bandwidth=max(max(wavelengths) - min(wavelengths), wavelength / 4),
        )
    ]
    monitors = [
        lumi.gf_monitor("wg", c, layer_stack, frame="wg"),
        lumi.monitor("fiber", fiber_origin, fiber_basis, frame="fiber"),
    ]

    views = [lumi.view("Ey", y=0, z=h_slab+h_teeth / 2)]
    lumi.mesh(path, c, layer_stack, materials, zmin, zmax)
    lumi.setup(
        path=path,
        wl=wavelength,
        wls=wavelengths,
        materials=materials,
        background=background,
        boundaries=boundaries,
        bbox=bbox,
        frames=frames,
        sources=sources,
        monitors=monitors,
        nres=nres,
        stop_decay=stop_decay,
        PML_decay=PML_decay,
        views=views,
        shot_every=shot_every,
        device="gpu",
    )

params_history=[]
efficiency_history=[]
params=[n_periods, period, fill_factor, taper_length, taper_angle]
efficiency=0
goal=.4 # uniform grating coupler can reach only 0.4 efficiency at best. for better, use apodization or inverse design
maxiter = 10
for iter in range(maxiter):
    path = os.path.join("runs", f"GC_opt_{iter}")
    setup(path, wavelength, wavelengths, *params)
    lumi.solve(path)
    
    sol = lumi.load(path)
    fiber_mode = sol["monitors"]["fiber"]["modes"][0]["-"]
    wg_mode = sol["monitors"]["wg"]["modes"][0]["+"]
    S = np.abs(wg_mode /fiber_mode) ** 2
    efficiency = S[nwls // 2]

    params_history.append(params.copy())
    efficiency_history.append(efficiency)
    print(f'{iter}: params={params}, efficiency={efficiency:.4f}')

    # todo for reader: use a black-box optimizer to update params
    break

    if efficiency > goal:
        break

lumi.movie(path)