# Elmer vs FEniCS - Part 1

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.