Version:
This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.

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
There should be something here! The notebook has been saved without widget output.
<salvus.flow.simple_config.simulation.Waveform object at 0x7fe42e987090>
During the simulation, the discrete Fourier transform is applied to the primary field phi and will record real and imaginary part at the specified frequencies of 500 and 1000 kHz.
api.run(
    input_file=w,
    site_name="local",
    ranks=4,
    output_folder="out",
    overwrite=True,
)
Job `job_2106022313362312_7467dbcc41` running on `local` with 4 rank(s).
Site information:
  * Salvus version: 0.11.34
  * Floating point size: 32
                                          
There should be something here! The notebook has been saved without widget output.
* Downloaded 20.4 MB of results to `out`.
* Total run time: 12.32 seconds.
* Pure simulation time: 11.98 seconds.
<salvus.flow.sites.salvus_job.SalvusJob at 0x7fe4307e8390>
The output is written to a file with the same format as the mesh, so we can either load it into an UnstructuredMesh object or open the associated xdmf file with Paraview.
For the validation we cut out the areas close to the center and the boundaries to avoid near-source effects and artifacts from the absorbing boundaries. Furthermore, we compute the amplitude of the complex-valued output and add it to the mesh object.
path = pathlib.Path("out", "fft.h5")
m = UnstructuredMesh.from_h5(path)

re = m.elemental_fields[f"FFT-Re_0"]
im = m.elemental_fields[f"FFT-Im_0"]
m.attach_field(f"p_amp", np.sqrt(re * re + im * im))
re = m.elemental_fields[f"FFT-Re_1"]
im = m.elemental_fields[f"FFT-Im_1"]
m.attach_field(f"p_amp_1", np.sqrt(re * re + im * im))
m = m.apply_element_mask(m.get_element_centroid_radius() > 0.005)
m = m.apply_element_mask(m.get_element_centroid_radius() < 0.02)

m
There should be something here! The notebook has been saved without widget output.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fe480e6fa50>
To validate the FFT output, we check the radial amplitude decay of the complex field.
# normalize fields
m.elemental_fields["p_amp_1"] /= np.max(m.elemental_fields["p_amp"])
m.elemental_fields["p_amp"] /= np.max(m.elemental_fields["p_amp"])
m.attach_field(
    "distance", 1 / np.sqrt(np.linalg.norm(m.points[m.connectivity], axis=2))
)
m.elemental_fields["distance"] /= np.max(m.elemental_fields["distance"])

np.testing.assert_allclose(
    m.elemental_fields["p_amp"], m.elemental_fields["distance"], rtol=1e-3
)
Furthermore, because the injected sine function produces a single-frequency sound wave with ω=500\omega = 500 kHz, the discrete Fourier transform evaluated at twice this frequency should be (close to) zero everywhere.
assert (m.elemental_fields["p_amp_1"] < 1e-3).all()
PAGE CONTENTS