Adaptive computation of vector leaky modes in a Bragg fiber

2.3. Adaptive computation of vector leaky modes in a Bragg fiber#

This notebook demonstrates the adaptive finite element workflow for computing leaky vector (Maxwell) eigenmodes of a Bragg fiber. It extends the computation in Notebook 2.2 by replacing the fixed mesh with an adaptively refined one driven by a dual-weighted residual (DWR) error estimator.

The adaptive loop follows the classical Solve → Estimate → Mark → Refine cycle. At each step, a Maxwell leaky-mode eigenproblem is solved on the current mesh using a non-selfadjoint FEAST contour integral eigensolver, with a mixed Nédélec–Lagrange discretization and a smooth PML. The DWR estimator then identifies which elements contribute most to the eigenvalue error; only those are refined.

The method and the construction of the DWR error estimators are described in detail in [1].

import numpy as np
import ngsolve as ng
from ngsolve.webgui import Draw
from scipy.optimize import newton

from fibermode.bragg import BraggExactVector, Bragg

Fiber geometry#

The fiber consists of three concentric regions — an air core, a glass ring, and an air cladding — surrounded by a PML absorbing layer. All refractive indices are specified at the operating wavelength \(\lambda = 2.45\,\mu\text{m}\).

ts    = [4.0775e-05, 1e-5, 1e-5]       # layer thicknesses (m): core | glass | air
mats  = ['air', 'glass', 'air']
ns    = [1.00027717, 1.4388164768221814, 1.00027717]   # refractive indices
wl    = 2.45e-6                        # wavelength (m)
scale = 15e-6                          # characteristic length L (m) for nondimensionalization
maxhs = [.1, .1, .1]                   # initial mesh sizes (fraction of layer radius)

# PML additions appending an absorbing outer layer
ts_pml   = ts   + [5e-5]
mats_pml = mats + ['Outer']
ns_pml   = ns   + [ns[0]]
maxhs_pml = maxhs + [.1]

Numerical setup#

Bragg fiber geometry#

The Bragg class builds the NGSolve mesh on the truncated domain including the smooth PML layer (see Notebook 2.2 for a detailed description of its parameters). The refractive index field is drawn to confirm the geometry.

bragg_n = Bragg(ts=ts_pml, scale=scale, mats=mats_pml, maxhs=maxhs_pml, ns=ns_pml, wl=wl)

Draw(bragg_n.index, bragg_n.mesh, 'Refractive index');

Boundary names not provided, using default names.
ModeSolver: Checking if mesh has required regions
Mesh has  1814  elements,  938  points,  and  2751  edges.

Search region in the \(Z^2\)-plane#

The eigensolver searches in the nondimensional complex \(Z^2\)-plane, which is related to the physical propagation constant \(\beta\) by

\[ Z^2 = L^2\bigl(k_0^2 n_0^2 - \beta^2\bigr), \]

where \(k_0\) is the operating wavenumber and \(n_0\) is the refractive index at infinity (see the \(Z^2\)-plane section of Notebook 2.2 for further discussion). We center the search disk at a \(Z^2\) value near the expected location of the fundamental mode; the exact value is computed below using BraggExactVector to assess the numerical error.

# Disk in the Z² plane to search for eigenvalues:  
# (The exact value computed afterward semianalytically is within this disk.)

center = 0.78
radius = 0.1

Adaptive loop#

Each iteration performs:

  1. Solve — assemble the smooth-PML Maxwell system and run FEAST inside the \(Z^2\) contour.

  2. Estimate — compute the DWR error estimator \(\eta_T\) element-by-element.

  3. Mark — flag elements where \(\eta_T > \theta \cdot \max_T \eta_T\) (Dörfler marking, \(\theta = 0.1\)).

  4. Refine — bisect marked elements and re-curve the mesh.

These steps are implemented in the generator leakyvecmodes_adapt_gen, which yields the current state dictionary at each iteration. We use the generator directly here (rather than the fully automatic leakyvecmodes_adapt) so that we can inspect each intermediate result before triggering the next refinement. Both are methods of ModeSolver, the base class from which Bragg inherits. The autoupdate=True argument enables automatic extension of finite element spaces, memory reallocation, and prolongation of grid function data upon each mesh refinement.

# Initialize the stepper / generator for the adaptive loop. 

stepper = bragg_n.leakyvecmodes_adapt_gen(
    p=3,
    radius=radius,
    center=center,
    alpha=2,
    maxndofs=200000,
    autoupdate=True,
    verbose=False,
    npts=4,
    nspan=4,
    niterations=100,
    nrestarts=0)

This generator yields a dictionary containing the current state of the adaptive iterations, including the current computed eigenvalue and eigenfunction approximations, and the DWR error estimators.

Iteration 1#

state = next(stepper)

# Plot intensity of electric field & the DWR error estimator on the current mesh:

Draw(ng.Norm(state['uR'].gridfun(i=0).components[0])**2, bragg_n.mesh);
Draw(state['eevis']);  # plot the DWR error estimator on the current mesh

WARNING: kwarg 'nested' is an undocumented flags option for class <class 'ngsolve.comp.L2'>, maybe there is a typo?

ADAPTIVITY at  47405  ndofs:
  Assembling system...

  Computed eigenvalues: [0.81043041-0.00034947j 0.82212337-0.0005752j ]
   Normalizing by M-norms: max=7.5e+02, min=5.2e+02
   Normalizing by M-norms: max=7.4e+02, min=5.4e+02

  Error estimator: 0.0014224442121348168

  • The eigensolver produced two distinct but close-by eigenvalues (both printed out above). We will see as the iteration proceeds that they appear to merge into a near-degenerate pair representing the two polarizations of the HE₁₁ mode.

  • The computed (right) eigenfunctions corresponding to these two eigenvalues, stored in state['uR'], are returned in an object that wraps an NGSolve MultiVector object with further facilities for FEAST iterations.

  • From this “multieigenfunction” object uR, the ith eigenmode can be made into a GridFunction by state['uR'].gridfun(i). Here it gives the two eigenmodes (each containing three electric field components) for i=0 and i=1. For each i, the grid function state['uR'].gridfun(i) has two components:

    • the first gives the transverse electric field in the Nedelec space (state['uR'].gridfun(i).components[0])

    • the second gives the scaled longitudinal electric field component in the Lagrange space (state['uR'].gridfun(i).components[1])

  • Only the intensity of the transverse electric field of the first (i=0) mode is plotted above, as second looks similar (as you can verify by changing i=0 to i=1 in the above code).

  • Perhaps the most interesting plot is the DWR error estimator (plotted on the same mesh). It seems to indicate that the mesh is too coarse in the high-index glass layer. This may be surprising since the intensity of transverse electric field is almost negligible there. But this finding is consistent with the findings of [1, 2].

Iteration 2#

state = next(stepper)

# Plot intensity of electric field & the DWR error estimator on the current mesh:

Draw(ng.Norm(state['uR'].gridfun(i=0).components[0])**2, bragg_n.mesh);
Draw(state['eevis']);  # plot the DWR error estimator on the current mesh

  Marked  243  elements for refinement

ADAPTIVITY at  79576  ndofs:
  Assembling system...

  Computed eigenvalues: [0.79571044-0.00155303j 0.8368579 -0.00271257j]
   Normalizing by M-norms: max=1.5e+04, min=3.5e+02
   Normalizing by M-norms: max=2.3e+03, min=3.7e+02

  Error estimator: 0.00041007144795758056

Iteration 3#

state = next(stepper)

# Plot intensity of electric field &  DWR error estimator on the current mesh:

Draw(ng.Norm(state['uR'].gridfun(i=0).components[0])**2, bragg_n.mesh);
Draw(state['eevis']);  # plot the DWR error estimator on the current mesh

  Marked  582  elements for refinement

ADAPTIVITY at  157645  ndofs:
  Assembling system...

  Computed eigenvalues: [0.78726837-0.00089917j 0.78722105-0.00088223j]
   Normalizing by M-norms: max=1.0e+02, min=9.3e+01
   Normalizing by M-norms: max=1.2e+02, min=7.9e+01

  Error estimator: 0.00029651698989286725

The reason for the refinements#

Examining the longitudinal component closely reveals fine-scale oscillations in the glass layer — exactly what the DWR estimator was flagging for refinement. These oscillations are easy to miss if we look only at the total intensity of the electric field. They are the same physical feature noted in Notebook 2.2 and studied in detail in [2].

Draw(ng.Norm(state['uR'].gridfun(i=0).components[1])**2, bragg_n.mesh, 
     settings={"Objects": {"Wireframe": False}, "Colormap":{"autoscale": True, "ncolors": 16}});

The importance of capturing these fine scale oscillations when computing confinement losses were clearly shown in the paper [2].

Exact eigenvalue and numerical errors#

We use fibermode’s semi-analytical BraggExactVector class (see Notebook 2.1) to compute the exact propagation constant \(\beta\) of the fundamental leaky mode. We then map this \(\beta\) to the nondimensional exact \(Z^2\) eigenvalue and compare it with the adaptive iterates above.

bragg_e = BraggExactVector(ts=ts, scale=scale, mats=mats, ns=ns, wl=wl)

nu    = 1      # azimuthal mode number for vector fundamental mode (HE₁₁)
outer = 'h1'   # outgoing solution is the Hankel function of the first kind
k_low = bragg_e.k0 * bragg_e.ns[0] * bragg_e.scale   # lower scaled wavenumber bound

beta_exact = newton(bragg_e.determinant, np.array(.9999 * k_low),
                    args=(nu, outer), tol=1e-15)

print(f'Exact β (scaled)    = {beta_exact}')
print(f'Exact β physical    = {beta_exact/bragg_e.scale}')
print(f'|det| residual      = {abs(bragg_e.determinant(beta_exact, nu, outer)):.2e}')

exact_z2 = bragg_n.sqrZfrom(beta_exact / bragg_n.L)
print(f'Exact Z² = {exact_z2}')

Exact β (scaled)    = (38.46892217725964+1.1394349840994064e-05j)
Exact β physical    = (2564594.81181731+0.7596233227329376j)
|det| residual      = 1.25e-10
Exact Z² = (0.78653268033122-0.000876656714587343j)

Note that solving for the zero of the transfer-matrix determinant can be rather ill-conditioned, as seen by the determinant residual at the computed root, which is only 1.25e-10. Yet the result we obtained is sufficient for studying convergence as our discretization errors are larger. The table below shows the \(Z^2\) error against the exact transfer-matrix value and the DWR estimator at each refinement level.

Zsqrs        = state['Zsqrs']
errestimates = state['errestimates']
ndofs        = state['ndofs']

print(f'{"ndofs":>10}  {"Z² error":>14}  {"DWR estimator":>14}')
for n, z, (eta, _) in zip(ndofs[1:], Zsqrs, errestimates):
    err = abs(z[0] - exact_z2)
    print(f'{n:10d}  {err:14.6e}  {eta:14.6e}')

beta_num  = bragg_n.betafrom(Zsqrs[-1])[0]
beta_exact_phys = beta_exact / bragg_e.scale
print(f'\nNumerical β = {beta_num:.6e} m⁻¹')
print(f'Exact     β = {beta_exact_phys:.6e} m⁻¹')
print(f'|Δβ| / |β|  = {abs(beta_num - beta_exact_phys) / abs(beta_exact_phys):.2e}')

     ndofs        Z² error   DWR estimator
     47405    2.390354e-02    1.422444e-03
     79576    9.202652e-03    4.100714e-04
    157645    7.360375e-04    2.965170e-04

Numerical β = 2.564594e+06+7.791276e-01j m⁻¹
Exact     β = 2.564595e+06+7.596233e-01j m⁻¹
|Δβ| / |β|  = 2.49e-07

References#

[1] J. Gopalakrishnan, J. Grosek, G. Pinochet-Soto, and P. Vandenberge, “Adaptive resolution of fine scales in modes of microstructured optical fibers,” SIAM Journal on Scientific Computing (2025). Open Access, DOI: 10.1137/24M1651605

[2] P. Vandenberge, J. Gopalakrishnan, and J. Grosek, “Sensitivity of confinement losses in optical fibers to modeling approach,” Optics Express 31(16), 26735–26760 (2023). Open Access, DOI: 10.1364/OE.495467

Contents