Elmer vs FEniCS - Part 2
Stefano Tronci
In the Elmer vs FEniCS - Part 1 episode we solved the Pulsating Sphere problem over a fairly coarse mesh. We solved the problem both with Elmer and acoupy_helmholtz, our FEniCS based Helmholtz solver. We found that Elmer and FEniCS produce the same solution. This solution compares well to the exact solution, but has some appreciable magnitude error. In this episode we will study how this error changes with mesh size, and whether Elmer and FEniCS keep on providing the same solution.
Project Files
All the files used for this project are available at the repositories below:
The Problem
We solve exactly the same problem outlined in Elmer vs FEniCS - Part 1 | The Problem, but with meshes of different resolution.
We use a finer mesh with $10$ $\text{mm}$ maximum size and a fine mesh with $5$ $\text{mm}$ maximum size. We will refer to these as 10 mm and 5 mm throughout.
Again, we solve the problem on all meshes with both Elmer and acoupy_helmholtz. Refer to simulate_and_compare.py from acoupy_helmholtz for reference on how to implement the simulation. Remember that acoupy_helmholtz (as well as acoupy_meshutil) has a documentation with tutorials and examples accessible from their README.md file.
Results
Solution
As done in Elmer vs FEniCS - Part 1 | Results, we take each mesh node and evaluate the solution at the node. We compute the distance of the node from the source centre and plot magnitude and phase of the solutions as functions of this distance. We do this because the Pulsating Sphere field depends only on distance. The solutions are shown below.
10 mm
Click Here for Larger Version
5 mm
Click Here for Larger Version
Discussion
We can see that both solutions agree fairly well with the expected solution. The agreement in phase is good in both cases. However, the solutions underestimate the field slightly, the 10 mm having the higher error as expected. We will do a more quantitative error analysis in the next section.
Error
As done in the Mesh Order and Accuracy episode, for each mesh and solver we compute the quotient between the solution and the exact field. By taking the decibel value of the magnitude of this quotient we will have the magnitude error in $\text{dB}$. The angle of this quotient is the phase error. The results are shown below.
10 mm
Click Here for Larger Version
5 mm
Click Here for Larger Version
Discussion
We can see that the errors for the phase are fairly low for both meshes, being within $\pm 0.017$ $\text{rad}$ ($\pm 1$ $\text{deg}$) for all solutions. However, it is also evident that the 5 mm solution has way lower errors.
Magnitude errors on the other hand are somewhat significant for the 10 mm mesh. Again, the errors breach the $\pm 0.1$ $\text{dB}$ boundary. In fact, in absolute terms, they reach $0.3$ $\text{dB}$. We observed a similar total spread of error in Elmer vs FEniCS - Part 1 | Results for a $28$ $\text{mm}$ mesh. However, the 10 mm mesh error is actually better behaved. This because the error pattern is relatively flat around a lower mean value. The magnitude errors for the 5 mm solution are instead very low, well within the $\pm 0.1$ $\text{dB}$ threshold.
As observed in in Elmer vs FEniCS - Part 1 | Results the error patterns are almost identical. This shows that the solvers are providing the same solution and the inaccuracies are rather due to the mesh itself.
Since there is no substantial difference between the Elmer error and the FEniCS error we can pick any of the solutions to further characterise the error as a function of mesh size. One good way to analyse error is by computing its statistics. As evident from the plots above the error is not really random. In this case the errors appear larger the closer we are to the source. However, error statistics are still useful to quickly gauge the overall trends of the error. In the plots below we report the Kernel Density Estimates of the error as computed with KDEpy. Kernel Density Estimation (KDE) is a technique to estimate the probability distribution of a random process. It is mathematically equivalent to a histogram with uniform bins passed through a smoothing filter. Wikipedia has a good overall introduction. We will not discuss the details of KDE in this episode. Rather, we use a KDE as an instrument to visualise spread and mean. We chosen Gaussian kernels and narrow bandwidth to show the trends without excessive bias.
Click Here for Larger Version
The plots show how, by reducing mesh size, the magnitude error reduces in spread as well as gaining an average value closer to $0$. On the other hand, the phase error appears to be centred more closely around $0$, and mainly reducing in spread.
A Note on Execution Time
Elmer solutions were very fast to compute due to having used the BiCGStabl iterative method. With FEniCS we used the MUMPS direct method, which is slower. FEniCS allows to select multiple solution methods as well. We will investigate those in the next episodes. Many of the methods are also available with Elmer. Hence, you can refer to the Elmer Model and Solver Parameters episode for an overall introduction.
Conclusion
Elmer and FEniCS output very well matched solutions for meshes of different density when the element parameters are kept consistent between the two solvers. This shows that, in terms of accuracy, acoupy_helmholtz essentially matches Elmer’s Helmholtz solver. acoupy_helmholtz has the same features as the Elmer Helmholtz solver, plus two additional material parameters to describe damping. This shows how versatile writing our own software is: we are now able to expand the acoupy_helmholtz solver further with source terms and adiabatic absorbers (for example). All by not affecting accuracy. However, acoupy_helmholtz solutions were computed with a direct MUMPS solver, which is slower with respect iterative solvers. We will explore iterative FEniCS solvers in the next episodes.
License Information
This work is licensed under a Creative Commons Attribution 4.0 International License.