2D acoustic problem

Overview

Solving the acoustic wave equation is a core technology in fields such as medical ultrasound and geological exploration. Large-scale acoustic wave equation solvers face challenges in terms of computational power and storage. Solvers for the wave equation generally use either frequency domain algorithms or time domain algorithms. The representative time domain algorithm is the Time Domain Finite Difference (TDFD) method, while frequency domain algorithms include Frequency Domain Finite Difference (FDFD), Finite Element Method (FEM), and Convergent Born Series (CBS) iterative method. The CBS method, due to its low memory requirement and absence of dispersion error, has gained widespread attention in engineering and academia. In particular, Osnabrugge et al. (2016) have addressed the convergence issue of this method, expanding the application prospects of the CBS method.

This case study will demonstrate how to invoke the CBS API provided by MindFlow to solve the two-dimensional acoustic wave equation.

Problem Description

In solving the acoustic wave equation, the input parameters are the velocity field and source information, and the output is the spatiotemporal distribution of the wave field.

The expression of the two-dimensional acoustic wave equation is as follows

Time Domain Expression	Frequency Domain Expression
$\frac{\partial^{2} u}{\partial t^{2}} - c^{2} Δ u = f$	$- ω^{2} \hat{u} - c^{2} Δ \hat{u} = \hat{f}$

Where

$u (x, t) [L]$ is the deformation displacement (pressure divided by density), a scalar.
$c (x) [L / T]$ is the wave velocity, a scalar.
$f (x, t) [L / T^{2}]$ is the excitation source (volume distributed force), a scalar.

In practical solving, in order to reduce the parameter dimension, the parameters are generally made dimensionless first, and then the dimensionless equations and parameters are solved, and finally the dimensional solutions are restored. By selecting $ω$ , $\hat{f}$ , and $d$ (grid spacing, which requires equal spacing in all directions) to nondimensionalize the frequency domain equation, we can obtain the dimensionless frequency domain equation:

u^{*} + c^{* 2} \tilde{Δ} + f^{*} = 0

Where

$u^{*} = \hat{u} ω^{2} / \hat{f}$ is the dimensionless deformation displacement.
$c^{*} = c / (ω d)$ is the dimensionless wave velocity.
$\tilde{Δ}$ is the normalized Laplace operator, which is the Laplace operator when the grid spacing is 1.
$f^{*}$ the mask that marks the source position, with a value of 1 at the source and 0 at other positions.

The src package in this case can be downloaded at src.

[1]:

import os
import numpy as np
import scipy
import pandas as pd
import mindspore as ms
from mindspore import Tensor

from mindflow.utils import load_yaml_config

from cbs.cbs import CBS
from src import visual
from solve_acoustic import solve_cbs

Define input parameters and output sampling method

The required inputs for this case are dimensional 2D velocity field, source location list, and source waveform. The input file name is specified in the config.yaml file. For user convenience, pre-set inputs are provided here. Please download the data and put them in ./dataset in the case directory. The data include the velocity field velocity.npy, source location list srclocs.csv, and source waveform srcwaves.csv. Users can modify the input parameters based on the input file format.

The output is a spatiotemporal distribution of the wavefield. To specify how the output is sampled in time and frequency, parameters such as dt and nt need to be specified in the config.yaml file.

Since the sampling rate of the input source waveform in time may differ from the required output, interpolation needs to be performed.

[2]:

ms.set_context(device_target='Ascend', device_id=0, mode=ms.GRAPH_MODE)

config = load_yaml_config('config.yaml')

data_config = config['data']
solve_config = config['solve']
summary_config = config['summary']

# read time & frequency points
dt = solve_config['dt']
nt = solve_config['nt']
ts = np.arange(nt) * dt
omegas_all = np.fft.rfftfreq(nt) * (2 * np.pi / dt)

# read source locations
df = pd.read_csv(os.path.join(data_config['root_dir'], data_config['source_locations']), index_col=0)
slocs = df[['y', 'x']].values # shape (ns, 2)

# read & interp source wave
df = pd.read_csv(os.path.join(data_config['root_dir'], data_config['source_wave']))
inter_func = scipy.interpolate.interp1d(df.t, df.f, bounds_error=False, fill_value=0)
src_waves = inter_func(ts) # shape (nt)
src_amplitudes = np.fft.rfft(src_waves) # shape (nt//2+1)

# read velocity array
velo = np.load(os.path.join(data_config['root_dir'], data_config['velocity_field']))
nz, nx = velo.shape
dx = data_config['velocity_dx']

Select desired frequency points

With the output sampling method determined, all the desired frequency points are in turn determined. However, in order to reduce computational load, it is also possible to select only a portion of the frequency points for calculation, while obtaining the remaining frequency points through interpolation. The specific frequency point downsampling method is specified by the downsample_mode and downsample_rate in the config.yaml file. The default is no downsampling, which means solving all frequency points except $ω = 0$ .

[3]:

# select omegas
no = len(omegas_all) // solve_config['downsample_rate']

if solve_config['downsample_mode'] == 'exp':
    omegas_sel = np.exp(np.linspace(np.log(omegas_all[1]), np.log(omegas_all[-1]), no))
elif solve_config['downsample_mode'] == 'square':
    omegas_sel = np.linspace(omegas_all[1]**.5, omegas_all[-1]**.5, no)**2
else:
    omegas_sel = np.linspace(omegas_all[1], omegas_all[-1], no)

Perform simulation

Define the relevant arrays as Tensors, call solve_cbs(), and execute the solution on the NPU. Due to memory limitations, the solution process is executed in batches in the frequency domain. The number of batches is specified by the user in config.yaml and does not need to be divisible by the number of frequency points (allowing the size of the last batch to be different from the other batches). After the solution is completed, the frequency domain solution results will be saved to the file u_star.npy.

[ ]:

# send to NPU and perform computation
os.makedirs(summary_config['root_dir'], exist_ok=True)
velo = Tensor(velo, dtype=ms.float32, const_arg=True)
cbs = CBS((nz, nx), remove_pml=False)

ur, ui = solve_cbs(cbs, velo, slocs, omegas_sel, dx=dx, n_batches=solve_config['n_batches']) # shape (ns, no, len(receiver_zs), nx)

u_star = np.squeeze(ur.numpy() + 1j * ui.numpy()) # shape (ns, no, len(krs), nx)
np.save(os.path.join(summary_config['root_dir'], 'u_star.npy'), u_star)

Post-processing

CBS solves the dimensionless frequency domain equation, but downstream tasks often require observing the evolution process of dimensional wavefields in the time domain. Therefore, the final solution is restored to dimensional and converted back to the time domain. The restoration method is given by $\hat{u} = u^{*} h a t f / ω^{2}$ . If downsampling is performed on the frequency points in the “Select desired frequency points” step, interpolation along the frequency direction is required here to restore the solutions for all frequency points. Then, perform a Fourier inverse transform on the dimensional frequency domain wavefield $\hat{u}$ to obtain the time domain wavefield $u$ . Save the time domain wavefield to the file u_time.npy.

[ ]:

# recover dimension and interpolate to full frequency domain
u_star /= omegas_sel.reshape(-1, 1, 1)**2
u_star = scipy.interpolate.interp1d(omegas_sel, u_star, axis=1, kind='cubic', bounds_error=False, fill_value=0)(omegas_all)
u_star *= src_amplitudes.reshape(-1, 1, 1)

# transform to time domain
u_time = np.fft.irfft(u_star, axis=1)
np.save(os.path.join(summary_config['root_dir'], 'u_time.npy'), u_time)

Finally, read the time-domain wave field and visualize.

[ ]:

# visualize the result
u_time = np.load(os.path.join(summary_config['root_dir'], 'u_time.npy'))
visual.anim(velo.numpy(), u_time, ts, os.path.join(summary_config['root_dir'], 'wave.gif'))