Add symmetry planes for simulation to filter TE and TM for optics, or even and odd modes for RF

[duarte-jfs]

Adding the possibility to ]have symmetry planes to save on mesh cells to increase computational speed would be a great add-on when fine meshes are required. Perhaps I missed it on the documentation, but I could only spot examples with full 2D mesh even when symmetry is present.

[HelgeGehring]

Hey, @simbilod is currently working on an example for this.
You can just end the mesh at the symmetry plane and either use there a PMC(standard setting) or PEC (using metalic_boundaries) to have either even or odd modes.
Currently, you can already see it in the benchmark (yeah, we should add there more text, but they explain in the paper which sides are metalic :) ).

[duarte-jfs]

Indeed, you're right. That is a very simple solution and it seems to work. However, I tried to benchmark it with the rectangular waveguide and I get different results for the fundamental mode, which I don't know if I missed something or something else is going on. The code below gives an effective index of 1.871 for the fundamental TE mode when considering the full mesh and 1.899 when considering the symmetry plane, which is quite a significant difference. Did I miss something?

EDIT: Indeed I missed something. When I used the full mesh I define 'core = box(0,0,width/2, 1)', which is wrong. By replacing it with 'core = box(-width/2,0,width/2, 1)' it is solved and I now get the neff matching up to 1.899.

from collections import OrderedDict

import matplotlib.pyplot as plt
import numpy as np
import scipy.constants
import shapely
from shapely.ops import clip_by_rect
from shapely.geometry import LineString, box
from skfem import Basis, ElementTriP0, ElementTriP1
from skfem.io.meshio import from_meshio
from tqdm import tqdm

from femwell.maxwell.waveguide import compute_modes
from femwell.mesh import mesh_from_OrderedDict
from femwell.visualization import plot_domains


wavelength = 1.55 #um
num_modes = 8
width = 4 #um

#Find modes with full mesh
core = box(0,0,width/2, 1)
polygons = OrderedDict(
    core = core,
    box = clip_by_rect(core.buffer(2, resolution = 5),
                       -np.inf,
                      -np.inf,
                      np.inf,
                      0),
    clad = clip_by_rect(core.buffer(2, resolution = 5),
                       -np.inf,
                      0,
                      np.inf,
                      np.inf),

)

resolutions = {"core": {"resolution": 0.1, "distance": 1}}

mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.6))

basis0 = Basis(mesh, ElementTriP0())
epsilon = basis0.zeros(dtype = complex)

for subdomain, n in {"core": 1.9963, "box": 1.444, "clad": 1}.items():
    epsilon[basis0.get_dofs(elements = subdomain)] = n**2
    
modes = compute_modes(basis0, 
                      epsilon, 
                      wavelength=wavelength, 
                      num_modes=num_modes,
                      metallic_boundaries=True)

modes_full = modes
### Find modes with symmetry plane
polygons = OrderedDict(
    core = core,
    box = clip_by_rect(core.buffer(2, resolution = 5),
                       0,
                      -np.inf,
                      np.inf,
                      0),
    clad = clip_by_rect(core.buffer(2, resolution = 5),
                       0,
                      0,
                      np.inf,
                      np.inf),

)

resolutions = {"core": {"resolution": 0.1, "distance": 1}}

mesh = from_meshio(mesh_from_OrderedDict(polygons, resolutions, default_resolution_max=0.6))

basis0 = Basis(mesh, ElementTriP0())
epsilon = basis0.zeros(dtype = complex)

for subdomain, n in {"core": 1.9963, "box": 1.444, "clad": 1}.items():
    epsilon[basis0.get_dofs(elements = subdomain)] = n**2
    
modes = compute_modes(basis0, 
                      epsilon, 
                      wavelength=wavelength, 
                      num_modes=num_modes,
                      metallic_boundaries=True)

modes_sym = modes

print(modes_full.n_effs[0])
print(modes_sym.n_effs[0])

modes_full[0].plot_component('E', 'x')
modes_sym[0].plot_component('E', 'x')

Output:

(1.8716330817857192-6.7088318229796605e-18j)
(1.8991017992315486-3.272880415020128e-17j)

[simbilod]

Nice!

I am making an example with this, but it is a bit trivial, so I am also trying to add a the functionality to reflect the half-mode into a full mode so it is easier to take overlap with other modes