TIDY3D
LEARNING CENTER

Lectures
Tutorials
Applications

Example Library / Mode Solver / Mode solver benchmark against Femwell

Mode solver benchmark against Femwell

Sign up now to simulate .ipynb

Tidy3D includes a fast local mode solver, ModeSimulation, that computes waveguide modes on the FDTD Yee grid with subpixel-averaged permittivity. Here we benchmark it against Femwell, an open-source finite element (FEM) mode solver, on four silicon-on-insulator structures at 1.55 µm.

We compare accuracy against the number of unknowns each solver assembles (the degrees of freedom): for Tidy3D, the number of grid cells in the mode plane times the two transverse field components; for Femwell, the modal FEM basis size (order-2 elements on the narrow edge coupler tip, order 1 elsewhere). This is only an approximate, machine-independent proxy for cost, it does not reflect actual runtime.

Note. The local mode solver uses subpixel averaging only when the tidy3d-extras package is installed, which this benchmark requires. Install it (matching your tidy3d version) with:

pip install "tidy3d[extras]"
from importlib.metadata import version

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.interpolate import UnivariateSpline

import tidy3d as td
from tidy3d.constants import C_0

assert version("tidy3d-extras") == version("tidy3d"), "tidy3d-extras must match tidy3d exactly"
td.config.simulation.use_local_subpixel = True

Tidy3D Models

Four silicon-on-insulator structures (silicon core n = 3.48, SiO₂ cladding n = 1.44), single mode at 1.55 µm.

WAVELENGTH_UM = 1.55
FREQ0 = C_0 / WAVELENGTH_UM

silicon = td.Medium(permittivity=3.48**2)
cladding = td.Medium(permittivity=1.44**2)

MODELS = {
    "simple_strip": dict(title="Simple Strip", width=0.50, thickness=0.22,
        domain=(4.3, 3.02), mode_plane=(4.3, 3.02), target_neff=2.6),
    "sidewall_angle": dict(title="Sidewall Angle", width=0.50, thickness=0.22, sidewall_angle_deg=10.0,
        domain=(4.3, 3.02), mode_plane=(4.3, 3.02), target_neff=2.6),
    "narrow_edge_coupler_tip": dict(title="Narrow Edge Coupler Tip", width=0.19, thickness=0.22,
        domain=(7.0, 5.0), mode_plane=(7.0, 5.0), target_neff=1.8),
    "waveguide_bend": dict(title="Waveguide Bend", width=0.50, thickness=0.22, bend_radius=5.0,
        domain=(7.22, 4.94), mode_plane=(7.22, 4.94), target_neff=2.45),
}

STEPS = list(range(10, 101, 10))


def build_geometry(model_key):
    m = MODELS[model_key]
    if "sidewall_angle_deg" in m:
        half_width = m["width"] / 2
        return td.PolySlab(
            vertices=[(-half_width, -2.0), (half_width, -2.0), (half_width, 2.0), (-half_width, 2.0)],
            axis=2, slab_bounds=(-m["thickness"] / 2, m["thickness"] / 2),
            sidewall_angle=np.deg2rad(m["sidewall_angle_deg"]),
        )
    return td.Box(center=(0, 0, 0), size=(m["width"], td.inf, m["thickness"]))


def build_mode_simulation(model_key, steps_per_wvl):
    m = MODELS[model_key]
    sim = td.Simulation(
        center=(0, 0, 0), size=(m["domain"][0], 4.0, m["domain"][1]), medium=cladding,
        structures=[td.Structure(geometry=build_geometry(model_key), medium=silicon)],
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=steps_per_wvl, wavelength=WAVELENGTH_UM),
        boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()),
        run_time=1e-12, subpixel=td.SubpixelSpec(dielectric=td.PolarizedAveraging()),
    )
    mode_spec_kwargs = dict(num_modes=1, target_neff=m["target_neff"])
    if "bend_radius" in m:
        mode_spec_kwargs.update(bend_radius=m["bend_radius"], bend_axis=1, num_pml=(12, 12), precision="double")
    return td.ModeSimulation.from_simulation(
        simulation=sim, plane=td.Box(center=(0, 0, 0), size=(m["mode_plane"][0], 0, m["mode_plane"][1])),
        mode_spec=td.ModeSpec(**mode_spec_kwargs), freqs=[FREQ0],
    )


# Narrow edge coupler tip: the weakly confined tip mode converges far faster with a mesh
# override refining a small box around the core plus snapping points on the tip corners,
# instead of a uniform grid. Sweep the override cell size dl (not min_steps_per_wvl).
NEC_TIP_HALF = (0.19 / 2, 0.22 / 2)
NEC_SNAP_POINTS = [(sx * NEC_TIP_HALF[0], 0.0, sz * NEC_TIP_HALF[1]) for sx in (-1, 1) for sz in (-1, 1)]
NEC_OVERRIDE_BOX = 0.5
NEC_BACKGROUND_STEPS = 10
NEC_DL_NM = [45, 30, 20, 14, 10, 7, 5, 3.6, 2.7, 2.0]


def build_nec_mode_simulation(dl_um):
    m = MODELS["narrow_edge_coupler_tip"]
    override = td.MeshOverrideStructure(
        geometry=td.Box(center=(0, 0, 0), size=(NEC_OVERRIDE_BOX, td.inf, NEC_OVERRIDE_BOX)),
        dl=(dl_um, None, dl_um),
    )
    sim = td.Simulation(
        center=(0, 0, 0), size=(m["domain"][0], 4.0, m["domain"][1]), medium=cladding,
        structures=[td.Structure(geometry=build_geometry("narrow_edge_coupler_tip"), medium=silicon)],
        grid_spec=td.GridSpec.auto(min_steps_per_wvl=NEC_BACKGROUND_STEPS, wavelength=WAVELENGTH_UM,
                                   override_structures=[override], snapping_points=NEC_SNAP_POINTS),
        boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()),
        run_time=1e-12, subpixel=td.SubpixelSpec(dielectric=td.PolarizedAveraging()),
    )
    return td.ModeSimulation.from_simulation(
        simulation=sim, plane=td.Box(center=(0, 0, 0), size=(m["mode_plane"][0], 0, m["mode_plane"][1])),
        mode_spec=td.ModeSpec(num_modes=1, target_neff=m["target_neff"]), freqs=[FREQ0],
    )

Geometries

Cross-section of the four structures at the mode plane, before running the sweeps.

fig, axes = plt.subplots(1, 4, figsize=(16, 3.5))
for ax, model_key in zip(axes, MODELS):
    if model_key == "narrow_edge_coupler_tip":
        mode_sim = build_nec_mode_simulation(5e-3)
    else:
        mode_sim = build_mode_simulation(model_key, 20)
    mode_sim.plot(y=0, ax=ax)
    ax.set_title(MODELS[model_key]["title"], fontsize=10)
fig.tight_layout()
plt.show()

Run the Tidy3D Sweeps

The number of unknowns for each point is the count of grid cells in the mode plane times the two transverse electric field components.

rows = []
for model_key in MODELS:
    if model_key == "narrow_edge_coupler_tip":
        for dl_nm in NEC_DL_NM:
            mode_sim = build_nec_mode_simulation(dl_nm * 1e-3)
            dofs = 2 * mode_sim.grid.num_cells[0] * mode_sim.grid.num_cells[2]
            n_eff = float(mode_sim.run_local().modes.n_eff.values.squeeze())
            rows.append(dict(model_key=model_key, dofs=dofs, n_eff=n_eff))
            print(f"{model_key:26s} dl={dl_nm:5.1f}nm  dofs={dofs:8d}  n_eff={n_eff:.6f}", flush=True)
    else:
        for steps_per_wvl in STEPS:
            mode_sim = build_mode_simulation(model_key, steps_per_wvl)
            dofs = 2 * mode_sim.grid.num_cells[0] * mode_sim.grid.num_cells[2]
            n_eff = float(mode_sim.run_local().modes.n_eff.values.squeeze())
            rows.append(dict(model_key=model_key, dofs=dofs, n_eff=n_eff))
            print(f"{model_key:26s} steps={steps_per_wvl:3d}  dofs={dofs:8d}  n_eff={n_eff:.6f}", flush=True)

tidy3d_df = pd.DataFrame(rows)
simple_strip               steps= 10  dofs=    3500  n_eff=2.438241
simple_strip               steps= 20  dofs=   12804  n_eff=2.444471
simple_strip               steps= 30  dofs=   27936  n_eff=2.446241
simple_strip               steps= 40  dofs=   48384  n_eff=2.447098
simple_strip               steps= 50  dofs=   75366  n_eff=2.447542
simple_strip               steps= 60  dofs=  107920  n_eff=2.447835
simple_strip               steps= 70  dofs=  145418  n_eff=2.448034
simple_strip               steps= 80  dofs=  189504  n_eff=2.448166
simple_strip               steps= 90  dofs=  239984  n_eff=2.448263
simple_strip               steps=100  dofs=  294532  n_eff=2.448339
sidewall_angle             steps= 10  dofs=    3570  n_eff=2.436385
sidewall_angle             steps= 20  dofs=   13068  n_eff=2.443754
sidewall_angle             steps= 30  dofs=   28130  n_eff=2.445848
sidewall_angle             steps= 40  dofs=   49408  n_eff=2.446876
sidewall_angle             steps= 50  dofs=   76002  n_eff=2.447447
sidewall_angle             steps= 60  dofs=  109060  n_eff=2.447765
sidewall_angle             steps= 70  dofs=  147186  n_eff=2.447968
sidewall_angle             steps= 80  dofs=  191016  n_eff=2.448100
sidewall_angle             steps= 90  dofs=  241682  n_eff=2.448186
sidewall_angle             steps=100  dofs=  297044  n_eff=2.448258
narrow_edge_coupler_tip    dl= 45.0nm  dofs=    9086  n_eff=1.491508
narrow_edge_coupler_tip    dl= 30.0nm  dofs=   10880  n_eff=1.487710
narrow_edge_coupler_tip    dl= 20.0nm  dofs=   14100  n_eff=1.486264
narrow_edge_coupler_tip    dl= 14.0nm  dofs=   18576  n_eff=1.485420
narrow_edge_coupler_tip    dl= 10.0nm  dofs=   25168  n_eff=1.484978
narrow_edge_coupler_tip    dl=  7.0nm  dofs=   36792  n_eff=1.484720
narrow_edge_coupler_tip    dl=  5.0nm  dofs=   54984  n_eff=1.484588
narrow_edge_coupler_tip    dl=  3.6nm  dofs=   85932  n_eff=1.484511
narrow_edge_coupler_tip    dl=  2.7nm  dofs=  130380  n_eff=1.484476
narrow_edge_coupler_tip    dl=  2.0nm  dofs=  206544  n_eff=1.484455
waveguide_bend             steps= 10  dofs=    8268  n_eff=2.432795
waveguide_bend             steps= 20  dofs=   30804  n_eff=2.442550
waveguide_bend             steps= 30  dofs=   66752  n_eff=2.444833
waveguide_bend             steps= 40  dofs=  118404  n_eff=2.446517
waveguide_bend             steps= 50  dofs=  184262  n_eff=2.447476
waveguide_bend             steps= 60  dofs=  264032  n_eff=2.448048
waveguide_bend             steps= 70  dofs=  358110  n_eff=2.448425
waveguide_bend             steps= 80  dofs=  466496  n_eff=2.448690
waveguide_bend             steps= 90  dofs=  591848  n_eff=2.448886
waveguide_bend             steps=100  dofs=  729144  n_eff=2.449033

Femwell Data

Femwell results, loaded from a single CSV (misc/femwell.csv): the effective index versus problem size for the four structures (dataset = neff_vs_dof) and the effective index versus wavelength for the group-index comparison (dataset = group_index). Femwell uses second-order elements on the narrow edge coupler tip and first order elsewhere.

Note on the comparison. This compares accuracy against problem size (the number of unknowns each solver assembles). It is an approximate, indicative comparison: the matrix dimension is not strictly equivalent between a structured FDFD grid (Tidy3D) and an unstructured FEM mesh (Femwell), and it does not correspond to solver runtime, which further depends on the machine and on sparse-solver fill-in. The curves show convergence trends per unknown, not an exact head-to-head cost.

femwell = pd.read_csv("misc/femwell.csv")
FEMWELL_DOF = {
    k: dict(dofs=list(g.dofs), n_eff=list(g.n_eff))
    for k, g in femwell[femwell.dataset == "neff_vs_dof"].groupby("model_key")
}

Effective Index Versus Problem Size (Degrees of Freedom)

fig, axes = plt.subplots(2, 2, figsize=(13, 9))

for ax, model_key in zip(axes.flatten(), MODELS):
    m = MODELS[model_key]
    td_df = tidy3d_df[tidy3d_df.model_key == model_key].sort_values("dofs")
    fw = FEMWELL_DOF[model_key]

    ax.plot(td_df.dofs, td_df.n_eff, marker="o", color="#1f77b4",
            linewidth=1.5, markersize=6, label="Tidy3D")
    ax.plot(fw["dofs"], fw["n_eff"], marker="s", color="#d62728",
            linewidth=1.5, markersize=5, alpha=0.85, label="Femwell")

    ax.set_xscale("log")
    ax.set_title(m["title"], fontsize=11, fontweight="bold")
    ax.set_xlabel("Degrees of freedom (matrix dimension)")
    ax.set_ylabel("$n_\\mathrm{eff}$")
    ax.legend(fontsize=8)
    ax.grid(True, which="major", linestyle="--", alpha=0.5)
    ax.grid(True, which="minor", linestyle=":", alpha=0.3)

fig.suptitle("n_eff vs problem size (matrix dimension), Tidy3D vs Femwell", fontsize=13)
fig.tight_layout()
plt.show()

Note on the waveguide bend

The residual offset between the two curves in the waveguide-bend panel is not a convergence issue: both solvers are converged well beyond it. It reflects a difference in how the curvature is modeled. Femwell uses the widely adopted first-order equivalent-index transformation, \(n_{eq}(x) = n(x)\,(1 + x/R)\), as described in its documentation. The Tidy3D mode solver applies the full cylindrical-coordinate transformation of Maxwell’s equations, which accounts for the curvature to all orders in \(x/R\). The two treatments agree in the straight limit and approach each other as the bend radius grows (the difference scales as \((w/R)^2\)); at \(R = 5\) µm it amounts to a few \(10^{-3}\) in \(n_\mathrm{eff}\). The Tidy3D result is therefore the more accurate of the two, since it is a higher-order treatment of the curvature, retaining all orders in \(x/R\) where the equivalent-index transformation truncates at first order.

Group Index and GVD

We now compare the two solvers on a \(\text{Si}_3\text{N}_4\) strip waveguide (0.88 × 0.69 µm) embedded in \(\text{SiO}_2\), over 50 wavelengths from 500 nm to 2500 nm. Both the core and the cladding carry Sellmeier material dispersion, so the group index and GVD reported below are the total waveguide values: the material contribution of the two media plus the geometric contribution from the confinement, which is what the mode solver returns and what a real device sees.

  • Tidy3D: a single ModeSimulation with all 50 frequencies and group_index_step=True, so the group index and dispersion come natively from the solver. The ModeSortSpec with fill_fraction_box keeps the solver locked on the fundamental mode across the whole band.
  • Femwell: effective index versus wavelength from the CSV file; the group index and GVD follow from analytic derivatives of an interpolating quintic spline, which is the approach recommended in the Femwell documentation.
fw_gvd = femwell[femwell.dataset == "group_index"].sort_values("wavelength_nm")

gvd_width, gvd_height = 0.88, 0.69
gvd_domain_x, gvd_domain_z = 4 * gvd_width, 5 * gvd_height

si3n4 = td.Sellmeier(coeffs=[(3.0249, 0.1353406**2), (40314, 1239.842**2)])
sio2 = td.Sellmeier(
    coeffs=[
        (0.6961663, 0.0684043**2),
        (0.4079426, 0.1162414**2),
        (0.8974794, 9.896161**2),
    ]
)

wavelength_nm = np.linspace(500, 2500, 50)
gvd_freqs = list(C_0 / (wavelength_nm / 1000))

gvd_core_box = td.Box(center=(0, 0, 0), size=(gvd_width, 0, gvd_height))

gvd_sim = td.Simulation(
    center=(0, 0, 0),
    size=(gvd_domain_x, 4.0, gvd_domain_z),
    medium=sio2,
    structures=[
        td.Structure(
            geometry=td.Box(center=(0, 0, 0), size=(gvd_width, td.inf, gvd_height)),
            medium=si3n4,
        )
    ],
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=25, wavelength=0.8),
    boundary_spec=td.BoundarySpec.all_sides(boundary=td.Periodic()),
    run_time=1e-12,
)

gvd_mode_sim = td.ModeSimulation.from_simulation(
    simulation=gvd_sim,
    plane=td.Box(center=(0, 0, 0), size=(gvd_domain_x - 0.2, 0, gvd_domain_z - 0.2)),
    mode_spec=td.ModeSpec(
        num_modes=1,
        target_neff=2.0,
        sort_spec=td.ModeSortSpec(
            sort_key="fill_fraction_box",
            sort_order="descending",
            bounding_box=gvd_core_box,
            track_freq="central",
        ),
        group_index_step=True,
    ),
    freqs=gvd_freqs,
)

gvd_data = gvd_mode_sim.run_local()
18:22:30 UTC WARNING: Mode '0' appears to undergo a discontinuous change between
             frequencies '557104890347547.06' and '512855483852280.7' (overlap: 
             '0.00').                                                           

n_eff_td = np.real(gvd_data.modes.n_eff.isel(mode_index=0).values)
n_group_td = np.real(gvd_data.modes.n_group.isel(mode_index=0).values)
dispersion_td = np.real(gvd_data.modes.dispersion.isel(mode_index=0).values)

spl = UnivariateSpline(fw_gvd.wavelength_nm, fw_gvd.n_eff, s=0, k=5)
wl_fine = np.linspace(wavelength_nm[0], wavelength_nm[-1], 1000)
n_group_fw = spl(wl_fine) - wl_fine * spl.derivative(1)(wl_fine)
gvd_fw = -wl_fine / 2.99792e-7 * spl.derivative(2)(wl_fine)
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

axes[0].plot(wavelength_nm, n_eff_td, label="Tidy3D")
axes[0].plot(fw_gvd.wavelength_nm, fw_gvd.n_eff, label="Femwell", linestyle="--")
axes[0].set_ylabel("$n_\\mathrm{eff}$")
axes[0].set_title("Effective Index")

axes[1].plot(wavelength_nm, n_group_td, label="Tidy3D")
axes[1].plot(wl_fine, n_group_fw, label="Femwell", linestyle="--")
axes[1].set_ylabel("$n_g$")
axes[1].set_title("Group Index")

axes[2].plot(wavelength_nm, dispersion_td, label="Tidy3D")
axes[2].plot(wl_fine, gvd_fw, label="Femwell", linestyle="--")
axes[2].axhline(0, color="k", linestyle="--", linewidth=0.8)
axes[2].set_ylabel("D (ps / (nm $\\cdot$ km))")
axes[2].set_title("Group Velocity Dispersion")

for ax in axes:
    ax.set_xlabel("Wavelength (nm)")
    ax.legend()

plt.tight_layout()
plt.show()

Example library
Previous Tutorial Next Tutorial