kspaceFirstOrder3D-CUDA  1.1
The CUDA/C++ implementation of the k-wave toolbox for the time-domain simulation of acoustic wave fields in 3D
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SolverCUDAKernels.cuh
Go to the documentation of this file.
1 /**
2  * @file SolverCUDAKernels.cuh
3  *
4  * @author Jiri Jaros \n
5  * Faculty of Information Technology \n
6  * Brno University of Technology \n
7  * jarosjir@fit.vutbr.cz
8  *
9  * @brief Name space for all CUDA kernels used in the 3D solver
10  *
11  * @version kspaceFirstOrder3D 3.4
12  *
13  * @date 11 March 2013, 13:10 (created) \n
14  * 27 July 2016, 15:09 (revised)
15  *
16  * @section License
17  * This file is part of the C++ extension of the k-Wave Toolbox
18  * (http://www.k-wave.org).\n Copyright (C) 2016 Jiri Jaros and Bradley Treeby.
19  *
20  * This file is part of the k-Wave. k-Wave is free software: you can redistribute it and/or modify
21  * it under the terms of the GNU Lesser General Public License as published by the Free Software
22  * Foundation, either version 3 of the License, or (at your option) any later version.
23  *
24  * k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even
25  * the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
26  * General Public License for more details.
27  *
28  * You should have received a copy of the GNU Lesser General Public License along with k-Wave.
29  * If not, see http://www.gnu.org/licenses/.
30  */
31 
32 #ifndef SOLVER_CUDA_KERNELS_CUH
33 #define SOLVER_CUDA_KERNELS_CUH
34 
39 
40 #include <Utils/DimensionSizes.h>
41 
42 #include <Parameters/Parameters.h>
44 
45 /**
46  * @namespace SolverCUDAKernels
47  * @brief List of cuda kernels used k-space first order 3D solver
48  * @details List of cuda kernels used k-space first order 3D solver
49  *
50  */
51 namespace SolverCUDAKernels
52 {
53  /// Get the CUDA architecture and GPU code version the code was compiled with
54  int GetCUDACodeVersion();
55 
56  //----------------------------------- velocity operations --------------------------------------//
57  /// Compute acoustic velocity for default case (heterogeneous).
58  void ComputeVelocity(TRealMatrix& ux_sgx,
59  TRealMatrix& uy_sgy,
60  TRealMatrix& uz_sgz,
61  const TRealMatrix& ifft_x,
62  const TRealMatrix& ifft_y,
63  const TRealMatrix& ifft_z,
64  const TRealMatrix& dt_rho0_sgx,
65  const TRealMatrix& dt_rho0_sgy,
66  const TRealMatrix& dt_rho0_sgz,
67  const TRealMatrix& pml_x,
68  const TRealMatrix& pml_y,
69  const TRealMatrix& pml_z);
70 
71 
72  /// Compute acoustic velocity, scalar and uniform case.
74  TRealMatrix& uy_sgy,
75  TRealMatrix& uz_sgz,
76  const TRealMatrix& ifft_x,
77  const TRealMatrix& ifft_y,
78  const TRealMatrix& ifft_z,
79  const TRealMatrix& pml_x,
80  const TRealMatrix& pml_y,
81  const TRealMatrix& pml_z);
82 
83  /// Compute acoustic velocity, scalar, non-uniform case.
85  TRealMatrix& uy_sgy,
86  TRealMatrix& uz_sgz,
87  const TRealMatrix& ifft_x,
88  const TRealMatrix& ifft_y,
89  const TRealMatrix& ifft_z,
90  const TRealMatrix& dxudxn_sgx,
91  const TRealMatrix& dyudyn_sgy,
92  const TRealMatrix& dzudzn_sgz,
93  const TRealMatrix& pml_x,
94  const TRealMatrix& pml_y,
95  const TRealMatrix& pml_z);
96 
97  //----------------------------------------- Sources --------------------------------------------//
98  /// Add transducer data source to X component.
99  void AddTransducerSource(TRealMatrix& ux_sgx,
100  const TIndexMatrix& u_source_index,
101  TIndexMatrix& delay_mask,
102  const TRealMatrix & transducer_signal);
103 
104  /// Add in velocity source terms.
105  void AddVelocitySource(TRealMatrix& uxyz_sgxyz,
106  const TRealMatrix& u_source_input,
107  const TIndexMatrix& u_source_index,
108  const size_t t_index);
109 
110  /// Add in pressure source term
111  void AddPressureSource(TRealMatrix& rhox,
112  TRealMatrix& rhoy,
113  TRealMatrix& rhoz,
114  const TRealMatrix& p_source_input,
115  const TIndexMatrix& p_source_index,
116  const size_t t_index);
117 
118  /// Compute velocity for the initial pressure problem.
119  void Compute_p0_Velocity(TRealMatrix& ux_sgx,
120  TRealMatrix& uy_sgy,
121  TRealMatrix& uz_sgz,
122  const TRealMatrix& dt_rho0_sgx,
123  const TRealMatrix& dt_rho0_sgy,
124  const TRealMatrix& dt_rho0_sgz);
125 
126  /// Compute acoustic velocity for initial pressure problem, if rho0_sgx is scalar, uniform grid.
127  void Compute_p0_Velocity(TRealMatrix& ux_sgx,
128  TRealMatrix& uy_sgy,
129  TRealMatrix& uz_sgz);
130 
131 
132  /// Compute acoustic velocity for initial pressure problem, if rho0_sgx is scalar, non uniform grid, x component.
134  TRealMatrix& uy_sgy,
135  TRealMatrix& uz_sgz,
136  const TRealMatrix& dxudxn_sgx,
137  const TRealMatrix& dyudyn_sgy,
138  const TRealMatrix& dzudzn_sgz);
139 
140 
141  //------------------------------------- pressure kernels ---------------------------------------//
142  /// Compute part of the new velocity - gradient of p.
144  TCUFFTComplexMatrix& fft_y,
145  TCUFFTComplexMatrix& fft_z,
146  const TRealMatrix& kappa,
147  const TComplexMatrix& ddx,
148  const TComplexMatrix& ddy,
149  const TComplexMatrix& ddz);
150 
151  /// Compute gradient of acoustic velocity on uniform grid.
153  TCUFFTComplexMatrix& fft_y,
154  TCUFFTComplexMatrix& fft_z,
155  const TRealMatrix& kappa,
156  const TComplexMatrix& ddx_k_shift_neg,
157  const TComplexMatrix& ddy_k_shift_neg,
158  const TComplexMatrix& ddz_k_shift_neg);
159 
160  /// Shift gradient of acoustic velocity on non-uniform grid.
162  TRealMatrix& duydy,
163  TRealMatrix& duzdz,
164  const TRealMatrix& dxudxn,
165  const TRealMatrix& dyudyn,
166  const TRealMatrix& dzudzn);
167 
168  /// Add initial pressure to p0 (as p0 source).
170  TRealMatrix& rhox,
171  TRealMatrix& rhoy,
172  TRealMatrix& rhoz,
173  const TRealMatrix& p0,
174  const bool Is_c2_scalar,
175  const float* c2);
176 
177  //------------------------------------- density kernels ----------------------------------------//
178  /// Calculate acoustic density for non-linear case, homogenous case.
180  TRealMatrix& rhoy,
181  TRealMatrix& rhoz,
182  const TRealMatrix& pml_x,
183  const TRealMatrix& pml_y,
184  const TRealMatrix& pml_z,
185  const TRealMatrix& duxdx,
186  const TRealMatrix& duydy,
187  const TRealMatrix& duzdz);
188 
189  /// Calculate acoustic density for non-linear case, heterogenous case.
191  TRealMatrix& rhoy,
192  TRealMatrix& rhoz,
193  const TRealMatrix& pml_x,
194  const TRealMatrix& pml_y,
195  const TRealMatrix& pml_z,
196  const TRealMatrix& duxdx,
197  const TRealMatrix& duydy,
198  const TRealMatrix& duzdz,
199  const TRealMatrix& rho0);
200 
201  /// Calculate acoustic density for linear case, homogenous case.
203  TRealMatrix& rhoy,
204  TRealMatrix& rhoz,
205  const TRealMatrix& pml_x,
206  const TRealMatrix& pml_y,
207  const TRealMatrix& pml_z,
208  const TRealMatrix& duxdx,
209  const TRealMatrix& duydy,
210  const TRealMatrix& duzdz);
211 
212  /// Calculate acoustic density for linear case, heterogeneous case.
214  TRealMatrix& rhoy,
215  TRealMatrix& rhoz,
216  const TRealMatrix& pml_x,
217  const TRealMatrix& pml_y,
218  const TRealMatrix& pml_z,
219  const TRealMatrix& duxdx,
220  const TRealMatrix& duydy,
221  const TRealMatrix& duzdz,
222  const TRealMatrix& rho0);
223 
224  //---------------------------------- new value of pressure -------------------------------------//
225  /// Calculate three temporary sums in the new pressure formula, non-linear absorbing case.
227  TRealMatrix& BonA_sum,
228  TRealMatrix& du_sum,
229  const TRealMatrix& rhox,
230  const TRealMatrix& rhoy,
231  const TRealMatrix& rhoz,
232  const TRealMatrix& duxdx,
233  const TRealMatrix& duydy,
234  const TRealMatrix& duzdz,
235  const bool is_BonA_scalar,
236  const float* BonA_matrix,
237  const bool is_rho0_scalar,
238  const float* rho0_matrix);
239 
240  /// Compute absorbing term with abosrb_nabla1 and absorb_nabla2.
242  TCUFFTComplexMatrix& fft2,
243  const TRealMatrix& absorb_nabla1,
244  const TRealMatrix& absorb_nabla2);
245 
246  /// Sum sub-terms to calculate new pressure, non-linear case.
248  const TRealMatrix& BonA_temp,
249  const bool is_c2_scalar,
250  const float* c2_matrix,
251  const bool is_tau_eta_scalar,
252  const float* absorb_tau,
253  const float* tau_matrix,
254  const float* absorb_eta,
255  const float* eta_matrix);
256 
257  /// Sum sub-terms to calculate new pressure, linear case.
259  const TRealMatrix& absorb_tau_temp,
260  const TRealMatrix& absorb_eta_temp,
261  const TRealMatrix& sum_rhoxyz,
262  const bool is_c2_scalar,
263  const float* c2_matrix,
264  const bool is_tau_eta_scalar,
265  const float* tau_matrix,
266  const float* eta_matrix);
267 
268  /// Sum sub-terms for new p, linear lossless case.
270  const TRealMatrix& rhox,
271  const TRealMatrix& rhoy,
272  const TRealMatrix& rhoz,
273  const bool is_c2_scalar,
274  const float* c2_matrix,
275  const bool is_BonA_scalar,
276  const float* BonA_matrix,
277  const bool is_rho0_scalar,
278  const float* rho0_matrix);
279 
280  /// Calculate two temporary sums in the new pressure formula, linear absorbing case.
281  void ComputePressurePartsLinear(TRealMatrix& sum_rhoxyz,
282  TRealMatrix& sum_rho0_du,
283  const TRealMatrix& rhox,
284  const TRealMatrix& rhoy,
285  const TRealMatrix& rhoz,
286  const TRealMatrix& duxdx,
287  const TRealMatrix& duydy,
288  const TRealMatrix& duzdz,
289  const bool is_rho0_scalar,
290  const float* rho0_matrix);
291 
292  /// Sum sub-terms for new p, linear lossless case.
294  const TRealMatrix& rhox,
295  const TRealMatrix& rhoy,
296  const TRealMatrix& rhoz,
297  const bool is_c2_scalar,
298  const float* c2_matrix);
299 
300  //----------------------------------- unstaggered velocity -------------------------------------//
301 
302  /// Transpose a real 3D matrix in the X-Y direction
303  void TrasposeReal3DMatrixXY(float* outputMatrix,
304  const float* inputMatrix,
305  const dim3& dimSizes);
306 
307  /// Transpose a real 3D matrix in the X-Y direction
308  void TrasposeReal3DMatrixXZ(float* outputMatrix,
309  const float* inputMatrix,
310  const dim3& dimSizes);
311 
312  /// Compute the velocity shift in Fourier space over the X axis
313  void ComputeVelocityShiftInX(TCUFFTComplexMatrix& cufft_shift_temp,
314  const TComplexMatrix& x_shift_neg_r);
315  /// Compute the velocity shift in Fourier space over the Y axis
316  void ComputeVelocityShiftInY(TCUFFTComplexMatrix& cufft_shift_temp,
317  const TComplexMatrix& y_shift_neg_r);
318  /// Compute the velocity shift in Fourier space over the Z axis
319  void ComputeVelocityShiftInZ(TCUFFTComplexMatrix& cufft_shift_temp,
320  const TComplexMatrix& z_shift_neg_r);
321 }
322 
323 #endif /*SOLVER_CUDA_KERNELS_CUH*/
void AddVelocitySource(TRealMatrix &uxyz_sgxyz, const TRealMatrix &u_source_input, const TIndexMatrix &u_source_index, const size_t t_index)
Add in velocity source terms.
int GetCUDACodeVersion()
Get the CUDA architecture and GPU code version the code was compiled with.
The header file for the class for setting CUDA kernel parameters.
void SumPressureLinearLossless(TRealMatrix &p, const TRealMatrix &rhox, const TRealMatrix &rhoy, const TRealMatrix &rhoz, const bool is_c2_scalar, const float *c2_matrix)
Sum sub-terms for new p, linear lossless case.
The header file containing the class for real matrices.
void ComputeVelocityGradientNonuniform(TRealMatrix &duxdx, TRealMatrix &duydy, TRealMatrix &duzdz, const TRealMatrix &dxudxn, const TRealMatrix &dyudyn, const TRealMatrix &dzudzn)
Shift gradient of acoustic velocity on non-uniform grid.
void ComputePressurePartsNonLinear(TRealMatrix &rho_sum, TRealMatrix &BonA_sum, TRealMatrix &du_sum, const TRealMatrix &rhox, const TRealMatrix &rhoy, const TRealMatrix &rhoz, const TRealMatrix &duxdx, const TRealMatrix &duydy, const TRealMatrix &duzdz, const bool is_BonA_scalar, const float *BonA_matrix, const bool is_rho0_scalar, const float *rho0_matrix)
Calculate three temporary sums in the new pressure formula, non-linear absorbing case.
void AddPressureSource(TRealMatrix &rhox, TRealMatrix &rhoy, TRealMatrix &rhoz, const TRealMatrix &p_source_input, const TIndexMatrix &p_source_index, const size_t t_index)
Add in pressure source term.
void ComputeVelocityScalarNonuniform(TRealMatrix &ux_sgx, TRealMatrix &uy_sgy, TRealMatrix &uz_sgz, const TRealMatrix &ifft_x, const TRealMatrix &ifft_y, const TRealMatrix &ifft_z, const TRealMatrix &dxudxn_sgx, const TRealMatrix &dyudyn_sgy, const TRealMatrix &dzudzn_sgz, const TRealMatrix &pml_x, const TRealMatrix &pml_y, const TRealMatrix &pml_z)
Compute acoustic velocity, scalar, non-uniform case.
The header file containing the parameters of the simulation.
void AddTransducerSource(TRealMatrix &ux_sgx, const TIndexMatrix &u_source_index, TIndexMatrix &delay_mask, const TRealMatrix &transducer_signal)
Add transducer data source to X component.
void SumPressureTermsNonlinear(TRealMatrix &p, const TRealMatrix &BonA_temp, const bool is_c2_scalar, const float *c2_matrix, const bool is_tau_eta_scalar, const float *absorb_tau, const float *tau_matrix, const float *absorb_eta, const float *eta_matrix)
Sum sub-terms to calculate new pressure, non-linear case.
void Compute_p0_AddInitialPressure(TRealMatrix &p, TRealMatrix &rhox, TRealMatrix &rhoy, TRealMatrix &rhoz, const TRealMatrix &p0, const bool Is_c2_scalar, const float *c2)
Add initial pressure to p0 (as p0 source).
void ComputePressurelGradient(TCUFFTComplexMatrix &fft_x, TCUFFTComplexMatrix &fft_y, TCUFFTComplexMatrix &fft_z, const TRealMatrix &kappa, const TComplexMatrix &ddx, const TComplexMatrix &ddy, const TComplexMatrix &ddz)
Compute part of the new velocity - gradient of p.
void ComputeDensityNonlinearHeterogeneous(TRealMatrix &rhox, TRealMatrix &rhoy, TRealMatrix &rhoz, const TRealMatrix &pml_x, const TRealMatrix &pml_y, const TRealMatrix &pml_z, const TRealMatrix &duxdx, const TRealMatrix &duydy, const TRealMatrix &duzdz, const TRealMatrix &rho0)
Calculate acoustic density for non-linear case, heterogenous case.
void ComputeVelocityGradient(TCUFFTComplexMatrix &fft_x, TCUFFTComplexMatrix &fft_y, TCUFFTComplexMatrix &fft_z, const TRealMatrix &kappa, const TComplexMatrix &ddx_k_shift_neg, const TComplexMatrix &ddy_k_shift_neg, const TComplexMatrix &ddz_k_shift_neg)
Compute gradient of acoustic velocity on uniform grid.
The header file containing the class that implements 3D FFT using the cuFFT interface.
void ComputeDensityLinearHomogeneous(TRealMatrix &rhox, TRealMatrix &rhoy, TRealMatrix &rhoz, const TRealMatrix &pml_x, const TRealMatrix &pml_y, const TRealMatrix &pml_z, const TRealMatrix &duxdx, const TRealMatrix &duydy, const TRealMatrix &duzdz)
Calculate acoustic density for linear case, homogenous case.
The header file containing the class for 64b integer matrices.
The header file containing the structure with 3D dimension sizes.
void ComputeVelocity(TRealMatrix &ux_sgx, TRealMatrix &uy_sgy, TRealMatrix &uz_sgz, const TRealMatrix &ifft_x, const TRealMatrix &ifft_y, const TRealMatrix &ifft_z, const TRealMatrix &dt_rho0_sgx, const TRealMatrix &dt_rho0_sgy, const TRealMatrix &dt_rho0_sgz, const TRealMatrix &pml_x, const TRealMatrix &pml_y, const TRealMatrix &pml_z)
Compute acoustic velocity for default case (heterogeneous).
void ComputeDensityNonlinearHomogeneous(TRealMatrix &rhox, TRealMatrix &rhoy, TRealMatrix &rhoz, const TRealMatrix &pml_x, const TRealMatrix &pml_y, const TRealMatrix &pml_z, const TRealMatrix &duxdx, const TRealMatrix &duydy, const TRealMatrix &duzdz)
Calculate acoustic density for non-linear case, homogenous case.
The header file with the class for complex matrices.
The class for real matrices.
Definition: RealMatrix.h:45
void ComputePressurePartsLinear(TRealMatrix &sum_rhoxyz, TRealMatrix &sum_rho0_du, const TRealMatrix &rhox, const TRealMatrix &rhoy, const TRealMatrix &rhoz, const TRealMatrix &duxdx, const TRealMatrix &duydy, const TRealMatrix &duzdz, const bool is_rho0_scalar, const float *rho0_matrix)
Calculate two temporary sums in the new pressure formula, linear absorbing case.
void ComputeVelocityShiftInX(TCUFFTComplexMatrix &cufft_shift_temp, const TComplexMatrix &x_shift_neg_r)
Compute the velocity shift in Fourier space over the X axis.
void SumPressureTermsLinear(TRealMatrix &p, const TRealMatrix &absorb_tau_temp, const TRealMatrix &absorb_eta_temp, const TRealMatrix &sum_rhoxyz, const bool is_c2_scalar, const float *c2_matrix, const bool is_tau_eta_scalar, const float *tau_matrix, const float *eta_matrix)
Sum sub-terms to calculate new pressure, linear case.
The class for 64b unsigned integers (indices). It is used for sensor_mask_index or sensor_corners_mas...
Definition: IndexMatrix.h:46
void ComputeDensityLinearHeterogeneous(TRealMatrix &rhox, TRealMatrix &rhoy, TRealMatrix &rhoz, const TRealMatrix &pml_x, const TRealMatrix &pml_y, const TRealMatrix &pml_z, const TRealMatrix &duxdx, const TRealMatrix &duydy, const TRealMatrix &duzdz, const TRealMatrix &rho0)
Calculate acoustic density for linear case, heterogeneous case.
void ComputeAbsorbtionTerm(TCUFFTComplexMatrix &fft1, TCUFFTComplexMatrix &fft2, const TRealMatrix &absorb_nabla1, const TRealMatrix &absorb_nabla2)
Compute absorbing term with abosrb_nabla1 and absorb_nabla2.
void Compute_p0_VelocityScalarNonUniform(TRealMatrix &ux_sgx, TRealMatrix &uy_sgy, TRealMatrix &uz_sgz, const TRealMatrix &dxudxn_sgx, const TRealMatrix &dyudyn_sgy, const TRealMatrix &dzudzn_sgz)
Compute acoustic velocity for initial pressure problem, if rho0_sgx is scalar, non uniform grid...
void Compute_p0_Velocity(TRealMatrix &ux_sgx, TRealMatrix &uy_sgy, TRealMatrix &uz_sgz, const TRealMatrix &dt_rho0_sgx, const TRealMatrix &dt_rho0_sgy, const TRealMatrix &dt_rho0_sgz)
Compute velocity for the initial pressure problem.
void TrasposeReal3DMatrixXY(float *outputMatrix, const float *inputMatrix, const dim3 &dimSizes)
Transpose a real 3D matrix in the X-Y direction.
void ComputeVelocityScalarUniform(TRealMatrix &ux_sgx, TRealMatrix &uy_sgy, TRealMatrix &uz_sgz, const TRealMatrix &ifft_x, const TRealMatrix &ifft_y, const TRealMatrix &ifft_z, const TRealMatrix &pml_x, const TRealMatrix &pml_y, const TRealMatrix &pml_z)
Compute acoustic velocity, scalar and uniform case.
void ComputeVelocityShiftInZ(TCUFFTComplexMatrix &cufft_shift_temp, const TComplexMatrix &z_shift_neg_r)
Compute the velocity shift in Fourier space over the Z axis.
The class for complex matrices.
Definition: ComplexMatrix.h:54
void TrasposeReal3DMatrixXZ(float *outputMatrix, const float *inputMatrix, const dim3 &dimSizes)
Transpose a real 3D matrix in the X-Y direction.
Class implementing 3D Real-To-Complex and Complex-To-Real transforms using CUDA FFT interface...
void SumPressureNonlinearLossless(TRealMatrix &p, const TRealMatrix &rhox, const TRealMatrix &rhoy, const TRealMatrix &rhoz, const bool is_c2_scalar, const float *c2_matrix, const bool is_BonA_scalar, const float *BonA_matrix, const bool is_rho0_scalar, const float *rho0_matrix)
Sum sub-terms for new p, linear lossless case.
void ComputeVelocityShiftInY(TCUFFTComplexMatrix &cufft_shift_temp, const TComplexMatrix &y_shift_neg_r)
Compute the velocity shift in Fourier space over the Y axis.