I am using k-Wave to simulate wave propagation in strongly contrasting media. A strange effect shows up during propagation when there is strong absorption. The attached figure illustrates the effect: Shortly after starting the excitation, and before the pulse arrives, a wave is generated inside the neighbouring layer.
The magnitude of this effect is apparently influenced by several factors:
* The level of absorption — more absorption _increases_ the effect, contrary to my intuition.
* The distance from the source to the medium interface.
* The number of layers in the total medium (four layers is worse than three).
* When the wave reaches the PML it gets even worse.
* Using a power law exponent of 2 seems to remove the problem (this was pointed out to me by Dr. Treeby.)
I can reduce the effect by increasing the temporal and spatial sampling rates, but some times even using 40 points per wavelength and a timestep less than a 1ns is not good enough for a transmit frequency of 500kHz.
I have enclosed an example script that illustrates the effect. It seems that the effect gets worse at some ’sweet spot’ of combined impedances and layer thicknesses, which I cannot quite re-create.
My questions:
* Do you see why this is happening? Could there be a bug in the absorption implementation, or is this more or less expected behaviour?
* What can I do to reduce this effect? For certain geometries, cranking up the sampling rates leads to 20h+ long 1D simulation.
Here's an example script:
%%
% This example shows a strange effect where a spurious wave is generated inside a
% medium layer before the excitation wave arrives in that layer.
%
% This happens when there is strong absorption inside one of the layers (e.g. 10dB/cm/MHz), and the
% effect is smaller when the absorption is reduced (e.g. to 1dB/cm/MHz). Tweaking the
% geometry in various ways also seems to change the effect, having four layers seems
% to be worse than having three layers, and the thickness of the various layers also
% make a difference.
%
% The effect can be reduced by reducing the timestep, but sometimes to extreme levels.
%% Geometry
padding = 5e-3;
Lx = [6,2,10,20]*1e-3;
%% Material properties
c = [600, 3000, 6000, 3000];
rho = [200, 2000, 8000, 2000];
attn = [0.2, 10, 0.05, 0.1];
% attn = [0.2, 1.0, 0.05, 0.1]; % reducing the absorption reduces the effect
attn_power = 1.0; % changing it to 2.0 removes the problem
useattn = true;
%% Transmit pulse
Npulse = 2;
fc = 0.5e6;
%% Computational grid
ppw = 20; % points per wavelength
cfl = 0.1; % Courant-Friedrich-Lewy
% grid size
dx = min(c)/(ppw*fc);
Nx = ceil(sum([Lx,padding])/dx);
kgrid = makeGrid(Nx,dx);
% index positions of interfaces
ind = round(cumsum([padding,Lx(1:end-1)])/dx);
medium = struct();
% sound velocity
medium.sound_speed = c(1) + zeros(Nx,1);
for jj = 2:length(ind),
medium.sound_speed(ind(jj):end) = c(jj);
end
% density
medium.density = rho(1) + zeros(Nx,1);
for jj = 2:length(ind),
medium.density(ind(jj):end) = rho(2);
end
% attenuation
if useattn,
medium.alpha_coeff = attn(1) + zeros(Nx,1);
for jj = 2:length(ind),
medium.alpha_coeff(ind(jj):end) = attn(jj);
end
medium.alpha_power = attn_power;
medium.alpha_mode = 'no_dispersion';
end
[kgrid.t_array,dt] = makeTime(kgrid,medium.sound_speed,cfl);
%% Transducer
Amp = 1e-6;
source = struct();
source.u_mask = zeros(Nx,1);
source.u_mask(ind(1)) = 1;
source.ux = -Amp*toneBurst(1/kgrid.dt,fc,Npulse,'Envelope','Gaussian');
sensor = struct();
sensorind = ind(1):10:Nx;
sensor.mask = zeros(Nx,1);
sensor.mask(sensorind) = 1;
sensor.record = {'p'};
% display
display_mask = zeros(Nx,1);
display_mask(ind) = 1;
%% 1D sim
args = {'DataCast','double', ...
'LogScale',true, ...
'DisplayMask',display_mask, ...
'PMLSize',ppw*2, ... % 2 wavelengths
'PMLInside',false, ...
'PlotPML',false, ...
'PlotSim',true};
data = kspaceFirstOrder1D(kgrid,medium,source,sensor,args{:});