In the Dealing with Convergence Issues we made use of first order meshes in order to ease convergence of our simulation at high frequency. However, the accuracy of FEM solutions is higher the higher the order of the mesh, so doing so will come at the expenses of accuracy. Still, we argued that the accuracy is mostly controlled by the mesh size, so as long as we have more than ten elements per wavelength the solution should be reasonably accurate. In this episode we will use one of our benchmark models, The Pulsating Sphere, to quantify the impact of reducing mesh order on the accuracy of a FEM solution. We will also look at the effect of parallelisation. We will discover that in reality, although it is possible to reach good accuracy with first order meshes, second order meshes allow the Helmholtz Solver to provide significantly more accurate results.

Project Files

All the files used for this episode are available at the repositories below:

• Pulsating Sphere 1.
• Pulsating Sphere 2.
• Pulsating Sphere 3.
• Pulsating Sphere 4.
• Pulsating Sphere 5.
• Pulsating Sphere Study, which contains most of the code to compare the simulations set in the repositories listed above.

The Method

As we seen in The FEM Pipeline episode it is customary, when developing numerical solutions, to use a benchmark models to assess the accuracy of our numerical modelling setup. Once our modelling setup is well tuned on the benchmark model, which we can do as we know the exact solution for it, then we are safe to apply the setup to more complex problems. Here we will solve The Pulsating Sphere model with various mesh configurations and compare with the exact solution in order to assess how accuracy varies with mesh parameters.

Comparing in Practice

In The Pulsating Sphere we compared the field only along a radial line. In this episode we will take it further and produce comparison fields defined all over the domain.

Elmer outputs vtu files for each of its solutions. These files contain the Cartesian coordinates of all of our mesh nodes and the solution field(s) at those points. As we introduced in the FEM in a Nutshell this is all is needed, as the field in any other position in our domain will be computed (interpolated) in terms of that at the nodes by using a basis of functions whose complexity is related to the mesh order (typically polynomials).

So, in order to build a comparison field, all we have to do is take those coordinates, use them to compute the exact field, and produce a “comparison field” between the exact solution and the FEM solution. We can then save this comparison field as a vtu file itself, so we can visualize it with ParaView.

In order to load and manipulate the vtu files we can use the python package meshio. Thanks to PyCall it is very easy to use this package through Julia. For more details, you can read the code used to implement this, which is available here. Additionally, an intro to meshio is presented in the Intro to meshio.

As for the comparison field definition, let’s denote with $P_{\text{FEM}}\left(x,y,z\right)$ the FEM field and with $P_{\text{Exact}}\left(x,y,z\right)$ the exact field, with $x$, $y$ and $y$ the Cartesian coordinates defined over our domain (remember that these fields are solutions of the The Helmholtz Model so they are complex fields dependent upon the spatial variables only). Then, a natural way to define the comparison field $C\left(x,y,z\right)$ is as follows:

$$C\left(x,y,z\right) = \frac{P_{\text{FEM}}\left(x,y,z\right)}{P_{\text{Exact}}\left(x,y,z\right)}$$

In the ideal case in which the FEM solution is the same as the exact one this field would be $1$ everywhere (note that the comparison field is adimensional). Otherwise, it will be a different complex number, varying with position. We can check both magnitude and phase of this complex field, this will tell us the magnitude error (as a ratio) and the phase error (as a difference). In fact we can write:

$$P_{\text{FEM}}\left(x,y,z\right)=M_{\text{FEM}}\left(x,y,z\right)\exp\left(j\Phi_{\text{FEM}}\left(x,y,z\right)\right)$$ $$P_{\text{Exact}}\left(x,y,z\right)=M_{\text{Exact}}\left(x,y,z\right)\exp\left(j\Phi_{\text{Exact}}\left(x,y,z\right)\right)$$

where $M$ and $\Phi$ denote magnitude and phase respectively. Hence, thanks to the properties of the $\exp$ function:

$$C\left(x,y,z\right) = \frac{P_{\text{FEM}}\left(x,y,z\right)}{P_{\text{Exact}}\left(x,y,z\right)} = \frac{M_{\text{FEM}}\left(x,y,z\right)}{M_{\text{Exact}}\left(x,y,z\right)}\exp\left(j\left[\Phi_{\text{FEM}}\left(x,y,z\right) - \Phi_{\text{Exact}}\left(x,y,z\right) \right]\right)$$

Tested Meshes

The table below reports the mesh parameters used for the test.

Mesh ID	Algorithm	Min. Size $\left[\text{mm}\right]$	Max. Size $\left[\text{mm}\right]$	Element Order	Parallel Jobs	Project Files
1	NETGEN 1D-2D-3D	$1$	$29$	$2$	$1$	click here
2	NETGEN 1D-2D-3D	$1$	$5$	$2$	$1$	click here
3	NETGEN 1D-2D-3D	$1$	$29$	$1$	$1$	click here
4	NETGEN 1D-2D-3D	$1$	$5$	$1$	$1$	click here
5	NETGEN 1D-2D-3D	$1$	$5$	$2$	$8$	click here

It is possible to see that all meshes differ only for size and order. This way it will be possible to study the effect of both mesh size and order on the accuracy. In addition, the fine second order mesh was also used to solve the problem by using parallelisation (see the Dealing with Convergence Issues episode for details on how to setup a parallel study with Elmer). We can use the parallel solution to investigate whether the parallelisation itself introduces additional errors.

Results

Comparison Field Magnitude

The figure above shows a comparison of the various error fields for the different values of mesh size and order. Going left to right we refine the mesh, while going top to bottom we increase the order. The figures show isosurfaces of the magnitude of the comparison field expressed in $\text{dB}$. $10$ isosurfaces are drawn for each comparison field, at $10$ equally spaced values between the minimum and maximum of the comparison field magnitude. The closer the comparison field is to $0$ $\text{dB}$ the higher the agreement of the magnitude of the FEM field is to the magnitude of the exact field. The isosurfaces are semi-transparent so to be able to see the various layers of them. The colour scheme used for them is divergent, centred on $0$ $\text{dB}$ (white). The positive values are in shades of red, while the negative values are in shades of blue. The outer surface of the domain, a sphere of $0.1$ $\text{m}$ radius, is shown as a wireframe for reference (the nodes of the frame are that of the mesh). To create similar plots with ParaView, you can refer to the Intro to Paraview episode.

It is possible to see that the agreement between the FEM solution magnitude and the exact one increases the higher the order and the finer the mesh, as expected. For low order and high mesh size, the error can be found all over the domain while it gets progressively confined close to the field source the higher the order and the smaller the mesh size. It is interesting to note that the comparison field tends to assume values closer to $0$ $\text{dB}$ more significantly when the order is raised to $2$. Also interesting is that the FEM solution tends to underestimate the magnitude of the field, given that most error values are negative $\text{dB}$.

The figure below shows the magnitude of the comparison field $C$ in $\text{dB}$ as a function of distance from the centre of the pulsating sphere. The data is obtained from all the nodes. Note that all the plots use the same axes. These plots largely tell the same story as the 3D plots above, but allow to observe few additional things.

Click Here for Larger Version

First of all, linear meshes solutions appear to be affected by a “level offset”, that is, the FEM solution is a fraction of a $\text{dB}$ lower then the exact one. This offset seems to to be reduced the finer the mesh, although it is still fairly significant for the fine first order mesh (more than $-0.2$ $\text{dB}$). On the other hand, the second order meshes values of the comparison field appear centred on $0$ $\text{dB}$. Higher errors are seen closer to the source.

But how good are these values? To understand this, we can refer to acoustic measurements. Let’s suppose that we get hold of a real sphere source (which doesn’t quite exist, but let’s suppose it does for the sake of argument). Then, we could try to setup a very good measurement system to figure out the pressure field that it is emitting. Even by attempting to control the experimental conditions as much as possible (temperature, humidity, position of microphones and probes, etc…) it will be extremely hard to get to an accuracy and precision of $1$ $\text{dB}$. Typically, $0.1$ $\text{dB}$ is about how accurate (and precise) state of the art measurements can get, and that will require a very good effort. We should treat then $0.1$ $\text{dB}$ as a sort of “threshold” for the performance of our FEM solution.

Referring to the plots above, we can see that the first order meshes have an overall shift (mean value) which is well above $0.1$ $\text{dB}$. Similarly, the spread of the values is above $0.1$ $\text{dB}$ as well. The opposite happens to second order meshes, where the overall error (and the spread around it) is well below $0.1$ $\text{dB}$. We can therefore conclude that second order meshes significantly outperform first order meshes in terms of accuracy.

Still, we should be careful not to take the concept of accuracy rigidly. In fact, the level of accuracy required from a simulation depends on the application. Note how the errors from the first order meshes are mostly within $0.6$ $\text{dB}$. In many applications this can be an acceptable error (think, for example, of the mass production variation of loudspeakers which can easily be $\pm4$ $\text{dB}$). As always, the best thing to do is to perform an accuracy study on a simplified model first, and then apply the settings that produce the accuracy we aim for to the more complex problem.

Additionally: punctilious quantification and analysis of error is more important than its blind minimisation. For example, our plots above suggest that errors are higher the closer we are to the field sources. A way to increase accuracy of the solution could then be local refinement of the mesh only close to the sources, which in turn will propagate lesser errors in the far field. By keeping this open mindset we can devise strategies to reach higher accuracy without resorting to second order meshes and the associated computational costs and harder to manage convergence, in the case those concerns matter.

Comparison Field Phase

The figure above shows the phase of the comparison field in a similar fashion to Figure 1. In this case we can still see that accuracy increases left to right - top to bottom, just as for the magnitude. However, in this case the error is more uniformly distributed in the domain, with a less pronounced “retreat” towards the source with respect the magnitude case.

The comparison field phase as a function of distance from the source centre is shown below.

Click Here for Larger Version

This figure does a better job at showing that the errors affecting the first order meshes are higher closer to the source. We can draw parallel conclusions as those that we draw from the magnitude analysis, observing how significantly the error is reduced by mesh order. In this case, though, only the coarsest first order mesh is affected by truly significant errors. If radiants are not intuitive for you, you can keep in mind that $0.02$ $\text{rad}$ are approximatively $1.15\text{°}$. We can see that even the most inaccurate solution rarely exceeds an error of $1\text{°}$, with the exception of the volume up to $1$ $\text{cm}$ away from the source centre. In a real measurement accuracy and, especially, precision up to $1\text{°}$ can be very challenging to get (for measurements of sources in free field) so we can say that all of our simulations are as accurate as a state of the art measurement of a real system even for the largest errors we observe, and even more so for the more refined meshes. As far as phase error is concerned, a denser first order meshes already provides essentially negligible error.

Effect of Parallelisation

The figure above shows a comparison of the solution for the $5$ $\text{mm}$ second order mesh (left) results with the same study solved with parallelisation (right). When parallelising execution the mesh is sliced in as many sections as there are processes. The boundaries of these sections can be seen in the right hand plots, shown as wireframes. The plot shows the magnitude (top) and phase (bottom) of the comparison fields for the two studies respectively. It is possible to see that the two solutions are essentially exactly the same, with the exception that the parallel solution has a phase shift of $\frac{\pi}{2}$ with respect the exact solution. This is more evident in the figure below, where $\frac{\pi}{2}$ is subtracted to the phase of the comparison field.

This shows that parallelisation did not introduce any error beside an overall $\frac{\pi}{2}$ phase shift. This should not be considered problematic as overall phase shifts are easily accounted for. Spatial patterns in the comparison field are not observed, suggesting that the partitioning of the domain did not results in additional errors.

Conclusion

In most of the previous episodes we used second order meshes when solving Helmholtz problems. Already in The Pulsating Sphere episode we seen how, when using second order meshes, the results are very accurate as soon as the mesh maximum size is below one tenth of a wavelength. However, in the Dealing with Convergence Issues we resorted to first order meshes in order to ease convergence. We argued that accuracy to to the underlying PDE real solution is strongly controlled by mesh size, so it should not be a problem to “downgrade” the mesh, so to speak. In this episode we took a quantitative look at the effects of mesh size and order, so to understand more properly what happens when we reduce mesh order. To do so, we compared the results from different FEM studies, with different meshes, to the exact solution of the PDE being studied.

In terms of magnitude of the FEM field we found that use of first order meshes produces higher errors. If high magnitude accuracy is needed, sticking to second order meshes could be necessary. Alternatively, first order meshes could be prepared with smaller size than one tenth of a wavelength, eventually only close to the field sources where the error is largest. It should be noted that, however, our results for a tenth-wavelength first order mesh agreed with the exact solution well within $1$ $\text{dB}$. This is normally as accurate and precise as ordinary free field acoustic measurements get. In other words we can say that, if we could measure a real sphere source in free field, the measurement results might very well be affected by errors as large, or larger, than those observed in the worst performing mesh.

In terms of phase similar considerations hold, but errors are in general much less significant, only seldom exceeding $1\text{°}$ in the far field.

Finally, the effect of parallelisation on the accuracy of a FEM solution was investigated for the finest second order mesh. It was found that parallelisation only introduced a $\frac{\pi}{2}$ phase shift. It is reasonable to expect that the presence of this phase shift, and its eventual value, might differ between solver versions.

Putting all of these observation together, we can say:

1. Second order meshes with maximum size below one tenth of a wavelength provide very accurate results and should be preferred.
2. First order meshes with maximum size below one tenth of a wavelength result in higher errors with respect their second order counterparts. Magnitude errors of up to $1$ $\text{dB}$ should be expected, while phase errors up to $2\text{°}$ should be expected. Errors are normally more significant the closer we are to the field sources, suggesting that a way to increase accuracy of first order meshes could be local mesh refinement close to the field sources.
3. Parallelisation can introduce an overall phase shift to our FEM solution.

With reference to the above, we can finally conclude that we should expect the final results from the Dealing with Convergence Issues to be accurate up to $1$ $\text{dB}$ and $2\text{°}$, with potentially a phase shift coming from the parallel solver procedure. To be noted, however, that our benchmark study (pulsating sphere in free field) differs significantly from a realistic room, with is the context of the previous 2 episodes. Whilst it is still a valuable system to understand the effect of the various parameters controlling the accuracy of FEM solutions, we would be able to trust our analysis more if we had a benchmark model closer to the real deal. As always in physics, problems should be dealt with in increasing levels of complexity. In the next episodes we will look into designing a better benchmark model for the room problem.

License Information

This work is licensed under a Creative Commons Attribution 4.0 International License.

• Physics
• Elmer
• FEM
• Linear System
• Single Physics