Intro to FEniCS - Part 2

December 11, 2021

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.

Prerequisites

Dealing with PDEs requires fairly advanced mathematical knowledge. The material presented below will make the most sense if you have some knowledge of:

Topology;
Complex algebra;
Mathematical analysis;
Functional analysis.

Still, we will discuss each passage step by step. You are welcome to attempt reading this article even though you are not an expert in the topics listed above. It is very likely there might be mistakes in this article. Please reach out if you find them!

We recommend to read the first two chapters of Solving PDEs in Python before diving into this. If this material doesn’t make too much sense, do not worry about it. The important takeaway is that any PDE can be turned into a weak form. FEniCS will take the weak form to solve the PDE.

The Approach

Our aim is to develop a weak form for the Helmholtz equation. This equation is crucial in acoustics, as it governs steady state wave propagation at a given frequency. In the previous episodes we stressed the importance of increasing complexity gradually. Adding complexity a step at a time is a very good way to handle a problem in physics. It allows to build confidence in the results by ensuring the details do not contradict the previous coarser estimates. Yet, when building a mathematical formalism, another approach might be better.

In mathematics many problems turn out to be a particular form of a more general problem. Going from general to particular is easier. This is the case for the material presented in this episode. We will have a more general version of the Helmholtz equation. We will complement it with multiple boundary conditions. This PDE version is more complex. But we will simply need to zero some coefficients to recover the simpler version. This is why a more general PDE will work well for us. It is like we have written a lot of Helmholtz PDEs at once. This means we will not have to derive many weak forms over and over.

We will still use a simple to complex approach when writing a solver for the PDE. We will zero out most terms to begin with, and add them back little by little. This will allow us to validate “features” in baby steps. In turn, this enables us to write a robust solver.

In other words, we can use both approaches. We just need to make sure we use them where they shine.

Notation

In the following we will use the following notation.

Lower case letters to denote scalar and scalar functions, for example $f$;
Bold lower case letters to denote Euclidean vectors end vector functions, for example $\mathbf{v}$;
A dot to denote the inner product in an Euclidean space, for example $\mathbf{v} \cdot \mathbf{u}$;
Angled brackets to denote the inner product in a function space, for example $\langle f,g \rangle$;
We use $\nabla^2$ for the Laplace operator;

The Governing PDE

Let’s start from visualising our problem. We have a domain $\Omega \subset \mathbb{R}^3$. This set represents a body of fluid, so it is not particularly funky. We can require it to be an open set without sacrificing any valuable physical system. An open set in $\mathbb{R}^3$ can be defined as follows (based upon Analisi Matematica di Base, even though Gilardi’s definition is more general):

$$ \forall \mathbf{x} \in \Omega \space \exists \space r \in \left(0, +\infty\right) : \left\{ \mathbf{y} \in \mathbb{R}^{3} : \left\| \mathbf{y} - \mathbf{x} \right\| <\ r \right\} \subseteq \Omega $$

In essence, for every point in $\Omega$ we can find a ball (without surface), maybe very small, completely enclosed in $\Omega$. It is intuitive to see how any realistic body of fluid satisfies this requirement. We provided a definition of open set to clarify what kinds of bodies we are dealing with. However, below we will not provide the definition of other topological concepts such as boundary. If you need a proper definition for those, you can consult any mathematics textbook. If you know Italian, Analisi Matematica di Base is a good one. Alternatively, Wolfram MathWorld can provide quick facts to you. Wikipedia can also be useful, and it appears to be correct on mathematical topics. In fact, the definition of open set it provides is the same as that above.

We denote the boundary of $\Omega$ with $\partial \Omega$. In this domain (fluid) we want to determine a steady state acoustic field at some frequency. To do so we need a governing equation.

Let’s start from the wave equation as formulated in the Elmer Models Manual. The equation is reported below with minor notation changes:

$$ \nabla^2 p - \frac{1}{c^2} \frac{\partial^2 p}{\partial t^2} + \frac{\eta}{c^2} \nabla^2 \left( \frac{\partial p}{\partial t} \right) - \frac{\alpha}{c^2} \frac{\partial p}{\partial t} = f $$

The terms of this equation are detailed below.

$p$ is the pressure field. It is a complex function of space and time. By denoting the closure of $\Omega$ with $\overline{\Omega}$ we have:

$$ \overline{\Omega} \doteq \Omega \cup \partial \Omega $$ $$ p : \overline{\Omega} \times \mathbb{R} \rightarrow \mathbb{C} $$ $$ p = p \left( \mathbf{x}, t \right) $$

Above we denoted with $\mathbf{x}$ the generic point in $\overline{\Omega}$. $\mathbf{x}$ is a 3D vector:

$$ \mathbf{x} = \begin{bmatrix} x_{0} \\ x_{1} \\ x_{2} \end{bmatrix} $$

$t$ is the generic time instant in $\mathbb{R}$. $p$ is complex to capture, in one single quantity, all information about amplitude and phase of the field. The physical solution is found by simply taking the real part. Of course, for the PDE above to make sense, we need $p$ to be differentiable with respect time at least one time. Also, we need it to be differentiable with respect space at least two times. Additionally, we should expect pressure fields to be “well behaved”. That is, we do not expect them to diverge, as that would be nonphysical. We do not expect the derivatives of $p$ to diverge either. We can say that, for each time instant $t$, we expect the integral of $\left| p \left( \mathbf{x}, t \right) \right|^2$ over $\Omega$ to be finite, and the same for its derivatives. These constrains create a well defined class of functions $p$ is part of. Or, $p$ lives in a function space. We will deal with this after we switch over to the Helmholtz problem.

$f$ is a source term. This is actually a more complex object. It will be easier to deal with it once we transition to the Helmholtz equation. For now, let’s just think of it as an object dependent upon $\mathbf{x} \in \overline{\Omega}$ and $t \in \mathbb{R}$.

The coefficients are $c$, $\eta$ and $\alpha$. These are real numbers and describe the fluid. $c$ is the speed of sound in the material. $\eta$ and $\alpha$ introduce damping.

We now impose the unknown pressure field and the source term to be steady state harmonic. This means that they are factorised in a spatial part and temporal part. The temporal part is an harmonic function. Hence:

$$ p \left( \mathbf{x}, t \right) = u\left( \mathbf{x} \right) \exp \left( j \omega t \right) $$

$$ f \left( \mathbf{x}, t \right) = q\left( \mathbf{x} \right) \exp \left( j \omega t \right) $$

Where $\omega$ is the angular frequency of the harmonic source, end hence that of the field. $j$ is the imaginary unit. The following holds:

$$ \omega = 2 \pi \nu = k c$$

With $\nu$ the frequency and $k$ the wavenumber respectively.

Based on what we said above, $u$ is a complex function defined in $\overline{\Omega}$. $u$ is differentiable two times with respect each spatial direction. $\left|u\right|^2$ has a finite integral in $\Omega$ and so do all its first and second partial derivatives. The function space of functions having these properties is a Sobolev space. Let $\mathcal{L}_{\Omega}^2$ denote the space of complex functions in $\Omega$ whose square module has finite integral. Then, we can call our functional space $V$ and define it as follows:

$$ V \doteq \left\{ \phi \in \mathcal{L}_{\Omega}^2 \space : \space \forall \xi, \gamma \in \left\{ 0, 1 , 2\right\} \space \frac{\partial \phi}{\partial x_{\xi} \partial\ x_{\gamma}} \in \mathcal{L}_{\Omega}^2 \right\} $$

$V$, the home of $u$, is equipped with the following inner product (courtesy of $\mathcal{L}_{\Omega}^2$, but other inner products are also possible):

$$ \langle\phi_1, \phi_2 \rangle = \int_{\Omega} \phi_1 \phi_2^{\star} \space d\mathbf{x} $$

For any pair of functions $\phi_1$ and $\phi_2$ in $V$. Note that $\langle\phi_1, \phi_2 \rangle$ is complex in general. $^\star$ denotes complex conjugation.

By simple differentiation we have:

$$ \frac{\partial p}{\partial t} = j \omega \space u\left( \mathbf{x} \right) \exp \left( j \omega t \right) = j \omega p $$

$$ \Rightarrow \frac{\partial^2 p}{\partial t ^2} = \frac{\partial}{\partial t} \frac{\partial p}{\partial t} = j \omega \frac{\partial p}{\partial t} = - \omega^2 p$$

As we mentioned above, the source term is slightly different. If we had a simple point source in $\mathbf{x}_s \in \Omega$ we would have:

$$ q = \delta_{\mathbf{x}_s}$$

Where $\delta_{\mathbf{x}_s}$ is the three-dimensional Delta function centred in $\mathbf{x}_s$. Despite its name, this is not a function, but a linear functional. This objects takes functions as inputs and produces complex numbers. That is, if we take any function $\phi \in V$ we have, by definition:

$$ \delta_{\mathbf{x}_s} \left(\phi\right) = \int_\Omega \delta_{\mathbf{x}_s} \phi \left( \mathbf{x} \right) d\mathbf{x} = \phi \left(\mathbf{x}_s \right) $$

So, if we want to be able to describe all sources, including the point-like sources, $q$ cannot be an ordinary function. $q$ must be a linear functional over $V$. Or, in other words, an element of the dual space of $V$, denoted with $V^\prime$. Note that we can have an ordinary function defining a functional through the inner product. For example, we can have $\phi \in V$ inducing $q$ as follows:

$$ q\left(\varphi\right) = \langle \phi, \varphi \rangle \space \forall \varphi \in V $$

We will use these facts to derive the weak form. Before proceeding, we need to substitute our harmonic quantities in the wave equation. This leads to the Helmholtz equation:

$$ \begin{aligned} \nabla^2 \left[ u\left( \mathbf{x} \right) \exp \left( j \omega t \right) \right] - \frac{1}{c^2} \left[ - \omega^2 u\left( \mathbf{x} \right) \exp \left( j \omega t \right) \right] + \\ + \frac{\eta}{c^2} \nabla^2 \left[ j \omega \space u\left( \mathbf{x} \right) \exp \left( j \omega t \right) \right] - \frac{\alpha}{c^2} \left[ j \omega \space u\left( \mathbf{x} \right) \exp \left( j \omega t \right) \right] = \\ q\left( \mathbf{x} \right) \exp \left( j \omega t \right) \end{aligned} $$

The Laplace operator $ \nabla^2 $ is a differential operator acting only with respect the spatial variables:

$$ \nabla^2 \doteq \frac{\partial^2}{\partial x_{0}^2} + \frac{\partial^2}{\partial x_{1}^2} + \frac{\partial^2}{\partial x_{2}^2} $$

Hence $ \exp \left( j \omega t \right) $ will come out its argument. The harmonic function then appears on all terms, and can be removed. To remove some “fumo negli occhi” we will simplify the notation by removing the explicit dependency on $\mathbf{x}$ from the notation. By arranging the terms and using the equations linking $\omega$ and $k$ we arrive to:

$$ \gamma \nabla^2 u + \beta u = q $$

Where $\gamma$ and $\beta$ are material parameters defined as follows:

$$ \gamma \doteq 1 + j k \frac{\eta}{c} $$ $$ \beta \doteq k^2 - jk \frac{\alpha}{c} $$

By comparison with the Elmer Helmholtz Model we can see that Elmer is solving a version of this equation with $q = 0$, $\eta = 0$ and $D = \frac{\alpha}{c}$

Boundary Conditions

We formulate three types of boundary conditions, Dirichlet, Neumann and Robin. Each will be applied to a set of “slices” of the boundary $\partial \Omega$. We consider these slices as a partition of $\partial\Omega$. Hence, all slices are disjointed (empty mutual intersection) and their union provides $\partial\Omega$.

To define the boundary conditions we will make use of the normal vector to the boundary. For any point $\mathbf{x} \in \partial \Omega$ we denote with $\hat{\mathbf{n}}$ the unit vector normal to the boundary. We consider this normal unit vector an outward unit vector, meaning that $\hat{\mathbf{n}}$ points into the complement of $\Omega$, which is given by $\mathbb{R}^{3} \setminus \Omega$.

Dirichlet Boundary Conditions

This condition simply prescribes the value of $u$ on selected parts of the boundary. We want to apply a condition of this kind to $N_{D}$ slices of $\partial \Omega$. Hence, we impose:

$$ u \left( \mathbf{x} \right) = u_{D}^i \left( \mathbf{x} \right) \space \text{for} \space \mathbf{x} \in \partial \Omega_{D}^i \\ \space i = 0, 1, \dots, N_{D} - 1 $$

Where the complex functions $u_{D}^i$ defined in each of the slices $\Omega_{D}^i$ are given.

Neumann Boundary Conditions

When we think about acoustics we think about the pressure field. But there is also another field, the particle velocity field. We define it as follows:

$$ \mathbf{v} : \overline{\Omega} \times \mathbb{R} \rightarrow \mathbb{C}^3 $$ $$ \mathbf{v} = \mathbf{v} \left( \mathbf{x}, t \right) $$

Our problem is steady state harmonic, so we write:

$$ \mathbf{v} \left( \mathbf{x}, t \right) = \mathbf{w}\left( \mathbf{x} \right) \exp \left( j \omega t \right) $$ $$ \mathbf{w} : \overline{\Omega} \rightarrow \mathbb{C}^3 $$

With $\mathbf{w}$ the spatial part of the particle velocity vector field.

Thanks to the Euler equation in frequency domain we can relate the normal component of $\mathbf{w}$ at the boundary with the the normal partial derivative of $u$. In fact, it is even possible to think of the wave equation as two coupled equations, one for pressure and one for particle velocity (or flow velocity). This latter approach is followed by the DeltaEC Software. We have:

$$ \mathbf{w} \cdot \hat{\mathbf{n}} \doteq w = -\frac{1}{j\omega \rho} \nabla u \cdot \hat{\mathbf{n}} $$ $$ \Rightarrow \nabla u \cdot \hat{\mathbf{n}} = -j \omega \rho w $$

Where $\rho$ is the at rest density of the fluid.

Hence, this means that we can prescribe the normal particle velocity at a set of boundary slices. If we do so, we are applying a Neumann boundary condition. Hence, to apply a condition of this kind to $N_{N}$ slices of $\partial \Omega$ we impose:

$$ \nabla u \cdot \hat{\mathbf{n}} = -j \omega \rho w_{l} \space \text{for} \space \mathbf{x} \in \partial \Omega_{N}^l \\ \space l = 0, 1, \dots, N_{N} - 1 $$

Note that, for each slice $\Omega_{N}^l$, $w_{l}$ is a complex function defined over the whole slice. This ties in with the Elmer Helmholtz Model flux boundary condition. We can follow the same approach as the Elmer manual to put the equations above in terms of displacement. If we denote with $\mathbf{d}$ the displacement:

$$ \mathbf{d} : \overline{\Omega} \times \mathbb{R} \rightarrow \mathbb{C}^3 $$ $$ \mathbf{d} = \mathbf{d} \left( \mathbf{x}, t \right) $$

And we remember we are working with steady state harmonic quantities:

$$ \mathbf{d} \left( \mathbf{x}, t \right) = \mathbf{s}\left( \mathbf{x} \right) \exp \left( j \omega t \right) $$ $$ \mathbf{s} : \overline{\Omega} \rightarrow \mathbb{C}^3 $$

Since the velocity is the time derivative of displacement we have:

$$ \begin{aligned}
\mathbf{v} \left( \mathbf{x}, t \right) = \mathbf{w}\left( \mathbf{x} \right) \exp \left( j \omega t \right) = \\
= \frac{\partial}{\partial t} \mathbf{d} \left( \mathbf{x}, t \right) = \\
= \frac{\partial}{\partial t} \mathbf{s}\left( \mathbf{x} \right) \exp \left( j \omega t \right) = \\
= j \omega \space \mathbf{s}\left( \mathbf{x} \right) \exp \left( j \omega t \right) \end{aligned} $$

So:

$$ \Rightarrow \mathbf{w}\left( \mathbf{x} \right) = j \omega \space \mathbf{s}\left( \mathbf{x} \right) $$

$$ \Rightarrow \nabla u \cdot \hat{\mathbf{n}} = \omega^2 \rho s, \space \mathbf{s} \cdot \hat{\mathbf{n}} \doteq s $$

Reformulating this equation on each of the $\Omega_{N}^l$ slices is of course trivial. In the following we will keep a velocity based formulation for brevity, but keep in mind that displacement can be prescribed just as easily. In fact, one could prescribe acceleration following an analogous derivation.

A Note on Neumann Boundary Conditions

As we observed in other episodes, for example in the Rigid Walled Room episode, it often happens that FEM solutions have opposite sign with respect what is expected. This is equivalent to a $\pm \pi$ phase shift in the solution:

$$ -1 = \exp \left(j \pi \right) = \exp \left( -j \pi \right) $$

This happens fairly often when Neumann boundary conditions are applied. If this is the case for you, simply apply a $\pi$ phase shift to your velocity boundary condition. Or, which is the same, flip the sign of the Neumann terms.

On a related note, it is possible to find in literature Neumann terms formulated with the opposite sign (for example, see Computational Acoustics of Noise Propagation in Fluids, First Edition, page 8). Instead, Fundamentals of Acoustics reports the same sign as that used in this article (see page 119, Fourth Edition). The sign of this term will only affect the phase of the final solution. Hence, it doesn’t influence the solution accuracy in most cases. If phase is important for your application, solve an initial simple problem with known solution and select the sign that matches the FEM phase to the expected one.

Robin Boundary Conditions

The pressure field and the velocity field are not only conjugated through the Euler equation, but also through the specific acoustic impedance $z$. The specific acoustic impedance is defined as the ratio of the pressure field and the normal component of the particle velocity field at the boundary. Hence:

$$ z \doteq \frac{u}{\mathbf{w} \cdot \hat{\mathbf{n}}} = \frac{u}{w} \Rightarrow zw = u $$

Clearly, just like $w$ and $u$, $z$ is a complex function over the boundary. We can introduce an additional velocity component $w_{s}$ as well. This adds velocity to the field at the boundary. Hence, we can have a velocity applied while the boundary also offers an impedance. Of course, $w_{s}$ can be set to $0$ to simulate a simple impedance boundary. By using the equations in the previous section we have:

$$ z \left( w - w_s \right) = z \left( -\frac{\nabla u \cdot \hat{\mathbf{n}}}{j\omega \rho} - w_s \right) = u $$

$$ \Rightarrow \nabla u \cdot \hat{\mathbf{n}} = -j \omega \rho \left(\frac{u}{z} + w_s \right) $$

The last form of the equation shows that we are directly adding $w_s$ to the natural velocity that the field assumes due to the impedance.

Note that if we set the function $w_s$ to $0$ and rearrange we obtain:

$$ \nabla u \cdot \hat{\mathbf{n}} + \frac{j \omega \rho}{z}u = 0 $$

Which is the far field boundary condition of the Elmer solver, where Elmer’s $Z$ equals $\frac{z}{\rho}$, as observed in the previous episodes.

On the other hand, if $\left|z\right| \rightarrow +\infty$, and $w_s \neq 0$, we have a Neumann condition. We can see that a Neumann condition is that given by an ideally rigid boundary oscillating at a certain velocity. A Robin condition instead is that of an impedance boundary oscillating at a certain (eventually $0$) velocity.

As a conclusion, we can prescribe the specific acoustic impedance and the velocity at a set of boundary slices. If we do so, we are applying a Robin boundary condition. Hence, to apply a condition of this kind to $N_{R}$ slices of $\partial \Omega$ we impose:

$$ \nabla u \cdot \hat{\mathbf{n}} = -j \omega \rho \left(\frac{u}{z_r} + w_s^r \right) \\ \space r = 0, 1, \dots, N_{R} - 1 $$

More information can be found in the first chapter of Computational Acoustics of Noise Propagation in Fluids.

Model Summary

Before diving into the weak form, it is useful to collect all the equation we derived so far.

Domain and Boundaries

Domain $\Omega \in \mathbb{R}^3$ with boundary $\partial \Omega$.
$N_{D}$ Dirichlet subsets of $\partial \Omega$ denoted with $\partial \Omega_{D}^i$, $i=0,1,\dots,N_{D}-1$.
$N_{N}$ Neumann subsets of $\partial \Omega$ denoted with $\partial \Omega_{N}^l$, $l=0,1,\dots,N_{N}-1$.
$N_{R}$ Robin subsets of $\partial \Omega$ denoted with $\partial \Omega_{R}^r$, $r=0,1,\dots,N_{R}-1$.

All the subsets (slices) of $\partial\Omega$ are mutually disjointed. Their union provides $\partial\Omega$.

Material Parameters

Symbol	Name	Unit	Set
$c$	Speed of Sound	$\frac{\text{m}}{\text{s}}$	$\mathbb{R}^+$
$\rho$	Density (at rest)	$\frac{\text{kg}}{\text{m}^3} $	$\mathbb{R}^+$
$\eta$	Damping Parameter	$\frac{\text{rad}}{\text{s}} \text{m}^2 $	$\mathbb{R}$
$\alpha$	Damping Parameter	$\frac{\text{rad}}{\text{s}} $	$\mathbb{R}$

Model Parameters

Symbol	Name	Unit	Set
$\nu$	Frequency	$\text{Hz}$	$\mathbb{R}^+$
$q$	Source Function	$\frac{\text{Pa}}{\text{m}^2}$	$V^\prime$

Boundary Conditions

Type	Prescribed Quantity Symbol	Prescribed Quantity Name	Unit	Boundary Subsets	Equation (Over the Subsets)
Dirichlet	$u_{D}^i$	Pressure Field	$\text{Pa}$	$u_{D}^i : \partial \Omega_{D}^i \rightarrow \mathbb{C}$	$u \left( \mathbf{x} \right) = u_{D}^i \left( \mathbf{x} \right)$
Neumann	$w_l$	Normal Velocity Field	$\frac{\text{m}}{\text{s}}$	$w_l : \partial \Omega_{N}^l \rightarrow \mathbb{C}$	$\nabla u \cdot \hat{\mathbf{n}} = -j \omega \rho w_{l}$
Robin	$z_r$	Specific Acoustic Impedance	$\text{Rayl}_{\text{MKS}}$	$z_r : \partial \Omega_{R}^r \rightarrow \mathbb{C}$	$ \nabla u \cdot \hat{\mathbf{n}} = -j \omega \rho \left(\frac{u}{z_r} + w_s^r \right)$
Robin	$w_s^r$	External Normal Velocity Field	$\frac{\text{m}}{\text{s}}$	$w_s^r : \partial \Omega_{R}^r \rightarrow \mathbb{C}$	$ \nabla u \cdot \hat{\mathbf{n}} = -j \omega \rho \left(\frac{u}{z_r} + w_s^r \right) $

Governing PDE

$$ \begin{align} \gamma \nabla^2 u + \beta u = q \\ \gamma \doteq 1 + j k \frac{\eta}{c} \\ \beta \doteq k^2 - jk \frac{\alpha}{c} \\ k = \frac{2 \pi \nu}{c} = \frac{\omega}{c} \end{align} $$

When the equation is solved and $u$ is found, the physical solution is given by:

$$ \Re \left[ u \exp(j\omega t) \right] $$

Weak Form

The final step is to translate the PDE into its weak form. The weak formulation is what FEniCS uses as an input to define the model. A weak formulation is useful as it rewrites the PDE in terms of integrals. Domain discretization will cause these integrals to become sums. Sums can be expressed as matrix-vector operations. Hence, in a form suitable for numerical evaluation with a computer. We will leave discretization out of this episode, and focus on deriving the weak form.

To understand the weak form, we first recognise that our equation is written in these terms:

$$ Au = q $$

Where $A$ is a differential operator defined by:

$$ A \doteq \gamma \nabla^2 + \beta $$

This definition is complete only with domain and codomain for $A$.

As we discussed in the previous sections, $u$ belongs to the $V$ functional space, while $q$ belongs to its dual $V^\prime$. Hence:

$$ A : V \rightarrow V^\prime $$

Now, $A$ differentiates the ordinary function $u$ with respect space, and multiplies and adds real numbers to it. This has some consequences.

First, we can consider $V$ the space of possible solutions for our problem. In fact, we are searching for $u$ in $V$ and we know for a fact there are no solutions outside $V$, given that only functions in $V$ satisfy the basic requirements to be a solution. Now, let’s fix all the boundary conditions and, for a source term $q_1$, we find $u_1$. Then, keeping all the same, we change the source term to $q_2$ and find the solution $u_2$. The space $V$ is a vector space. So:

$$ u_1 + u_2 \in V $$

Plus, the Helmholtz PDE is linear, so the solution for $q = q_1 + q_2$ must be $u = u_1 + u_2$. This is also known as the superposition principle. But in all of this reasoning we kept all the same boundary conditions. This includes the Dirichlet boundary conditions. How can $u_1$, $u_2$ and $u = u_1 + u_2$ have all the same boundary value on the Dirichlet boundaries? The only solution can be that the boundary value must be $0$. So, we must add a requirement to $V$: all the functions in $V$ must be $0$ on the Dirichlet boundaries. How we can deal with the Dirichlet boundary condition then? We can find $u_b$ such that $Au_b = 0$ which satisfies the Dirichlet boundary conditions. Then, the final solution will be $u + u_b$. We do not need to do this ourselves, this step is automated by FEniCS. All we need to know is that the functions in $V$ are $0$ on the Dirichlet boundaries.

Second, $Au$ results in an ordinary function. So, the linear functional given by $Au$ is induced by an ordinary function. As we seen above, induction is done through inner product. We can consider also $q$ as an ordinary function inducing a functional. Even if $q$ is not induced by an ordinary function we can still use the inner product notation to apply it to any function in $V$. So, let’s apply the linear functional on both sides of $Au = q$ to an arbitrary function $\phi$ in $V$. We have:

$$ \langle{Au, \phi \rangle} = \langle{q, \phi \rangle} $$

Let’s fully expand this expression:

$$ \begin{aligned} \langle{Au, \phi \rangle} = \\ = \int_{\Omega} \left(\gamma \nabla^2 u + \beta u \right) \phi^{\star} d\mathbf{x} = \\ = \gamma \int_{\Omega} \left(\nabla^2 u \right) \phi^{\star} d\mathbf{x} + \beta \int_{\Omega} u \phi^{\star} d\mathbf{x} = \\ \langle{q, \phi \rangle} = \\ \int_{\Omega} q \phi^{\star} d\mathbf{x} \end{aligned} $$

In one line:

$$ \gamma \int_{\Omega} \left(\nabla^2 u \right) \phi^{\star} d\mathbf{x} + \beta \int_{\Omega} u \phi^{\star} d\mathbf{x} = \int_{\Omega} q \phi^{\star} d\mathbf{x} $$

The first integral has second order derivatives. We can get rid of those using integration by parts. Multidimensional integration by parts reads as follows (from Solving PDEs in Python):

$$ \int_{\Omega} \left( \nabla^2 u \right) \phi^\star d\mathbf{x} = \int_{\partial \Omega} \left( \nabla u \cdot \hat{\mathbf{n}} \right) \phi^\star ds - \int_{\Omega} \nabla u \cdot \nabla \phi^\star d\mathbf{x} $$

With $ds$ the surface element on $\partial \Omega$. Note how we formulated integration by part for $u$ and $\phi^\star$. This allows us to substitute the equation easily. Of course, the formula holds for any pair of functions in $V$, thanks to $V$ being a Sobolev space. By substitution we obtain:

$$ \gamma \int_{\partial \Omega} \left( \nabla u \cdot \hat{\mathbf{n}} \right) \phi^\star ds - \gamma \int_{\Omega} \nabla u \cdot \nabla \phi^\star d\mathbf{x} + \beta \int_{\Omega} u \phi^{\star} d\mathbf{x} = \int_{\Omega} q \phi^{\star} d\mathbf{x} $$

The first term is on the boundary. We can rewrite it by splitting it on the various boundaries we used for the boundary conditions.

$$ \begin{aligned}
\int_{\partial \Omega} \left( \nabla u \cdot \hat{\mathbf{n}} \right) \phi^\star ds = \\
= \sum_{i=0}^{N_D-1} \int_{\partial \Omega_{D}^i} \left( \nabla u \cdot \hat{\mathbf{n}} \right) \phi^\star ds + \\
+ \sum_{l=0}^{N_N-1} \int_{\partial \Omega_{N}^l} \left( \nabla u \cdot \hat{\mathbf{n}} \right) \phi^\star ds + \\
+ \sum_{r=0}^{N_R-1} \int_{\partial \Omega_{R}^r} \left( \nabla u \cdot \hat{\mathbf{n}} \right) \phi^\star ds
\end{aligned} $$

As we said, $\phi$ is $0$ on the Dirichlet boundaries. So, the first term above disappears. For the other terms, we use the equations for $\nabla u \cdot \hat{\mathbf{n}}$ we found for the Neumann and Robin conditions. We obtain:

$$ \begin{aligned}
\int_{\partial \Omega} \left( \nabla u \cdot \hat{\mathbf{n}} \right) \phi^\star ds = \\
-j \omega \rho \sum_{l=0}^{N_N-1} \int_{\partial \Omega_{N}^l} w_l \phi^\star ds + \\
- j \omega \rho \sum_{r=0}^{N_R-1} \int_{\partial \Omega_{R}^r} \left( \frac{u}{z_r} + w_s^r \right) \phi^\star ds = \\
= -j \omega \rho \sum_{l=0}^{N_N-1} \int_{\partial \Omega_{N}^l} w_l \phi^\star ds + \\
- j \omega \rho \sum_{r=0}^{N_R-1} \int_{\partial \Omega_{R}^r} \frac{u}{z_r} \phi^\star ds + \\
- j \omega \rho \sum_{r=0}^{N_R-1} \int_{\partial \Omega_{R}^r} w_s^r \phi^\star ds \end{aligned} $$

By rearranging and changing the sings we reach the final form below:

$$ \begin{aligned}
\gamma \int_{\Omega} \nabla u \cdot \nabla \phi^\star d\mathbf{x} + \\
- \beta \int_{\Omega} u \phi^{\star} d\mathbf{x} + \\
j \omega \rho \gamma \sum_{r=0}^{N_R-1} \int_{\partial \Omega_{R}^r} \frac{u}{z_r} \phi^\star ds = \\ = -\int_{\Omega} q \phi^{\star} d\mathbf{x} + \\
- j \omega \rho \gamma \sum_{l=0}^{N_N-1} \int_{\partial \Omega_{N}^l} w_l \phi^\star ds + \\
- j \omega \rho \gamma \sum_{r=0}^{N_R-1} \int_{\partial \Omega_{R}^r} w_s^r \phi^\star ds \end{aligned} $$

The equation above has the following form:

$$ a \left(u, \phi \right) = L \left(\phi \right) $$

$a$ and $L$ are called the bilinear and linear forms respectively. FEniCS will take the variational formulation in this form as an input. Additionally, the Dirichlet boundary conditions will be supplied, so that the final solution can be built.

This formulation allow us to appreciate that the Helmholtz equation by “default” is formulated for a rigid box. I.e., putting all the boundary terms to zero (zero velocity or magnitude infinite impedance) is the same as not having those terms at all in the formalism.

Note how all the integrands are of the type $gh^\star$, with $g$ and $h$ some functions. Expanding in real and imaginary parts we have:

$$ \begin{aligned}
gh^\star = \left[\Re\left(g\right) + j\Im \left(g\right)\right] \left[ \Re\left(h\right) + j\Im \left(h\right) \right]^\star = \\
= \left[ \Re\left(g\right)\Re\left(h\right) + \Im \left(g\right)\Im \left(h\right) \right] + \\ + j \left[ \Im \left(g\right) \Re\left(h\right) - \Re\left(g\right) \Im \left(h\right) \right]
\end{aligned}$$

We will need to expand all integrals in this form when solving the PDE with FEniCS, as at the time of writing FEniCS does not support complex fields. So, we will need to use mixed-spaces. FEniCSx does support complex fields, but at the moment it is a bit hard to install it.

Conclusion

This was a long episode focusing on fairly advanced topics. We selected a generalised wave equation and found a generalised Helmholtz equation. We constructed a variety of boundary conditions for it. Finally, we found the weak form. This is the initial process behind building a PDE solver with FEniCS. As we mentioned in the Intro to FEniCS - Part 1 episodes, FEniCS requires higher expertise with respect Elmer to get started. This can be thought of a con, but in reality it is also a pro. The pro being in the learning opportunities this approach offers.

In the following episode we will build a simple Helmholtz solver based on the weak form we derived, but without damping and with only Neumann boundary conditions.

References

License Information

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