examples/ivp_photoacoustic_waveforms/ivp_photoacoustic_waveforms.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "11cd3967",
   "metadata": {},
   "source": [
    "# Initial Value Problem: Photoacoustic Waveforms\n",
    "\n",
    "First load the modules needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c11efc20",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "!pip install k-wave-python"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35bbf8ce",
   "metadata": {},
   "source": [
    "Now import the libraries we need"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fe4153f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from kwave.data import Vector\n",
    "from kwave.kgrid import kWaveGrid\n",
    "from kwave.kmedium import kWaveMedium\n",
    "from kwave.ksource import kSource\n",
    "from kwave.ksensor import kSensor\n",
    "from kwave.utils.data import scale_SI\n",
    "from kwave.utils.mapgen import make_disc, make_ball\n",
    "from kwave.kspaceFirstOrder2D import kspaceFirstOrder2D\n",
    "from kwave.kspaceFirstOrder3D import kspaceFirstOrder3D\n",
    "from kwave.options.simulation_execution_options import SimulationExecutionOptions\n",
    "from kwave.options.simulation_options import SimulationOptions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96d23a4d",
   "metadata": {},
   "source": [
    "Next define the problem in terms of the domain dimension, size, source position, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7e5f5b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "Nx = 64      # number of grid points in the x (row) direction\n",
    "\n",
    "x = 1e-3     # size of the domain in the x direction [m]\n",
    "\n",
    "dx = x / Nx  # grid point spacing in the x direction [m]\n",
    "\n",
    "sound_speed = 1500\n",
    "\n",
    "# size of the initial pressure distribution\n",
    "source_radius = 2              # [grid points]\n",
    "\n",
    "# distance between the centre of the source and the sensor\n",
    "source_sensor_distance = 10    # [grid points]\n",
    "\n",
    "# time array\n",
    "dt = 2e-9                      # [s]\n",
    "t_end = 300e-9                 # [s]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87aa5a24",
   "metadata": {},
   "source": [
    "Code for 2d simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "964309b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "medium2 = kWaveMedium(sound_speed=sound_speed)\n",
    "\n",
    "# create the k-space grid\n",
    "kgrid2 = kWaveGrid([Nx, Nx], [dx, dx])\n",
    "\n",
    "# create the time array using an integer number of points per period\n",
    "kgrid2.setTime(int(np.round(t_end / dt)), dt)\n",
    "\n",
    "# create instance of a sensor\n",
    "sensor2 = kSensor()\n",
    "\n",
    "# set sensor mask: the mask says at which points data should be recorded\n",
    "sensor2.mask = np.zeros((Nx, Nx), dtype=bool)\n",
    "\n",
    "# define a single sensor point\n",
    "sensor2.mask[Nx // 2 + source_sensor_distance, Nx // 2] = True\n",
    "\n",
    "# set the record type: record the pressure waveform\n",
    "sensor2.record = ['p']\n",
    "\n",
    "# make a source object\n",
    "source2 = kSource()\n",
    "source2.p0 = make_disc(Vector([Nx, Nx]), Vector([Nx // 2, Nx // 2]), source_radius, plot_disc=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a53f8a88",
   "metadata": {},
   "source": [
    "simulation and execution options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37de7ce4",
   "metadata": {},
   "outputs": [],
   "source": [
    "simulation_options = SimulationOptions(\n",
    "    data_cast='single',\n",
    "    save_to_disk=True)\n",
    "\n",
    "execution_options = SimulationExecutionOptions(\n",
    "    is_gpu_simulation=True,\n",
    "    delete_data=False,\n",
    "    verbose_level=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "315e699c",
   "metadata": {},
   "source": [
    "Now run the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4e9d359",
   "metadata": {},
   "outputs": [],
   "source": [
    "# run the simulation\n",
    "sensor_data_2D = kspaceFirstOrder2D(\n",
    "    medium=medium2,\n",
    "    kgrid=kgrid2,\n",
    "    source=source2,\n",
    "    sensor=sensor2,\n",
    "    simulation_options=simulation_options,\n",
    "    execution_options=execution_options)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f25433e0",
   "metadata": {},
   "source": [
    "# 3D Simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76905e44",
   "metadata": {},
   "outputs": [],
   "source": [
    "medium3 = kWaveMedium(sound_speed=1500)\n",
    "\n",
    "# create the k-space grid\n",
    "kgrid3 = kWaveGrid([Nx, Nx, Nx], [dx, dx, dx])\n",
    "\n",
    "# create the time array using an integer number of points per period\n",
    "kgrid3.setTime(int(np.round(t_end / dt)), dt)\n",
    "\n",
    "# create instance of a sensor\n",
    "sensor3 = kSensor()\n",
    "\n",
    "# set sensor mask: the mask says at which points data should be recorded\n",
    "sensor3.mask = np.zeros((Nx, Nx, Nx), dtype=bool)\n",
    "\n",
    "# define a single sensor point\n",
    "sensor3.mask[Nx // 2 + source_sensor_distance, Nx // 2, Nx // 2] = True\n",
    "\n",
    "# set the record type: record the pressure waveform\n",
    "sensor3.record = ['p']\n",
    "\n",
    "# make a source object\n",
    "source3 = kSource()\n",
    "source3.p0 = make_ball(Vector([Nx, Nx, Nx]), Vector([Nx // 2, Nx // 2, Nx // 2]), source_radius)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfce9e85",
   "metadata": {},
   "outputs": [],
   "source": [
    "sensor_data_3D = kspaceFirstOrder3D(\n",
    "    medium=medium3,\n",
    "    kgrid=kgrid3,\n",
    "    source=source3,\n",
    "    sensor=sensor3,\n",
    "    simulation_options=simulation_options,\n",
    "    execution_options=execution_options)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1626ec9b",
   "metadata": {},
   "source": [
    "Plot the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "507443fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "t_sc, t_scale, t_prefix, _ = scale_SI(t_end)\n",
    "_, ax1 = plt.subplots()\n",
    "ax1.plot(np.squeeze(kgrid2.t_array * t_scale), sensor_data_2D['p'] / np.max(np.abs(sensor_data_2D['p'])), 'r-')\n",
    "ax1.plot(np.squeeze(kgrid3.t_array * t_scale), sensor_data_3D[\"p\"] / np.max(np.abs(sensor_data_3D[\"p\"])), \"k-\", label=\"3D\")\n",
    "ax1.set(xlabel= f\"Time [{t_prefix}s]\", ylabel='Recorded Pressure [au]')\n",
    "ax1.grid(True)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "gpuType": "T4",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}