Helmholtz problem with PML

helmholtz_pml.py

# Copyright (C) 2023 Igor A. Baratta
# SPDX-License-Identifier:    MIT

from functools import partial

import dolfinx
import numpy as np
import scipy.sparse
from dolfinx.fem import Function, FunctionSpace
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
from ufl import Measure, TestFunction, TrialFunction, dx, grad, inner

from mesh_generation import generate_mesh

scalar_type = np.complex128


k0 = 10                     # wavenumber
c = 340                     # wavespeed m/s
omega = c * k0              # angular frequency
wavelength = 2*np.pi / k0  
radius = 2 * wavelength     # scatterer radius
ref_ind = 1.2               # refractive index of scatterer
dim_x = 20*wavelength       # width of computational domain
degree = 2                  # polynomial degree
ne = 5                      # number of mesh elements per wavelength
w = 2 * wavelength          # depth of absorber
h_elem = wavelength / ne    # target mesh size
a_star = dim_x/2


# MPI communicator
mesh_order = 1
comm = MPI.COMM_WORLD
file_name = "domain.msh"

# generate mesh and read from file
generate_mesh(file_name, radius, dim_x, w, h_elem, mesh_order)
mesh, cell_tags, _ = gmshio.read_from_msh(file_name, comm, rank=0, gdim=2)

# Define function space
V = FunctionSpace(mesh, ("Lagrange", degree))


def gamma(component, x):
    eps = 2.22e-16
    sigma_i = c/(a_star - np.abs(x[component]) + eps)
    return 1j/omega * c/sigma_i


corners = cell_tags.find(1)
horizontal = cell_tags.find(2)
vertical = cell_tags.find(3)

W0 = FunctionSpace(mesh, ("DG", 0))
gamma_x = Function(W0, dtype=scalar_type)
gamma_x.interpolate(partial(gamma, 0), corners)
gamma_x.interpolate(partial(gamma, 0), vertical)
gamma_x.x.array[:] += 1


gamma_y = Function(W0, dtype=scalar_type)
gamma_y.interpolate(partial(gamma, 1), corners)
gamma_y.interpolate(partial(gamma, 1), horizontal)
gamma_y.x.array[:] += 1


# Interpolate wavenumber k onto V
W = FunctionSpace(mesh, ("DG", 0))
k = Function(W, dtype=scalar_type)

k.x.array[:] = k0
# wavenumber inside the scatterer = ref_ind * k0 
k.x.array[cell_tags.find(4)] = ref_ind * k0 


# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

dx = Measure("dx", subdomain_data=cell_tags, domain=mesh)
dx_omega = dx(4) + dx(5)
dx_pml = dx(1) + dx(2) + dx(3)

a_omega = (inner(grad(u), grad(v)) - k**2 * inner(u, v)) * dx_omega
a_pml = -gamma_x * gamma_y * k**2 * inner(u, v) * dx_pml
a_pml += gamma_y/gamma_x * inner(u.dx(0), v.dx(0)) * dx_pml
a_pml += gamma_x/gamma_y * inner(u.dx(1), v.dx(1)) * dx_pml

a = a_omega + a_pml
a = dolfinx.fem.form(a, dtype=scalar_type)


# Interpolate incident wave field onto V
ui = Function(V, dtype=scalar_type)
ui.interpolate(lambda x: np.exp(1.0j * k0 * x[0]))

L = - inner(grad(ui), grad(v)) * dx + k**2 * inner(ui, v) * dx
L = dolfinx.fem.form(L, dtype=scalar_type)


# Assemble matrix and vector and set up direct solver
A = dolfinx.fem.assemble_matrix(a)
A.finalize()

# Create a scipy sparse matrix that shares data with A
As = scipy.sparse.csr_matrix((A.data, A.indices, A.indptr))


# Create RHS vector
b = dolfinx.fem.assemble_vector(L)

# Solve the variational problem and return the solution
uh = dolfinx.fem.Function(V, dtype=scalar_type)
uh.x.array[:] = scipy.sparse.linalg.spsolve(As, b.array)

# U[Total] = U[scattered] + U[incident]
uh.x.array[:] += ui.x.array[:]

with XDMFFile(comm, "out/uh.xdmf", "w") as file:
    file.write_mesh(mesh)
    file.write_function(uh)

mesh_generation.py

from mpi4py import MPI

try:
    import gmsh
except ModuleNotFoundError:
    import sys
    print("This demo requires gmsh to be installed")
    sys.exit(0)


def generate_mesh(filename: str, radius: float, length: float, w: float, h_elem: float, order: int, verbose: bool = False):
    if MPI.COMM_WORLD.rank == 0:

        gmsh.initialize()
        gmsh.model.add("helmholtz_domain")
        gmsh.option.setNumber("General.Terminal", verbose)
        # Set the mesh size
        gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.75*h_elem)


        # Add domain
        r0 = gmsh.model.occ.addRectangle(-length/2, -length/2, 0.0, length, length)
        lw = length - w
        r1 = gmsh.model.occ.addRectangle(-length/2, length/2 - w, 0.0, length, w)
        r2 = gmsh.model.occ.addRectangle(-length/2, -length/2, 0.0, length, w)
        r3 = gmsh.model.occ.addRectangle(-length/2, -length/2, 0.0, w, length)
        r4 = gmsh.model.occ.addRectangle(length/2 - w, -length/2, 0.0, w, length)

        rfrag = gmsh.model.occ.fragment([(2, r1), (2, r2), (2, r3), (2, r4)], [(2, r0)])

        # Add scatterers
        c1 = gmsh.model.occ.addCircle(0.0, 0.0, 0.0, radius)
        gmsh.model.occ.addCurveLoop([c1], tag=c1)
        gmsh.model.occ.addPlaneSurface([c1], tag=c1)

        c2 = gmsh.model.occ.fragment([(2, c1)], [rfrag[0][-1]])

        gmsh.model.occ.synchronize()

        # Add physical groups
        gmsh.model.addPhysicalGroup(2, [rfrag[0][0][1], rfrag[0][2][1], rfrag[0][3][1], rfrag[0][5][1]], tag=1)
        gmsh.model.addPhysicalGroup(2, [rfrag[0][1][1], rfrag[0][4][1]], tag=2)
        gmsh.model.addPhysicalGroup(2, [rfrag[0][6][1], rfrag[0][7][1]], tag=3)
        gmsh.model.addPhysicalGroup(2, [c2[0][0][1]], tag=4)
        gmsh.model.addPhysicalGroup(2, [c2[0][1][1]], tag=5)

        # Generate mesh
        gmsh.model.mesh.generate(2)
        gmsh.model.mesh.setOrder(order)
        gmsh.model.mesh.optimize("HighOrder")
        gmsh.write(filename)

        gmsh.finalize()