Hello everyone. I am new to Fenicsx. I practically copied the Modal Analysis of a Dielectric Waveguide into vscode and changed the domain functions to include absorbing pml, vacuum, and dielectric. I then deleted the conducting boundary condition. There happens to be something very wrong with the “eps” function which defines the value of the dielectric permittivity over the domain. As a test, I made all of the permittivities equal to 1, and still, the program seems to solve for eigenvalues that reside on the outline of my domains which look like nested squares. Needless to say there is some serious error here. When I instead said eps.x.array[:] = 1.0 # Assign eps = 1 for all cells in the domain, it worked. But when I set it equal to 1 in the code below, it does not work at all. The simulational results of the malfunctioning code will be also included as an image. What you are seeing is two squares separating three dielectric domains. The innermost domain is dielectric. The middle domain is vacuum. The outermost and thickest domain is absorbing pml. Of course my code has their permittivities set to 1 just for the demonstrative purpose that the error is not in the physics but with something else. Thanks for your help - Frank.
The tutorial I am referencing: Electromagnetic modal analysis for a waveguide
# imports
import sys
from mpi4py import MPI
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type, fem, io, plot
from dolfinx.fem.petsc import assemble_matrix
from dolfinx.mesh import CellType, create_rectangle, exterior_facet_indices, locate_entities
try:
import pyvista
have_pyvista = True
except ModuleNotFoundError:
print("pyvista and pyvistaqt are required to visualise the solution")
have_pyvista = False
try:
from slepc4py import SLEPc
except ModuleNotFoundError:
print("slepc4py is required for this demo")
sys.exit(0)
# geometry
l0 = 1 # free space wavelength
p = 4 # thickness of PML region
g = 1 # gap thickness
d = 2 # dielectric thickness
w = p + g + d # total width/height
h = w
nx = 300
ny = nx
msh = create_rectangle(
MPI.COMM_WORLD, np.array([[0, 0], [w, h]]), np.array([nx, ny]), CellType.quadrilateral
)
msh.topology.create_connectivity(msh.topology.dim - 1, msh.topology.dim)
eps_p = 1 # complex permittivity in PML
eps_v = 1 # vacuum permittivity
eps_d = 1 # dielectric permittivity
# Center coordinate
c = w / 2.0
def Omega_d(x):
"""
Dielectric region:
A central square of side d, centered at (c, c).
"""
return (
(x[0] >= c - d/2) & (x[0] <= c + d/2) &
(x[1] >= c - d/2) & (x[1] <= c + d/2)
)
def Omega_v(x):
"""
Vacuum region:
A square of side (d + g), centered at (c, c),
but excluding the smaller dielectric square.
"""
in_big_square = (
(x[0] >= c - (d + g)/2) & (x[0] <= c + (d + g)/2) &
(x[1] >= c - (d + g)/2) & (x[1] <= c + (d + g)/2)
)
# Vacuum region is the big square minus the dielectric
return in_big_square & (~Omega_d(x))
def Omega_p(x):
"""
PML region:
Everything in the domain [0, w] x [0, w] not in Omega_d or Omega_v.
"""
in_domain = (
(x[0] >= 0) & (x[0] <= w) &
(x[1] >= 0) & (x[1] <= h)
)
return in_domain & (~Omega_d(x)) & (~Omega_v(x))
D = fem.functionspace(msh, ("DQ", 0))
eps = fem.Function(D)
cells_p = locate_entities(msh, msh.topology.dim, Omega_p)
cells_v = locate_entities(msh, msh.topology.dim, Omega_v)
cells_d = locate_entities(msh, msh.topology.dim, Omega_d)
eps.x.array[cells_p] = np.full_like(cells_p, eps_p, dtype=default_scalar_type)
eps.x.array[cells_v] = np.full_like(cells_v, eps_v, dtype=default_scalar_type)
eps.x.array[cells_d] = np.full_like(cells_d, eps_d, dtype=default_scalar_type)
# definition
degree = 1
RTCE = element("RTCE", msh.basix_cell(), degree)
Q = element("Lagrange", msh.basix_cell(), degree)
V = fem.functionspace(msh, mixed_element([RTCE, Q]))
lmbd0 = 5*l0 # If this value is too large, you might not have supported modes!
k0 = 2 * np.pi / lmbd0
et, ez = ufl.TrialFunctions(V)
vt, vz = ufl.TestFunctions(V)
a_tt = (ufl.inner(ufl.curl(et), ufl.curl(vt)) - (k0**2) * eps * ufl.inner(et, vt)) * ufl.dx
b_tt = ufl.inner(et, vt) * ufl.dx
b_tz = ufl.inner(et, ufl.grad(vz)) * ufl.dx
b_zt = ufl.inner(ufl.grad(ez), vt) * ufl.dx
b_zz = (ufl.inner(ufl.grad(ez), ufl.grad(vz)) - (k0**2) * eps * ufl.inner(ez, vz)) * ufl.dx
a = fem.form(a_tt)
b = fem.form(b_tt + b_tz + b_zt + b_zz)
# solving with SLEPc
A = assemble_matrix(a)
A.assemble()
B = assemble_matrix(b)
B.assemble()
eps = SLEPc.EPS().create(msh.comm)
eps.setOperators(A, B)
eps.setProblemType(SLEPc.EPS.ProblemType.GNHEP)
tol = 1e-9
eps.setTolerances(tol=tol)
eps.setType(SLEPc.EPS.Type.KRYLOVSCHUR)
# Get ST context from eps
st = eps.getST()
# Set shift-and-invert transformation
st.setType(SLEPc.ST.Type.SINVERT)
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_REAL)
eps.setTarget(-((0.5 * k0) ** 2))
eps.setDimensions(nev=10)
eps.solve()
eps.view()
eps.errorView()
# Save the kz
vals = [(i, np.sqrt(-eps.getEigenvalue(i))) for i in range(eps.getConverged())]
# Sort kz by real part
vals.sort(key=lambda x: x[1].real)
eh = fem.Function(V)
kz_list = []
for i, kz in vals:
# Save eigenvector in eh
eps.getEigenpair(i, eh.x.petsc_vec)
# Compute error for i-th eigenvalue
error = eps.computeError(i, SLEPc.EPS.ErrorType.RELATIVE)
# Verify, save and visualize solution
if error < tol and np.isclose(kz.imag, 0, atol=tol): # This chooses the eigenvalues which are close to zero
kz_list.append(kz)
# I deleted the assert statement which checks against the analytical solutions.
print(f"eigenvalue: {-kz**2}")
print(f"kz: {kz}")
print(f"kz/k0: {kz / k0}")
eh.x.scatter_forward()
eth, ezh = eh.split()
eth = eh.sub(0).collapse()
ez = eh.sub(1).collapse()
# Transform eth, ezh into Et and Ez
eth.x.array[:] = eth.x.array[:] / kz
ezh.x.array[:] = ezh.x.array[:] * 1j
gdim = msh.geometry.dim
V_dg = fem.functionspace(msh, ("DQ", degree, (gdim,)))
Et_dg = fem.Function(V_dg)
Et_dg.interpolate(eth)
# Save solutions
with io.VTXWriter(msh.comm, f"sols/Et_{i}.bp", Et_dg) as f:
f.write(0.0)
with io.VTXWriter(msh.comm, f"sols/Ez_{i}.bp", ezh) as f:
f.write(0.0)
# Visualize solutions with Pyvista
if have_pyvista:
V_cells, V_types, V_x = plot.vtk_mesh(V_dg)
V_grid = pyvista.UnstructuredGrid(V_cells, V_types, V_x)
Et_values = np.zeros((V_x.shape[0], 3), dtype=np.float64)
Et_values[:, : msh.topology.dim] = Et_dg.x.array.reshape(
V_x.shape[0], msh.topology.dim
).real
V_grid.point_data["u"] = Et_values
plotter = pyvista.Plotter()
plotter.add_mesh(V_grid.copy(), show_edges=False)
plotter.view_xy()
plotter.link_views()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
pyvista.start_xvfb()
plotter.screenshot("Et.png", window_size=[400, 400])
