Elmer vs FEniCS - Part 1
Stefano Tronci
In the Writing Your Packages - Part 1 and Writing Your Packages - Part 2 episodes we introduced the most important concepts in writing your own simulation packages. We chosen Python as a language and implemented two simple packages, acoupy_meshutil and acoupy_helmholtz. In this episode we will use the Pulsating Sphere benchmark problem to compare acoupy_helmholtz to our old trustworthy Elmer. Since our software is little more than syntactic sugar around FEniCS this will allow us to compare with FEniCS directly.
Project Files
All the files used for this project are available at the repositories below:
The Software
As we said, we implemented acoupy_meshutil and acoupy_helmholtz. We did not cover the details of the development of these packages in a tutorial form. Instead, this effort was put into the documentation of the software itself. The documentation, accessible through the README, has examples and tutorials as well as minimal theoretical explanations. The code was also written to be self explanatory and you are encouraged to read it with the documentation. The simulate_and_compare.py script from the FEniCS VS Elmer 1 repository servers also as a good example on how to use the packages.
The Problem
As we explained in the The FEM Pipeline episode benchmark problems are very useful to understand the accuracy of numerical solutions. In this episode we will solve the Pulsating Sphere problem over the same mesh with both Elmer and FEniCS. This will allow us to compare them to each other as well to the analytical solution, which is known. The domain that we setup for this episode is very similar to that of The Pulsating Sphere episode. It is a spheric domain with a small sphere cut-off in the middle used to apply a Neumann boundary condition. All the parameters are in the tables below.
Domain and Mesh
| Parameter | Value |
|---|---|
| Domain Radius | $0.1$ $\text{m}$ |
| Source Radius | $5$ $\text{mm}$ |
| Mesh Max Size | $28$ $\text{mm}$ |
| Mesh Min Size | $1$ $\text{mm}$ |
| Meshing Algorithm | NETGEN |
| Mesh Order | $1$ |
| Element Type | Second Order $P$-element * |
* See Elmer Solver Manual | Appendix E and the FEniCS’s Periodic Table of the Finite Elements for details.
Material
| Parameter | Symbol | Value |
|---|---|---|
| Speed of Sound | $c$ | $343$ $\frac{\text{m}}{\text{s}}$ |
| Density | $\rho$ | $1.205$ $\frac{\text{kg}}{\text{m}^3}$ |
Source
| Parameter | Value |
|---|---|
| Radius | $5$ $\text{mm}$ |
| Normal Velocity | $\frac{3}{4} \exp \left( j \frac{\pi}{4} \right)$ $\frac{\text{m}}{\text{s}}$ |
Simulation
| Parameter | Value |
|---|---|
| Frequency | $1$ $\text{kHz}$ |
| Termination | Matched Impedance |
By matched impedance we mean a Robin boundary condition where the specific impedance $z$ of the boundary is matched to the outgoing spherical wave impedance.
Results
The Sound Pressure Level (SPL) in units of $\text{dBSPL}$ (i.e. $\text{dB re 20} \mu\text{Pa}$) and the phase of all solutions are plotted below. We plot SPL and phase as functions of distance from the centre of the source. This because distance is the only parameter of interest for a pulsating sphere field. This information is obtained by the solution at all the mesh nodes: we plot the field at a node as a function of the distance of that node from the source centre.
Click Here for Larger Version
We can see that the the FEM solutions both underestimate the SPL slightly, while there is very good agreement with phase. Elmer and FEniCS solutions are essentially identical.
We can make the analysis of agreement more quantitative by computing the error. As done in the Mesh Order and Accuracy episode, we compute the ratio between the FEM solution and the exact solution. Then, we study its decibel value as well as its phase. This will tell us the error in $\text{dB}$ as well as in phase for both the Elmer and FEniCS solutions. The results are shown below.
Click Here for Larger Version
It is possible to see that the SPL error is not as low as we typically like. In absolute terms, it surpasses $0.3$ $\text{dB}$ while we prefer to have all errors within $0.1$ $\text{dB}$.
Phase errors are within $0.003$ $\text{rad}$ ($0.17$ $\text{deg}$) within most of the domain. Errors reach $0.007$ $\text{rad}$ ($0.40$ $\text{deg}$) at the very end of the domain. Overall, the phase accuracy is very good, within fractions of a degree.
It can be noted that errors generally increase with distance. This is most likely due to error in the radius introduced by the meshing. This error is higher with distance as the mesh gets coarser (from which $1$ $\text{mm}$ near the source to $28$ $\text{mm}$ near the termination). This makes the impedance value used for the boundary condition Slightly inaccurate, affecting the whole field.
Overall, it is evident that Elmer and FEniCS have strikingly similar errors. In fact, the errors are almost identical, showing that errors are mainly driven by the mesh rather than the solvers.
Conclusion
We implemented packages to perform FEM simulations for the Helmholtz equation. We solved the Pulsating Sphere benchmark problem for both the FEniCS base package and Elmer. We did so by using the same mesh. The mesh size was on the coarser side ($\frac{1}{12}$ of the wavelength) and this resulted in relatively high SPL errors (in excess of $0.3$ $\text{dB}$). However, phase accuracy was good and, most importantly, the FEniCS and Elmer solutions were essentially identical. This shows that inaccuracy stems mainly from mesh size.
In the next episode we will study how well accuracy increases with mesh size reduction for both solvers.
License Information
This work is licensed under a Creative Commons Attribution 4.0 International License.