Example 3
ΒΆ
The example shows how to prepare k-Wave input file describing simulation in linear homogeneous medium. The pressure source is specified in the form of two sinusoidal sources (planes) perpendicular to X-axis. Here each source point uses its own separate driving signal. While all signals are of \(1\) MHz with amplitude of \(10\) Pa, the points belonging to second plane source use inverted series. Therefore plane sources should be always in the opposite phase.
The simulation has spatial resolution of \(10^{-3}\) m and the simulation domain of \(256 \times 128 \times 64\) grid points. While temporal resolution is computed so that CFL condition number is 0.3.
Finally, the sensor mask is set up to allow collection of pressure evolution in time at a plane perpendicular to Z axis placed in the middle of the domain (in Z-direction).
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 | import sys
import numpy as np
import math
sys.path.append('..')
from kwave_input_file import KWaveInputFile
from kwave_data_filters import SpectralDataFilter
filename = 'example_003_p_source_many_py.h5'
grid_res = (1.0e-4, 1.0e-4, 1.0e-4)
grid_size = (256, 128, 64)
dt = 0.3 * grid_res[0] / 1500.0
end_time = ((grid_size[0] / 2) * grid_res[0]) / 1500.0
steps = int(end_time / dt)
with KWaveInputFile(filename, grid_size, steps, grid_res, dt) as file:
file.write_medium_sound_speed(1500.0)
file.write_medium_density(1000.0)
# Define source shape as two planes in the domain
source_mask = np.zeros(grid_size)
source_mask[ int(grid_size[0] / 4),
int(grid_size[1] / 2 - 16): int(grid_size[1] / 2 + 16),
int(grid_size[2] / 2 - 16): int(grid_size[2] / 2 + 16)] = 1
source_mask[-int(grid_size[0] / 4),
int(grid_size[1] / 2 - 16): int(grid_size[1] / 2 + 16),
int(grid_size[2] / 2 - 16): int(grid_size[2] / 2 + 16)] = 2
# Define sinusoidal source signal of 1 MHz and amplitude of 10 Pa
amplitude = 10
frequency = 1.0e6
source_signal = amplitude * np.sin(frequency * (2*math.pi) * np.arange(0.0, steps*dt, dt))
# Retrieve list of non-zero values in the source mask in the order of indexes that will be used to define the source
source_mask_non_zero = file.domain_mask_values_in_index_order(source_mask)
# Allocate space for source series for each source grid-point
source_signals = np.zeros((source_mask_non_zero.size, source_signal.size))
# Using non-zero value positions we can identify each plane and assign appropriate signal series to its grid-points
source_signals[source_mask_non_zero == 1, :] = source_signal
source_signals[source_mask_non_zero == 2, :] = - source_signal
# Finally write the source signals and source indexes to the input file
file.write_source_input_p(file.domain_mask_to_index(source_mask), source_signals, KWaveInputFile.SourceMode.ADDITIVE, 1500.0)
sensor_mask = np.zeros(grid_size)
sensor_mask[:, :, int(grid_size[2] / 2)] = 1
file.write_sensor_mask_index(file.domain_mask_to_index(sensor_mask))
|