Thermal¶

# !uv add femwell  # uncomment to install
#! uvx inc run gmsh

from collections import OrderedDict
from functools import partial

import gdsfactory as gf
import matplotlib.pyplot as plt
import numpy as np
from femwell.maxwell.waveguide import compute_modes
from femwell.mesh import mesh_from_OrderedDict
from femwell.thermal import solve_thermal
from gdsfactory.add_padding import add_padding, add_padding_container
from gdsfactory.generic_tech import get_generic_pdk
from gdsfactory.generic_tech.layer_map import LAYER
from gplugins.common.utils.geometry import region_to_shapely_polygons
from shapely.geometry import LineString, Polygon
from skfem import Basis, ElementTriP0
from skfem.io import from_meshio
from tqdm import tqdm

gf.config.rich_output()
PDK = get_generic_pdk()
PDK.activate()

LAYER_STACK = PDK.layer_stack

LAYER_STACK.layers["heater"].thickness = 0.13
LAYER_STACK.layers["heater"].zmin = 2.2

heater = gf.components.straight_heater_metal(length=50)
heater.plot()

s = heater.to_3d()
s.show()

w_sim = 16
h_substrate = 0.5
layer_names = ["core", "heater", "clad", "box", "substrate"]

heater_padded = add_padding_container(
    component=heater,
    function=partial(add_padding, layers=(LAYER.WAFER,), default=w_sim / 2),
)

x_center = (heater.dxmin + heater.dxmax) / 2
xsection_line = LineString([(x_center, -w_sim / 2), (x_center, w_sim / 2)])

substrate_bottom = LAYER_STACK.layers["box"].zmin - h_substrate
polygons = OrderedDict(
    bottom=LineString([(-w_sim / 2, substrate_bottom), (w_sim / 2, substrate_bottom)])
)

for layer_name in layer_names:
    layer_level = LAYER_STACK.layers[layer_name]
    region = layer_level.layer.get_shapes(heater_padded)
    if region.is_empty():
        continue

    shapely_polys = region_to_shapely_polygons(region)
    if not shapely_polys:
        continue

    zmin = layer_level.zmin
    thickness = layer_level.thickness
    if layer_name == "substrate":
        zmin = substrate_bottom
        thickness = h_substrate
    zmax = zmin + thickness

    for poly in (
        shapely_polys.geoms if hasattr(shapely_polys, "geoms") else [shapely_polys]
    ):
        intersection = poly.intersection(xsection_line)
        if intersection.is_empty:
            continue
        for geom in (
            intersection.geoms if hasattr(intersection, "geoms") else [intersection]
        ):
            y_min, y_max = geom.bounds[1], geom.bounds[3]
            rect = Polygon([(y_min, zmin), (y_max, zmin), (y_max, zmax), (y_min, zmax)])
            if layer_name in polygons:
                polygons[layer_name] = polygons[layer_name].union(rect)
            else:
                polygons[layer_name] = rect

resolutions = dict(
    core={"resolution": 0.04, "distance": 1},
    clad={"resolution": 0.6, "distance": 1},
    box={"resolution": 0.6, "distance": 1},
    heater={"resolution": 0.1, "distance": 1},
)

mesh = from_meshio(
    mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.6)
)
mesh.draw().show()

currents = np.linspace(0.0, 7.4e-3, 10)
current_densities = currents / polygons["heater"].area
neffs = []

for current_density in tqdm(current_densities):
    basis0 = Basis(mesh, ElementTriP0(), intorder=4)
    thermal_conductivity_p0 = basis0.zeros()
    for domain, value in {
        "core": 90,
        "box": 1.38,
        "clad": 1.38,
        "heater": 28,
        "substrate": 148,
    }.items():
        thermal_conductivity_p0[basis0.get_dofs(elements=domain)] = value
    thermal_conductivity_p0 *= 1e-12  # 1e-12 -> conversion from 1/m^2 -> 1/um^2

    basis, temperature = solve_thermal(
        basis0,
        thermal_conductivity_p0,
        specific_conductivity={"heater": 2.3e6},
        current_densities={"heater": current_density},
        fixed_boundaries={"bottom": 0},
    )

    if current_density == current_densities[-1]:
        basis.plot(temperature, shading="gouraud", colorbar=True)
        plt.show()

    temperature0 = basis0.project(basis.interpolate(temperature))
    epsilon = basis0.zeros() + (1.444 + 1.00e-5 * temperature0) ** 2
    epsilon[basis0.get_dofs(elements="core")] = (
        3.4777 + 1.86e-4 * temperature0[basis0.get_dofs(elements="core")]
    ) ** 2

    modes = compute_modes(basis0, epsilon, wavelength=1.55, num_modes=1, solver="scipy")

    if current_density == current_densities[-1]:
        modes[0].show(modes[0].E.real)

    neffs.append(np.real(modes[0].n_eff))

length = 320  # um
print(f"Phase shift: {2 * np.pi / 1.55 * (neffs[-1] - neffs[0]) * length}")
plt.xlabel("Current / mA")
plt.ylabel("Effective refractive index $n_{eff}$")
plt.plot(currents * 1e3, neffs)
plt.show()