ivp_1D_simulation.py

"""
Simulations In One Dimension Example

This example provides a simple demonstration of using k-Wave for the
simulation and detection of the pressure field generated by an initial
pressure distribution within a one-dimensional heterogeneous propagation
medium. It builds on the Homogeneous Propagation Medium and Heterogeneous
Propagation Medium examples.
"""

import matplotlib.pyplot as plt
import numpy as np

from kwave.data import Vector
from kwave.kgrid import kWaveGrid
from kwave.kmedium import kWaveMedium
from kwave.ksensor import kSensor
from kwave.ksource import kSource
from kwave.kspaceFirstOrder1D import kspace_first_order_1D
from kwave.options.simulation_execution_options import SimulationExecutionOptions
from kwave.options.simulation_options import SimulationOptions

# =========================================================================
# SIMULATION
# =========================================================================

# create the computational grid
Nx: int = 512         # number of grid points in the x (row) direction
dx: float = 0.05e-3   # grid point spacing in the x direction [m]

grid_size_points = Vector([Nx, ])
grid_spacing_meters = Vector([dx, ])

# create the k-space grid
kgrid = kWaveGrid(grid_size_points, grid_spacing_meters)

# define the properties of the propagation medium
sound_speed = 1500.0 * np.ones((Nx, 1))                    # [m/s]
sound_speed[:np.round(Nx / 3).astype(int) - 1] = 2000.0	   # [m/s]
density = 1000.0 * np.ones((Nx, 1))                        # [kg/m^3]
density[np.round(4 * Nx / 5).astype(int) - 1:] = 1500.0    # [kg/m^3]
medium = kWaveMedium(sound_speed=sound_speed, density=density)

# Create the source object
source = kSource()

# create initial pressure distribution using a smoothly shaped sinusoid
x_pos: int = 280    # [grid points]
width: int = 100    # [grid points]
height: int = 1     # [au]
p0 = np.linspace(0.0, 2.0 * np.pi, width + 1)

part1 = np.zeros(x_pos).astype(float)
part2 = (height / 2.0) * np.sin(p0 - np.pi / 2.0) + (height / 2.0)
part3 = np.zeros(Nx - x_pos - width - 1).astype(float)
source.p0 = np.concatenate([part1, part2, part3])

# create a Cartesian sensor mask recording the pressure
sensor = kSensor()
sensor.record = ["p"]

# this hack is needed to ensure that the sensor is in [1,2] dimensions
mask = np.array([-10e-3, 10e-3]) # [m]
mask = mask[:, np.newaxis].T
sensor.mask = mask 

# set the simulation time to capture the reflections
c_max = np.max(medium.sound_speed.flatten()) # [m/s]
t_end = 2.5 * kgrid.x_size / c_max           # [s]

# define the time array
kgrid.makeTime(c_max, t_end=t_end)

# define the simulation options
simulation_options = SimulationOptions(data_cast="off", save_to_disk=False)
execution_options = SimulationExecutionOptions(is_gpu_simulation=True)

# run the simulation
sensor_data = kspace_first_order_1D(kgrid, source, sensor, medium, simulation_options=simulation_options, execution_options=execution_options)

# =========================================================================
# VISUALISATION
# =========================================================================

# plot the recorded time signals
_, ax1 = plt.subplots()
ax1.plot(sensor_data['p'][0, :], 'b-')
ax1.plot(sensor_data['p'][1, :], 'r-')
ax1.grid(True)
ax1.set_ylim(-0.1, 0.7)
ax1.set_ylabel('Pressure')
ax1.set_xlabel('Time Step')
ax1.legend(['Sensor Position 1', 'Sensor Position 2'])
plt.show()