Hey guys,
I decided to paste this code in the forum.
I would be very happy if someone could help me out.
This is the code:
close all;
clearvars;
clc;
%% Grid definition
% create the computational grid
dx = 1e-4; % [m/gridpoint]
dy = dx; % [m/gridpoint]
dz = dx; % [m/gridpoint]
Nx = 500; % [gridpoints]
Ny = 60; % [gridpoints]
PML_size = 10; % [grid points]
kgrid2d = kWaveGrid(Nx, dx, Ny, dy);
%% Sensor Definition
sensor2d.record = {'u_max_all','p_max', 'p_final','p','u','p_max_all','p_min_all'};
sensor2d.mask = logical(zeros(Nx,Ny));
sensor2d.mask(1,:) = 1;
date = datestr(datetime('now'));
% ABS
m1.cp = 2000;
m1.rho = 1000;
m1.ac = 2.6;
% Water
m2.cp = 1500;
m2.rho = 1000;
m2.ac = 0.05;
% Silicon -- unused
m3.cp = 1000;
m3.rho = 1200;
m3.ac = 3;
%% Time definition
% create the time array
cfl = 0.1;
t_end = 40e-6;
kgrid2d.makeTime(m1.cp, cfl, t_end);
%% Medium definition
mf2d.sound_speed = m2.cp*ones(Nx, Ny);
mf2d.density = m2.rho*ones(Nx, Ny);
mf2d.alpha_coeff = m2.ac*ones(Nx, Ny);
% Alpha Power
mf2d.alpha_power = 2;
mf2d.sound_speed(251:end,:) = m1.cp;
mf2d.density(251:end,:) = m1.rho;
mf2d.alpha_coeff(251:end,:) = m1.ac;
%% Source definition
source_freq = 2000e3; % [Hz]
source_strength = 1e6; % [Pa]
source_cycles = 5; % number of tone burst cycles
source2d.p_mask = logical(zeros(Nx,Ny));
source2d.p_mask(1,:) = 1;
source2d.p = source_strength*toneBurst(1/kgrid2d.dt, source_freq, source_cycles, 'Plot', false, 'Envelope', 'Rectangular');
%% run the fluid simulation
input_args = { 'PMLSize', PML_size,'PlotLayout', true, 'PMLAlpha', 2, 'PlotPML', false, ...
'PlotScale', [-0.5, 0.5]*source_strength , 'DataCast', 'gpuArray-single'...
'PMLInside', false};
%% alpha_power = 2
sensor_data_fluid2d1 = kspaceFirstOrder2DG(kgrid2d, mf2d, source2d, sensor2d, input_args{:});
%% alpha_power = 2 & 'no_dispersion'
mf2d.alpha_mode = 'no_absorption';
sensor_data_fluid2d2 = kspaceFirstOrder2DG(kgrid2d, mf2d, source2d, sensor2d, input_args{:});
%% alpha_power = 1.01
clear mf2d
%% Correction of the absorption factors at frequency of sender for lower alpha_power
m1.ac = m1.ac*4;
m2.ac = m2.ac*4;
mf2d.sound_speed = m2.cp*ones(Nx, Ny);
mf2d.density = m2.rho*ones(Nx, Ny);
mf2d.alpha_coeff = m2.ac*ones(Nx, Ny);
mf2d.alpha_power = 1.01;
mf2d.sound_speed(251:end,:) = m1.cp;
mf2d.density(251:end,:) = m1.rho;
mf2d.alpha_coeff(251:end,:) = m1.ac;
sensor_data_fluid2d3 = kspaceFirstOrder2DG(kgrid2d, mf2d, source2d, sensor2d, input_args{:});
%% alpha_power = 1.01 & 'no_dispersion'
mf2d.alpha_mode = 'no_absorption';
sensor_data_fluid2d4 = kspaceFirstOrder2DG(kgrid2d, mf2d, source2d, sensor2d, input_args{:});
%% Visualization
figure
y_vec = (kgrid2d.y_vec)*1e3;
x_vec = (kgrid2d.x_vec)*1e3;
t = [1:size(sensor_data_fluid2d1.p,2)]*kgrid2d.dt*1e6;
plot(t,mean(sensor_data_fluid2d1.p),t,0.1E6+mean(sensor_data_fluid2d2.p),...
t,0.2E6+mean(sensor_data_fluid2d3.p),t,0.3E6+mean(sensor_data_fluid2d4.p))
xline((0.025/m2.cp*2)*1e6);
legend(["alpha power = 2" "alpha power = 2 & 'no dispersion'"...
"alpha power = 1.01" "alpha power = 1.01 & 'no dispersion'" ...
"expected first echo"]);
This code runs 4 simulations of the same impuls echo setup with 25mm water and 25mm ABS. The main parameters are alpha_power and alpha_mode. If you try out the code, you will see that alpha_mode has no influence in this example and that alpha_power has a huge influence on the time of travel.
I've seen the same effect in my other simulations and I am hard stuck, because I need proper results for time of travel and even the simplest examples like the one above yield wrong results. I am still hoping that it is a simple coding error of mine and not a limitation of k-wave.
Best regards,
Denis