Basic setting

Helmholtz equations are important for simulating wave propagation in seismic exploration, medical imaging. Usually, the Helmholtz equation has the following form:
$\begin{align} \left\{\begin{aligned} & \Delta u(x) + \kappa^2(1+q(x))u(x) = f(x) \quad \text{in } \Omega, \\ & \partial_{\mathbf{n}}u = Bu \quad \text{on }\partial\Omega, \end{aligned}\right. \end{align}$
where $B$ is an operator stands for some boundary conditions. If we try to solve problems in seismic exploration, we need to simulate waves in infinite domain, so we need perfect match layer (PML) technique to absorb waves.

On Page 9 of the following reference:
[1] G.Gao, S-N. Chow, P. Li, and H. Zhou, Numerical solution of an inverse medium scattering problem with a stochastic source, Inverse Problems, 26, 2010, 074014.
We can find some details of uniaxial PML technique for 2-D domain.

Let $D$ be the rectangle with contains $\Omega$ and let $d_1$ and $d_2$ be the thickness of the PML layers along $x$ and $y$ , respectively. Denote by $\partial D$ the boundary of the domain $D$ . Let $s_1(x)=1+i\sigma_1(x)$ and $s_{2}(y)=1+i\sigma_2(x)$ be the model medium property, where $\sigma_j$ are the positive continuous even functions and satisfy $\sigma_j(x) = 0$ in $\Omega$ .

Following the general idea in designing PML absorbing layers, we may deduce the truncated PML problem: find the PML solution, still denoted as $u$ , to the following system:
$\begin{align} \left\{\begin{aligned} & \nabla\cdot (s\nabla u) + s_1 s_2 \kappa^2 (1+q)u = f \quad \text{in }D, \\ & u = 0 \quad \text{on }\partial D, \end{aligned}\right. \tag{2.18} \end{align}$
where $s=\text{diag}(s_{2}(y)/s_{1}(x), s_{1}(x)/s_{2}(y))$ is a diagonal matrix. The variational problem can be formulated to find $u\in H_{0}^{1}(D)$ such that
$\begin{align} a_{PML}(u,v) = b(v) \quad \text{for all }v\in H_{0}^{1}(D), \\ \tag{2.19} \end{align}$
where the bilinear form $a_{PML}: H_0^1 (D)\times H_0^1(D) \rightarrow \mathbb{C}$
$\begin{align} a_{PML}(a,v)=(s\nabla u, \nabla v) - \kappa^2 (s_1 s_2 (1+qu, v)), \tag{2.20} \end{align}$
and the linear functional $b(v) = (f,v)$ .

Denote the physical domain $\Omega$ by the rectangle $[x_1, x_2]\times[y_1, y_2]$ . The computational domain is then $D=[x_1-d_1,x_2+d_1]\times [y_1-d_2, y_2+d_2]$ . The model medium property is usually taken as a power function:
$\begin{align} \sigma_1(x) = \left\{\begin{aligned} & \sigma_0 \left( \frac{x-x_2}{d_1} \right)^p \quad \text{for }x_2<x<x_2+d_1 \\ & 0 \quad \text{for }x_1 \leq x \leq x_2 \\ & \sigma_0\left( \frac{x_1-x}{d_1} \right)^p \quad \text{for }x_1-d_1 < x< x_1, \end{aligned}\right. \end{align}$
and
$\begin{align} \sigma_2(y) = \left\{\begin{aligned} & \sigma_0 \left( \frac{y-y_2}{d_2} \right)^p \quad \text{for }y_2<y<y_2+d_2 \\ & 0 \quad \text{for }y_1 \leq y \leq y_2 \\ & \sigma_0\left( \frac{y_1-y}{d_2} \right)^p \quad \text{for }y_1-d_2 < y< y_1, \end{aligned}\right. \end{align}$
where the constant $\sigma_0 > 1$ and the integer $p \geq 2$ .

We can see that the Helmholtz equation with PML has complex coefficients. However, FEniCS cannot solve PDE with complex variables directly. Hence, we need to split the complex coefficient system into real and imaginary parts.

Transformed system

Denote $x = x^R + ix^{I}$ , where $x$ can be $u$ , $s$ and so on. Relying on
$\begin{align} \frac{s_2(y)}{s_1(x)}=\frac{1+i\sigma_2(y)}{1+i\sigma_1(x)}=\frac{1+\sigma_1(x)\sigma_2(y)}{1+\sigma_1^2(x)}+i\frac{\sigma_2(y)-\sigma_1(x)}{1+\sigma_1^2(x)}, \end{align}$
and
$\begin{align} \frac{s_1(x)}{s_2(y)}=\frac{1+i\sigma_1(x)}{1+i\sigma_2(y)}=\frac{1+\sigma_1(x)\sigma_2(y)}{1+\sigma_2^2(y)}+i\frac{\sigma_1(x)-\sigma_2(y)}{1+\sigma_2^2(y)}, \end{align}$
we find that
$\begin{align} s = \text{diag}\left( \frac{s_2(y)}{s_1(x)}, \frac{s_1(x)}{s_2(y)} \right) = s^R+is^I, \end{align}$
with
$\begin{align} s^R = \text{diag}\left( \frac{1+\sigma_1(x)\sigma_2(y)}{1+\sigma_1^2(x)}, \frac{1+\sigma_1(x)\sigma_2(y)}{1+\sigma_2^2(y)} \right), \end{align}$
and
$\begin{align} s^I = \text{diag}\left( \frac{\sigma_2(y)-\sigma_1(x)}{1+\sigma_1^2(x)}, \frac{\sigma_1(x)-\sigma_2(y)}{1+\sigma_2^2(y)} \right). \end{align}$
Denote $c(x,y)=c^{R}(x,y) + ic^{I}(x,y)=s_1(x)s_2(y)=(1-\sigma_1(x)\sigma_2(y))+i(\sigma_1(x)+\sigma_2(y))$ , we know that
$\begin{align} c^R = 1-\sigma_1(x)\sigma_2(y), \quad c^{I}=\sigma_1(x)+\sigma_2(y). \end{align}$
Plugging $s$ and $c$ in the following equations
$\begin{align} \left\{\begin{aligned} & \nabla\cdot (s\nabla u) + s_1 s_2 \kappa^2 (1+q)u = f \quad \text{in }D, \\ & u = 0 \quad \text{on }\partial D, \end{aligned}\right. \end{align}$
we obtain the following system
$\begin{align} \left\{\begin{aligned} & \nabla\cdot(s^R\nabla u^R - s^I\nabla u^I)+\kappa^2(1+q)(c^R u^R - c^I u^I) = f^R \quad \text{in }D, \\ & \nabla\cdot(s^I\nabla u^R + s^R\nabla u^I)+\kappa^2(1+q)(c^I u^R + c^R u^I) = f^I \quad \text{in }D, \\ & u^R=u^I = 0 \quad \text{on }\partial D, \end{aligned}\right. \end{align}$
which can be solved by FEniCS efficiently.

FEniCS code

FEniCS is a software which can solve PDE efficiently. (https://fenicsproject.org/)
Environments: Python 3.6, FEniCS 2018.1.0
The following codes contain two class: Domain, Helmholtz

import numpy as np 
import matplotlib.pyplot as plt 
import fenics as fe

"""
Solve Helmholtz equation
"""

class Domain(object):
    def __init__(self, para = {'nx': 100, 'ny': 100, 'dPML': 0.15, \
                               'xx': 2.0, 'yy': 2.0, 'sig0': 1.5, 'p': 2.3}):
        self.nx, self.ny = para['nx'], para['ny']
        self.dPML = para['dPML']
        self.sig0, self.p = para['sig0'], para['p']
        self.xx, self.yy = para['xx'], para['yy']
        self.haveMesh = False
        
    def geneMesh(self):
        dPML, xx, yy, Nx, Ny = self.dPML, self.xx, self.yy, self.nx, self.ny
        self.mesh = fe.RectangleMesh(fe.Point(-dPML, -dPML), fe.Point(xx+dPML, yy+dPML), Nx, Ny)
        self.haveMesh = True
        
    def modifyNxNy(self, nx, ny):
        self.nx, self.ny = nx, ny
        self.haveMesh = False
        
    def numberOfMesh(self):
        if self.haveMesh == False:
            print('Mesh has not been generated!')
            return 0
        else:
            return self.nx*self.ny*2
    
    def sizeOfMesh(self):
        if self.haveMesh == False:
            print('Mesh has not been generated!')
            return 0
        else:
            xlen, ylen = self.xx/self.nx, self.yy/self.ny
            return [xlen, ylen, np.sqrt(xlen**2 + ylen**2), 0.5*xlen*ylen]
        

class Helmholtz(object):
    def __init__(self, domain, para={'kappa': 5.0}):
        self.domain = domain
        self.kappa = para['kappa']
        self.haveFunctionSpace = False
        
    def modifyDomain(self, domain):
        self.domain = domain
        self.haveFunctionSpace = False
        
    def geneFunctionSpace(self):
        P2 = fe.FiniteElement('P', fe.triangle, 2)
        element = fe.MixedElement([P2, P2])
        if self.domain.haveMesh == False:
            self.domain.geneMesh()
        self.V = fe.FunctionSpace(self.domain.mesh, element)
        self.haveFunctionSpace = True
    
    def geneForwardMatrix(self, q_fun=fe.Constant(0.0), fR=fe.Constant(0.0), \
                          fI=fe.Constant(0.0)):
        if self.haveFunctionSpace == False:
            self.geneFunctionSpace()
            
        xx, yy, dPML, sig0_, p_ = self.domain.xx, self.domain.yy, self.domain.dPML,\
                                  self.domain.sig0, self.domain.p
        # define the coefficents induced by PML
        sig1 = fe.Expression('x[0] > x1 && x[0] < x1 + dd ? sig0*pow((x[0]-x1)/dd, p) : (x[0] < 0 && x[0] > -dd ? sig0*pow((-x[0])/dd, p) : 0)', 
                     degree=3, x1=xx, dd=dPML, sig0=sig0_, p=p_)
        sig2 = fe.Expression('x[1] > x2 && x[1] < x2 + dd ? sig0*pow((x[1]-x2)/dd, p) : (x[1] < 0 && x[1] > -dd ? sig0*pow((-x[1])/dd, p) : 0)', 
                     degree=3, x2=yy, dd=dPML, sig0=sig0_, p=p_)
        
        sR = fe.as_matrix([[(1+sig1*sig2)/(1+sig1*sig1), 0.0], [0.0, (1+sig1*sig2)/(1+sig2*sig2)]])
        sI = fe.as_matrix([[(sig2-sig1)/(1+sig1*sig1), 0.0], [0.0, (sig1-sig2)/(1+sig2*sig2)]])
        cR = 1 - sig1*sig2
        cI = sig1 + sig2
        
        # define the coefficients with physical meaning
        angl_fre = self.kappa*np.pi
        angl_fre2 = fe.Constant(angl_fre*angl_fre)
        
        # define equations
        u_ = fe.TestFunction(self.V)
        du = fe.TrialFunction(self.V)
        
        u_R, u_I = fe.split(u_)
        duR, duI = fe.split(du)
        
        def sigR(v):
            return fe.dot(sR, fe.nabla_grad(v))
        def sigI(v):
            return fe.dot(sI, fe.nabla_grad(v))
        
        F1 = - fe.inner(sigR(duR)-sigI(duI), fe.nabla_grad(u_R))*(fe.dx) \
            - fe.inner(sigR(duI)+sigI(duR), fe.nabla_grad(u_I))*(fe.dx) \
            - fR*u_R*(fe.dx) - fI*u_I*(fe.dx)
        
        a2 = fe.inner(angl_fre2*q_fun*(cR*duR-cI*duI), u_R)*(fe.dx) \
             + fe.inner(angl_fre2*q_fun*(cR*duI+cI*duR), u_I)*(fe.dx) \
        
        # define boundary conditions
        def boundary(x, on_boundary):
            return on_boundary
        
        bc = [fe.DirichletBC(self.V.sub(0), fe.Constant(0.0), boundary), \
              fe.DirichletBC(self.V.sub(1), fe.Constant(0.0), boundary)]
        
        a1, L1 = fe.lhs(F1), fe.rhs(F1)
        self.u = fe.Function(self.V)
        self.A1 = fe.assemble(a1)
        self.b1 = fe.assemble(L1)
        self.A2 = fe.assemble(a2)
        bc[0].apply(self.A1, self.b1)
        bc[1].apply(self.A1, self.b1)
        bc[0].apply(self.A2)
        bc[1].apply(self.A2)
        self.A = self.A1 + self.A2
        
    def addPointSource(self, points=[(1, 2)], magnitude=[1]):
        if len(points) != len(magnitude):
            print('The length of points and magnitude must be equal!')
            
        for i in range(len(magnitude)):
            fe.PointSource(self.V.sub(0), fe.Point(points[i]), magnitude[i]).apply(self.b1)
    
    def solve(self):
        self.u = fe.Function(self.V)
        fe.solve(self.A, self.u.vector(), self.b1, 'mumps')
        self.uReal, self.uImag = self.u.split()

    def drawSolution(self):
        plt.figure(figsize=(10, 5))
        plt.subplot(1, 2, 1)
        imuR = fe.plot(self.uReal, title='Real part of the solution')
        plt.colorbar(imuR)
        plt.subplot(1, 2, 2)
        imuI = fe.plot(self.uImag, title='Imaginary part of the solution')
        plt.colorbar(imuI)
        plt.show()

The following code provide an example for using the above classes

import numpy as np 
import matplotlib.pyplot as plt 
import fenics as fe 

from HelmholtzPro import *

plt.close()

q_fun = fe.Expression('1.0+((0.5 <= x[0] && x[0] <= 1.5 && 0.5 <= x[1] && x[1] <= 1.5) ? 1 : 0)',
                      degree=3)

domain_para = {'nx': 300, 'ny': 300, 'dPML': 0.15, 'xx': 2.0, 'yy': 2.0, 'sig0': 1.5, 'p': 2.3}
equ_para = {'kappa': 20}

do = Domain(domain_para)
hel = Helmholtz(do, equ_para)
hel.geneForwardMatrix(q_fun)
hel.addPointSource()
hel.solve()

uR1 = hel.uReal
uI1 = hel.uImag

plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
fe.plot(uR1)
plt.subplot(1, 2, 2)
fe.plot(ul1)
plt.show()

Solve Helmholtz equations with PML by FEniCS

Solve Helmholtz equations with PML by FEniCS

Basic setting

Transformed system

FEniCS code

相关阅读更多精彩内容

友情链接更多精彩内容