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))