I am trying to find the eigenmodes of a cavity. In chapter 34 of the FEniCS book (by Lezar and Davidson) there's an example on how to calculate the cutoff-frequencies of a waveguide.
At first this seemed very similar to me, but it seems to be very different for several reasons:
In Ch.34 they consider purely transversal modes, i.e. the e-field is inplane and the out-of-plane component is entirely in the h-field. This makes the separation of the DE very simple.
My case is somewhat different:
I am looking at a cavity of cylindrical symmetry, i.e. rotational symmetry in phi but arbitrary geometry in r and z. I won't dive into special symmetry directly however. The code below is supposed to use a 3D cube instead.
My approach is to define a combined space:
VFS = VectorFunctionSpace(mesh, "CG", 2) # in-plane (r,z)
FS = FunctionSpace(mesh, "DG", 2) # phi
combined_space = VFS * FS
The weak form of the Maxwell equations (assuming PEC boundaries and a closed domain) is
$$ \frac{1}{\mu_0} \int_\Omega (\nabla \times \vec E_u) \cdot (\nabla \times \vec E_v) dx = \omega^2 \int_\Omega \varepsilon \vec E_u \vec E_v dx $$
($E_u$ being the trial function and $E_v$ the test function) and now I run into several problems:
1.) If I run the program for a box (3D but cartesian), it produces incredibly rough and unphysical modes ... and a lot of them. I suspect this is because I'm not using a curl-conforming VectorFunctionSpace. Is it?
2.) Defining boundary conditions is not necessary, since they are already implicitly included in the weak formulation. Adding them in the case described above (CG and DG Function Spaces) doesn't change anything (just as I expected).
Minimal code example
from dolfin import *
import numpy as np
import scipy.constants as co
(a,b,c) = (5e-3, 6e-3, 7e-3) # edge lengths of resonator box
(nx,ny,nz) = (5,5,2)
mesh = BoxMesh(Point(0,0,0), Point(a,b,c), nx, ny, nz)
VFS = VectorFunctionSpace(mesh, "CG", vector_order) # in-plane (r,z)
FS = FunctionSpace(mesh, "DG", scalar_order) # phi
combined_space = VFS * FS
(E_rz_0, E_phi_0) = TrialFunctions(combined_space)
(E_rz_1, E_phi_1) = TestFunctions(combined_space)
S_ij = inner(curl(E0), curl(E1)) *dx # LHS * mu0
T_ij = eps_r*inner(E0, E1) *dx # RHS / eps0
# hence, eigenvalue is k**2 = omega**2/(eps0*mu0)
S = PETScMatrix()
T = PETScMatrix()
assemble(S_ij, tensor=S)
assemble(T_ij, tensor=T)
eigensolver = SLEPcEigenSolver(S, T)
eigensolver.solve()
if eigensolver.get_number_converged() > 0:
print("%d results obtained..." % eigensolver.get_number_converged())
for i in range(eigensolver.get_number_converged()):
r, c, rx, cx = eigensolver.get_eigenpair(i)
print("%4d: %10.10f" % (i, r))
else:
print("No results obtained, eigensolver.get_number_converged() returned 0. Exiting...")
sys.exit(0)
To look at the results (here, modes 15 and 19), I've got this thing:
for i in [15, 19]: # or range(395,402)
if i > eigensolver.get_number_converged():
continue
r, c, rx, cx = eigensolver.get_eigenpair(i)
# Initialize function and assign eigenvector
u = Function(VFS)
u.vector()[:] = rx
norm_u = project(sqrt(inner(u, u)))
# Plot eigenfunction
plot(norm_u)
interactive()
