k-Wave User Forum » Topic: C++/GPU version of kspaceFirstOrder2D

k-Wave User Forum » Topic: C++/GPU version of kspaceFirstOrder2D http://www.k-wave.org/forum/topic/c-version-of-kspacefirstorder2d Support for the k-Wave MATLAB toolbox en-US Wed, 13 May 2026 01:14:49 +0000 http://bbpress.org/?v=1.0.2 <![CDATA[Search]]> q http://www.k-wave.org/forum/search.php bramth on "C++/GPU version of kspaceFirstOrder2D" http://www.k-wave.org/forum/topic/c-version-of-kspacefirstorder2d#post-7268 Thu, 05 Mar 2020 10:20:45 +0000 bramth 7268@http://www.k-wave.org/forum/ Hi Brad, thank you indeed! I'll just use the interference effect then. gpr on "C++/GPU version of kspaceFirstOrder2D" http://www.k-wave.org/forum/topic/c-version-of-kspacefirstorder2d#post-7262 Wed, 04 Mar 2020 19:11:45 +0000 gpr 7262@http://www.k-wave.org/forum/ Hi Brad, thank you very much, I appreciate. Bradley Treeby on "C++/GPU version of kspaceFirstOrder2D" http://www.k-wave.org/forum/topic/c-version-of-kspacefirstorder2d#post-7261 Wed, 04 Mar 2020 14:07:27 +0000 Bradley Treeby 7261@http://www.k-wave.org/forum/ 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. bramth on "C++/GPU version of kspaceFirstOrder2D" http://www.k-wave.org/forum/topic/c-version-of-kspacefirstorder2d#post-7260 Wed, 04 Mar 2020 13:54:41 +0000 bramth 7260@http://www.k-wave.org/forum/ 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 Bradley Treeby on "C++/GPU version of kspaceFirstOrder2D" http://www.k-wave.org/forum/topic/c-version-of-kspacefirstorder2d#post-7255 Wed, 04 Mar 2020 11:32:06 +0000 Bradley Treeby 7255@http://www.k-wave.org/forum/ Hi gpr, The problems is that <code>medium.alpha_mode = no_dispersion</code> 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. gpr on "C++/GPU version of kspaceFirstOrder2D" http://www.k-wave.org/forum/topic/c-version-of-kspacefirstorder2d#post-7251 Tue, 03 Mar 2020 18:02:47 +0000 gpr 7251@http://www.k-wave.org/forum/ 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 Bradley Treeby on "C++/GPU version of kspaceFirstOrder2D" http://www.k-wave.org/forum/topic/c-version-of-kspacefirstorder2d#post-7250 Tue, 03 Mar 2020 17:33:39 +0000 Bradley Treeby 7250@http://www.k-wave.org/forum/ 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. gpr on "C++/GPU version of kspaceFirstOrder2D" http://www.k-wave.org/forum/topic/c-version-of-kspacefirstorder2d#post-7249 Tue, 03 Mar 2020 16:57:33 +0000 gpr 7249@http://www.k-wave.org/forum/ 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? Regards Bradley Treeby on "C++/GPU version of kspaceFirstOrder2D" http://www.k-wave.org/forum/topic/c-version-of-kspacefirstorder2d#post-7222 Wed, 12 Feb 2020 12:31:46 +0000 Bradley Treeby 7222@http://www.k-wave.org/forum/ Hi dyangsteven, Yes - it will be out in a few weeks with the next release. Brad. dyangsteven on "C++/GPU version of kspaceFirstOrder2D" http://www.k-wave.org/forum/topic/c-version-of-kspacefirstorder2d#post-7207 Fri, 07 Feb 2020 21:18:05 +0000 dyangsteven 7207@http://www.k-wave.org/forum/ 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.