![]() |
kspaceFirstOrder3D-OMP 1.0
The C++ implementation of the k-wave toolbox for the time-domain simulation of acoustic wave fields in 3D
|
00001 00033 #include <MatrixClasses/UXYZ_SGXYZMatrix.h> 00034 #include <MatrixClasses/FFTWComplexMatrix.h> 00035 00036 00037 using namespace std; 00038 //----------------------------------------------------------------------------// 00039 // Constants // 00040 //----------------------------------------------------------------------------// 00041 00042 00043 //----------------------------------------------------------------------------// 00044 // Implementation // 00045 // public methods // 00046 //----------------------------------------------------------------------------// 00047 00048 00049 00055 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2(TRealMatrix& dt_rho0_sg, TFFTWComplexMatrix& FFT){ 00056 00057 FFT.Compute_iFFT_3D_C2R(*this); 00058 00059 00060 const float Divider = 1.0f/(2.0f *pTotalElementCount); 00061 //dt_rho0_sgx .* real... 00062 00063 #ifndef __NO_OMP__ 00064 #pragma omp parallel for schedule (static) 00065 #endif 00066 for (size_t i = 0; i < pTotalElementCount; i++){ 00067 pMatrixData[i] = dt_rho0_sg[i] * (pMatrixData[i] * Divider); 00068 } 00069 00070 00071 }// end of Compute_dt_div_rho_sgx_mul 00072 //------------------------------------------------------------------------------ 00073 00074 00075 00081 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2(float dt_rho_0_sgx, TFFTWComplexMatrix& FFT){ 00082 FFT.Compute_iFFT_3D_C2R(*this); 00083 00084 00085 const float Divider = 1.0f/(2.0f *pTotalElementCount) * dt_rho_0_sgx; 00086 //dt_rho0_sgx .* real... 00087 00088 #ifndef __NO_OMP__ 00089 #pragma omp parallel for schedule (static) 00090 #endif 00091 for (size_t i = 0; i < pTotalElementCount; i++){ 00092 pMatrixData[i] = pMatrixData[i] * Divider; 00093 } 00094 00095 00096 }// end of Compute_dt_rho_sg_mul_ifft_div_2 00097 //------------------------------------------------------------------------------ 00098 00099 00106 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_x 00107 (float dt_rho_0_sgx, TRealMatrix & dxudxn_sgx, TFFTWComplexMatrix& FFT){ 00108 00109 00110 00111 FFT.Compute_iFFT_3D_C2R(*this); 00112 00113 const float Divider = 1.0f/(2.0f *pTotalElementCount) * dt_rho_0_sgx; 00114 //dt_rho0_sgx .* real... 00115 00116 00117 #ifndef __NO_OMP__ 00118 #pragma omp for schedule (static) 00119 #endif 00120 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00121 00122 register size_t i = z* pDimensionSizes.Y * pDimensionSizes.X; 00123 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00124 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00125 pMatrixData[i] = pMatrixData[i] * Divider * dxudxn_sgx[x]; 00126 i++; 00127 } // x 00128 } // y 00129 } // z 00130 00131 00132 00133 00134 }// end of Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_x 00135 //------------------------------------------------------------------------------ 00136 00137 00144 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_y 00145 (float dt_rho_0_sgy, TRealMatrix & dyudyn_sgy, TFFTWComplexMatrix& FFT){ 00146 00147 00148 00149 FFT.Compute_iFFT_3D_C2R(*this); 00150 00151 const float Divider = 1.0f/(2.0f *pTotalElementCount) * dt_rho_0_sgy; 00152 //dt_rho0_sgx .* real... 00153 00154 00155 #ifndef __NO_OMP__ 00156 #pragma omp for schedule (static) 00157 #endif 00158 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00159 00160 register size_t i = z* pDimensionSizes.Y * pDimensionSizes.X; 00161 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00162 const float dyudyn_sgy_data = dyudyn_sgy[y] * Divider; 00163 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00164 pMatrixData[i] = pMatrixData[i] * dyudyn_sgy_data; 00165 i++; 00166 } // x 00167 } // y 00168 } // z 00169 00170 00171 00172 00173 }// end of Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_y 00174 //------------------------------------------------------------------------------ 00175 00182 void Tuxyz_sgxyzMatrix::Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_z 00183 (float dt_rho_0_sgz, TRealMatrix & dzudzn_sgz, TFFTWComplexMatrix& FFT){ 00184 00185 00186 00187 00188 FFT.Compute_iFFT_3D_C2R(*this); 00189 00190 const float Divider = 1.0f/(2.0f *pTotalElementCount) * dt_rho_0_sgz; 00191 //dt_rho0_sgx .* real... 00192 00193 00194 #ifndef __NO_OMP__ 00195 #pragma omp for schedule (static) 00196 #endif 00197 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00198 00199 register size_t i = z* pDimensionSizes.Y * pDimensionSizes.X; 00200 const float dzudzn_sgz_data = dzudzn_sgz[z] * Divider; 00201 00202 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00203 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00204 pMatrixData[i] = pMatrixData[i] * dzudzn_sgz_data; 00205 i++; 00206 } // x 00207 } // y 00208 } // z 00209 00210 00211 00212 }// end of Compute_dt_rho_sg_mul_ifft_div_2_scalar_nonuniform_z 00213 //------------------------------------------------------------------------------ 00214 00215 00216 00224 void Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize(TRealMatrix& FFT_p, TRealMatrix& dt_rho0, TRealMatrix& pml){ 00225 00226 const float Divider = 1.0f / pTotalElementCount; 00227 00228 #ifndef __NO_OMP__ 00229 #pragma omp for schedule (static) 00230 #endif 00231 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00232 00233 register size_t i = z* p2DDataSliceSize; 00234 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00235 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00236 00237 register float pMatrixElement = pMatrixData[i]; 00238 00239 //FFT_p.ElementMultiplyMatrices(dt_rho0); 00240 const float FFT_p_el = Divider * FFT_p[i] * dt_rho0[i]; 00241 00242 //BSXElementRealMultiply_1D_X(abc); 00243 pMatrixElement *= pml[x]; 00244 00245 //ElementSubMatrices(FFT_p); 00246 pMatrixElement -= FFT_p_el; 00247 00248 //BSXElementRealMultiply_1D_X(abc); 00249 pMatrixData[i] = pMatrixElement * pml[x]; 00250 00251 i++; 00252 } // x 00253 } // y 00254 } // z 00255 00256 }// end of Compute_ux_sgx_normalize_Optimized 00257 //------------------------------------------------------------------------------ 00258 00265 void Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize_scalar_uniform(TRealMatrix& FFT_p, float dt_rho0, TRealMatrix& pml){ 00266 00267 const float Divider = dt_rho0 / pTotalElementCount; 00268 00269 #ifndef __NO_OMP__ 00270 #pragma omp for schedule (static) 00271 #endif 00272 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00273 00274 register size_t i = z* p2DDataSliceSize; 00275 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00276 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00277 00278 register float pMatrixElement = pMatrixData[i]; 00279 00280 //FFT_p.ElementMultiplyMatrices(dt_rho0); 00281 const float FFT_p_el = Divider * FFT_p[i]; 00282 00283 //BSXElementRealMultiply_1D_X(abc); 00284 pMatrixElement *= pml[x]; 00285 00286 //ElementSubMatrices(FFT_p); 00287 pMatrixElement -= FFT_p_el; 00288 00289 //BSXElementRealMultiply_1D_X(abc); 00290 pMatrixData[i] = pMatrixElement * pml[x]; 00291 00292 i++; 00293 } // x 00294 } // y 00295 } // z 00296 00297 00298 }// end of Compute_ux_sgx_normalize 00299 //------------------------------------------------------------------------------ 00300 00301 00309 void Tuxyz_sgxyzMatrix::Compute_ux_sgx_normalize_scalar_nonuniform(TRealMatrix& FFT_p, float dt_rho0, TRealMatrix & dxudxn_sgx, TRealMatrix& pml){ 00310 00311 const float Divider = dt_rho0 / pTotalElementCount; 00312 00313 #ifndef __NO_OMP__ 00314 #pragma omp for schedule (static) 00315 #endif 00316 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00317 00318 register size_t i = z* p2DDataSliceSize; 00319 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00320 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00321 00322 register float pMatrixElement = pMatrixData[i]; 00323 00324 //FFT_p.ElementMultiplyMatrices(dt_rho0); 00325 const float FFT_p_el = (Divider * dxudxn_sgx[x]) * FFT_p[i]; 00326 00327 //BSXElementRealMultiply_1D_X(abc); 00328 pMatrixElement *= pml[x]; 00329 00330 //ElementSubMatrices(FFT_p); 00331 pMatrixElement -= FFT_p_el; 00332 00333 //BSXElementRealMultiply_1D_X(abc); 00334 pMatrixData[i] = pMatrixElement * pml[x]; 00335 00336 i++; 00337 } // x 00338 } // y 00339 } // z 00340 00341 00342 }// end of Compute_ux_sgx_normalize_scalar_nonuniform 00343 //------------------------------------------------------------------------------ 00344 00345 00353 void Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize(TRealMatrix& FFT_p, TRealMatrix& dt_rho0, TRealMatrix& pml){ 00354 00355 00356 const float Divider = 1.0f / pTotalElementCount; 00357 #ifndef __NO_OMP__ 00358 #pragma omp for schedule (static) 00359 #endif 00360 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00361 00362 size_t i = z* p2DDataSliceSize; 00363 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00364 const float pml_y = pml[y]; 00365 00366 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00367 00368 register float pMatrixElement = pMatrixData[i]; 00369 00370 //FFT_p.ElementMultiplyMatrices(dt_rho0); 00371 const float FFT_p_el = Divider * FFT_p[i] * dt_rho0[i] ; 00372 00373 //BSXElementRealMultiply_1D_X(abc); 00374 pMatrixElement *= pml_y; 00375 00376 //ElementSubMatrices(FFT_p); 00377 pMatrixElement -= FFT_p_el; 00378 00379 //BSXElementRealMultiply_1D_X(abc); 00380 pMatrixData[i] = pMatrixElement * pml_y; 00381 00382 i++; 00383 00384 } // x 00385 } // y 00386 } // z 00387 00388 }// end of Compute_uy_sgy_normalize_Optimized 00389 //------------------------------------------------------------------------------ 00390 00391 00398 void Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize_scalar_uniform(TRealMatrix& FFT_p, float dt_rho0, TRealMatrix& pml){ 00399 00400 const float Divider = dt_rho0 / pTotalElementCount; 00401 #ifndef __NO_OMP__ 00402 #pragma omp for schedule (static) 00403 #endif 00404 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00405 00406 size_t i = z* p2DDataSliceSize; 00407 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00408 const float pml_y = pml[y]; 00409 00410 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00411 00412 register float pMatrixElement = pMatrixData[i]; 00413 00414 //FFT_p.ElementMultiplyMatrices(dt_rho0); 00415 const float FFT_p_el = Divider * FFT_p[i] ; 00416 00417 //BSXElementRealMultiply_1D_X(abc); 00418 pMatrixElement *= pml_y; 00419 00420 //ElementSubMatrices(FFT_p); 00421 pMatrixElement -= FFT_p_el; 00422 00423 //BSXElementRealMultiply_1D_X(abc); 00424 pMatrixData[i] = pMatrixElement * pml_y; 00425 00426 i++; 00427 00428 } // x 00429 } // y 00430 } // z 00431 00432 00433 }// end of Compute_ux_sgx_normalize 00434 //------------------------------------------------------------------------------ 00435 00436 00444 void Tuxyz_sgxyzMatrix::Compute_uy_sgy_normalize_scalar_nonuniform(TRealMatrix& FFT_p, float dt_rho0, TRealMatrix & dyudyn_sgy, TRealMatrix& pml){ 00445 00446 const float Divider = dt_rho0 / pTotalElementCount; 00447 #ifndef __NO_OMP__ 00448 #pragma omp for schedule (static) 00449 #endif 00450 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00451 00452 size_t i = z* p2DDataSliceSize; 00453 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00454 const float pml_y = pml[y]; 00455 const float dyudyn_sgy_data = dyudyn_sgy[y]; 00456 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00457 00458 register float pMatrixElement = pMatrixData[i]; 00459 00460 //FFT_p.ElementMultiplyMatrices(dt_rho0); 00461 const float FFT_p_el = (Divider * dyudyn_sgy_data) * FFT_p[i] ; 00462 00463 //BSXElementRealMultiply_1D_X(abc); 00464 pMatrixElement *= pml_y; 00465 00466 //ElementSubMatrices(FFT_p); 00467 pMatrixElement -= FFT_p_el; 00468 00469 //BSXElementRealMultiply_1D_X(abc); 00470 pMatrixData[i] = pMatrixElement * pml_y; 00471 00472 i++; 00473 00474 } // x 00475 } // y 00476 } // z 00477 00478 00479 }// end of Compute_uy_sgy_normalize_scalar_nonuniform 00480 //------------------------------------------------------------------------------ 00481 00482 00490 void Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize(TRealMatrix& FFT_p, TRealMatrix& dt_rho0, TRealMatrix& pml){ 00491 00492 const float Divider = 1.0f / pTotalElementCount; 00493 00494 #ifndef __NO_OMP__ 00495 #pragma omp for schedule (static) 00496 #endif 00497 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00498 size_t i = z* p2DDataSliceSize; 00499 const float pml_z = pml[z]; 00500 00501 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00502 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00503 00504 register float pMatrixElement = pMatrixData[i]; 00505 00506 //FFT_p.ElementMultiplyMatrices(dt_rho0); 00507 const float FFT_p_el = Divider * FFT_p[i] * dt_rho0[i]; 00508 00509 //BSXElementRealMultiply_1D_X(abc); 00510 pMatrixElement *= pml_z; 00511 00512 //ElementSubMatrices(FFT_p); 00513 pMatrixElement -= FFT_p_el; 00514 00515 //BSXElementRealMultiply_1D_X(abc); 00516 00517 pMatrixData[i] = pMatrixElement * pml_z; 00518 00519 i++; 00520 00521 } // x 00522 } // y 00523 } // z 00524 00525 }// end of Compute_uz_sgz_normalize_fast 00526 //------------------------------------------------------------------------------ 00527 00534 void Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize_scalar_uniform(TRealMatrix& FFT_p, float& dt_rho0, TRealMatrix& pml){ 00535 00536 const float Divider = dt_rho0 / pTotalElementCount; 00537 00538 #ifndef __NO_OMP__ 00539 #pragma omp for schedule (static) 00540 #endif 00541 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00542 size_t i = z* p2DDataSliceSize; 00543 const float pml_z = pml[z]; 00544 00545 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00546 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00547 00548 register float pMatrixElement = pMatrixData[i]; 00549 00550 //FFT_p.ElementMultiplyMatrices(dt_rho0); 00551 const float FFT_p_el = Divider * FFT_p[i]; 00552 00553 //BSXElementRealMultiply_1D_X(abc); 00554 pMatrixElement *= pml_z; 00555 00556 //ElementSubMatrices(FFT_p); 00557 pMatrixElement -= FFT_p_el; 00558 00559 //BSXElementRealMultiply_1D_X(abc); 00560 00561 pMatrixData[i] = pMatrixElement * pml_z; 00562 00563 i++; 00564 00565 } // x 00566 } // y 00567 } // z 00568 00569 00570 00571 }// end of Compute_uz_sgz_normalize_scalar_uniform 00572 //------------------------------------------------------------------------------ 00573 00574 00582 void Tuxyz_sgxyzMatrix::Compute_uz_sgz_normalize_scalar_nonuniform(TRealMatrix& FFT_p, float& dt_rho0,TRealMatrix & dzudzn_sgz, TRealMatrix& pml){ 00583 00584 const float Divider = dt_rho0 / pTotalElementCount; 00585 00586 #ifndef __NO_OMP__ 00587 #pragma omp for schedule (static) 00588 #endif 00589 for (size_t z = 0; z < pDimensionSizes.Z; z++){ 00590 size_t i = z* p2DDataSliceSize; 00591 const float pml_z = pml[z]; 00592 const float dzudzn_sgz_data = dzudzn_sgz[z]; 00593 00594 for (size_t y = 0; y< pDimensionSizes.Y; y++){ 00595 for (size_t x = 0; x < pDimensionSizes.X; x++){ 00596 00597 register float pMatrixElement = pMatrixData[i]; 00598 00599 //FFT_p.ElementMultiplyMatrices(dt_rho0); 00600 const float FFT_p_el = (Divider * dzudzn_sgz_data)* FFT_p[i]; 00601 00602 //BSXElementRealMultiply_1D_X(abc); 00603 pMatrixElement *= pml_z; 00604 00605 //ElementSubMatrices(FFT_p); 00606 pMatrixElement -= FFT_p_el; 00607 00608 //BSXElementRealMultiply_1D_X(abc); 00609 00610 pMatrixData[i] = pMatrixElement * pml_z; 00611 00612 i++; 00613 00614 } // x 00615 } // y 00616 } // z 00617 00618 00619 }// end of Compute_uz_sgz_normalize_scalar_nonuniform 00620 //------------------------------------------------------------------------------ 00621 00622 00629 void Tuxyz_sgxyzMatrix::AddTransducerSource(TLongMatrix& u_source_index, TLongMatrix& delay_mask, TRealMatrix& transducer_signal){ 00630 00631 //#ifndef __NO_OMP__ 00632 // #pragma omp parallel for schedule (static) 00633 //#endif 00634 for (size_t i = 0; i < u_source_index.GetTotalElementCount(); i++){ 00635 pMatrixData[u_source_index[i]] += transducer_signal[delay_mask[i]]; 00636 delay_mask[i]++; 00637 } 00638 00639 00640 }// end of AddTransducerSource 00641 //------------------------------------------------------------------------------ 00642 00643 00644 00645 00656 void Tuxyz_sgxyzMatrix::Add_u_source(TRealMatrix &u_source_input, TLongMatrix & u_source_index, int t_index, long u_source_mode, long u_source_many){ 00657 00658 00659 size_t index2D = t_index; 00660 if (u_source_many != 0) { // is 2D 00661 index2D = t_index * u_source_index.GetTotalElementCount(); 00662 } 00663 00664 00665 00666 if (u_source_mode == 0){ 00667 for (size_t i = 0; i < u_source_index.GetTotalElementCount(); i++){ 00668 00669 pMatrixData[u_source_index[i]] = u_source_input[index2D]; 00670 00671 if (u_source_many != 0) index2D ++; 00672 } 00673 }// end of dirichlet 00674 00675 if (u_source_mode == 1){ 00676 for (size_t i = 0; i < u_source_index.GetTotalElementCount(); i++){ 00677 00678 pMatrixData[u_source_index[i]] += u_source_input[index2D]; 00679 00680 if (u_source_many != 0) index2D ++; 00681 } 00682 }// end of add 00683 00684 }// end of Add_u_source 00685 //------------------------------------------------------------------------------ 00686 00687 00688 00689 00690 00691 00692 00693 00694 //----------------------------------------------------------------------------// 00695 // Implementation // 00696 // private methods // 00697 //----------------------------------------------------------------------------//