Galbrun's equation: $H^1$-conforming discretizations¶

This notebooks contains an implementation of the $H^1$ conforming discretization of Galbrun's equation from [1].

We start by importing the basic functionality of the finite element library $\texttt{ngsolve}$

In [1]:
import sys
from math import pi


#NGSolve related imports 
from ngsolve import *
from netgen.geom2d import unit_square, SplineGeometry
from netgen.geom2d import *
from ngsolve.meshes import *
from ngsolve import *

#Setting ngsolve options
ngsglobals.msg_level=2
SetNumThreads(6)
In [2]:
#Custom imports
from barycenter_ref import barycenter_ref
from Galbrun import *

The following code sets general parameters for our example. Specifially, we have choosen:

density $\rho$ $1.5+0.2 \cos(\frac{\pi}{4}x)\cdot\sin(\frac{\pi}{2}y)$
wave number $\omega$ $0.78\cdot2\cdot\pi$
angular velocity of the frame $\Omega$ $\begin{pmatrix} 0 \\ 0 \end{pmatrix}$
sound speed $c_s$ $\sqrt{1.44+0.16\cdot\rho}$
pressure $p$ $1.44 \cdot \rho + 0.08\cdot \rho^2$
damping coefficient $\gamma$ $0.1$
In [3]:
#Density
rho = 1.5 + 0.2*cos(pi*x/4)*sin(pi*y/2)
#wave number
omega = 0.78*2*pi
#Angular velocity
Omega = CF((0,0))
#Sound speed
cs = sqrt(1.44+0.16*rho)
#Pressure
pres=1.44*rho+0.08*rho**2
#Damping coefficient
gamma = 0.1

Now, we set the geometry for our example. Here, we have choosen the square $[-4,4]^2 \subset \mathbb{R}^2$.

In [4]:
#Geometry description
geo = SplineGeometry()
pnts = [ (-4,-4), (4,-4), (4,4), (-4,4) ]
pnums = [geo.AppendPoint(*p) for p in pnts]
lbot = geo.Append ( ["line", pnums[0], pnums[1]],bc="bnd")
lright = geo.Append ( ["line", pnums[1], pnums[2]], bc="bnd")
geo.Append ( ["line", pnums[2], pnums[3]], bc="bnd")
geo.Append ( ["line", pnums[3], pnums[0]], bc="bnd");

TThe following function solves Galbrun's equation using a $H^1$ conforming discretization as introduced in [1]. To simplify the presentation in this notebook, the definitions of the bilinear forms $a(\cdot,\cdot)$ and the linear form $f(\cdot)$ are outsourced and can be found in the script $\texttt{Galbrun.py}$, which is called by the $\texttt{H1Assemble}$ function. This function allows you to specify the order of approximation $k$ and the parameter $\alpha$, which scales the background flow $b$ defined by

\begin{equation} b = \frac{\alpha}{\rho} \begin{pmatrix} \sin(\pi x) \cos(\pi y) \\ -\cos(\pi x) \sin(\pi y) \end{pmatrix}. \end{equation}
Further, we calculate the source term is this correct in the script? such that the solution is given by

\begin{equation} \frac{1}{\rho} \begin{pmatrix} (1+i)g \\ -(1+i)g \end{pmatrix}, \end{equation}
where $g$ is the Gaussian defined by

\begin{equation} g(x,y) = \sqrt{\frac{\log(10^6)}{\pi}} \exp(-\log(10^6)(x^2+y^2)). \end{equation} Finally, the function allows the user to specify whether a barycentric refinement should be applied to the mesh.

In [5]:
def Solve_Galbrun(alpha, order, bary):
    mesh = Mesh(geo.GenerateMesh(maxh=0.5))
    # Generating the mesh
    if bary: 
        mesh = barycenter_ref(mesh) 

    #Gaussian
    ga = sqrt(log(10**6)/pi)*exp( -log(10**6)*(x**2+y**2) )
    b = alpha/rho * CF(( sin(pi*x)*cos(pi*y), -cos(pi*x)*sin(pi*y) ))
    
    # For calculating the source term
    convcf = lambda b, f, omega, Omega: (b[0] * f.Diff(x) + b[1] * f.Diff(y)) - 1j*omega*f
    
    #Source term 
    s = convcf( b, CF(( 1 , 0 )) * ga, omega, Omega )
    
    ## Assembling the linear system
    fes = VectorH1(mesh,order=order,dirichlet="",complex=True)
    a,f,c = H1Assemble(fes,b,s,rho,cs,omega,Omega,gamma,p=pres,p_reduction=False,nitsche=True,precond=None)

    ## Solve the linear system
    gfu = GridFunction(fes)
    gfu.vec.data = a.mat.Inverse() * f.vec
    
    # Drawing the mesh and the solution
    from ngsolve.webgui import Draw
    print("------------------------------------- Barycentric refinement: {} ----------------------------------".format(bary))
    Draw(mesh)
    print("------------------------------------- Solution with k = {} ----------------------------------------".format(order,alpha))
    Draw(gfu,min=-0.03, max=0.03, autoscale=True, settings={"mesh":False})

Now it's your turn to experiment with the parameters of the function $\texttt{Solve_Galbrun}$ using the interactive widget below. Each time you adjust a parameter, two images will be displayed. The first image shows the mesh type, while the second image displays the numerical solution. Turn off the barycentric mesh refinement and observe the difference! You can also increase the order of the approximation and the value of $\alpha$ to observe their effects. Below, you will find some explanations of the differences that you might observe.

In [6]:
import ipywidgets as widgets
from ipywidgets import HBox, VBox
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display
%matplotlib inline

@widgets.interact(
    alpha=(0.1,1),order=range(2,5))
def plot(alpha=.1, order=2, bary=True):
    Solve_Galbrun(alpha, order, bary)

By experimenting with the interactive widget above, you may have observed the following phenomena:

  • Disabling the barycentric mesh refinement results in a decrease in the quality of the approximation. However, increasing the order of approximation to either $k=3$ or $k=4$ eliminates this difference.
  • The solution's quality degrades as the value of $\alpha$ increases.

[1] Martin Halla, Christoph Lehrenfeld, and Paul Stocker. A new T-compatibility condition and its application to the discretization of the damped time-harmonic Galbrun’s equation. 2022. ArXiv Preprint.

Code by Martin Halla, Christoph Lehrenfeld, and Paul Stocker. Written by Tim van Beeck. Published under the MIT license.

In [ ]: