Acoustic Helmholtz Solver

MATLAB implementation of an iterative solver for the acoustic Helmholtz equation with heterogeneous sound speed, density, and absorption. The solver uses a convergent Born-like method based on the Universal Split Preconditioner formulation.

Features

• Solves the first-order linear acoustic equations for continuous wave sources
• Supports heterogeneous sound speed, density, and absorption
• 2D and 3D simulations
• Perfectly Matched Layer (PML) boundary conditions
• Optional GPU acceleration

Installation

1. Clone this repository:

   git clone https://github.com/ucl-bug/acoustic-helmholtz-solver.git

2. Add the src folder to your MATLAB path:

   addpath('path/to/acoustic-helmholtz-solver/src')

Dependencies

This solver requires the following packages to be installed and on the MATLAB path:

• k-Wave Toolbox - Used for helper functions (makeDisc, makeGrid, etc.) and validation
• Anderson Sphere Scattering - For analytical sphere scattering solutions (only needed for validation examples)

Quick Start

% Define grid
N = 256;
dx = 0.5e-3;  % grid spacing [m]
f0 = 500e3;   % frequency [Hz]

% Define sources (pressure and velocity)
source_p = zeros(N, N);
source_p(21, N/4:3*N/4) = 1.0;  % line source
source_u = zeros(N, N, 1, 2);   % no velocity source

% Define medium properties
c0 = 1500 * ones(N, N);         % sound speed [m/s]
rho0 = 1000 * ones(N, N);       % density [kg/m^3]
alpha = zeros(N, N);            % absorption [Np/m]

% Run solver
[pressure, velocity] = acousticHelmholtzSolver(...
   dx, source_p, source_u, f0, c0, rho0, alpha);

% Plot result
imagesc(abs(pressure));
axis image;
colorbar;

Required Inputs

Parameter	Description
dx	Grid spacing in all directions [m]
source_p	Complex pressure source matrix [Pa]
source_u	Complex velocity source (4D array, components in dim 4) [m/s]
f0	Source frequency [Hz]
c0	Sound speed (scalar or matrix) [m/s]
rho0	Mass density (scalar or matrix) [kg/m^3]
alpha	Absorption coefficient (scalar or matrix) [Np/m]

Optional Parameters

Parameter	Default	Description
number_iterations	200	Maximum iterations
atol	1e-3	Absolute tolerance
rtol	1e-3	Relative tolerance
pml_size	50	PML thickness in grid points
pml_alpha	0.6	PML absorption coefficient
pml_inside	false	Include PML inside defined grid
use_gpu	false	Enable GPU acceleration
plot_sim	false	Plot field during iteration
verbose	true	Print status messages

Examples

See the examples/ directory for demonstration scripts:

• heterogeneous_all_2D.m - 2D simulation with heterogeneous density, sound speed, and absorption
• sphere_scattering.m - Scattering from a sphere
• convergence.m - Convergence analysis

License

This project is licensed under the GNU Lesser General Public License v3.0 - see the LICENSE file for details.

Acknowledgments

• Some utility functions are adapted from the k-Wave Toolbox
• The split preconditioner formulation is based on Vettenburg & Vellekoop (2022)