Time Domain: Driven Trapping Square

We now revisit the trapping square in the time domain. A localized source drives the cavity for a short time and is then switched off. The animation lets us watch part of the wave leave through the hyperboloidal layer while another part rings inside the cavity.

The frequency-domain form \(k^2M+ikC+K_{\rm Helmholtz}\) becomes, with our convention \(e^{-ikt}\),

\[ M\ddot q+C\dot q+Kq=f(t),\qquad K=-K_{\rm Helmholtz}. \]

The matrix \(C\) contains the hyperboloidal transport and outflow terms. No absorbing boundary condition is added at infinity. We evolve this system with the average-acceleration Newmark method.

Geometry, resolution, and forcing

We use the same shifted cavity as in the DG notebook. Standard physical coordinates extend to R=1.5; the hyperboloidal layer starts there. The source is supported in this unchanged region and oscillates with frequency \(\omega=10.8\). We turn it off at drive_until so that the remaining cavity response can be seen clearly. With fourth-order elements, the nominal physical-region resolution is \(\omega h/p=0.27\).

from ngsolve import *
from netgen.occ import *
from ngsolve.webgui import Draw
import numpy as np
import matplotlib.pyplot as plt
import time
notebook_start = time.perf_counter()

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

# Discretization parameters
maxh = 0.1  
order = 4

# Geometry and forcing parameters
xOffset = 0.5
drive_frequency = 10.8
drive_until = 6.0
t_end = 10.0

def create_geo(Rout = Rout, RScat = RScat, DScat = DScat, DHole = DHole, R = R, xOffset = xOffset,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+xOffset,-RScat/2).Rectangle(RScat,RScat).Face()
            -MoveTo(-RScat/2+DScat/2+xOffset,-RScat/2+DScat/2).Rectangle(RScat-DScat,RScat-DScat).Face()
            -MoveTo(-RScat/2-DScat+xOffset,-DHole/2).Rectangle(2*DScat,DHole).Face()
        )
    scatterer.edges.name = 'scatterer'

    scattererdom = MoveTo(-RScat/2+xOffset,-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')

Semidiscrete hyperboloidal forms

The function below assembles the three forms from Weak Form Used in NGSolve without inserting a wavenumber. The mass weight is \(M=(1-H^2)/G\), C_ contains the mixed time-space and outflow terms, and K_ is the Helmholtz stiffness form. For time evolution we reverse its sign, \(K=-K_{\rm Helmholtz}\), and obtain \(M\ddot q+C\dot q+Kq=f\).

def get_bilinear_forms(maxh = maxh, order = order):
    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, dirichlet = "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
    Mu = 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_ = Mu*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
        )


    return M_, C_, K_,fes

Physical-region diagnostics

Inside R the transformed field is the physical field, so we can measure it directly. We record \(\left(\|\dot q\|_{L^2}^2+\|\nabla q\|_{L^2}^2\right)^{1/2}\) both in the full physical region and in the cavity alone. Comparing the two curves shows how much of the response remains trapped. These are local diagnostics, not a global energy balance for the compactified system. The four matrices used to evaluate them are assembled only once.

def assemble_energy_matrices(fes):
    u, v = fes.TnT()
    return {
        "mass_cavity": BilinearForm(
            u * v * dx("cavity"), check_unused=False
        ).Assemble().mat,
        "mass_physical": BilinearForm(
            u * v * dx("inner|cavity"), check_unused=False
        ).Assemble().mat,
        "stiffness_cavity": BilinearForm(
            grad(u) * grad(v) * dx("cavity"), check_unused=False
        ).Assemble().mat,
        "stiffness_physical": BilinearForm(
            grad(u) * grad(v) * dx("inner|cavity"), check_unused=False
        ).Assemble().mat,
    }


def compute_energy_norms(displacement, velocity, matrices):
    norm_physical_sq = (
        velocity.vec.InnerProduct(matrices["mass_physical"] * velocity.vec)
        + displacement.vec.InnerProduct(matrices["stiffness_physical"] * displacement.vec)
    ).real
    norm_cavity_sq = (
        velocity.vec.InnerProduct(matrices["mass_cavity"] * velocity.vec)
        + displacement.vec.InnerProduct(matrices["stiffness_cavity"] * displacement.vec)
    ).real
    return np.sqrt(max(0, norm_physical_sq)), np.sqrt(max(0, norm_cavity_sq))

Average-acceleration Newmark system

We use the average-acceleration choice \(\beta=1/4\) and \(\gamma=1/2\). Every time step reuses the same effective inverse

\[ S^{-1}=\left(M+\gamma\Delta t C+\beta\Delta t^2K\right)^{-1}. \]

The initial displacement, velocity, acceleration, and load all vanish. We therefore don’t need a separate mass-matrix inverse. In the code we keep the positive time-domain stiffness K explicit so that its sign relative to the Helmholtz form remains visible.

def get_newmark_beta_matrices(mass_form, transport_form, helmholtz_stiffness_form, dt, gamma, beta, freedofs):
    time_stiffness_form = (-1) * helmholtz_stiffness_form
    C = BilinearForm(transport_form).Assemble().mat
    K = BilinearForm(time_stiffness_form).Assemble().mat
    effective_form = mass_form + gamma * dt * transport_form + beta * dt**2 * time_stiffness_form
    Sinv = BilinearForm(effective_form).Assemble().mat.Inverse(freedofs=freedofs)
    return C, K, Sinv

dt = 2e-2

beta = 1/4
gamma = 1/2
assembly_start = time.perf_counter()
mass_form, transport_form, helmholtz_stiffness_form, fes = get_bilinear_forms(maxh=maxh, order=order)
C, K, Sinv = get_newmark_beta_matrices(
    mass_form, transport_form, helmholtz_stiffness_form, dt, gamma, beta, fes.FreeDofs()
)
energy_matrices = assemble_energy_matrices(fes)
assembly_seconds = time.perf_counter() - assembly_start
print(f"degrees of freedom: {fes.ndof}")
print(f"omega L = {drive_frequency * RScat:.1f}, nominal physical-region omega h/p = {drive_frequency * maxh / order:.3f}")
print(f"matrix assembly and factorization: {assembly_seconds:.3f} s")

degrees of freedom: 21898
omega L = 10.8, nominal physical-region omega h/p = 0.270
matrix assembly and factorization: 0.483 s

Driven evolution and ringdown

The Newmark loop first predicts the displacement and velocity. It then needs one multiplication by K, one by C, and one application of the effective inverse per step. The time step is \(\Delta t=0.02\). We evaluate the diagnostic curves every \(0.5\) time units and save an animation frame every \(1.0\) time unit; this keeps the visualization inexpensive without changing the evolution.

def run_simulation(omega, drive_until=drive_until, tend=t_end):
    source_profile = exp(-30 * ((x + 0.8)**2 + (y + 0.8)**2))
    source_amplitude = 100.0
    force = LinearForm(
        source_amplitude * source_profile * fes.TestFunction() * dx("inner|cavity")
    ).Assemble().vec
    effective_force = Sinv * force

    displacement = GridFunction(fes, name="u")
    velocity = GridFunction(fes, name="ut")
    acceleration = GridFunction(fes, name="utt")  # zero for the present initial data
    old_acceleration = acceleration.vec.CreateVector()
    predicted_displacement = acceleration.vec.CreateVector()
    predicted_velocity = acceleration.vec.CreateVector()
    acceleration_rhs = acceleration.vec.CreateVector()
    animation = GridFunction(fes, multidim=0, name="u_history")

    diagnostic_stride = max(1, int(round(0.5 / dt)))
    frame_stride = max(1, int(round(1.0 / dt)))
    sample_times = []
    physical_norms = []
    cavity_norms = []

    # Average-acceleration Newmark method
    nsteps = int(round(tend / dt))
    for step in range(nsteps):
        t = step * dt
        if step % diagnostic_stride == 0:
            norm_inner, norm_scat = compute_energy_norms(
                displacement, velocity, energy_matrices
            )
            sample_times.append(t)
            physical_norms.append(norm_inner)
            cavity_norms.append(norm_scat)
            print(
                f"t={t:5.2f}, physical norm={norm_inner:.4e}, cavity norm={norm_scat:.4e}",
                end="\r",
            )
        if step % frame_stride == 0:
            animation.AddMultiDimComponent(displacement.vec)

        old_acceleration.data = acceleration.vec
        predicted_displacement.data = displacement.vec
        predicted_displacement.data += dt * velocity.vec
        predicted_displacement.data += (0.5 - beta) * dt**2 * old_acceleration
        predicted_velocity.data = velocity.vec
        predicted_velocity.data += (1 - gamma) * dt * old_acceleration

        acceleration_rhs.data = -K * predicted_displacement
        acceleration_rhs.data += -C * predicted_velocity
        acceleration.vec.data = Sinv * acceleration_rhs
        if t + dt < drive_until:
            acceleration.vec.data += np.sin(omega * (t + dt)) * effective_force

        displacement.vec.data = predicted_displacement
        displacement.vec.data += beta * dt**2 * acceleration.vec
        velocity.vec.data = predicted_velocity
        velocity.vec.data += gamma * dt * acceleration.vec

    norm_inner, norm_scat = compute_energy_norms(
        displacement, velocity, energy_matrices
    )
    sample_times.append(tend)
    physical_norms.append(norm_inner)
    cavity_norms.append(norm_scat)
    animation.AddMultiDimComponent(displacement.vec)
    print(f"t={tend:5.2f}, physical norm={norm_inner:.4e}, cavity norm={norm_scat:.4e}")
    print("simulation done")
    return animation, np.array(sample_times), np.array(physical_norms), np.array(cavity_norms)

evolution_start = time.perf_counter()
gf_anim, times, physical_norms, cavity_norms = run_simulation(drive_frequency)
evolution_seconds = time.perf_counter() - evolution_start
print(f"Newmark evolution ({int(round(t_end/dt))} steps): {evolution_seconds:.3f} s")

t= 0.00, physical norm=0.0000e+00, cavity norm=0.0000e+00

t= 0.50, physical norm=4.8036e+00, cavity norm=5.5217e-07

t= 1.00, physical norm=5.4721e+00, cavity norm=4.9064e-02

t= 1.50, physical norm=6.4727e+00, cavity norm=2.9197e-01

t= 2.00, physical norm=7.1669e+00, cavity norm=3.1320e-01

t= 2.50, physical norm=7.0191e+00, cavity norm=3.8893e-01

t= 3.00, physical norm=6.9423e+00, cavity norm=4.9165e-01

t= 3.50, physical norm=7.3249e+00, cavity norm=5.5310e-01

t= 4.00, physical norm=7.2621e+00, cavity norm=6.3141e-01

t= 4.50, physical norm=6.9134e+00, cavity norm=7.3821e-01

t= 5.00, physical norm=7.1242e+00, cavity norm=8.3078e-01

t= 5.50, physical norm=7.3981e+00, cavity norm=9.0133e-01

t= 6.00, physical norm=7.0277e+00, cavity norm=1.0038e+00

t= 6.50, physical norm=5.6422e+00, cavity norm=1.1167e+00

t= 7.00, physical norm=4.4677e+00, cavity norm=1.1901e+00

t= 7.50, physical norm=3.3639e+00, cavity norm=1.2093e+00

t= 8.00, physical norm=2.2067e+00, cavity norm=1.2001e+00

t= 8.50, physical norm=1.4274e+00, cavity norm=1.1890e+00

t= 9.00, physical norm=1.2149e+00, cavity norm=1.1713e+00

t= 9.50, physical norm=1.2025e+00, cavity norm=1.1574e+00

t=10.00, physical norm=1.1904e+00, cavity norm=1.1476e+00
simulation done
Newmark evolution (500 steps): 10.405 s

Draw(gf_anim,min=-0.1,max=0.1,animate=True)

fig, ax = plt.subplots(figsize=(7, 3.5))
ax.semilogy(times, physical_norms, label="physical region")
ax.semilogy(times, cavity_norms, label="cavity")
ax.axvline(drive_until, color="0.4", ls="--", label="source switched off")
ax.set(xlabel="time", ylabel="energy norm")
ax.grid(True, which="both", alpha=0.3)
ax.legend()
plt.show()
print(f"notebook compute time after imports: {time.perf_counter() - notebook_start:.3f} s")

notebook compute time after imports: 11.879 s

Interpretation

The animation shows the compactified time-domain field. It’s the physical field inside the layer interface and remains regular all the way to null infinity. After the source is switched off, most of the radiation leaves through the compactified exterior. The resonant component continues to ring inside the cavity and dies slowly as shown by the two norm curves.