Finite-Element Mode Solver#
Using femwell, you can mesh any component cross-section and solve PDEs with its powerful mode solver.
Unlike other solvers that rely on predefined geometries, femwell works directly with the actual component geometry. You can compute the modes of a GDSFactory cross-section, which internally defines a “uz” mesh perpendicular to a straight component using the provided cross-section.
Additionally, you can downsample layers from the LayerStack and modify both the cross-section and LayerStack before running the simulation to adjust the geometry. You can also define refractive indices based on the active PDK.
import logging
import sys
import gdsfactory as gf
import matplotlib.pyplot as plt
from femwell.maxwell.waveguide import compute_modes
from femwell.visualization import plot_domains
from gdsfactory.generic_tech import LAYER_STACK, get_generic_pdk
from gdsfactory.technology import LayerStack
from gplugins.gmsh import get_mesh
from rich.logging import RichHandler
from skfem import Basis, ElementTriP0
from skfem.io.meshio import from_meshio
gf.config.rich_output()
PDK = get_generic_pdk()
PDK.activate()
logger = logging.getLogger()
logger.removeHandler(sys.stderr)
logging.basicConfig(level="WARNING", datefmt="[%X]", handlers=[RichHandler()])
First we choose a component to simulate. Here, a straight strip waveguide:
xs = gf.cross_section.strip(width=1)
c = gf.components.straight(cross_section=xs)
c
Then we choose a Layer Stack. Here, we simply downsample the generic stack:
filtered_layer_stack = LayerStack(
layers={
k: LAYER_STACK.layers[k]
for k in (
"core",
"clad",
"slab90",
"box",
)
}
)
We can also change some of the values:
filtered_layer_stack.layers[
"core"
].thickness = 0.22 # Perturb the layer_stack before simulating
filtered_layer_stack.layers[
"slab90"
].thickness = 0.09 # Perturb the layer_stack before simulating
# When selecting resolutions, the names must match the keys of the layerstack
# Here, choose a finer mesh inside and close to the core
resolutions = {
"core": {"resolution": 0.02, "DistMax": 2, "SizeMax": 0.2},
}
Using gplugins, we quickly generate a cross-sectional mesh:
mesh_gmsh = get_mesh(
component=c,
layer_stack=filtered_layer_stack,
type="uz", # we want a cross-section
xsection_bounds=((1, -3), (1, 3)), # the line from which we take a cross-section
wafer_padding=3, # pad simulation domain 3 microns around the component
filename="mesh.msh",
resolutions=resolutions,
default_characteristic_length=0.5,
)
Info : [ 0%] Difference
Info : [ 10%] Difference
Info : [ 20%] Difference
Info : [ 30%] Difference
Info : [ 40%] Difference
Info : [ 50%] Difference
Info : [ 70%] Difference
Info : [ 80%] Difference - Splitting faces
Info : [ 0%] Fragments
Info : [ 10%] Fragments
Info : [ 20%] Fragments
Info : [ 30%] Fragments
Info : [ 40%] Fragments
Info : [ 50%] Fragments
Info : [ 70%] Fragments
Info : [ 80%] Fragments - Splitting faces
Info : [ 0%] Difference
Info : [ 10%] Difference
Info : [ 20%] Difference
Info : [ 30%] Difference
Info : [ 40%] Difference
Info : [ 50%] Difference
Info : [ 60%] Difference
Info : [ 70%] Difference
Info : [ 80%] Difference - Making faces
Info : [ 90%] Difference
Info : [ 0%] Fragments
Info : [ 10%] Fragments
Info : [ 20%] Fragments
Info : [ 30%] Fragments
Info : [ 40%] Fragments
Info : [ 50%] Fragments
Info : [ 60%] Fragments
Info : [ 70%] Fragments
Info : [ 80%] Fragments - Splitting faces
Info : Meshing 1D...
Info : [ 0%] Meshing curve 1 (Line)
Info : [ 0%] Meshing curve 4 (Line)
Info : [ 0%] Meshing curve 2 (Line)
Info : [ 0%] Meshing curve 3 (Line)
Info : [ 10%] Meshing curve 5 (Line)
Info : [ 10%] Meshing curve 8 (Line)
Info : [ 40%] Meshing curve 9 (Line)
Info : [ 10%] Meshing curve 6 (Line)
Info : [ 10%] Meshing curve 7 (Line)
Info : [ 40%] Meshing curve 10 (Line)
Info : [ 40%] Meshing curve 11 (Line)
Info : [ 40%] Meshing curve 12 (Line)
Info : [ 40%] Meshing curve 13 (Line)
Info : [ 40%] Meshing curve 14 (Line)
Info : Done meshing 1D (Wall 0.00828209s, CPU 0.01615s)
Info : Meshing 2D...
Info : [ 0%] Meshing surface 3 (Plane, Frontal-Delaunay)
Info : [ 0%] Meshing surface 2 (Plane, Frontal-Delaunay)
Info : [ 0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info : Done meshing 2D (Wall 0.0423437s, CPU 0.125454s)
Info : 3708 nodes 7607 elements
Info : Writing 'mesh.msh'...
Info : Done writing 'mesh.msh'
Info : Writing '/tmp/tmpx2fvxge7/mesh.msh'...
Info : Done writing '/tmp/tmpx2fvxge7/mesh.msh'
We can now throw this mesh into FEMWELL directly:
mesh = from_meshio(mesh_gmsh)
mesh.draw().show()
plot_domains(mesh)
plt.show()
Assign material values
basis0 = Basis(mesh, ElementTriP0())
epsilon = basis0.zeros()
for subdomain, n in {"core": 3.45, "box": 1.444, "clad": 1.444}.items():
epsilon[basis0.get_dofs(elements=subdomain)] = n**2
basis0.plot(epsilon, colorbar=True).show()
Solve for the modes:
wavelength = 1.55
modes = compute_modes(basis0, epsilon, wavelength=wavelength, num_modes=2, order=1)
You can use them as inputs to other femwell mode solver functions to inspect or analyze the modes:
print(modes[0].te_fraction)
0.9978478799141409
modes[0].show("E", part="real")
dir(modes[0])
[
'E',
'H',
'Sx',
'Sy',
'Sz',
'__annotations__',
'__class__',
'__dataclass_fields__',
'__dataclass_params__',
'__delattr__',
'__dict__',
'__dir__',
'__doc__',
'__eq__',
'__format__',
'__ge__',
'__getattribute__',
'__getstate__',
'__gt__',
'__hash__',
'__init__',
'__init_subclass__',
'__le__',
'__lt__',
'__match_args__',
'__module__',
'__ne__',
'__new__',
'__reduce__',
'__reduce_ex__',
'__repr__',
'__setattr__',
'__sizeof__',
'__str__',
'__subclasshook__',
'__weakref__',
'basis',
'basis_epsilon_r',
'calculate_confinement_factor',
'calculate_coupling_coefficient',
'calculate_effective_area',
'calculate_intensity',
'calculate_overlap',
'calculate_pertubated_neff',
'calculate_power',
'calculate_propagation_loss',
'epsilon_r',
'frequency',
'k',
'k0',
'n_eff',
'omega',
'plot',
'plot_component',
'plot_intensity',
'poynting',
'show',
'te_fraction',
'tm_fraction',
'transversality',
'wavelength'
]
modes[0].plot_component?
modes[0].plot_component("E", component="x", part="real", colorbar=True)
<Axes: title={'center': 'Ex (real. part)'}>
modes[1].plot_component("E", component="x", part="real", colorbar=True)
<Axes: title={'center': 'Ex (real. part)'}>
modes[1].show("E", part="real")
