Example 2ΒΆ

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. Both sources use single driving signal of \(1\) MHz with amplitude of \(10\) Pa.

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
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_002_p_source_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)] = 1

    # 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))
    file.write_source_input_p(file.domain_mask_to_index(source_mask), source_signal,
                              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))