It has been quite some time since I contributed a new tutorial to this website. In addition to the general lack of time, the main reason is that I was busy studying FEniCSx. My intention was to create a piece of software to solve the Helmholtz equation in heterogeneous media, with support for adiabatic absorbers boundary conditions. Whilst I was able to draft something up (see the acoupy_helmholtz fenicsx branch) I did not make it very far into workable code yet. This also because most of my time has been spent dealing with the Helmholtz equation itself, and the modifications required to support such simulations. Eventually all of this work will find its way into various articles in this website. But for the time being I decided to stop for a second. There are many more Elmer solvers that I have not explored yet. I decided that, while I carry on this “research” work, I can explore these solvers and write a few new episodes while doing so. This article is about a first model of the ear canal. The model is not particularly sophisticated, as it uses Elmer’s HelmholtzSolve solver that we have learned and love. However, it deals with a realistic geometry. This allows me to show how to manipulate complicated geometries, produce initial simple results, assess whether the results are reasonable, and lay a foundation for more complicated simulations. So, let’s dive in!

In the Elmer vs FEniCS - Part 2 episode we studied the effect of mesh size on accuracy. We largely found the same results observed in the Mesh Order and Accuracy episode. We also confirmed that the agreement between FEniCS and Elmer is very good. However, we noticed that FEniCS was significantly slower than Elmer. This because we used a MUMPS direct method on FEniCS and a BiCGStabl iterative method on Elmer. In this episode we will explore the iterative solver options for FEniCS and attempt to optimize execution time for our simulation.

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.

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.

In the last episode we explained why we might want to make our own packages. We selected Python as a language for this series as we will be writing software to facilitated FEniCS simulations. In this episode we will take a look at the outline of the acoupy_meshutil and acoupy_helmholtz packages. We will not go into the details of creating Python packages. For more information about that you can see Hosting Python Packages on Git Repositories. Instead, we will explain the basic outline of the packages.

In the last three episodes about FEniCS we introduced the software and installed it. We also derived the weak form for the Helmholtz equation and developed some code to solve it. This is all nice and good, but writing scripts from scratch every time can be tedious and error prone. We can help ourselves by writing reusable code packages to serve as tools. For example, we could write a software package to solve the Helmholtz equation with FEniCS. Or any other package we might find useful. The idea is that we can create an entire toolbox to deal with common tasks when creating simulations. We started the development of some tools under the acoupy project. This series of articles will serve as a sort of “development” diary, illustrating how the packages are being developed and why. This will allow us to illustrate how to develop scientific software.

In the Intro to FEniCS - Part 2 episode we completed the first step in developing a FEniCS model. We characterised the domain. Then, we selected a governing PDE. We then provided suitable boundary conditions. Finally, we expressed the PDE in weak form. We chosen a generalised Helmholtz equation as our PDE. In this episode we will start using the resulting weak form to develop a model.

In the Intro to FEniCS - Part 1 episode we introduced FEniCS. We compared it to Elmer and covered how to install it on an Ubuntu workstation. This episode will cover the process of deriving the weak form for the Helmholtz equation. We will focus on the mathematics, which are advanced. Still, we will attempt to present the formalism in an approachable manner.

In the previous episodes we used Elmer to develop a few different acoustics models. We were able to verify the accuracy of Elmer with few benchmark problems. Additionally, we studied the effect of mesh order and size on the accuracy. Finally, we started probing vibro-acoustic problems. We only scratched the surface of Elmer capabilities in this field. Still, it is worth to take a step back and reconsider these problems under a different light. The FEniCS project is perfect for this. We will be introducing FEniCS in this episode.

In many of the previous episodes we produced animations of steady state Elmer FEM solutions. For simplicity we did not discuss the details of how to prepare them. The process is actually very simple, but there are few tips and tricks to keep in mind. This episode aims to provide a concise introduction to animations with ParaView.

In the Tuning Fork - Part 3 episode we seen how to setup a weak vibroacoustic coupling problem. We applied an harmonic body force to our tuning fork, and used its displacement field to wake the acoustic field within a coupled volume of air. In this episode we will review the results.

In the Tuning Fork - Part 1 and Tuning Fork - Part 2 we set up and solved an Elmer study for the modes of vibration of a tuning fork. We learned that the principal mode of the fork, the one related to the nominal frequency of the fork, is very marginally affected by boundary conditions affecting the handle of the fork. In this episode we will attempt studying the acoustic field emitted by the fork in the surrounding air. In other words, we will start exploring the coupling of vibration and acoustics.

In the Tuning Fork - Part 1 episode we setup a new study for the determination of the elastic modes of vibration of a tuning fork. Our study made use a few meshing and solver setting ideas we learned along the way, including the use of $p$-elements. In this episode we will examine the results from the Elmer simulation.

The very first episode in which we introduced Elmer was the Elastic Modes of a Metal Bar episode. In that episode we introduced Elmer by solving a linear elasticity eigenproblem, and mentioned how vibration is an integral part of acoustics, being vibrating bodies one of the principal causes of airborne sound radiation. In this new series of episode we will explore vibration further, integrating in it what we learned so far, and we will explore vibro-acoustic coupling with Elmer.

This project has now quite a few articles and simulations going. Unfortunately, the way the Julia code repositories were set up is not very functional and is not scaling well to many repositories. For this reason, this site will undergo a period of refactoring. The Status Quo At the moment there are few additional Julia code repositories that contain helper code, mainly for postprocessing and to create the plots in this website.

In the Home Studio - Part 3 and Home Studio - Part 4 episodes we referred to the results of our frequency sweep as frequency responses. In those episodes we calculated the steady state pressure disturbance in a room at different frequencies and sampled the field at a couple of locations so to be able to plot the steady state magnitude at those locations as a function of frequency. It is intuitive to refer to this result as a frequency response, in the sense that the field magnitude as a function of frequency gives us an understanding of how strong the steady state disturbance is at any given frequency in our study.

In the Home Studio - Part 3 episode we simulated the steady state pressure field for a realistically shaped room governed by the linear wave equation. The boundary condition for the room was that of ideally rigid walls. We seen that the effect of such a boundary condition is an extremely uneven frequency response. This due to how sharp the resonances of the room are. In this episode we will introduce slightly more realistic boundary conditions and see what the effect on the solution is.

In the Home Studio - Part 1 episode we gave a first look at the response of a realistically shaped room with rigid walls. We understood how to develop and run a model with an uniform velocity sphere source placed somewhere in the room. We run the model up to $400$ $\text{Hz}$ with the help of the convergence considerations outlined in the Dealing with Convergence Issues episode. In the Home Studio - Part 2 episode instead we used Elmer’s WaveSolver to compute the eigenmodes of the room.

In the previous episodes we manipulated some of the result fields from Elmer with meshio, a very nice Python package that allows us to do many useful operations with meshes and fields. However, in the episodes we did not focus into the details of how to use the package. Whilst this webiste is perhaps not the best place to do so, as the best place to understand how to use a package is, of course, its official documentation, it is useful to collect here few examples and remarks about using it with Julia, the programming language we are using for this project.

In the Home Studio - Part 1 episode we computed the steady state field of a realistic, but still rigid walled, room. We put a source in the room, modelled as a flux boundary condition, and run the study at few different frequencies. The solutions that we found were reminiscent of modal patterns, which is expected as the low frequency response of a room is dominated by its resonances. However, that kind of study does not inform us on the actual resonance frequencies of the room, which are properties of great interest, as well as how the associated modal shapes (eigenmodes, or eigenfunctions) look like.

In the previous episodes we obtained various numerical FEM solutions for Helmholtz problems. Although we discussed their meaning in the previous episodes it is useful to discuss that again more deeply in a single page, which can act as a quicker reference. This episode will focus on the Helmholtz solver solution fields, while other episodes will focus on the output of other solvers. The Physical Meaning of a Solution Field The meaning of a solution field largely depends on the solver itself.

In the Rigid Walled Room Revisited - Part 2 episode we reviewed the solution provided by the Elmer Wave Equation solver. In that study we solved for the eigenfrequencies and eigenmodes of a rectangular room by using a first order mesh. The results we found were very accurate already, with eigenfrequencies within $1$ $\text{Hz}$ from the exact value and eigenmodes accurate within $1.4\text{%}$. However, we seen in the Mesh Order and Accuracy episode that we can significantly increase the accuracy of results by acting on the mesh fineness and, more effectively, the order.

In the Rigid Walled Room Revisited - Part 1 episode we seen how to setup an eigen system problem with Elmer by making use of the Wave Equation solver. In this episode we will review the simulation results and check the agreement with the analytical solutions in the Acoustic Modes of a Rectangular Room episode. Project Files All the files used for this project are available at the repositories below:

We covered the rigid walled rectangular room previously in the Rigid Walled Room episode. In that episode we solved for the steady state field in a rectangular rigid walled room as driven by a source placed somewhere in the room. This allowed us to see how the steady state field is sustained by a modal superposition, the purer the lower the driving frequency (assuming that this driving frequency matched an eigenfrequency of the room).

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 the Home Studio - Part 1 episode we faced convergence issues when dealing with the highest driving frequency for our room ($400$ $\text{Hz}$). This meant that we could not quite trust the solution, and so we discarded it. In this episode we will look at what to do in this cases, and how to reach convergence. Rather that dealing with the issues in abstract and general terms (which would require writing an entire book about it) we will use the Home Studio - Part 1 episode to introduce the problem and figure out how to deal with it practically.

In the Rigid Walled Room episode we seen how to model a rectangular room with rigid walls. We driven the room at the modal frequencies and compared the solution field with the theoretical modal shapes, finding that the results matched single modal shapes real well until, at a frequency high enough, the contribution of multiple modes (in addition to that related to the driving modal frequency) became important. In this episode we will look at making the model more realistic.

In the previous episodes we often made use of the ParaView postprocessor to visualise our solution field from the Elmer solver. ParaView can do all sorts of cool visualisations and animations, as well as providing the way of doing quantitative analysis. It is by far the best option to visualise and postprocess results from Elmer. It also allows to export data in various formats, such as CSV, that allow us to do any additional kind of postprocessing or verification, by either using Julia, Python or any other language, or even spreadsheet software if you fancy that (for whatever reason…).

In the Acoustic Modes of a Rectangular Room episode we explored the analytical model of a rigid walled room with some Julia code. We focused on finding the resonance frequencies (or eigenfrequencies) of the room and calculating the related modal patterns (eigenfunctions). Now that, thanks to The Pulsating Sphere episode, we know how to setup Helmholtz problems with Elmer we can approach the problem with the FEM method. In this episode we will solve for the modal superposition in a rectangular rigid walled room and use the results from the Acoustic Modes of a Rectangular Room episode to check the accuracy.

In the previous episodes we solved a few equations with Elmer. We did some choices when we setup the solver parameters. What those parameters do, and how should we set them? This is perhaps the trickiest part in FEM (beside making the mesh right). In this episode we will step back and look at those solver options more closely. This post is really not meant to be an exhaustive explanation. For that, refer to the Elmer documentation.

In this episode we will build a model of a pulsating sphere source. The pulsating sphere source is an ideal source which forms the base for the development of point sources. In essence, a point source is a pulsating sphere in the limit of $a$, the radius of the sphere, approaching $0$. For this reason, although abstract, the pulsating sphere is a very powerful theoretical tool that enables the study of point sources which in turn, through integration and wave propagation principles, enable to study of any arbitrary acoustic field source.

In the last episode we seen that our model solution for the modes of vibration of a metal bar wasn’t looking particularly good for the highest mode. In fact, when we clipped into the bar, we seen few bubbles of discontinuity in the displacement field that are not expected for linearly elastic homogeneous and isotropic bodies. To understand this we have to step back a little and think about FEM some more.

In the last episode we examined the analytical solution of the acoustic modes of a rectangular room and we are now ready to take steps into moving in the world of FEM modelling. Elmer is a powerful package, but it is not extremely user friendly. So, it is best to have a gentle introduction to it first before dwelling into the intricacies of FEM modelling of acoustic fields. One of the simplest problems to solve with Elmer is that of the elastic vibration modes of solids.

In this episode we will look at how to make a simple Julia model of one of the simplest systems in acoustics, a rectangular room with rigid walls, assuming adiabatic wave propagation. Even if this system is among the simplest in acoustics it is actually already very complicated. As such, we will focus only on the modes, one part of the problem, without attempting impulse response simulation or other fancy things like that, for now.

In the previous article an overview of the open-source ecosystem for scientific and technical computing was presented. In this article the focus will be on FEM. We will first look at how we build a FEM problem, and then select among the packages we listed the software that is best suited for implementing FEM studies. The Basic Pipeline The diagram above shows the process of creating a FEM study in its simplest form.

In the previous episodes we seen that physical systems are described by PDEs. However, PDEs are seldom solved analytically, unless for certain simple cases, the reason being that finding the solution of a PDE on a complex, realistic, domain with realistic boundary and initial conditions is extremely complicated, or downright impossible, even when the PDE is known to have analytical solutions from theoretical considerations. For this reason, numerical methods have arisen in order to handle these problems.

Welcome to the first actual episode of the series about acoustic modelling with open-source software. We will first try to understand what modelling acoustics means. In reality it doesn’t mean just one thing, as many phenomena of acoustic wave production and propagation can be modelled and simulated in various different ways, with higher or lower degree of accuracy. However, the core of the modelling problem resides in partial differential equations. This post will be a very, very, brief, intuitive and nonrigorous introduction to the topic, mainly to give context to those that are not accustomed to the concept.

In the first episode of this series the open-source ecosystem for the modelling of acoustics was covered. As we will see in the next episodes, we will make extensive use of FEM, at least initially, to model acoustic phenomena. This episode focus on building a simple simulation workstation based around Elmer, a powerful multiphysics FEM suite with excellent acoustics capabilities. Why an Ubuntu Based Workstation? Ubuntu Linux is a very convenient operating system for this purpose as all the needed packages are straightforward to install, especially Elmer and Salome, which can be quite complicated to get up and running in other distributions.

One of the backbones of the scientific paradigms is the repeatability of results. Ideally, all software used for scientific research should be open-source, so to allow complete repeatability of data analysis and models. And in fact, a lot of scientific software is indeed open-source, and a lot of this software can be used to model acoustic phenomena. What follows below is a panoramic overview of the scientific open-source software ecosystem. This overview is far from comprehensive.