Hi Brad, Ben, and the k-Wave community,
As I mentioned at Photonics West, I am applying k-Wave to model scattering of planar shear waves in heterogeneous tissue, and so far the toolbox has been very easy to use and is working quite well for elastic scattering problems. I recently started trying to model shear wave absorption using the medium.alpha_coeff_shear factor to account for viscoelasticity. The results look quite reasonable for certain cases, but when I increase the absorption factor enough, the result seems to become unstable. I am using a value on the order of 3e6 dB/(MHz^2 cm), since this is around what we measure for shear wave absorption in tissue-mimicking phantoms with wave frequency on the order of 600 Hz, i.e. 1.08 dB/cm. You can see the instability forming in the example code below, where the edges of the plane wave start to increase with an amplitude much larger than the initial wave amplitude.
I saw in the Modelling Power Law Absorption Example that you discuss numerical errors (http://www.k-wave.org/documentation/example_na_modelling_absorption.php#heading3), but there it is mentioned that the errors can be minimized by reducing the time step. I found the instability I see seems to worsen with a decreased time step (which I implemented by reducing the CFL number). The instability is also worse when I reduce the grid spacing (which in turn also reduced the time step).
See the comments in my code below for the factors that influence the instability that I investigated. All the parameters used in the example code are typical for our experiments, except maybe the inital velocity u_0, but there doesn't seem to be any dependence of the instability on this value. This leads me to the following questions:
1. Is this instability related to the numerical errors discussed in the Modelling Power Law Absorption Example?
2. Is there some relationship that tells me for which conditions (up to which absorption) the code works for before becoming unstable?
3. If I have a result that appears stable (i.e. with a low enough medium.alpha_coeff_shear value), how can I check whether the results are reliable, since when I decrease the grid size to check convergence the simulation then blows up?
I would be grateful if someone had any insight to any of these questions.
Thanks,
Genny
% -------------------------------------------------------------------------------------
% Modified from Shear Waves And Critical Angle Reflection Example
% Simplified code that shows instability during absorption of planar shear wave
% Heterogeneities removed to focus on instability
% Grid size cropped for quicker simulation
% The parameters coded below result in an amplitude that appears to blow up
% May 2019
%
% If one of the factors is changed to as decribed below (with all other parameters constant):
% * decreasing alpha0_s1 to 2e6 the solution appears stable
% * increasing alpha0_s1 to 4e6 the solution blows up more
% * increasing t_end to 5000 the solution blows up more
% * decreasing t_end to 2000 the solution appears stable
% * decreasing cfl to 0.1 (decreases dt by a factor of 2) the solution blows up more
% * increasing scale to 1.25 (decreases dx & dt by a factor of 0.8) the solution blows up much more
% * decreasing scale to 0.75 (increases dx & dt by a factor of 1.3) the solution appears stable
% * increasing cp1 to 3000 (decreases dt by a factor of 2) the solution blows up more
% * decreasing cp1 to 750 (increases dt by a factor of 2) the solution appears stable
% * increasing cs1 to 4.5 the solution blows up more
% * decreasing cs1 to 3.5 the solution appears stable
% * instability appears insensitive to frequency (300-1800 Hz tested)
% * instability appears insensitive to u_0 (0.5 to 1000 tested)
close all
clearvars
% define the medium properties for the outise material
cp1 = 1500; % compressional wave speed [m/s]
cs1 = 4.35; % shear wave speed [m/s]
rho1 = 1000; % density [kg/m^3]
alpha0_p1 = 0; % compressional absorption [dB/(MHz^2 cm)]
alpha0_s1 = 3e6; % shear absorption [dB/(MHz^2 cm)] (measured approx 1.08 dB/cm at 600 Hz)
% create the computational grid
scale = 1; % grid spacing factor
PML_size = 10; % [grid points]
Nx = 40*scale ; % [grid points]
Ny = 20*scale ; % [grid points]
dx = 0.4e-3/scale; % [m]
dy = 0.4e-3/scale; % [m]
kgrid = kWaveGrid(Nx, dx, Ny, dy);
% create the time array to run until wave passes through entire grid
cfl = 0.2;
t_end = 4500*1e-6;
kgrid.makeTime(max([cp1 cs1]), cfl, t_end);
% define the source
source_freq = 600; % [Hz]
u_0 = 500; % initial velocity
source_mask = zeros(Nx, Ny);
source_mask(:,1) = 1;
source.u_mask = source_mask;
source.ux = u_0 * sin(2 * pi * source_freq * kgrid.t_array);
source.uy = 0 * sin(2 * pi * source_freq * kgrid.t_array);
% define sensor
sensor.record = {'u_max_all','u_min_all','u_final'};
% define the medium properties
clear medium
medium.sound_speed_compression = cp1*ones(Nx, Ny);
medium.sound_speed_shear = cs1*ones(Nx, Ny);
medium.density = rho1*ones(Nx, Ny);
medium.alpha_coeff_compression = alpha0_p1*ones(Nx, Ny);
medium.alpha_coeff_shear = alpha0_s1*ones(Nx, Ny);
% set the input arguments
input_args = {'PMLSize', PML_size, 'PMLAlpha', 2, 'PlotPML', false, ...
'PMLInside', false, 'PlotScale', 'auto', ...
'DataCast', 'single'};
% run the elastic simulation
sensor_data_elastic = pstdElastic2D(kgrid, medium, source, sensor, input_args{:});
% visualization
ui = u_0*cp1/cs1; % initial velocity
figure(1)
subplot(2,2,1);
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.ux_max_all/ui); caxis([0 1.5])
xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_x');title('ux max all')
subplot(2,2,2);
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.ux_final/ui); caxis([-1.5 1.5])
xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_x');title('ux final')
subplot(2,2,3);
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.uy_max_all/ui); caxis([0 1.5])
xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_y');title('uy max all')
subplot(2,2,4);
imagesc(kgrid.y_vec*1e3, kgrid.x_vec*1e3, sensor_data_elastic.uy_final/ui); caxis([-1.5 1.5])
xlabel('y mm');ylabel('x mm');axis image;h = colorbar;ylabel(h, 'Normalized u_y');title('uy final')