k-Wave 工具箱学习-2定义计算网格

3.2 定义计算网格

第一个输入 kgrid 定义计算网格的属性。这决定了如何将连续介质划分为均匀分布的网格点。一维、二维、三维网格定义代码示例如下：

kgrid = kWaveGrid(Nx, dx);
kgrid = kWaveGrid(Nx, dx, Ny, dy);
kgrid = kWaveGrid(Nx, dx, Ny, dy, Nz, dz);
x_size = Nx*dx; % x_size为实际x方向长度，单位为米

其中，Nx为网格点数，dx为网格大小（建议dx = dy = dz）。matlab没有单位，k-Wave单位大多为国际单位。

生成后的kgrid结构体包含函数使用的属性，如下所示（属性以k开头的基本不用理解）

kgrid.k	标量波数的格子 ND 网格（不太用理解）
kgrid.k_max	网格支持的最大空间频率
kgrid.t_array	等间隔的时间值数组（即时间步，时域）
kgrid.Nt	时间步数
kgrid.dt	时间步长
kgrid.dim	空间维度数(1、2或3)
kgrid.total_grid_points	网格点的总数
kgrid.Nx	x方向网格点数
kgrid.dx	x方向网格点间距[m]
kgrid.x	x坐标以0 为中心的格子ND网格 [m]
kgrid.x_vec	x坐标的一维向量 [m]
kgrid.x_size	x方向上的实际总长度 [m]
kgrid.kx	x方向波数的平铺ND网格
kgrid.kx_vec	x方向波数的一维矢量
kgrid.kx_max	x方向的最大空间频率

网格数量和网格大小的确定

创建新的模拟时，需确定Nx和dx。

dx越小越精确，计算量越大，求dx的公式如下

% 根据所需的x_size和f_max计算dx和Nx

points_per_wavelength = 3;

dx = c0_min/(points_per_wavelength*f_max);

Nx = round(x_size/dx);

这里，c0_min是介质中的最小声速（声速越小，dx越小，越精确），f_max是最大频率（频率越大，dx越小，越精确），points_per_wavelength每个超声波长的期望网格点数。对于均匀介质中的线性模拟，可以使用接近奈奎斯特限制的每波长两点(建议每波长至少三点以使PML有效工作)来运行模拟。但是，对于非均匀介质，则建议每个波长使用四个或更多点数。对于非线性仿真，最大频率应设置为具有大量能量的最高谐波的频率。

k-Wave 中使用的空间梯度计算大量使用快速傅立叶变换 (FFT)。根据模拟的复杂性，每个时间步最多计算 14 个 FFT。通过将每个方向（包括 PML）上的网格点总数选择为 2 的幂或具有较小的质因数的数，可以最大限度地减少计算每个 FFT 的时间。应避免使用大素数（例如 149）的网格总数。在许多情况下，可以通过稍微修改 kgrid.Nx、kgrid.Ny 等的值以使最大素因子变小来提高 k-Wave 的性能。可以使用函数 checkFactors 获得为任何给定范围选择的适当值。这将返回指定范围内最大质因数为 7 或更少的数字。

checkFactors(100, 150) % 寻找 100 到 150 之间合适网格数量的示例

% 输出结果展示
Numbers with a maximum prime factor of 2
128
Numbers with a maximum prime factor of 3
108 144
Numbers with a maximum prime factor of 5
100 120 125 135 150
Numbers with a maximum prime factor of 7
105 112 126 140 147


Nx = 128； %选择2的幂

完美匹配层（PML）

当声波到达计算域的边缘时，它们被一种特殊类型的各向异性吸收边界层吸收，称为完美匹配层或 PML。可以通过运行工具箱中包含的示例观察传播波接近计算域边缘时会发生什么，可以看到该层的效果。默认情况下，PML 占据域边缘周围的 20 个网格点（3D 中为 10 个网格点）。值得注意的是，PML 放置在使用 makeGrid 指定的网格大小内。

这意味着用户必须小心，不要将源或传感器点放置在该层内。

或者，通过将可选输入参数“PMLInside”设置为 false，可以将 PML 设置为用户定义的网格大小之外（有关如何使用可选输入参数的概述，请参阅第 3.6 节）。在这种情况下，网格大小和介质输入会在模拟开始之前自动放大。（我也不填明白，没用过）

创建时间步的方法。包括时间步数量和大小

创建 kgrid 对象后，用户可以修改的唯一参数是 kgrid.t_array。这描述了运行模拟的时间值数组，默认设置为“自动”。在这种情况下，时间数组是使用 makeTime 在模拟函数中自动计算的。此函数将总时间设置为声波以最小声速穿过最长网格对角线所需的时间。时间步长基于 0.3 的 Courant-Friedrichs-Lewy (CFL) 数和介质中的最大声速，其中 kgrid.dt = CFL*dx_min/c0_max（请参阅第 2.7 节中的讨论）。还可以通过调用 makeTime 或显式设置时间数组来手动设置时间数组。下面给出几个例子。时间数组必须均匀分布且单调递增。创建后，可以使用kgrid.Nt和kgrid.dt查询时间点的数量和时间步的大小。

% 使用 makeTime 创建时间数组，设置 CFL 和结束时间

kgrid.t_array = makeTime(kgrid, medium.sound_speed, CFL, t_end);

% 使用 makeTime 创建时间数组并设置结束时间

kgrid.t_array = makeTime(kgrid, medium.sound_speed, [], t_end);

% 显式创建时间数组

kgrid.t_array = 0:1e-9:1e-6;

 References

[1] B. E. Treeby and B. T. Cox, "k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave-fields," J. Biomed. Opt., vol. 15, no. 2, p. 021314, 2010. download pdf

[2] B. E. Treeby, J. Jaros, A. P. Rendell, and B. T. Cox, "Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a k-space pseudospectral method," J. Acoust. Soc. Am., vol. 131, no. 6, pp. 4324-4336, 2012. download pdf

[3] B. E. Treeby and B. T. Cox, "Modeling power law absorption and dispersion for acoustic propagation using the fractional Laplacian," J. Acoust. Soc. Am., vol. 127, no. 5, pp. 2741-2748, 2010. download pdf (see also erratum)