MondaicThis 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.%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 apiWe 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_buffermesh = 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 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
<salvus.flow.simple_config.simulation.waveform.Waveform object at 0x7a53adb1f810>
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,
)SalvusJob `job_2411201442041174_f8c1d405d3` running on `local` with 4 rank(s). Site information: * Salvus version: 2024.1.2 * Floating point size: 32
* Downloaded 28.3 MB of results to `out`. * Total run time: 19.12 seconds. * Pure simulation time: 18.66 seconds.
<salvus.flow.executors.salvus_job.SalvusJob at 0x7a53ac774190>
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
<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7a53acd7d690>
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 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