Version:

This documentation is not for the latest stable Salvus 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,
max_frequency=omega,
elements_per_wavelength=6,
tensor_order=4,
).create_mesh()

## Source time function and steady state

We use a single scalar point source in the center of the domain, continuously emitting a sine wave with a frequency of $500$ 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,
)
]

## Simulation

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.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 0x7f7348273c10>
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_2105271706858307_04b915d514 running on local with 4 rank(s).
Site information:
* Salvus version: 0.11.33
* 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: 13.59 seconds.
* Pure simulation time: 13.04 seconds.

<salvus.flow.sites.salvus_job.SalvusJob at 0x7f73430115d0>
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

There should be something here! The notebook has been saved without widget output.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f73482e82d0>

## Validation

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 $\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