C++/GPU version of kspaceFirstOrder2D

Hi Brad,
yes the masks are binary.

Absolutely, here is a code to paste in a matlab script. I just change the variable "simvis" at the beginning, which only changes the run mode (or so I thought). This, however, provides a different result.

Regards

__________________________________________________________________________________________
__________________________________________________________________________________________

clear all
clc
close all

gpuDeviceCount

d = gpuDevice

simvis = 1;

% an_lim = 20;%26.5651;
% angle = linspace(-an_lim,an_lim,11);
angle = 0;

ref_speed_of_sound = 1500;
centerfreq = 3.5e6;
lamb = ref_speed_of_sound/centerfreq;
press = 1e5; %;%5e6;

afety_co = 0.1e-6;
afety_co_2 = 0.5e-6;

run_on_gpu = 1;
if run_on_gpu == 0
p = gcp;
end

%% create the computational grid

trans_foc = 30/1e3;

sub_lamb_coef = 4;
pix_size = lamb/sub_lamb_coef;

diam_trans_m = 40e-3;
diam_trans_pi = ceil(diam_trans_m/pix_size);

depth_im_m = trans_foc*1.1 ; % 3e-3; % 100; %120; % in wavelengths
depth_im_pi = depth_im_m./pix_size;

PML_size = 32; % size of the PML in grid points
PMLAlpha = 50;

x = 2;
a = 1;
Nx = x^a;
while Nx*pix_size <= diam_trans_m + 2 * PML_size*pix_size
display(Nx*pix_size)
a = a+1;
Nx = x^a;
end
Nx = Nx - 2 * PML_size

x = 2;
a = 1;
Nz = x^a;
while Nz*pix_size <= depth_im_m + 2 * PML_size*pix_size
a = a+1;
Nz = x^a;
end
Nz = Nz - 2 * PML_size

display('image width (mm)');
display(Nx*pix_size*1e3);
display('image depth (mm)');
display(Nz*pix_size*1e3);
display('transducer focus (mm)');
display(trans_foc*1e3);
display('transducer diameter (mm)');
display(diam_trans_m*1e3);

kgrid = kWaveGrid(Nx, pix_size, Nz, pix_size);

%% material properties

crod = 2000;
rhorod = 5000;
BonArod = 0;
alpha_coeffrod = 1;

Cliver = 1559;
BonAliver = 8;
rholiver = 1020;
alpha_coeff_liver = 0.8;

%% map

% dx = 0.02;
dx = 0.002;
R = random('Uniform',1-dx,1+dx,size(kgrid.x,1),size(kgrid.x,2),size(kgrid.x,3));

sound_speed_lo = Cliver.*R;
density_lo = rholiver.*R;
BonA_lo = BonAliver.*R;
alpha_coeff_lo = alpha_coeff_liver.*R;
c_med = Cliver;

medium.sound_speed = sound_speed_lo;
medium.density = density_lo;
medium.BonA = BonA_lo; % 0; %
medium.alpha_coeff = alpha_coeff_lo;

medium.alpha_mode = 'no_dispersion';
medium.alpha_power = 1.01; %

%% make rod

disc = zeros(Nx,Nz);

cx = round(Nx/2);
cy = round(Nz/2);
disc = disc+makeDisc(Nx, Nz, cx, cy, sub_lamb_coef/2);

dx = 0.001;
R = random('Uniform',1-dx,1+dx,size(kgrid.x,1),size(kgrid.x,2),size(kgrid.x,3));

medium.sound_speed = (medium.sound_speed.*(1-disc))+disc.*crod.*R;
medium.density = (medium.density.*(1-disc))+disc.*rhorod.*R;
medium.BonA = (medium.BonA.*(1-disc))+disc.*BonArod.*R;
medium.alpha_coeff = (medium.alpha_coeff.*(1-disc))+disc.*alpha_coeffrod.*R;

ftsize = 20;

figure(33)
colormap('hot');
sos_map = medium.sound_speed( :,:,round(size(medium.sound_speed,3)/2));
imagesc(kgrid.y_vec.*1e3,(kgrid.x_vec-min(kgrid.x_vec)).*1e3,sos_map)
xlabel('lateral position (mm)','interpreter', 'latex','fontsize',ftsize)
ylabel('axial position (mm)','interpreter', 'latex','fontsize',ftsize)
c = colorbar;
colormap('hot');

ylabel(c,'speed of sound (m/s)');
drawnow

%% to drive the transducer

deltt = 1.5/centerfreq;
tau =1/centerfreq/4;

%% create the time array and input signal

coti = 1;
kgrid.makeTime(medium.sound_speed);
ti = 0:mean(diff(kgrid.t_array))/coti:(Nz*pix_size/c_med+deltt)+(sqrt(Nz^2+Nx^2)*pix_size/c_med+deltt)+5e-6;

input_signal = sin(2*pi*centerfreq.*ti);

po = find(ti>deltt,1,'first');

ti_sig = ti;
ti_sig(po:end) = [];
input_signal(po:end) = [];
input_signal = input_signal./max(input_signal);

figure(11)
plot(ti_sig.*1e6,input_signal);
grid on

kgrid.t_array = ti;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% simulation %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

total_depth = Nz*pix_size;
iter = 1;

sensor_data_c = cell(1,length(angle));

for p = 1:length(angle)

%% sources
d = Nx*pix_size*sind(angle(p));

x = 1:Nx;
xtau = 30;
env = (tanh(x./xtau)+fliplr(tanh(x./xtau)))-1;

source.p = zeros(Nx,length(ti));

n_t = (linspace(0,(d/c_med/mean(diff(ti))),Nx));
n_t = round(n_t);

if d<0
n_t = n_t-min(n_t);
end

for i = 1:Nx
if n_t(i)>0
temp = [zeros(1,n_t(i)), input_signal];
temp(end+1 : length(ti)) = 0;
source.p(i,:) = temp.*press.*env(i);
else
temp = input_signal;
temp(end+1 : length(ti)) = 0;
source.p(i,:) = temp.*press.*env(i);
end
end

source.p_mask = zeros(Nx,Nz);
center = [Nx/2,trans_foc/pix_size];

delt = 0;
source.p_mask = zeros(Nx,Nz);
source.p_mask(:,1+delt) = 1;

figure;
imagesc( source.p_mask)
drawnow

display(sum(source.p_mask(:)));

%% sensor

sensor.mask = zeros(Nx,Nz);

sensor.mask(:,1+delt) = 1;

figure;
imagesc( sensor.mask )
display(sum(sensor.mask(:)))

%% set the input arguements

if run_on_gpu == 1 && simvis == 0

tic
input_args = {'PMLSize', PML_size, 'PMLAlpha', PMLAlpha,'PMLInside', false, 'PlotPML', false, ...
'Smooth', false, 'DataCast', 'gpuArray-single', 'CartInterp', 'nearest','PlotSim' , false,...
'PlotScale',[-press/10 press/10],'RecordMovie', false }; % 'PlotScale',,[-1e5 1e5]te'auto','RecordMovie', true }; % [-1e5 1e5] %
% run the simulation
filename = 'test_gui.h5';
sensor_data = kspaceFirstOrder2DG(kgrid, medium, source, sensor,input_args{:});
if isa(sensor_data,'gpuArray')
sensor_data = gather(sensor_data);
end
toc
%
elseif run_on_gpu == 1 && simvis == 1

tic
input_args = {'PMLSize', PML_size, 'PMLAlpha', PMLAlpha,'PMLInside', false, 'PlotPML', false, ...
'Smooth', false, 'DataCast', 'gpuArray-single', 'CartInterp', 'nearest','PlotSim' , true,...
'PlotScale',[-press/2 press/2],'RecordMovie', false }; % 'PlotScale',,[-1e5 1e5]te'auto','RecordMovie', true }; % [-1e5 1e5] %
% run the simulation
filename = 'test_gui.h5';
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,input_args{:});
if isa(sensor_data,'gpuArray')
sensor_data = gather(sensor_data);
end
toc
%

else

tic
input_args = {'PMLSize', PML_size, 'PMLAlpha', PMLAlpha,'PMLInside', false, 'PlotPML', false, ...
'Smooth', false, 'CartInterp', 'nearest','PlotSim' , true,...
'PlotScale','auto','RecordMovie', false }; % 'PlotScale',,[-1e5 1e5]te'auto','RecordMovie', true }; % [-1e5 1e5] %
% run the simulation
sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,input_args{:});
toc

end

sensor_data_c{p} = sensor_data;

figure(1+p);
imagesc(ti.*1e6,kgrid.x_vec,sensor_data)
caxis([-3e3 3e3]);

end

k-Wave

A MATLAB toolbox for the time-domain
simulation of acoustic wave fields