In this note I present a series of benchmark problems to ensure that
the Maxwell solver in GkeyllZero (G0) works on non-Cartesian
geometries. G0 allows an arbitrary mapping from computational to
physical space. Here I assume that the physical domain is a coaxial
cylinder with inner and outer radii \(r_0\) and \(r_1\) and
height \(L_z\). In this geometry I study two sets of problems:
transverse magnetic (TM) modes and transverse electric (TE) modes. For
a derivation of these solutions see this document.
Warning
This note is incomplete: I still need to test the TE modes, and
convergence of the scheme has not been numerically verified. For
the latter, we need to store cell volumes so the appropriate
integrals can be computed.
Transverse Magnetic modes are modes in which the magnetic field is
perpendicular to the direction of propagation. In this case, the
electric field is in the \(Z\)-direction and is given by
where \(m\) is the azimuthal mode number, \(k_n = 2\pi/L_z\)
describes the axial mode, \(\omega\) is the frequency. The radial
dependence of the field is given by
\[E_z(r) = a J_m\big[r\sqrt{\omega^2-k_n^2}\big] + b Y_m\big[r\sqrt{\omega^2-k_n^2}\big]\]
where \(J_m\) and \(Y_m\) are Bessel functions of the first
and second kind, and \(a\) and \(b\) are constants. The
boundary conditions \(E_z(r_0) = E_z(r_1) = 0\) determine the
possible frequencies as roots of the equation
Once we compute the desired root (for propagating modes we must have
\(\omega>k_n\)) we can choose set \(a\) to whatever we like
and then compute \(b\) from
\[b = -a \frac{J_m\big[r_0\sqrt{\omega^2-k_n^2}\big]}{Y_m\big[r_0\sqrt{\omega^2-k_n^2}\big]}.\]
We can compute the frequencies in a few lines of Maxima. In general,
it is not possible to “guess” the interval in which the roots lie and
so it is best to first make a plot to determine the interval. Then
Maxima’s root finder can compute the bracketed root to machine
precision. The function we need to solve can be written in Maxima as
Plot of the nonlinear function whose roots (zero crossings) are the
allowed frequencies. Maxima root-finder requires we find the
interval in which the root is desired. We also need to ensure that
the function changes sign only once in the interval.
Using this figure we can choose the interval \([1,2]\) and find
the root as
w1 : find_root( F(2,0,2,5,w), w, 1, 2 )$
This will yield \(1.19318673737701\). We can also find
higher-frequency roots by passing other intervals to the above
command. Once we have the frequency we can determine \(a\) and
\(b\) as described above, thus completing the solution.
First consider the case in which \(k_n = 0\). This a 2D standing
mode inside an annular disk (i.e. there is no variation in the
\(Z\)-direction). We will choose \(r_0 = 2\) and \(r_1 =
5\) and \(m=4\). For this the first two roots are \(\omega =
1.557919724821651\) and \(\omega = 2.430327042902498\).
The following plot shows the solution at the constant radius
\(r=3.5\) of \(E_z\) at \(t=0\) and \(t=2 T_0\), where
\(T_0 = 2\pi/\omega\) on a \(36\times 144\) grid.
\(E_z\) at \(t=0\) (blue) and \(t=2 T_0\) (orange) for \(m=2\),
\(k_n=0\) first harmonic mode for a \(36\times 144\)
grid. See 2d-tm-w1-1 for simulation input.
The following plot shows the solution at \(t=0\) and \(t=2
T_0\) is the mode period, for the second harmonic mode.
\(E_z\) at \(t=0\) (left) and \(t=2 T_0\) (right) for \(m=2\),
\(k_n=0\) second harmonic mode for a \(64\times 256\)
grid. See 2d-tm-w2-2 for simulation input.
The \(m=0\) modes are axisymmetric. In Gkeyll, axisymmetric
simulations must be done on a 3D \((r,\phi,z)\) domain with
“wedge” periodic BCs in the \(\phi\) direction. These BCs ensure
that the solution remains axisymmetric and also that the geometric
effects of non-rectangular grid are taken into account. Due to the
present limitations of the Gkeyll solver, one needs to use at least 2
cells in the \(\phi\) direction.
The following figure shows the \(E_z\) field at \(t=0\) and
\(t=2 T_0\) for an axisymmetric mode with \(n=1\),
\(L_z=5\), and the second harmonic frequency. A coarse resolution
of \(16\times 32\) in the \(RZ\) plane is used, with \(3\)
cells in \(\phi\) in a narrow wedge of \(1/5\) radians.
\(E_z\) at \(t=0\) (left) and \(t=2 T_0\) (right) for \(m=0\),
\(k_n=2\pi/5\) second harmonic mode for a \(16\times 32\)
grid. See 2d-rz-tm-1 for simulation input.
Finally, I look at a 3D TM mode. I consider the case in which
\(k_n = 2\pi/5\), \(r_0 = 2\), \(r_1 = 5\) and
\(m=4\). The second harmonic frequency is
\(2.73598725136604\).
The following figures show the solution in the \(r,\phi\) and
\(r,z\) plane for \(E_z\).
\(E_z\) at \(t=0\) (left) and \(t=2 T_0\) (right) in the
\(r,\phi\) plane at \(z=2.5\) for
\(m=4\), \(k_n=2\pi/5\) second harmonic mode on a
\(48\times 96\times 96\) grid. See 3d-tm-1 for simulation input.
\(E_z\) at \(t=0\) (left) and \(t=2 T_0\) (right) in the
\(r,z\) plane at \(\phi=\pi\) for
for
\(m=4\), \(k_n=2\pi/5\) second harmonic mode on a
\(48\times 96\times 96\) grid. See 3d-tm-1 for simulation input.
The geometry implementation for Maxwell equation seems to be working
fine. At present, there are some limitations in performing RZ
simulations as at least 2 cells need to be used in the \(\phi\)
direction. Also, the TE modes are not yet tested, though there is no
reason why they should not just work.
