Solve Helmholtz equations with P

作者: jjx323 | 来源:发表于2019-02-11 22:52 被阅读0次

Solve Helmholtz equations with P
有趣的 Javascript 库
VR英语试译 - Developers Test Running
Burny Math
Quotes From Charlie Munger
[HFE] U5L4 Homework
[HFE] U2L4 Homework
[HFE] U3L4 Homework
[HFI] U5L4 Homework
[HFE] U4L4 Homework

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 P
Basic setting Helmholtz equations are important for simul...
有趣的 Javascript 库
algebra.js Build, display, and solve algebraic equations ...
VR英语试译 - Developers Test Running
Developers Test Running in Place To Solve VR’s Movement P...
Burny Math
We always encounter situations where we must solve math p...
Quotes From Charlie Munger
Don't drift into self-pity because it doesn't solve any p...
[HFE] U5L4 Homework
本次学习内容：U5L4 灰色练习册 Workbook - P63, Activity A. Solve the p...
[HFE] U2L4 Homework
本次作业：灰色练习册 Workbook - P27, Activity A. Solve the puzzle....
[HFE] U3L4 Homework
本次作业：灰色练习册 Workbook - P39, Activity A. Solve the puzzle....
[HFI] U5L4 Homework
本次作业：灰色练习册： Workbook - P63. activity A.Solve the puzzle....
[HFE] U4L4 Homework
本次作业：灰色练习册 Workbook - P51, Activity A. Solve the puzzle....

Solve Helmholtz equations with P

Basic setting

Transformed system

FEniCS code

相关文章

Solve Helmholtz equations with P

有趣的 Javascript 库

VR英语试译 - Developers Test Running

Burny Math

Quotes From Charlie Munger

[HFE] U5L4 Homework

[HFE] U2L4 Homework

[HFE] U3L4 Homework

[HFI] U5L4 Homework

[HFE] U4L4 Homework

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读