Hi,
I'm wondering is there a C++ or GPU version of kspaceFirstOrder2D code? I'm generating a large amount of 2D k-wave simulation data, and I'm looking for ways to speed up 2D simulation.
k-Wave
A MATLAB toolbox for the time-domain
simulation of acoustic wave fields
C++/GPU version of kspaceFirstOrder2D
(10 posts) (4 voices)-
Posted 4 years ago #
-
Hi dyangsteven,
Yes - it will be out in a few weeks with the next release.
Brad.
Posted 4 years ago # -
Hi Brad,
I encounter an issue when running the new kspaceFirstOrder2DG: with the same code, I obtain a time of flight that is ~10% shorter than with kspaceFirstOrder2D. kspaceFirstOrder2D is correct.
Did I miss something?
RegardsPosted 4 years ago # -
Hi gpr,
What sort of sensor mask are you using?
With a binary sensor mask, the two codes should agree to machine precision.
If you're happy to post an example, I'll take a look.
Brad.
Posted 4 years ago # -
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 allgpuDeviceCount
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_sizex = 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_sizedisplay('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);
endfor 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
endsource.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)
drawnowdisplay(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 == 1tic
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{:});
tocend
sensor_data_c{p} = sensor_data;
%%
figure(1+p);
imagesc(ti.*1e6,kgrid.x_vec,sensor_data)
caxis([-3e3 3e3]);end
Posted 4 years ago # -
Hi gpr,
The problems is that
medium.alpha_mode = no_dispersion
is not implemented in the C++/CUDA codes. If you remove this line, the two agree to machine precision.I'll add it to the bug list to add an error condition if the code is used in this way.
Brad.
Posted 4 years ago # -
Hi Bradley,
Adding to this, is the sensor directivity implemented in C++/CUDA versions?
sensor.directivity_angle = (any angle)
sensor.directivity_pattern = 'pressure'; % alternatively: 'gradient'
sensor.directivity_size = transducer.element_width;Seems to give a large difference between 2D and 2DG. (On Linux, own compiled source.)
Bram
Posted 4 years ago # -
Hi Bram,
No, this is not included in the C++/CUDA version either.
I'll add this to the list of error conditions also.
Brad.
Posted 4 years ago # -
Hi Brad, thank you very much, I appreciate.
Posted 4 years ago # -
Hi Brad, thank you indeed! I'll just use the interference effect then.
Posted 4 years ago #
Reply
You must log in to post.