Example 4
¶
The example shows how to use k-Wave Binary Driver to call k-Wave solver directly from Python code. First input file is setup in a similar way to other examples and k-Wave Output File is prepared as it has to be provided to the solver driver.
The solver driver is setup and later used with input file and output file to run the solver. Note that it is necessary to correctly which data should be stored by the solver (“store_*” calls) and then read only the appropriate data from the output file.
Finally output file is open, final pressure and pressure evolution at the sensor is read and plotted together with initial pressure using matplotlib. The figure is stored as an image:

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 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 | import sys
import numpy as np
import matplotlib.pyplot as plt
sys.path.append('..')
from kwave_input_file import KWaveInputFile
from kwave_output_file import KWaveOutputFile, DomainSamplingType, SensorSamplingType
from kwave_data_filters import SpectralDataFilter
from kwave_bin_driver import KWaveBinaryDriver
filename = 'example_004_bin_driver_in.h5'
grid_res = (1.0e-4, 1.0e-4, 1.0e-4)
grid_size = (64, 128, 256)
dt = 0.3 * grid_res[0] / 1500.0
end_time = 2 * ((grid_size[0] / 2) * grid_res[0]) / 1500.0
steps = int(end_time / dt)
# Define simulation input and output files
input_file = KWaveInputFile(filename, grid_size, steps, grid_res, dt)
output_file = KWaveOutputFile('example_004_bin_driver_out.h5')
# Prepare initial pressure distribution
p0 = np.zeros(grid_size)
p0[int(grid_size[0] / 2 - 16): int(grid_size[0] / 2 + 16),
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.0e3
p0 = SpectralDataFilter.smooth(p0)
# Open the simulation input file and fill it as usual
with input_file as file:
file.write_medium_sound_speed(1500.0)
file.write_medium_density(1000.0)
file.write_source_input_p0(p0)
sensor_mask = np.zeros(grid_size)
sensor_mask[:, :, int(grid_size[2] / 4)] = 1
file.write_sensor_mask_index(file.domain_mask_to_index(sensor_mask))
# Create k-Wave solver driver, which will call C++/CUDA k-Wave binary.
# It is usually necessary to specify path to the binary: "binary_path=..."
driver = KWaveBinaryDriver()
# Specify which data should be sampled during the simulation (final pressure in the domain and
# RAW pressure at the sensor mask
driver.store_pressure_everywhere([DomainSamplingType.FINAL])
driver.store_pressure_at_sensor([SensorSamplingType.RAW])
# Execute the solver with specified input and output files
driver.run(input_file, output_file)
# Open the output file and generate plots from the results
with output_file as file:
fig = plt.figure()
gs = fig.add_gridspec(2, 2)
ax0 = fig.add_subplot(gs[0])
ax1 = fig.add_subplot(gs[1])
ax2 = fig.add_subplot(gs[2:4])
# Read final pressure over whole domain and plot it as a heat map together with initial pressure
p_final = file.read_pressure_everywhere(DomainSamplingType.FINAL)
ax0.imshow(np.squeeze(p0[32, :, :]), cmap='hot', interpolation='nearest')
ax1.imshow(np.squeeze(p_final[32, :, :]), cmap='hot', interpolation='nearest')
# Read pressure series at each of sensor mask points
p = file.read_pressure_at_sensor(SensorSamplingType.RAW)
# Take absolute value and maximum of the pressure at the sensor and plot the result over time.
ax2.plot(np.arange(0, output_file.read_temporal_properties()[0]) * dt, np.squeeze(np.max(np.abs(p), 2)))
ax2.set(xlabel='time [s]', ylabel='Acoustic pressure [Pa]', title='Max. abs. pressure on sensor mask')
ax2.grid()
fig.tight_layout(pad=3.0)
# Save figure to an image
fig.savefig('example_004_bin_driver_out.png')
|