Solving Maxwell's Eigenvalue Problem via Isogeometric Boundary Elements and a Contour Integral Method
Numerical Examples
In the following we will discuss some numerical experiments showcasing the application of the contour integral method to the isogeometric boundary element method. The particular implementation used is Bembel, which is available open-source . A branch containing the code to recreate all of the presented numerical examples is also available .
Since no quasi-optimal preconditioners for the isogeometric discretisation of the electric field integral equation are known, iterative solvers yield suboptimal runtimes. Thus, the following numerical experiments will use a dense matrix assembly together with the partially pivoted LU decomposition of the linear algebra library Eigen to solve the arising systems. The utilised higher-order approaches yield systems small enough for this approach to be more than feasible.
As a first example, we investigate the first two eigenvalues of the unit cube, where an analytical solution is given by the eigenvalues $`\lambda_{\mathrm{ana},0}=\pi\sqrt{2}`$ of multiplicity three and $`\lambda_{\mathrm{ana},1}=\pi\sqrt{3}`$ of multiplicity two. The ellipse was defined as
\sin(t) + i\cdot0.05\cdot\cos(t)+5.0,\quad\text{ for }t=[0,2\pi),
where the contour integrals were evaluated, again, with $`N=25.`$ The error visualised corresponds to
\mathrm{error} = \frac 15\sum_{0\leq j<5}\min_{i=0,1}\vert\lambda_j-\lambda_{\mathrm{ana},i}\vert.
As one can see in Figure 1, the multiplicity of the eigenvalues is reflected perfectly by the non-zero singular values, i.e., all eigenvalues have been recognised. Moreover, we have the theoretically obtainable convergence order of $`{\mathcal O}(h^{2p+1})`$.
Comparison to IGA-FEM and Commercial Tools
As a second example, we compute the first eigenvalue for the sphere, cf. Figure 2 for the results. The contour was chosen as the curve
0.5\cdot\sin(t) + i\cdot0.05\cdot\cos(t)+2.5,\quad\text{ for }t=[0,2\pi).
The trapezoidal rule for lines 5 and 15 of the algorithm was chosen with $`N=25`$. A closed-form solution $`\lambda_{\mathrm{ana}}`$ is known and given as the first root of the spherical Bessel function of the first kind, cf. Figure 3. It is an eigenvalue of multiplicity three, thus three non-zero singular values in $`\bb \Sigma`$ are expected. This behaviour is reflected by the numerical examples perfectly. The error in Figure 2 refers to the average error of all three computed eigenvalues, i.e.,
\mathrm{error} = \frac 13\sum_{0\leq j < 3} \vert \lambda_j - \lambda_{\mathrm{ana}}\vert.
The convergence behaviour w.r.t. $`h`$ matches the orders $`{\mathcal O}(h^{2p+1})`$ predicted by [ErrorEstimate3] once more.
We have also computed approximations of the smallest eigenvalues of the unit sphere with the volume-based IGA software package GeoPDEs 3.0 . In Table 1 we showcase results for different polynomial degrees and refinements the number of degrees of freedom and the reached accuracy for the IGA-BEM and the volume-based IGA. One can observe that for the IGA-BEM a notably fewer (at least $`20\times`$) number of degrees of freedom are necessary for each polynomial degree in order to reach the same accuracy as for the volume-based IGA. The difference to commercial tools, in this case CST Microwave Studio 2018, is even more pronounced, cf. Table 2. However, the matrices for the IGA-BEM are dense and the eigenvalue problem non-linear.
| IGA-BEM | volume-based IGA | ||||||
| p | m | DOFs | error | p | SD | DOFs | error |
| 3 | 1 | 192 | 2.66e-05 | 3 | 4 | 4350 | 3.345e-05 |
| 3 | 2 | 432 | 4.63e-07 | 3 | 8 | 20450 | 3.45e-07 |
| 3 | 3 | 1200 | 2.12e-09 | 3 | 14 | 84560 | 1.10e-08 |
| 4 | 1 | 300 | 1.62e-06 | 4 | 4 | 6944 | 2.36e-06 |
| 4 | 2 | 588 | 3.05e-08 | 4 | 8 | 27280 | 3.94e-09 |
| 4 | 3 | 1452 | 3.89e-11 | 4 | 13 | 84560 | 7.88e-11 |
| 5 | 1 | 432 | 8.97e-08 | 5 | 4 | 10408 | 1.88e-07 |
| 5 | 2 | 768 | 2.04e-09 | 5 | 8 | 35484 | 7.33e-11 |
| 5 | 3 | 1728 | 7.17e-12 | 5 | 11 | 69600 | 2.1498e-12 |
| CST Microwave Studio 2018 | |||
| elements | SD | DOFs | error |
| 1089 | 6 | 21597 | 1.8638e-07 |
| 5292 | 10 | 101997 | 5.4204e-09 |
| 9191 | 12 | 175818 | 1.9524e-09 |
An Industrial Application: TESLA Cavity
As a third example we discuss an eigenvalue computation of the one-cell TESLA geometry as shown in Figures [fig::intro::tesla1] and [fig::intro::tesla2], with the results depicted in Figure 4. For this example no analytical solution is known, but experience dictates a resonant frequency around $`\kappa \approx 26.5`$. We choose the contour
0.5\cdot\sin(t) + i\cdot0.01\cdot\cos(t)+26.5,\quad\text{ for }t=[0,2\pi)
and $`N=8.`$ As a reference solution we utilise the result of a computation with $`p=5`$ and $`h=1/8.`$ This reference solution was compared to a computation with CST Microwave Studio 2018. The set solver order was 3rd (constant) and the mesh was generated with 200 771 curved tetrahedral elements of order 5. CST yields the solution of 1.27666401260 GHz. Our reference solution of $`\lambda_{\mathrm{cim}} = 26.75690023`$ corresponds to a frequency of 1.276664064 GHz. This results in a relative error of 4.018e-08. Thus our experiments are in good agreement with those of the CST software. However, note that in Figure 4 one can clearly see stagnation w.r.t. the CST Solution on higher-order computations and small $`h`$. This suggests that the Bembel reference solution is more accurate, provided the convergence of the contour integral approach behaves as observed in the previous numerical experiments.
The order of convergence for the cavity is not as pronounced as in the examples with analytical solution, which, in light of the estimate [ErrorEstimate3] and Remark [rem::nonsmoothGeom] is most likely due to the reduced regularity of the corresponding eigenfunction. Either way, one can clearly see an increased accuracy in higher-order approaches.