Frequency Domain: Trapping Square

We solve the sound-soft Exterior Scattering Problem problem for a nonconvex trapping geometry: a square cavity with a narrow aperture. At frequencies close to a cavity resonance, the total field is strongly amplified inside the square. For comparison, we observe that a nearby frequency is non-resonant. As a simple quantitative diagnostic we report the fraction of the total-field \(L^2\) mass inside the cavity.

The incident plane wave is \(U^{\rm inc}=e^{ik a\cdot x}\). We compare the hyperboloidal layer with a radial PML on the same physical mesh and with the same finite-element order. Both solvers return \(U\) in the unchanged physical region.

Note that compactification does not remove the difficulty caused by trapping itself. Near resonances, the Helmholtz system can become poorly conditioned, and in the time domain the trapped signal can persist for a long time. Hyperboloidal compactification is designed for treating the unbounded exterior; it doesn’t make the interior dynamics easier. But it does allow long-time simulations without artificial reflections from a finite outer boundary.

Geometry and numerical parameters

A square wall of thickness DScat surrounds an interior cavity except for an aperture of width DHole. The unchanged physical region extends to radius R; the annulus from R to Rout is treated either as a hyperboloidal layer or as a radial PML. At \(k=10.8\) the unit-width cavity has nondimensional frequency \(kL=10.8\). With \(p=4\) and maxh=0.1, the nominal physical-region indicator is \(kh/p=0.27\).

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
import matplotlib.pyplot as pl
import numpy as np

# domain parameters
Rout = 2
RScat = 1
DScat = 0.1
DHole = 0.2
R=1.5

# Discretization parameters
maxh = 0.1  
order = 4
alpha = 1j  # radial PML stretching parameter

# Problem parameters
k = 10.8        #10.8 resonant, 10.6 non-resonant
theta = 0

def create_geo(Rout = Rout, RScat = RScat, DScat = DScat, DHole = DHole, R = R, draw = True):
    domain = Circle((0, 0), Rout).Face()
    domain.edges.name = 'outer'
    domain.faces.name = 'outer'
    
    inner = Circle((0, 0), R).Face()
    inner.edges.name = 'interface'
    inner.faces.name = 'inner'

    scatterer = (
            MoveTo(-RScat/2,-RScat/2).Rectangle(RScat,RScat).Face()
            -MoveTo(-RScat/2+DScat/2,-RScat/2+DScat/2).Rectangle(RScat-DScat,RScat-DScat).Face()
            -MoveTo(-RScat/2-DScat,-DHole/2).Rectangle(2*DScat,DHole).Face()
        )
    scatterer.edges.name = 'scatterer'

    scattererdom = MoveTo(-RScat/2,-RScat/2).Rectangle(RScat,RScat).Face()-scatterer
    scattererdom.faces.name = 'cavity'
    


    if draw:
        Draw(Glue([domain-inner,inner-scattererdom-scatterer, scattererdom-scatterer]))
    geo = OCCGeometry(Glue([domain-inner,inner-scattererdom-scatterer, scattererdom]), dim=2)
    return geo

geo = create_geo(draw = False)
mesh = Mesh(geo.GenerateMesh(maxh=maxh))
mesh.Curve(order)
Draw(mesh)
print(mesh.GetMaterials(),mesh.GetBoundaries())

('outer', 'inner', 'cavity') ('outer', 'interface', 'default', 'scatterer', 'scatterer', 'scatterer', 'scatterer', 'scatterer', 'scatterer', 'scatterer', 'scatterer', 'scatterer', 'scatterer', 'scatterer', 'scatterer')

Hyperboloidal-layer discretization

The first solver uses the same \(k^2M+ikC+K_{\rm Helmholtz}\) grouping as Weak Form Used in NGSolve. In the unchanged region, \(\Omega=1\), \(H=0\), and the transformed scattered field is the physical scattered field. In the layer, \(H=1-G\) makes \(M=(1-H^2)/G=1+H\) bounded. The boundary term in C_ is the transformed outflow contribution at infinity; no essential condition is imposed on outer.

def solve_scattering(maxh = maxh, order = order, k = k, theta = theta, draw = True):
    geo = create_geo(Rout = Rout, RScat = RScat, DScat = DScat, DHole = DHole, R = R, draw = False)
    mesh = Mesh(geo.GenerateMesh(maxh=maxh))
    mesh.Curve(order)

    fes = H1(mesh,order = order,complex = True,dirichlet = mesh.Boundaries("scatterer"))

    w,w_ = fes.TnT()

    rho = sqrt(x**2+y**2)

    #Omega = CF([(Rout-rho)/(Rout-R),1,1])
    Omega = mesh.MaterialCF({"outer": (Rout-rho)/(Rout-R)},default=1)

    DOmega = Omega.Diff(rho)
    L = Omega-rho*DOmega
    G = Omega**2/L

    #H = CF([1 - G,0,0]) #characteristic preserving
    H = mesh.MaterialCF({"outer": 1-G},default=0) #characteristic preserving
    M = 1+H  # (1-H**2)/G = 1+H when H=1-G


    vx = CF((x,y))

    Q = OuterProduct(vx,vx)/rho**2 # radial projector
    P = Id(2)-Q                    # tangential projector

    ### k**2
    M_ = M*w*w_*dx

    ### 1j*k
    C_ = (
        -H/rho*w*InnerProduct(vx,grad(w_))*dx
        +H/rho*w_*InnerProduct(vx,grad(w))*dx
        +w_*w*ds('outer')
        )

    ## 1
    K_ = (
        -grad(w)*( (G*Q+L*P)*grad(w_))*dx
        -Omega/L/2/rho*DOmega*InnerProduct(grad(w),vx)*w_*dx
        -Omega/L/2/rho*DOmega*InnerProduct(grad(w_),vx)*w*dx
        -1/L/4*DOmega**2*w*w_*dx
        )

    A = BilinearForm(k**2*M_+1j*k*C_+K_).Assemble().mat
    Ainv = A.Inverse(freedofs=fes.FreeDofs())

    direction = CF((cos(theta),sin(theta)))

    incident = exp(1j*k*InnerProduct(direction,(x,y))/Norm(direction))
    scattered = GridFunction(fes)

    scattered.Set(-incident, definedon=mesh.Boundaries('scatterer'))
    res = - A * scattered.vec
    scattered.vec.data += Ainv*res
    total = GridFunction(fes)
    total.Set(incident)
    total.vec.data += scattered.vec
    if draw:
        Draw(scattered,mesh,"hyperboloidal transformed scattered field", min=-1, max=1)
    return scattered,total

Radial-PML comparison

For sanity check, we replace the hyperboloidal annulus by NGSolve’s radial PML on the same geometry, mesh scale, and polynomial order. The PML field is terminated by a homogeneous condition at outer. Its parameter alpha=1j is not tuned here, so the comparison tests agreement of the physical response rather than relative efficiency.

def solve_scattering_pml(maxh = maxh, order = order, k = k, theta = theta, draw = True, alpha = alpha):
    geo = create_geo(Rout = Rout, RScat = RScat, DScat = DScat, DHole = DHole, R = R, draw = False)
    mesh = Mesh(geo.GenerateMesh(maxh=maxh))
    mesh.Curve(order)
    
    mesh.SetPML(pml.Radial((0,0),R,alpha),'outer')

    fes = H1(mesh,order = order,complex = True,dirichlet = mesh.Boundaries("scatterer|outer"))

    w,w_ = fes.TnT()
    M_ = w*w_*dx(bonus_intorder=4)    #bonus_intorder necessary due to additional coefficients from trafos (hidden in PML mesh trafo)

    K_ = (
        -grad(w)*grad(w_)*dx(bonus_intorder=4)
        )

    A = BilinearForm(k**2*M_+K_).Assemble().mat
    
    Ainv = A.Inverse(freedofs=fes.FreeDofs())

    direction = CF((cos(theta),sin(theta)))

    incident = exp(1j*k*InnerProduct(direction,(x,y))/Norm(direction))
    scattered = GridFunction(fes)

    scattered.Set(-incident, definedon=mesh.Boundaries('scatterer'))
    res = - A * scattered.vec
    scattered.vec.data += Ainv*res
    total = GridFunction(fes)
    total.Set(incident)
    total.vec.data += scattered.vec
    if draw:
        Draw(scattered,mesh,"PML scattered field", min=-1, max=1)
    return scattered,total

Cavity response

For each exterior treatment we solve at \(k=10.8\) and \(k=10.6\) with incidence angle \(\pi/4\) and report the fraction of total-field \(L^2\) mass inside the cavity relative to the full unchanged physical region. Agreement of this scalar response confirms that both exterior formulations resolve the large resonant amplification and the nearby low-response case.

def cavity_mass_fraction(total):
    mesh = total.space.mesh
    cavity = Integrate(total * Conj(total), mesh, definedon=mesh.Materials("cavity")).real
    physical = Integrate(total * Conj(total), mesh, definedon=mesh.Materials("inner|cavity")).real
    return cavity / physical

hyp_res_s, hyp_res_total = solve_scattering(k=10.8, theta=pi/4, draw=True)
hyp_non_s, hyp_non_total = solve_scattering(k=10.6, theta=pi/4, draw=True)
print(f"hyperboloidal k=10.8: cavity mass fraction = {cavity_mass_fraction(hyp_res_total):.4f}")
print(f"hyperboloidal k=10.6: cavity mass fraction = {cavity_mass_fraction(hyp_non_total):.4f}")

hyperboloidal k=10.8: cavity mass fraction = 0.6869
hyperboloidal k=10.6: cavity mass fraction = 0.0312

pml_res_s, pml_res_total = solve_scattering_pml(k=10.8, theta=pi/4, draw=True)
pml_non_s, pml_non_total = solve_scattering_pml(k=10.6, theta=pi/4, draw=True)
print(f"PML k=10.8: cavity mass fraction = {cavity_mass_fraction(pml_res_total):.4f}")
print(f"PML k=10.6: cavity mass fraction = {cavity_mass_fraction(pml_non_total):.4f}")

PML k=10.8: cavity mass fraction = 0.6899
PML k=10.6: cavity mass fraction = 0.0317