Version:

On-the-fly Fourier Transform

This notebook validates the frequency domain output when transforming the time-dependent wavefield on the fly during the simulation, and evaluating the discrete Fourier transform (DFT) for a list of specific frequencies.
This type of output can be configured similarly to volume or surface output with salvus.flow.simple_config.
sim.output.frequency_domain.filename = "fft.h5"
sim.output.frequency_domain.format = "hdf5"
sim.output.frequency_domain.fields = ["phi"]
sim.output.frequency_domain.frequencies = [50.0, 60.0, 70.0, 100.0]
The real and imaginary parts of the wavefield are stored in the mesh file format. For each of the specified frequencies, this will generate fields FFT-Re_N, FFT-Im_N, where N is the index of the frequency.
Copy
%config Completer.use_jedi = False

import numpy as np
import pathlib

from salvus.mesh import simple_mesh
from salvus.mesh.unstructured_mesh import UnstructuredMesh
from salvus.flow import simple_config as sc
from salvus.flow import api

Setup

We consider a homogeneous 2D disk with a source centered in the domain. To compare with a delta pulse in frequency domain, we continuously inject a sine wave as source time function, let the simulation run for a while to fill up the spatial domain, and then Fourier-transform a single cycle.
vp = 1500.0
rho = 1000.0
omega = 500000.0
time_steps_per_cycle = 256
propagation_length_in_meters = 0.04
ncycles_buffer = 1

cycle_in_seconds = 1 / omega
dt = cycle_in_seconds / time_steps_per_cycle
ncycles = int(omega * propagation_length_in_meters / vp)
ncycles_total = ncycles + ncycles_buffer
mesh = simple_mesh.basic_mesh.SphericalHomogeneousAcoustic2D(
    vp=vp,
    rho=rho,
    radius=0.03,
    max_frequency=omega,
    elements_per_wavelength=6,
    tensor_order=4,
).create_mesh()
We use a single scalar point source in the center of the domain, continuously emitting a sine wave with a frequency of 500500 kHz as a source time function. To avoid numerical artifacts, we apply a tanh taper at the beginning to ramp up the amplitude. Furthermore, we let the simulation run long enough to fill up the entire domain, so that we can extract the steady-state solution in the next step.
def create_stf(ncycles_total, time_steps_per_cycle, dt, plot=False):
    nt = time_steps_per_cycle * ncycles_total
    nt_buffer = 3 * time_steps_per_cycle
    taper = np.ones([nt])
    t = np.linspace(0, 1, nt_buffer)
    taper[:nt_buffer] = np.tanh((t - 0.5) * 2.0 * np.pi) / 2.0 + 0.5

    stf = sc.stf.Custom.from_array(
        taper
        * np.sin(
            np.linspace(
                0,
                ncycles_total * 2 * np.pi,
                time_steps_per_cycle * ncycles_total,
            )
        ),
        sampling_rate_in_hertz=1 / dt,
        start_time_in_seconds=0.0,
    )

    if plot:
        stf.plot()

    return stf


stf = create_stf(ncycles_total, time_steps_per_cycle, dt, plot=True)
Note that the plotted spectrum contains frequency aliasing, because the signal is upsampled internally before the DFT is evaluated. The correct power spectrum would contain a single peak at 500 kHz and no contributions for any other frequency.
srcs = [
    sc.source.cartesian.ScalarPoint2D(
        x=0.0,
        y=0.0,
        f=1.0,
        source_time_function=stf,
    )
]
w = sc.simulation.Waveform(mesh=mesh, sources=srcs)
w.physics.wave_equation.start_time_in_seconds = 0.0
w.physics.wave_equation.end_time_in_seconds = (
    ncycles_total * time_steps_per_cycle * dt
)
w.physics.wave_equation.time_step_in_seconds = dt

bc = sc.boundary.Absorbing(
    side_sets=["r1"], width_in_meters=0.0, taper_amplitude=1.0
)
w.add_boundary_conditions(bc)

w.output.frequency_domain.fields = ["phi"]
w.output.frequency_domain.filename = "fft.h5"
w.output.frequency_domain.frequencies = [omega, 2 * omega]
w.output.frequency_domain.format = "hdf5"
w.output.frequency_domain.start_time_in_seconds = (
    ncycles * time_steps_per_cycle * dt
)
w.output.frequency_domain.end_time_in_seconds = (
    (ncycles + 1) * time_steps_per_cycle - 1
) * dt

w.validate()
w