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
KSpaceFirstOrder3DSolver.cpp
Go to the documentation of this file.
1 /**
2  * @file KSpaceFirstOrder3DSolver.cpp
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 The implementation file containing the main class of the project responsible for the
10  * entire the entire 3D fluid simulation..
11  *
12  * @version kspaceFirstOrder3D 3.4
13  *
14  * @date 12 July 2012, 10:27 (created)\n
15  * 29 July 2016, 16:41 (revised)
16  *
17  * @section License
18  * This file is part of the C++ extension of the k-Wave Toolbox
19  * (http://www.k-wave.org).\n Copyright (C) 2016 Jiri Jaros and Bradley Treeby.
20  *
21  * This file is part of the k-Wave. k-Wave is free software: you can redistribute it and/or modify
22  * it under the terms of the GNU Lesser General Public License as published by the Free Software
23  * Foundation, either version 3 of the License, or (at your option) any later version.
24  *
25  * k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even
26  * the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser
27  * General Public License for more details.
28  *
29  * You should have received a copy of the GNU Lesser General Public License along with k-Wave.
30  * If not, see http://www.gnu.org/licenses/.
31  */
32 
33 // Linux build
34 #ifdef __linux__
35  #include <sys/resource.h>
36  #include <cmath>
37 #endif
38 
39 // Windows build
40 #ifdef _WIN64
41  #define _USE_MATH_DEFINES
42  #include <cmath>
43  #include <Windows.h>
44  #include <Psapi.h>
45  #pragma comment(lib, "Psapi.lib")
46 #endif
47 
48 #ifdef _OPENMP
49  #include <omp.h>
50 #endif
51 
52 #include <limits>
53 
55 
56 #include <Logger/ErrorMessages.h>
57 #include <Logger/Logger.h>
58 
61 
62 using std::string;
63 using std::ios;
64 
65 
66 //------------------------------------------------------------------------------------------------//
67 //------------------------------------------ Constants -------------------------------------------//
68 //------------------------------------------------------------------------------------------------//
69 
70 //------------------------------------------------------------------------------------------------//
71 //--------------------------------------- Public methods -----------------------------------------//
72 //------------------------------------------------------------------------------------------------//
73 
74 /**
75  * Constructor of the class.
76  */
78  matrixContainer(), outputStreamContainer(), parameters(TParameters::GetInstance()),
79  actPercent(0), isTimestepRightAfterRestore(false),
80  totalTime(), preProcessingTime(), dataLoadTime(), simulationTime(),
81  postProcessingTime(), iterationTime()
82 {
83  totalTime.Start();
84 
85  //Switch off default HDF5 error messages
86  H5Eset_auto(H5E_DEFAULT, NULL, NULL);
87 }// end of TKSpaceFirstOrder3DSolver
88 //--------------------------------------------------------------------------------------------------
89 
90 /**
91  * Destructor of the class.
92  */
94 {
95  // Delete CUDA FFT plans and related data
97 
98  // Free memory
99  FreeMemory();
100 
101  //Reset device after the run - recommended by CUDA SDK
102  cudaDeviceReset();
103 }// end of ~TKSpaceFirstOrder3DSolver
104 //--------------------------------------------------------------------------------------------------
105 
106 
107 /**
108  * The method allocates the matrix container and create all matrices and creates all output streams.
109  */
111 {
114 
115  // create container, then all matrices
118 
119  // add output streams into container
121 
123 }// end of AllocateMemory
124 //--------------------------------------------------------------------------------------------------
125 
126 /*
127  * The method frees all memory allocated by the class.
128  */
130 {
133 }// end of FreeMemory
134 //--------------------------------------------------------------------------------------------------
135 
136 /**
137  * Load data from the input file provided by the parameter class and creates the output time
138  * series streams.
139  */
141 {
142  // Load data from disk
145 
147 
148  // open and load input file
149  THDF5_File& inputFile = parameters.GetInputFile(); // file is opened (in Parameters)
150  THDF5_File& outputFile = parameters.GetOutputFile();
151  THDF5_File& checkpointFile = parameters.GetCheckpointFile();
152 
153  // Load data from disk
157 
158  // load data from the input file
160 
161  // close the input file
162  inputFile.Close();
163 
165 
166  // The simulation does not use check pointing or this is the first turn
167  bool recoverFromCheckpoint = (parameters.IsCheckpointEnabled() &&
169 
170 
171  if (recoverFromCheckpoint)
172  {
173  //--------------------------- Read data from the checkpoint file -----------------------------//
176 
177  // Open checkpoint file
178  checkpointFile.Open(parameters.GetCheckpointFileName());
179 
180  // Check the checkpoint file
182 
183  // read the actual value of t_index
184  size_t new_t_index;
185  checkpointFile.ReadScalarValue(checkpointFile.GetRootGroup(), t_index_NAME, new_t_index);
186  parameters.Set_t_index(new_t_index);
187 
188  // Read necessary matrices from the checkpoint file
190 
191  checkpointFile.Close();
193 
194  //----------------------------- Read data from the output file -------------------------------//
195  // Reopen output file for RW access
198 
199  outputFile.Open(parameters.GetOutputFileName(), H5F_ACC_RDWR);
200 
201  //Read file header of the output file
203 
204  // Restore elapsed time
205  LoadElapsedTimeFromOutputFile(outputFile);
206 
207  // Reopen streams
210  }
211  else
212  {
213  //-------------------------- First round of multi-leg simulation -----------------------------//
214  // Create the output file
217 
218  outputFile.Create(parameters.GetOutputFileName());
220 
221  // Create the steams, link them with the sampled matrices
222  // however DO NOT allocate memory!
224  }
225 
226  dataLoadTime.Stop();
227 
229  {
231  }
232 }// end of LoadInputData
233 //--------------------------------------------------------------------------------------------------
234 
235 /**
236  * This method computes k-space First Order 3D simulation. It launches calculation on a given
237  * dataset going through FFT initialization, pre-processing, main loop and post-processing phases.
238  *
239  */
241 {
243 
246 
247  TCUDAParameters& cudaParameters = parameters.GetCudaParameters();
248 
249  // fft initialisation and preprocessing
250  try
251  {
252  // initilaise all cuda FFT plans
255 
258 
259  // preprocessing is done on CPU and must pretend the CUDA configuration
261 
263  // Set kernel configurations
264  cudaParameters.SetKernelConfiguration();
265 
266  // Set up constant memory - copy over to GPU
267  // Constant memory uses some variables calculated during preprocessing
268  cudaParameters.SetUpDeviceConstants();
269 
271  }
272  catch (const std::exception& e)
273  {
276 
278  }
279 
280  // Logger header for simulation
285 
286 
288  cudaParameters.GetSolverGridSize1D(),
289  cudaParameters.GetSolverBlockSize1D());
290 
292 
294  cudaParameters.GetSamplerGridSize1D(),
295  cudaParameters.GetSamplerBlockSize1D());
296 
298 
299  // Main loop
300  try
301  {
303 
304  ComputeMainLoop();
305 
307 
309  }
310  catch (const std::exception& e)
311  {
314  }
315 
316  // Post processing region
318 
319  try
320  {
322  { // Checkpoint
328 
330  {
332  }
333 
335 
337  {
339  }
340  }
341  else
342  { // Finish
343 
348 
349  PostProcessing();
350 
351  // if checkpointing is enabled and the checkpoint file was created in the past, delete it
353  {
354  std::remove(parameters.GetCheckpointFileName().c_str());
355  }
357  }
358  }
359  catch (const std::exception &e)
360  {
363 
365  }
367 
368  // Final data written
369  try
370  {
373 
375  }
376  catch (const std::exception &e)
377  {
380  }
381 }// end of Compute()
382 //--------------------------------------------------------------------------------------------------
383 
384 /**
385  * Get peak CPU memory usage.
386  *
387  * @return peak memory usage in MBs.
388  */
390 {
391  // Linux build
392  #ifdef __linux__
393  struct rusage memUsage;
394  getrusage(RUSAGE_SELF, &memUsage);
395 
396  return memUsage.ru_maxrss >> 10;
397  #endif
398 
399  // Windows build
400  #ifdef _WIN64
401  HANDLE hProcess;
402  PROCESS_MEMORY_COUNTERS pmc;
403 
404  hProcess = OpenProcess(PROCESS_QUERY_INFORMATION | PROCESS_VM_READ,
405  FALSE,
406  GetCurrentProcessId());
407 
408  GetProcessMemoryInfo(hProcess, &pmc, sizeof(pmc));
409  CloseHandle(hProcess);
410 
411  return pmc.PeakWorkingSetSize >> 20;
412  #endif
413 }// end of GetHostMemoryUsageInMB
414 //--------------------------------------------------------------------------------------------------
415 
416 /**
417  * Get peak GPU memory usage.
418  *
419  * @return Peak memory usage in MBs.
420  */
422 {
423  size_t free, total;
424  cudaMemGetInfo(&free,&total);
425 
426  return ((total-free) >> 20);
427 }// end of GetDeviceMemoryUsageInMB
428 //--------------------------------------------------------------------------------------------------
429 
430 /**
431  * Get release code version.
432  *
433  * @return core name
434  */
436 {
437  return string(OUT_FMT_KWAVE_VERSION);
438 }// end of GetCodeName
439 //--------------------------------------------------------------------------------------------------
440 
441 /**
442  * Print full code name and the license
443  */
445 {
448  10,11,__DATE__,
449  8,8,__TIME__);
450 
451  if (parameters.GetGitHash() != "")
452  {
454  }
455 
457 
458  // OS detection
459  #ifdef __linux__
461  #elif __APPLE__
463  #elif _WIN32
465  #endif
466 
467  // Compiler detections
468  #if (defined(__GNUC__) || defined(__GNUG__)) && !(defined(__clang__) || defined(__INTEL_COMPILER))
470  #endif
471  #ifdef __INTEL_COMPILER
473  #endif
474  #ifdef _MSC_VER
476  #endif
477 
478  // instruction set
479  #if (defined (__AVX2__))
481  #elif (defined (__AVX__))
483  #elif (defined (__SSE4_2__))
485  #elif (defined (__SSE4_1__))
487  #elif (defined (__SSE3__))
489  #elif (defined (__SSE2__))
491  #endif
492 
494 
495  // CUDA detection
496  int cudaRuntimeVersion;
497  if (cudaRuntimeGetVersion(&cudaRuntimeVersion) != cudaSuccess)
498  {
500  }
501  else
502  {
505  cudaRuntimeVersion/1000, (cudaRuntimeVersion%100)/10);
506  }
507 
508  int cudaDriverVersion;
509  cudaDriverGetVersion(&cudaDriverVersion);
512  cudaDriverVersion/1000, (cudaDriverVersion%100)/10);
513 
514  const TCUDAParameters& cudaParameters = parameters.GetCudaParameters();
515  // no GPU was found
516  if (cudaParameters.GetDeviceIdx() == TCUDAParameters::DEFAULT_DEVICE_IDX)
517  {
519  }
520  else
521  {
527 
528  int paddingLength = 65 - ( 22 + strlen(cudaParameters.GetDeviceName().c_str()));
529 
532  cudaParameters.GetDeviceName().c_str(),
533  paddingLength,
535 
538  cudaParameters.GetDeviceProperties().major,
539  cudaParameters.GetDeviceProperties().minor);
540  }
541 
543 }// end of GetFullCodeAndLincence
544 //--------------------------------------------------------------------------------------------------
545 
546  /**
547  * Get total simulation time.
548  *
549  * @return total simulation time.
550  */
552 {
553  return totalTime.GetElapsedTime();
554 }// end of GetTotalTime()
555 //--------------------------------------------------------------------------------------------------
556 
557 /**
558  * Get pre-processing time.
559  *
560  * @return pre-processing time
561  */
563 {
565 }// end of GetPreProcessingTime
566 //--------------------------------------------------------------------------------------------------
567 
568 /**
569  * Get data load time.
570  *
571  * @return time to load data
572  */
574 {
575  return dataLoadTime.GetElapsedTime();
576 }// end of GetDataLoadTime
577 //--------------------------------------------------------------------------------------------------
578 
579 
580 /**
581  * Get simulation time (time loop).
582  *
583  * @return simulation time
584  */
586 {
588 }// end of GetSimulationTime
589 //--------------------------------------------------------------------------------------------------
590 
591 /**
592  * Get post-processing time.
593  *
594  * @return post-processing time
595  */
597 {
599 }// end of GetPostProcessingTime
600 //--------------------------------------------------------------------------------------------------
601 
602 /**
603  * Get total simulation time cumulated over all legs.
604  * @return simulation time cumulated over multiple legs
605  */
607 {
609 }// end of GetCumulatedTotalTime
610 //--------------------------------------------------------------------------------------------------
611 
612 /**
613  * Get pre-processing time cumulated over all legs.
614  *
615  * @return pre-processing time cumulated over multiple legs
616  */
618 {
620 } // end of GetCumulatedPreProcessingTime
621 //--------------------------------------------------------------------------------------------------
622 
623 /**
624  * Get data load time cumulated over all legs.
625  *
626  * @return time to load data cumulated over multiple legs
627  */
629 {
631 }// end of GetCumulatedDataLoadTime
632 //--------------------------------------------------------------------------------------------------
633 
634 /**
635  * Get simulation time (time loop) cumulated over all legs.
636  *
637  * @return simulation time cumulated over multiple legs
638  */
640 {
642 }// end of GetCumulatedSimulationTime
643 //--------------------------------------------------------------------------------------------------
644 
645 /**
646  * Get post-processing time cumulated over all legs.
647  *
648  * @return cumulated time to do post-processing over multiple legs
649  */
651 {
653 }// end of GetCumulatedPostProcessingTime
654 //--------------------------------------------------------------------------------------------------
655 
656 
657 //------------------------------------------------------------------------------------------------//
658 //-------------------------------------- Protected methods ---------------------------------------//
659 //------------------------------------------------------------------------------------------------//
660 
661 
662 /**
663  * Initialize FFT plans.
664  */
666 {
667  // create real to complex plans
669 
670  // create complex to real plans
672 
673  // if necessary, create 1D shift plans.
674  // in this case, the matrix has a bit bigger dimensions to be able to store shifted matrices.
676  {
677  // X shifts
680 
681  // Y shifts
684 
685  // Z shifts
688  }// end u_non_staggered
689 }// end of InitializeFFTPlans
690 //--------------------------------------------------------------------------------------------------
691 
692 /*
693  * Compute pre-processing phase. \n
694  * Initialize all indices, pre-compute constants such as c^2, rho0_sg* x dt and create kappa,
695  * absorb_eta, absorb_tau, absorb_nabla1, absorb_nabla2 matrices. Calculate this on the CPU side.
696  */
698 {
699  // get the correct sensor mask and recompute indices
700  if (parameters.Get_sensor_mask_type() == TParameters::INDEX)
701  {
703  }
704 
705  if (parameters.Get_sensor_mask_type() == TParameters::CORNERS)
706  {
708  }
709 
710  if ((parameters.Get_transducer_source_flag() != 0) ||
711  (parameters.Get_ux_source_flag() != 0) ||
712  (parameters.Get_uy_source_flag() != 0) ||
714  )
715  {
717  }
718 
720  {
722  }
723 
724  if (parameters.Get_p_source_flag() != 0)
725  {
727  }
728 
729  // compute dt / rho0_sg...
731  { // rho is scalar
735  }
736  else
737  { // non-uniform grid cannot be pre-calculated :-(
738  // rho is matrix
740  {
742  }
743  else
744  {
748  }
749  }
750 
751  // generate different matrices
752  if (parameters.Get_absorbing_flag() != 0)
753  {
756  }
757  else
758  {
759  GenerateKappa();
760  }
761 
762  // calculate c^2. It has to be after kappa gen... because of c modification
763  Compute_c2();
764 
765 }// end of PreProcessingPhase
766 //--------------------------------------------------------------------------------------------------
767 
768 /**
769  * Compute the main time loop of KSpaceFirstOrder3D.
770  */
772 {
773  actPercent = 0;
774 
775  // if resuming from a checkpoint,
776  // set ActPercent to correspond the t_index after recovery
777  if (parameters.Get_t_index() > 0)
778  {
779  // We're restarting after checkpoint
782  }
783 
784  // Progress header
786 
787  // Initial copy of data to the GPU
789 
791 
792  // execute main loop
794  {
795  const size_t t_index = parameters.Get_t_index();
796 
797  // compute velocity
798  ComputeVelocity();
799 
800  // add in the velocity u source term
802 
803  // add in the transducer source term (t = t1) to ux
804  if (parameters.Get_transducer_source_flag() > t_index)
805  {
808  Get_delay_mask(),
810  }
811 
812  // compute gradient of velocity
814 
815  // compute density
817  {
819  }
820  else
821  {
823  }
824 
825  // add in the source pressure term
827 
829  {
831  }
832  else
833  {
835  }
836 
837  //-- calculate initial pressure
838  if ((t_index == 0) && (parameters.Get_p0_source_flag() == 1)) Calculate_p0_source();
839 
840  StoreSensorData();
841  PrintStatistics();
842 
845  }
846 
847  // Since disk operations are one step delayed, we have to do the last one here.
848  // However we need to check if the loop wasn't skipped due to very short checkpoint interval
850  {
852  }
853 }// end of ComputeMainLoop()
854 //--------------------------------------------------------------------------------------------------
855 
856 /*
857  * Post processing, and closing the output streams.
858  */
860 {
862  {
863  Get_p().CopyFromDevice();
865  p_final_NAME,
867  }// p_final
868 
870  {
874 
884  }// u_final
885 
886  // Apply post-processing, flush data on disk/
889 
890  // store sensor mask if wanted
892  {
893  if (parameters.Get_sensor_mask_type() == TParameters::INDEX)
894  {
898  }
899 
900  if (parameters.Get_sensor_mask_type() == TParameters::CORNERS)
901  {
905  }
906  }
907 }// end of PostProcessing
908 //--------------------------------------------------------------------------------------------------
909 
910 /**
911  * Store sensor data.
912  * This routine exploits asynchronous behavior. It first performs IO from the i-1th step while
913  * waiting for ith step to come to the point of sampling.
914  */
916 {
917  // Unless the time for sampling has come, exit.
919  {
920 
921  // Read event for t_index-1. If sampling did not occur by then, ignored it.
922  // if it did store data on disk (flush) - the GPU is running asynchronously.
923  // But be careful, flush has to be one step delayed to work correctly.
924  // when restoring from checkpoint we have to skip the first flush
926  {
928  }
929 
930  // if --u_non_staggered is switched on, calculate unstaggered velocity.
932  {
934  }
935 
936  // Sample data for step t (store event for sampling in next turn)
938 
939  // the last step (or data after) checkpoint are flushed in the main loop
940  }
941 }// end of StoreSensorData
942 //--------------------------------------------------------------------------------------------------
943 
944 /**
945  * Write statistics and the header into the output file.
946  */
948 {
949  // write t_index into the output file
951  t_index_NAME,
953 
954  // Write scalars
956  THDF5_FileHeader& fileHeader = parameters.GetFileHeader();
957 
958  // Write File header
959  fileHeader.SetCodeName(GetCodeName());
960  fileHeader.SetMajorFileVersion();
961  fileHeader.SetMinorFileVersion();
962  fileHeader.SetActualCreationTime();
963  fileHeader.SetFileType(THDF5_FileHeader::OUTPUT);
964  fileHeader.SetHostName();
965 
967 
968  // Stop total timer here
969  totalTime.Stop();
970  fileHeader.SetExecutionTimes(GetTotalTime(),
971  GetDataLoadTime(),
975 
976  fileHeader.SetNumberOfCores();
977 
979 }// end of WriteOutputDataInfo
980 //--------------------------------------------------------------------------------------------------
981 
982 /**
983  * Save checkpoint data into the checkpoint file, flush aggregated outputs into the output file.
984  */
986 {
987  // Create Checkpoint file
988  THDF5_File& checkpointFile = parameters.GetCheckpointFile();
989  // if it happens and the file is opened (from the recovery, close it)
990  if (checkpointFile.IsOpen()) checkpointFile.Close();
991 
994 
995  // Create the new file (overwrite the old one)
996  checkpointFile.Create(parameters.GetCheckpointFileName());
997 
998  //-------------------------------------- Store Matrices ----------------------------------------//
999 
1000  // Store all necessary matrices in Checkpoint file
1002  // Write t_index
1003  checkpointFile.WriteScalarValue(checkpointFile.GetRootGroup(),
1004  t_index_NAME,
1006 
1007  // store basic dimension sizes (Nx, Ny, Nz) - Nt is not necessary
1008  checkpointFile.WriteScalarValue(checkpointFile.GetRootGroup(),
1009  Nx_NAME,
1011  checkpointFile.WriteScalarValue(checkpointFile.GetRootGroup(),
1012  Ny_NAME,
1014  checkpointFile.WriteScalarValue(checkpointFile.GetRootGroup(),
1015  Nz_NAME,
1017 
1018  // Write checkpoint file header
1019  THDF5_FileHeader checkpointFileHeader = parameters.GetFileHeader();
1020 
1021  checkpointFileHeader.SetFileType(THDF5_FileHeader::CHECKPOINT);
1022  checkpointFileHeader.SetCodeName(GetCodeName());
1023  checkpointFileHeader.SetActualCreationTime();
1024 
1025  checkpointFileHeader.WriteHeaderToCheckpointFile(checkpointFile);
1026 
1027  checkpointFile.Close();
1029 
1030  // checkpoint only if necessary (t_index > start_index), we're here one step ahead!
1032  {
1035 
1037 
1039  }
1040 
1042 }// end of SaveCheckpointData()
1043 //--------------------------------------------------------------------------------------------------
1044 
1045 
1046 /**
1047  * Compute new values of ux_sgx, uy_sgy, uz_sgz (acoustic velocity).
1048  */
1050 {
1052 
1054  Get_cufft_y_temp(),
1055  Get_cufft_z_temp(),
1056  Get_kappa(),
1060 
1064 
1066  { // scalars
1068  {
1070  Get_uy_sgy(),
1071  Get_uz_sgz(),
1075  Get_dxudxn_sgx(),
1076  Get_dyudyn_sgy(),
1077  Get_dzudzn_sgz(),
1078  Get_pml_x_sgx(),
1079  Get_pml_y_sgy(),
1080  Get_pml_z_sgz());
1081  }
1082  else
1083  {
1085  Get_uy_sgy(),
1086  Get_uz_sgz(),
1090  Get_pml_x_sgx(),
1091  Get_pml_y_sgy(),
1092  Get_pml_z_sgz());
1093  }
1094  }
1095  else
1096  {// matrices
1098  Get_uy_sgy(),
1099  Get_uz_sgz(),
1103  Get_dt_rho0_sgx(),
1104  Get_dt_rho0_sgy(),
1105  Get_dt_rho0_sgz(),
1106  Get_pml_x_sgx(),
1107  Get_pml_y_sgy(),
1108  Get_pml_z_sgz());
1109  }
1110 }// end of ComputeVelocity
1111 //--------------------------------------------------------------------------------------------------
1112 
1113 
1114 /**
1115  * Compute new values for duxdx, duydy, duzdz, gradient of velocity.
1116  */
1118 {
1122 
1123  // calculate duxyz on uniform grid
1125  Get_cufft_y_temp(),
1126  Get_cufft_z_temp(),
1127  Get_kappa(),
1131 
1135 
1136  // Non-uniform grid
1138  {
1140  Get_duydy(),
1141  Get_duzdz(),
1142  Get_dxudxn(),
1143  Get_dyudyn(),
1144  Get_dzudzn());
1145  }// nonlinear
1146 }// end of ComputeGradientVelocity()
1147 //--------------------------------------------------------------------------------------------------
1148 
1149 /**
1150  * Calculate new values of rhox, rhoy and rhoz for non-linear case.
1151  */
1153 {
1154  // Scalar
1156  {
1158  Get_rhoy(),
1159  Get_rhoz(),
1160  Get_pml_x(),
1161  Get_pml_y(),
1162  Get_pml_z(),
1163  Get_duxdx(),
1164  Get_duydy(),
1165  Get_duzdz());
1166 }
1167 else
1168 {
1169  // rho0 is a matrix
1171  Get_rhoy(),
1172  Get_rhoz(),
1173  Get_pml_x(),
1174  Get_pml_y(),
1175  Get_pml_z(),
1176  Get_duxdx(),
1177  Get_duydy(),
1178  Get_duzdz(),
1179  Get_rho0());
1180  } // end matrix
1181 }// end of ComputeDensityNonliner
1182 //--------------------------------------------------------------------------------------------------
1183 
1184 /**
1185  * Calculate new values of rhox, rhoy and rhoz for linear case.
1186  *
1187  */
1189 {
1190  // Scalar
1192  {
1194  Get_rhoy(),
1195  Get_rhoz(),
1196  Get_pml_x(),
1197  Get_pml_y(),
1198  Get_pml_z(),
1199  Get_duxdx(),
1200  Get_duydy(),
1201  Get_duzdz());
1202  }
1203  else
1204  {
1205  // rho0 is a matrix
1207  Get_rhoy(),
1208  Get_rhoz(),
1209  Get_pml_x(),
1210  Get_pml_y(),
1211  Get_pml_z(),
1212  Get_duxdx(),
1213  Get_duydy(),
1214  Get_duzdz(),
1215  Get_rho0());
1216  } // end matrix
1217 }// end of ComputeDensityLinear
1218 //--------------------------------------------------------------------------------------------------
1219 
1220 
1221 /**
1222  * Compute acoustic pressure for non-linear case.
1223  */
1225 {
1227  { // absorbing case
1228  TRealMatrix& rhoxyz_part = Get_temp_1_real_3D();
1229  TRealMatrix& BonA_part = Get_temp_2_real_3D();
1230  TRealMatrix& du_part = Get_temp_3_real_3D();
1231 
1232  TRealMatrix& absorb_tau_temp = du_part;
1233  TRealMatrix& absorb_eta_temp = rhoxyz_part;
1234 
1235  ComputePressurePartsNonLinear(rhoxyz_part, BonA_part, du_part);
1236 
1238  Get_cufft_y_temp().Compute_FFT_3D_R2C(rhoxyz_part);
1239 
1241  Get_cufft_y_temp(),
1243  Get_absorb_nabla2());
1244 
1245  Get_cufft_x_temp().Compute_FFT_3D_C2R(absorb_tau_temp);
1246  Get_cufft_y_temp().Compute_FFT_3D_C2R(absorb_eta_temp);
1247 
1248  SumPressureTermsNonlinear(absorb_tau_temp, absorb_eta_temp, BonA_part);
1249  }
1250  else
1251  {
1253  }
1254 }// end of ComputePressureNonlinear
1255 //--------------------------------------------------------------------------------------------------
1256 
1257 /*
1258  * Compute new p for linear case.
1259  */
1261 {
1263  { // absorbing case
1264  TRealMatrix& sum_rhoxyz = Get_temp_1_real_3D();
1265  TRealMatrix& sum_rho0_du = Get_temp_2_real_3D();
1266 
1267  TRealMatrix& absorb_tau_temp = Get_temp_2_real_3D();
1268  TRealMatrix& absorb_eta_temp = Get_temp_3_real_3D();
1269 
1270  ComputePressurePartsLinear(sum_rhoxyz, sum_rho0_du);
1271 
1272  // ifftn ( absorb_nabla1 * fftn (rho0 * (duxdx+duydy+duzdz))
1273  Get_cufft_x_temp().Compute_FFT_3D_R2C(sum_rho0_du);
1274  Get_cufft_y_temp().Compute_FFT_3D_R2C(sum_rhoxyz);
1275 
1277  Get_cufft_y_temp(),
1279  Get_absorb_nabla2());
1280 
1281  Get_cufft_x_temp().Compute_FFT_3D_C2R(absorb_tau_temp);
1282  Get_cufft_y_temp().Compute_FFT_3D_C2R(absorb_eta_temp);
1283 
1284  SumPressureTermsLinear(absorb_tau_temp, absorb_eta_temp, sum_rhoxyz);
1285  }
1286  else
1287  {
1288  // lossless case
1290  }
1291 }// end of ComputePressureLinear()
1292 //--------------------------------------------------------------------------------------------------
1293 
1294 /**
1295  * Add velocity source to the particle velocity.
1296  */
1298 {
1299  size_t t_index = parameters.Get_t_index();
1300 
1301  if (parameters.Get_ux_source_flag() > t_index)
1302  {
1306  t_index);
1307  }
1308  if (parameters.Get_uy_source_flag() > t_index)
1309  {
1313  t_index);
1314  }
1315  if (parameters.Get_uz_source_flag() > t_index)
1316  {
1320  t_index);
1321  }
1322 }// end of AddVelocitySource
1323 //--------------------------------------------------------------------------------------------------
1324 
1325 /*
1326  * Add in pressure source.
1327  */
1329 {
1330  size_t t_index = parameters.Get_t_index();
1331 
1332  if (parameters.Get_p_source_flag() > t_index)
1333  {
1335  Get_rhoy(),
1336  Get_rhoz(),
1339  t_index);
1340  }//
1341 }// end of AddPressureSource
1342 //--------------------------------------------------------------------------------------------------
1343 
1344 /**
1345  * Calculate p0 source when necessary.
1346  */
1348 {
1349  // get over the scalar problem
1350  bool is_c2_scalar = parameters.Get_c0_scalar_flag();
1351  const float* c2 = (is_c2_scalar) ? nullptr : Get_c2().GetDeviceData();
1352 
1353  //-- add the initial pressure to rho as a mass source --//
1355  Get_rhox(),
1356  Get_rhoy(),
1357  Get_rhoz(),
1359  is_c2_scalar,
1360  c2);
1361 
1362  //-----------------------------------------------------------------------//
1363  //--compute u(t = t1 + dt/2) based on the assumption u(dt/2) = -u(-dt/2)-//
1364  //-- which forces u(t = t1) = 0 -//
1365  //-----------------------------------------------------------------------//
1367 
1369  Get_cufft_y_temp(),
1370  Get_cufft_z_temp(),
1371  Get_kappa(),
1375 
1379 
1381  {
1383  { // non uniform grid, homogeneous
1385  Get_uy_sgy(),
1386  Get_uz_sgz(),
1387  Get_dxudxn_sgx(),
1388  Get_dyudyn_sgy(),
1389  Get_dzudzn_sgz());
1390  }
1391  else
1392  { //uniform grid, homogeneous
1394  }
1395  }
1396  else
1397  {
1398  // heterogeneous, uniform grid
1399  // divide the matrix by 2 and multiply with st./rho0_sg
1401  Get_uy_sgy(),
1402  Get_uz_sgz(),
1403  Get_dt_rho0_sgx(),
1404  Get_dt_rho0_sgy(),
1405  Get_dt_rho0_sgz());
1406  }
1407 }// end of Calculate_p0_source
1408 //--------------------------------------------------------------------------------------------------
1409 
1410 
1411 /**
1412  * Generate kappa matrix for non-absorbing mode.
1413  */
1415 {
1416  #pragma omp parallel
1417  {
1418  const float dx_sq_rec = 1.0f / (parameters.Get_dx() * parameters.Get_dx());
1419  const float dy_sq_rec = 1.0f / (parameters.Get_dy() * parameters.Get_dy());
1420  const float dz_sq_rec = 1.0f / (parameters.Get_dz() * parameters.Get_dz());
1421 
1422  const float c_ref_dt_pi = parameters.Get_c_ref() * parameters.Get_dt() * float(M_PI);
1423 
1424  const float nx_rec = 1.0f / static_cast<float>(parameters.GetFullDimensionSizes().nx);
1425  const float ny_rec = 1.0f / static_cast<float>(parameters.GetFullDimensionSizes().ny);
1426  const float nz_rec = 1.0f / static_cast<float>(parameters.GetFullDimensionSizes().nz);
1427 
1428  const size_t nx = parameters.GetReducedDimensionSizes().nx;
1429  const size_t ny = parameters.GetReducedDimensionSizes().ny;
1430  const size_t nz = parameters.GetReducedDimensionSizes().nz;
1431 
1432  float* kappa = Get_kappa().GetHostData();
1433 
1434  #pragma omp for schedule (static)
1435  for (size_t z = 0; z < nz; z++)
1436  {
1437  const float z_f = (float) z;
1438  float z_part = 0.5f - fabs(0.5f - z_f * nz_rec );
1439  z_part = (z_part * z_part) * dz_sq_rec;
1440 
1441  for (size_t y = 0; y < ny; y++)
1442  {
1443  const float y_f = (float) y;
1444  float y_part = 0.5f - fabs(0.5f - y_f * ny_rec);
1445  y_part = (y_part * y_part) * dy_sq_rec;
1446 
1447  const float yz_part = z_part + y_part;
1448  for (size_t x = 0; x < nx; x++)
1449  {
1450  const float x_f = (float) x;
1451  float x_part = 0.5f - fabs(0.5f - x_f * nx_rec);
1452  x_part = (x_part * x_part) * dx_sq_rec;
1453 
1454  float k = c_ref_dt_pi * sqrt(x_part + yz_part);
1455 
1456  // kappa element
1457  kappa[(z*ny + y) * nx + x ] = (k == 0.0f) ? 1.0f : sin(k)/k;
1458  }//x
1459  }//y
1460  }// z
1461  }// parallel
1462 }// end of Generate_kappa
1463 //-------------------------------------------------------------------------------------------------
1464 
1465 /*
1466  * Generate kappa, absorb_nabla1, absorb_nabla2 for absorbing media.
1467  */
1469 {
1470  #pragma omp parallel
1471  {
1472  const float dx_sq_rec = 1.0f / (parameters.Get_dx() * parameters.Get_dx());
1473  const float dy_sq_rec = 1.0f / (parameters.Get_dy() * parameters.Get_dy());
1474  const float dz_sq_rec = 1.0f / (parameters.Get_dz() * parameters.Get_dz());
1475 
1476  const float c_ref_dt_2 = parameters.Get_c_ref() * parameters.Get_dt() * 0.5f;
1477  const float pi_2 = float(M_PI) * 2.0f;
1478 
1479  const size_t nx = parameters.GetFullDimensionSizes().nx;
1480  const size_t ny = parameters.GetFullDimensionSizes().ny;
1481  const size_t nz = parameters.GetFullDimensionSizes().nz;
1482 
1483  const float nx_rec = 1.0f / (float) nx;
1484  const float ny_rec = 1.0f / (float) ny;
1485  const float nz_rec = 1.0f / (float) nz;
1486 
1487  const size_t nxComplex = parameters.GetReducedDimensionSizes().nx;
1488  const size_t nyComplex = parameters.GetReducedDimensionSizes().ny;
1489  const size_t nzComplex = parameters.GetReducedDimensionSizes().nz;
1490 
1491  float* kappa = Get_kappa().GetHostData();
1492  float* absorb_nabla1 = Get_absorb_nabla1().GetHostData();
1493  float* absorb_nabla2 = Get_absorb_nabla2().GetHostData();
1494  const float alpha_power = parameters.Get_alpha_power();
1495 
1496  #pragma omp for schedule (static)
1497  for (size_t z = 0; z < nzComplex; z++)
1498  {
1499  const float z_f = (float) z;
1500  float z_part = 0.5f - fabs(0.5f - z_f * nz_rec );
1501  z_part = (z_part * z_part) * dz_sq_rec;
1502 
1503  for (size_t y = 0; y < nyComplex; y++)
1504  {
1505  const float y_f = (float) y;
1506  float y_part = 0.5f - fabs(0.5f - y_f * ny_rec);
1507  y_part = (y_part * y_part) * dy_sq_rec;
1508 
1509  const float yz_part = z_part + y_part;
1510 
1511  size_t i = (z * nyComplex + y) * nxComplex;
1512 
1513  for (size_t x = 0; x < nxComplex; x++)
1514  {
1515  const float x_f = (float) x;
1516  float x_part = 0.5f - fabs(0.5f - x_f * nx_rec);
1517  x_part = (x_part * x_part) * dx_sq_rec;
1518 
1519  float k = pi_2 * sqrt(x_part + yz_part);
1520  float c_ref_k = c_ref_dt_2 * k;
1521 
1522  absorb_nabla1[i] = pow(k, alpha_power - 2);
1523  absorb_nabla2[i] = pow(k, alpha_power - 1);
1524 
1525  kappa[i] = (c_ref_k == 0.0f) ? 1.0f : sin(c_ref_k)/c_ref_k;
1526 
1527  if (absorb_nabla1[i] == std::numeric_limits<float>::infinity()) absorb_nabla1[i] = 0.0f;
1528  if (absorb_nabla2[i] == std::numeric_limits<float>::infinity()) absorb_nabla2[i] = 0.0f;
1529 
1530  i++;
1531  }//x
1532  }//y
1533  }// z
1534  }// parallel
1535 }// end of GenerateKappaAndNablas
1536 //--------------------------------------------------------------------------------------------------
1537 
1538 /**
1539  * Generate absorb_tau and absorb_eta in for heterogenous media.
1540  */
1542 {
1543  // test for scalars
1545  {
1546  const float alpha_power = parameters.Get_alpha_power();
1547  const float tan_pi_y_2 = tan(static_cast<float> (M_PI_2) * alpha_power);
1548  const float alpha_db_neper_coeff = (100.0f * pow(1.0e-6f /(2.0f * static_cast<float>(M_PI)), alpha_power)) /
1549  (20.0f * static_cast<float>(M_LOG10E));
1550 
1551  const float alpha_coeff_2 = 2.0f * parameters.Get_alpha_coeff_scalar() * alpha_db_neper_coeff;
1552 
1553  parameters.Get_absorb_tau_scalar() = (-alpha_coeff_2) * pow(parameters.Get_c0_scalar(), alpha_power - 1);
1554  parameters.Get_absorb_eta_scalar() = alpha_coeff_2 * pow(parameters.Get_c0_scalar(), alpha_power) * tan_pi_y_2;
1555  }
1556  else
1557  { // matrix
1558  #pragma omp parallel
1559  {
1560  const size_t nx = parameters.GetFullDimensionSizes().nx;
1561  const size_t ny = parameters.GetFullDimensionSizes().ny;
1562  const size_t nz = parameters.GetFullDimensionSizes().nz;
1563 
1564  float* absorb_tau = Get_absorb_tau().GetHostData();
1565  float* absorb_eta = Get_absorb_eta().GetHostData();
1566 
1567  float* alpha_coeff;
1568  size_t alpha_shift;
1569 
1571  {
1572  alpha_coeff = &(parameters.Get_alpha_coeff_scalar());
1573  alpha_shift = 0;
1574  }
1575  else
1576  {
1577  alpha_coeff = Get_temp_1_real_3D().GetHostData();
1578  alpha_shift = 1;
1579  }
1580 
1581  float* c0;
1582  size_t c0_shift;
1584  {
1585  c0 = &(parameters.Get_c0_scalar());
1586  c0_shift = 0;
1587  }
1588  else
1589  {
1590  c0 = Get_c2().GetHostData();
1591  c0_shift = 1;
1592  }
1593 
1594  const float alpha_power = parameters.Get_alpha_power();
1595  const float tan_pi_y_2 = tan(static_cast<float>(M_PI_2) * alpha_power);
1596 
1597  //alpha = 100*alpha.*(1e-6/(2*pi)).^y./
1598  // (20*log10(exp(1)));
1599  const float alpha_db_neper_coeff = (100.0f * pow(1.0e-6f / (2.0f * static_cast<float>(M_PI)), alpha_power)) /
1600  (20.0f * static_cast<float>(M_LOG10E));
1601 
1602 
1603  #pragma omp for schedule (static)
1604  for (size_t z = 0; z < nz; z++)
1605  {
1606  for (size_t y = 0; y < ny; y++)
1607  {
1608  size_t i = (z * ny + y) * nx;
1609  for (size_t x = 0; x < nx; x++)
1610  {
1611  const float alpha_coeff_2 = 2.0f * alpha_coeff[i*alpha_shift] * alpha_db_neper_coeff;
1612 
1613  absorb_tau[i] = (-alpha_coeff_2) * pow(c0[i * c0_shift],alpha_power - 1);
1614  absorb_eta[i] = alpha_coeff_2 * pow(c0[i * c0_shift],alpha_power) * tan_pi_y_2;
1615 
1616  i++;
1617  }//x
1618  }//y
1619  }// z
1620  }// parallel
1621  } // absorb_tau and aborb_eta = matrics
1622 }// end of GenerateTauAndEta
1623 //--------------------------------------------------------------------------------------------------
1624 
1625 /**
1626  * Prepare dt./ rho0 for non-uniform grid.
1627  *
1628  */
1630 {
1631  #pragma omp parallel
1632  {
1633  float* dt_rho0_sgx = Get_dt_rho0_sgx().GetHostData();
1634  float* dt_rho0_sgy = Get_dt_rho0_sgy().GetHostData();
1635  float* dt_rho0_sgz = Get_dt_rho0_sgz().GetHostData();
1636 
1637  const float dt = parameters.Get_dt();
1638 
1639  const float* duxdxn_sgx = Get_dxudxn_sgx().GetHostData();
1640  const float* duydyn_sgy = Get_dyudyn_sgy().GetHostData();
1641  const float* duzdzn_sgz = Get_dzudzn_sgz().GetHostData();
1642 
1643  const size_t nz = Get_dt_rho0_sgx().GetDimensionSizes().nz;
1644  const size_t ny = Get_dt_rho0_sgx().GetDimensionSizes().ny;
1645  const size_t nx = Get_dt_rho0_sgx().GetDimensionSizes().nx;
1646 
1647  const size_t sliceSize = (nx * ny );
1648 
1649  #pragma omp for schedule (static)
1650  for (size_t z = 0; z < nz; z++)
1651  {
1652  register size_t i = z* sliceSize;
1653  for (size_t y = 0; y < ny; y++)
1654  {
1655  for (size_t x = 0; x < nx; x++)
1656  {
1657  dt_rho0_sgx[i] = (dt * duxdxn_sgx[x]) / dt_rho0_sgx[i];
1658  i++;
1659  } // x
1660  } // y
1661  } // z
1662 
1663  #pragma omp for schedule (static)
1664  for (size_t z = 0; z < nz; z++)
1665  {
1666  register size_t i = z* sliceSize;
1667  for (size_t y = 0; y < ny; y++)
1668  {
1669  const float duydyn_el = duydyn_sgy[y];
1670  for (size_t x = 0; x < nx; x++)
1671  {
1672  dt_rho0_sgy[i] = (dt * duydyn_el) / dt_rho0_sgy[i];
1673  i++;
1674  } // x
1675  } // y
1676  } // z
1677 
1678 
1679  #pragma omp for schedule (static)
1680  for (size_t z = 0; z < nz; z++)
1681  {
1682  register size_t i = z* sliceSize;
1683  const float duzdzn_el = duzdzn_sgz[z];
1684  for (size_t y = 0; y < ny; y++)
1685  {
1686  for (size_t x = 0; x < nx; x++)
1687  {
1688  dt_rho0_sgz[i] = (dt * duzdzn_el) / dt_rho0_sgz[i];
1689  i++;
1690  } // x
1691  } // y
1692  } // z
1693  } // parallel
1694 }// end of GenerateInitialDenisty
1695 //--------------------------------------------------------------------------------------------------
1696 
1697 
1698 /**
1699  * Compute c^2 on the CPU side.
1700  */
1702 {
1704  { // scalar
1705  float c = parameters.Get_c0_scalar();
1706  parameters.Get_c0_scalar() = c * c;
1707  }
1708  else
1709  { // matrix
1710  float* c2 = Get_c2().GetHostData();
1711 
1712  #pragma omp parallel for schedule (static)
1713  for (size_t i=0; i < Get_c2().GetElementCount(); i++)
1714  {
1715  c2[i] = c2[i] * c2[i];
1716  }
1717  }// matrix
1718 }// ComputeC2
1719 //--------------------------------------------------------------------------------------------------
1720 
1721 /**
1722  * Compute three temporary sums in the new pressure formula, non-linear absorbing case.
1723  *
1724  * @param [out] rho_part - rhox_sgx + rhoy_sgy + rhoz_sgz
1725  * @param [out] BonA_part - BonA + rho ^2 / 2 rho0 + (rhox_sgx + rhoy_sgy + rhoz_sgz)
1726  * @param [out] du_part - rho0* (duxdx + duydy + duzdz)
1727  *
1728  */
1730  TRealMatrix& BonA_part,
1731  TRealMatrix& du_part)
1732 {
1733  const bool is_BonA_scalar = parameters.Get_BonA_scalar_flag();
1734  const bool is_rho0_scalar = parameters.Get_rho0_scalar_flag();
1735 
1736  const float* BonA_data = (is_BonA_scalar) ? nullptr : Get_BonA().GetDeviceData();
1737  const float* rho0_data = (is_rho0_scalar) ? nullptr : Get_rho0().GetDeviceData();
1738 
1740  BonA_part,
1741  du_part,
1742  Get_rhox(),
1743  Get_rhoy(),
1744  Get_rhoz(),
1745  Get_duxdx(),
1746  Get_duydy(),
1747  Get_duzdz(),
1748  is_BonA_scalar,
1749  BonA_data,
1750  is_rho0_scalar,
1751  rho0_data);
1752 
1753 }// end of ComputePressurePartsNonLinear
1754 //--------------------------------------------------------------------------------------------------
1755 
1756 
1757 /**
1758  * Calculate two temporary sums in the new pressure formula, linear absorbing case.
1759  *
1760  * @param [out] rhoxyz_sum -rhox_sgx + rhoy_sgy + rhoz_sgz
1761  * @param [out] rho0_du_sum - rho0* (duxdx + duydy + duzdz);
1762  */
1764  TRealMatrix& rho0_du_sum)
1765 {
1766  const bool is_rho0_scalar = parameters.Get_rho0_scalar();
1767  const float* rho0_matrix = (is_rho0_scalar) ? nullptr : Get_rho0().GetDeviceData();
1768 
1770  rho0_du_sum,
1771  Get_rhox(),
1772  Get_rhoy(),
1773  Get_rhoz(),
1774  Get_duxdx(),
1775  Get_duydy(),
1776  Get_duzdz(),
1777  is_rho0_scalar,
1778  rho0_matrix);
1779 }// end of ComputePressurePartsLinear
1780 //--------------------------------------------------------------------------------------------------
1781 
1782 /**
1783  * Sum sub-terms to calculate new pressure, non-linear case.
1784 
1785  * @param [in] absorb_tau_temp - tau component
1786  * @param [in] absorb_eta_temp - BonA + rho ^2 / 2 rho0 +
1787  * (rhox_sgx + rhoy_sgy + rhoz_sgz)
1788  * @param [in] BonA_temp - rho0* (duxdx + duydy + duzdz)
1789  */
1791  TRealMatrix& absorb_eta_temp,
1792  TRealMatrix& BonA_temp)
1793 {
1794  const bool is_c2_scalar = parameters.Get_c0_scalar_flag();
1795  const bool is_tau_eta_scalar = parameters.Get_c0_scalar_flag() &&
1797 
1798  const float* c2_data_matrix = (is_c2_scalar) ? nullptr : Get_c2().GetDeviceData();
1799  const float* tau_data_matrix = (is_tau_eta_scalar) ? nullptr : Get_absorb_tau().GetDeviceData();
1800  const float* eta_data_matrix = (is_tau_eta_scalar) ? nullptr : Get_absorb_eta().GetDeviceData();
1801 
1802  const float* Absorb_tau_data = absorb_tau_temp.GetDeviceData();
1803  const float* Absorb_eta_data = absorb_eta_temp.GetDeviceData();
1804 
1806  BonA_temp,
1807  is_c2_scalar,
1808  c2_data_matrix,
1809  is_tau_eta_scalar,
1810  Absorb_tau_data,
1811  tau_data_matrix,
1812  Absorb_eta_data,
1813  eta_data_matrix);
1814 }// end of SumPressureTermsNonlinear
1815 //--------------------------------------------------------------------------------------------------
1816 
1817 /**
1818  * Sum sub-terms to calculate new pressure, linear case.
1819  *
1820  * @param [in] absorb_tau_temp - sub-term with absorb_tau
1821  * @param [in] absorb_eta_temp - sub-term with absorb_eta
1822  * @param [in] rhoxyz_sum - rhox_sgx + rhoy_sgy + rhoz_sgz
1823  */
1825  TRealMatrix& absorb_eta_temp,
1826  TRealMatrix& rhoxyz_sum)
1827 {
1828  const bool is_c2_scalar = parameters.Get_c0_scalar_flag();
1829  const bool is_tau_eta_scalar = parameters.Get_c0_scalar_flag() &&
1831 
1832  const float* c2_data_matrix = (is_c2_scalar) ? nullptr : Get_c2().GetDeviceData();
1833  const float* tau_data_matrix = (is_tau_eta_scalar) ? nullptr : Get_absorb_tau().GetDeviceData();
1834  const float* eta_data_matrix = (is_tau_eta_scalar) ? nullptr : Get_absorb_eta().GetDeviceData();
1835 
1837  absorb_tau_temp,
1838  absorb_eta_temp,
1839  rhoxyz_sum,
1840  is_c2_scalar,
1841  c2_data_matrix,
1842  is_tau_eta_scalar,
1843  tau_data_matrix,
1844  eta_data_matrix);
1845 }// end of SumPressureTermsLinear
1846 //--------------------------------------------------------------------------------------------------
1847 
1848 /**
1849  * Sum sub-terms for new p, non-linear lossless case.
1850  */
1852 {
1853  const bool is_c2_scalar = parameters.Get_c0_scalar_flag();
1854  const bool is_BonA_scalar = parameters.Get_BonA_scalar_flag();
1855  const bool is_rho0_scalar = parameters.Get_rho0_scalar_flag();
1856 
1857  const float* c2_data_matrix = (is_c2_scalar) ? nullptr : Get_c2().GetDeviceData();
1858  const float* BonA_data_matrix = (is_BonA_scalar) ? nullptr : Get_BonA().GetDeviceData();
1859  const float* rho0_data_matrix = (is_rho0_scalar) ? nullptr : Get_rho0().GetDeviceData();
1860 
1862  Get_rhox(),
1863  Get_rhoy(),
1864  Get_rhoz(),
1865  is_c2_scalar,
1866  c2_data_matrix,
1867  is_BonA_scalar,
1868  BonA_data_matrix,
1869  is_rho0_scalar,
1870  rho0_data_matrix);
1871 
1872 }// end of SumPressureNonlinearLossless
1873 //--------------------------------------------------------------------------------------------------
1874 
1875 /**
1876  * Sum sub-terms for new p, linear lossless case.
1877  */
1879 {
1880  const float is_c2_scalar = parameters.Get_c0_scalar();
1881  const float* c2_matrix = (is_c2_scalar) ? nullptr : Get_c2().GetDeviceData();
1882 
1884  Get_rhox(),
1885  Get_rhoy(),
1886  Get_rhoz(),
1887  is_c2_scalar,
1888  c2_matrix);
1889 
1890 }// end of SumPressureLinearLossless
1891 //--------------------------------------------------------------------------------------------------
1892 
1893 
1894 /**
1895  * Calculated shifted velocities.
1896  * \n
1897  * ux_shifted = real(ifft(bsxfun(\@times, x_shift_neg, fft(ux_sgx, [], 1)), [], 1)); \n
1898  * uy_shifted = real(ifft(bsxfun(\@times, y_shift_neg, fft(uy_sgy, [], 2)), [], 2)); \n
1899  * uz_shifted = real(ifft(bsxfun(\@times, z_shift_neg, fft(uz_sgz, [], 3)), [], 3)); \n
1900  */
1901 
1903 {
1904 
1905  // ux_shifted
1909 
1910  // uy_shifted
1914 
1915  // uz_shifted
1919 
1920 }// end of CalculateShiftedVelocity
1921 //--------------------------------------------------------------------------------------------------
1922 
1923 /**
1924  * Print progress statistics.
1925  */
1927 {
1928  const float nt = (float) parameters.Get_nt();
1929  const size_t t_index = parameters.Get_t_index();
1930 
1931  if (t_index > (actPercent * nt * 0.01f) )
1932  {
1934 
1935  iterationTime.Stop();
1936 
1937  const double elTime = iterationTime.GetElapsedTime();
1938  const double elTimeWithLegs = iterationTime.GetElapsedTime() +
1940  const double toGo = ((elTimeWithLegs / (float) (t_index + 1)) * nt) - elTimeWithLegs;
1941 
1942  struct tm* current;
1943  time_t now;
1944  time(&now);
1945  now += toGo;
1946  current = localtime(&now);
1947 
1950  (size_t) ((t_index) / (nt * 0.01f)),'%',
1951  elTime, toGo,
1952  current->tm_mday, current->tm_mon+1, current->tm_year-100,
1953  current->tm_hour, current->tm_min, current->tm_sec);
1955  }
1956 }// end of PrintStatistics
1957 //--------------------------------------------------------------------------------------------------
1958 
1959 /**
1960  * Is time to checkpoint?
1961  *
1962  * @return true if it is necessary to stop to checkpoint
1963  */
1965 {
1966  if (!parameters.IsCheckpointEnabled()) return false;
1967 
1968  totalTime.Stop();
1969 
1970  return (totalTime.GetElapsedTime() > static_cast<float>(parameters.GetCheckpointInterval()));
1971 }// end of IsTimeToCheckpoint
1972 //--------------------------------------------------------------------------------------------------
1973 
1974 
1975 /**
1976  * Was the loop interrupted to checkpoint?
1977  *
1978  * @return true if it is time to checkpoint
1979  */
1981 {
1982  return (parameters.Get_t_index() != parameters.Get_nt());
1983 }// end of IsCheckpointInterruption
1984 //--------------------------------------------------------------------------------------------------
1985 
1986 /**
1987  * Check the output file has the correct format and version.
1988  *
1989  * @throw ios::failure if an error happens
1990  */
1992 {
1993  // The header has already been read
1994  THDF5_FileHeader& fileHeader = parameters.GetFileHeader();
1995  THDF5_File& outputFile = parameters.GetOutputFile();
1996 
1997  // test file type
1998  if (fileHeader.GetFileType() != THDF5_FileHeader::OUTPUT)
1999  {
2001  parameters.GetOutputFileName().c_str()));
2002  }
2003 
2004  // test file major version
2005  if (!fileHeader.CheckMajorFileVersion())
2006  {
2009  fileHeader.GetCurrentHDF5_MajorVersion().c_str()));
2010  }
2011 
2012  // test file minor version
2013  if (!fileHeader.CheckMinorFileVersion())
2014  {
2017  fileHeader.GetCurrentHDF5_MinorVersion().c_str()));
2018  }
2019 
2020  // Check dimension sizes
2021  TDimensionSizes outputDimSizes;
2022  outputFile.ReadScalarValue(outputFile.GetRootGroup(), Nx_NAME, outputDimSizes.nx);
2023  outputFile.ReadScalarValue(outputFile.GetRootGroup(), Ny_NAME, outputDimSizes.ny);
2024  outputFile.ReadScalarValue(outputFile.GetRootGroup(), Nz_NAME, outputDimSizes.nz);
2025 
2026  if (parameters.GetFullDimensionSizes() != outputDimSizes)
2027  {
2029  outputDimSizes.nx,
2030  outputDimSizes.ny,
2031  outputDimSizes.nz,
2035  }
2036 }// end of CheckOutputFile
2037 //--------------------------------------------------------------------------------------------------
2038 
2039 
2040 /**
2041  * Check the file type and the version of the checkpoint file.
2042  *
2043  * @throw ios::failure if an error happens
2044  *
2045  */
2047 {
2048  // read the header and check the file version
2049  THDF5_FileHeader fileHeader;
2050  THDF5_File& checkpointFile = parameters.GetCheckpointFile();
2051 
2052  fileHeader.ReadHeaderFromCheckpointFile(checkpointFile);
2053 
2054  // test file type
2055  if (fileHeader.GetFileType() != THDF5_FileHeader::CHECKPOINT)
2056  {
2058  parameters.GetCheckpointFileName().c_str()));
2059  }
2060 
2061  // test file major version
2062  if (!fileHeader.CheckMajorFileVersion())
2063  {
2066  fileHeader.GetCurrentHDF5_MajorVersion().c_str()));
2067  }
2068 
2069  // test file minor version
2070  if (!fileHeader.CheckMinorFileVersion())
2071  {
2074  fileHeader.GetCurrentHDF5_MinorVersion().c_str()));
2075  }
2076 
2077  // Check dimension sizes
2078  TDimensionSizes checkpointDimSizes;
2079  checkpointFile.ReadScalarValue(checkpointFile.GetRootGroup(), Nx_NAME, checkpointDimSizes.nx);
2080  checkpointFile.ReadScalarValue(checkpointFile.GetRootGroup(), Ny_NAME, checkpointDimSizes.ny);
2081  checkpointFile.ReadScalarValue(checkpointFile.GetRootGroup(), Nz_NAME, checkpointDimSizes.nz);
2082 
2083  if (parameters.GetFullDimensionSizes() != checkpointDimSizes)
2084  {
2086  checkpointDimSizes.nx,
2087  checkpointDimSizes.ny,
2088  checkpointDimSizes.nz,
2092  }
2093 }// end of CheckCheckpointFile
2094 //--------------------------------------------------------------------------------------------------
2095 
2096 /**
2097  * Restore cumulated elapsed time from the output file. Open the header, read this and store
2098  * into TMPI_Time classes
2099  *
2100  * @param [in] outputFile - Output file
2101  */
2103 {
2105 
2106  // Get execution times stored in the output file header
2108  dataLoadTime,
2109  preProcessingTime,
2110  simulationTime,
2111  postProcessingTime);
2112 
2113  this->totalTime.SetCumulatedElapsedTimeOverPreviousLegs(totalTime);
2114  this->dataLoadTime.SetCumulatedElapsedTimeOverPreviousLegs(dataLoadTime);
2115  this->preProcessingTime.SetCumulatedElapsedTimeOverPreviousLegs(preProcessingTime);
2116  this->simulationTime.SetCumulatedElapsedTimeOverPreviousLegs(simulationTime);
2117  this->postProcessingTime.SetCumulatedElapsedTimeOverPreviousLegs(postProcessingTime);
2118 
2119 }// end of LoadElapsedTimeFromOutputFile
2120 //--------------------------------------------------------------------------------------------------
void Compute_FFT_1DY_R2C(TRealMatrix &inMatrix)
Compute 1D out-of-place Real-to-Complex FFT in the Y dimension.
Class for HDF5 header.
Definition: HDF5_File.h:730
size_t GetDeviceMemoryUsageInMB()
Get memory usage in MB on the GPU side.
size_t nx
number of elements in the x direction
TOutputMessage OUT_FMT_MEMORY_ALLOCATION
Output message.
TDimensionSizes GetReducedDimensionSizes() const
Reduced dimension sizes of the simulation (complex classes).
Definition: Parameters.h:97
TRealMatrix & Get_absorb_nabla1()
Get the absorb_nabla1 matrix from the container.
TOutputMessage OUT_FMT_LICENCE
Print version output message.
std::string GetCheckpointFileName() const
Get checkpoint filename.
Definition: Parameters.h:215
float & Get_c0_scalar()
Get c0_scalar value.
Definition: Parameters.h:187
TRealMatrix & Get_dt_rho0_sgy()
Get the dt.*rho0_sgy matrix from the container.
void ComputePressurePartsLinear(TRealMatrix &rhoxyz_sum, TRealMatrix &rho0_du_sum)
Calculate two temporary sums in the new pressure formula, linear absorbing case.
float Get_dz() const
Get dz value.
Definition: Parameters.h:115
void AddPressureSource()
Add in pressure source.
void AddMatrices()
Populate the container based on the simulation type.
void FreeMatrices()
Destroy and free all matrices.
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.
TOutputMessage OUT_FMT_STORING_SENSOR_DATA
Output message.
double GetTotalTime() const
Get total simulation time.
TRealMatrix & Get_ux_sgx()
Get the ux_sgx matrix from the container.
static std::string WordWrapString(const std::string &inputString, const std::string &delimiters, const int indentation=0, const int lineSize=65)
Wrap the line based on logger conventions.
Definition: Logger.cpp:130
int GetSamplerGridSize1D() const
Get Number of blocks for the 1D data sampling kernels.
size_t GetProgressPrintInterval() const
Get progress print interval.
Definition: Parameters.h:222
TRealMatrix & Get_ux_source_input()
Get the ux_source_input matrix from the container.
size_t GetStartTimeIndex() const
Get start time index for sensor recording.
Definition: Parameters.h:225
static void Create_FFT_Plan_1DZ_R2C(const TDimensionSizes &inMatrixDims)
Create static cuFFT plan for Real-to-Complex in the Z dimension.
int GetCUDACodeVersion()
Get the CUDA architecture and GPU code version the code was compiled with.
TRealMatrix & Get_absorb_tau()
Get the absorb_tau matrix from the container.
TComplexMatrix & Get_x_shift_neg_r()
Get the x_shift_neg_r matrix from the container.
void ComputeDensityLinear()
Compute new values of acoustic density for linear case.
void SetFileType(const THDF5_FileHeader::TFileType fileType)
Set file type.
Definition: HDF5_File.cpp:1491
static const int DEFAULT_DEVICE_IDX
Default Device Index - no default GPU.
static TLogLevel GetLevel()
Get the log level.
Definition: Logger.h:74
void LoadDataFromCheckpointFile(THDF5_File &checkpointFile)
Load all matrices from the output HDF5 file.
TErrorMessage ERR_FMT_PATH_DELIMITERS
delimiters for linux paths
Definition: ErrorMessages.h:52
void SumPressureLinearLossless()
Sum sub-terms for new p, linear lossless case.
TErrorMessage ERR_FMT_BAD_OUTPUT_FILE_FORMAT
KSpaceFirstOrder3DSolver error message.
void Compute_FFT_3D_C2R(TRealMatrix &outMatrix)
Compute 3D out-of-place Complex-to-Real FFT.
float Get_c_ref() const
Get c_ref value.
Definition: Parameters.h:118
void SetMinorFileVersion()
Set minor file version.
Definition: HDF5_File.h:849
void PrintStatistics()
Print progress statistics.
bool isTimestepRightAfterRestore
This variable is true when calculating first time step after restore from checkpoint (to allow asynch...
TRealMatrix & Get_duydy()
Get the duydy matrix from the container.
static void DestroyAllPlansAndStaticData()
Destroy all static plans and error messages.
TErrorMessage ERR_FMT_BAD_CHECKPOINT_FILE_FORMAT
KSpaceFirstOrder3DSolver error message.
TOutputMessage OUT_FMT_CHECKPOINT_HEADER
Output message.
TRealMatrix & Get_pml_x()
Get the pml_x matrix from the container.
bool IsCheckpointInterruption() const
Was the loop interrupted to checkpoint?
TRealMatrix & Get_uz_shifted()
Get the uz_shifted matrix from the container.
bool IsTimeToCheckpoint()
Is time to checkpoint (save actual state on disk).
void InitializeFFTPlans()
Initialize FFT plans.
void SetMajorFileVersion()
Set major file version.
Definition: HDF5_File.h:840
TRealMatrix & Get_dzudzn_sgz()
Get the dzudzn_sgz matrix from the container.
void Open(const std::string &fileName, unsigned int flags=H5F_ACC_RDONLY)
Open a file.
Definition: HDF5_File.cpp:127
void ComputeGradientVelocity()
Compute new values of acoustic velocity gradients.
void GetExecutionTimes(double &totalTime, double &loadTime, double &preProcessingTime, double &simulationTime, double &postprocessingTime)
Get execution times stored in the output file header.
Definition: HDF5_File.cpp:1575
void AddVelocitySource()
Add in velocity source.
static void Create_FFT_Plan_1DX_R2C(const TDimensionSizes &inMatrixDims)
Create static cuFFT plan for Real-to-Complex in the X dimension.
TComplexMatrix & Get_ddz_k_shift_neg()
Get the ddz_k_shift_neg matrix from the container.
TCUFFTComplexMatrix & Get_cufft_z_temp()
Get the FFT_Z_temp from the container.
TOutputStreamContainer outputStreamContainer
Output stream container.
TOutputMessage OUT_FMT_CURRENT_DEVICE_MEMORY
Output message.
TRealMatrix & Get_kappa()
Get the kappa matrix from the container.
TOutputMessage OUT_FMT_READING_INPUT_FILE
Output message.
virtual void FreeMemory()
Memory deallocation.
TTimeMeasure dataLoadTime
Data load time of the simulation.
TMatrixName uz_final_NAME
uz_final variable name
Definition: MatrixNames.h:319
TOutputMessage OUT_FMT_VISUAL_STUDIO_COMPILER
Print version output message.
TOutputMessage OUT_FMT_CURRENT_HOST_MEMORY
Output message.
TRealMatrix & Get_pml_z()
Get the pml_z matrix from the container.
TOutputMessage OUT_FMT_SIMULATION_FINAL_SEPARATOR
Output message.
TComplexMatrix & Get_ddz_k_shift_pos()
Get the ddz_k_shift_pos matrix from the container.
TRealMatrix & Get_dxudxn_sgx()
Get the dxudxn_sgx matrix from the container.
TOutputMessage OUT_FMT_READING_CHECKPOINT_FILE
Output message.
TOutputMessage OUT_FMT_CUDA_DEVICE_INFO_NA
Print version output message.
TMatrixName sensor_mask_index_NAME
sensor_mask_index variable name
Definition: MatrixNames.h:166
TOutputMessage OUT_FMT_SIMULATION_PROGRESS
Output message.
size_t Get_ux_source_flag() const
Get ux_source_flag value.
Definition: Parameters.h:137
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.
TRealMatrix & Get_pml_z_sgz()
Get the pml_z_sgz matrix from the container.
void Compute_FFT_1DZ_R2C(TRealMatrix &inMatrix)
Compute 1D out-of-place Real-to-Complex FFT in the Z dimension.
TOutputMessage OUT_FMT_CUDA_RUNTIME_NA
Print version output message.
TMatrixName p_final_NAME
p_final variable name
Definition: MatrixNames.h:279
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 GenerateKappa()
Generate kappa matrix for non-absorbing media.
static TParameters & GetInstance()
Get instance of the singleton class.
Definition: Parameters.cpp:70
TRealMatrix & Get_c2()
Get the c^2 matrix from the container.
void Set_t_index(const size_t new_t_index)
Set simulation time step - should be used only when recovering from checkpoint.
Definition: Parameters.h:104
void SetUpDeviceConstants() const
Upload useful simulation constants into device constant memory.
virtual void WriteDataToHDF5File(THDF5_File &file, TMatrixName &matrixName, const size_t compressionLevel)
Write data into the HDF5 file.
TOutputMessage OUT_FMT_GNU_COMPILER
Print version output message.
static std::string GetCurrentHDF5_MinorVersion()
Get string version of current Minor version.
Definition: HDF5_File.h:831
static void Create_FFT_Plan_1DY_C2R(const TDimensionSizes &outMatrixDims)
Create static cuFFT plan for Complex-to-Real in the Y dimension.
void LoadDataFromInputFile(THDF5_File &inputFile)
Load all matrices from the input HDF5 file.
TOutputMessage OUT_FMT_CHECKPOINT_TIME_STEPS
Output message.
TCUFFTComplexMatrix & Get_cufft_y_temp()
Get the FFT_Y_temp from the container.
TOutputMessage OUT_FMT_KWAVE_VERSION
Output message.
void Compute_FFT_1DY_C2R(TRealMatrix &outMatrix)
Compute 1D out-of-place Complex-to-Real FFT in the Y dimension.
TOutputMessage OUT_FMT_SSE2
Print version output message.
TRealMatrix & Get_ux_shifted()
Get the ux_shifted matrix from the container.
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 SetExecutionTimes(const double totalTime, const double loadTime, const double preProcessingTime, const double simulationTime, const double postprocessingTime)
Set execution times.
Definition: HDF5_File.cpp:1552
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.
double GetPreProcessingTime() const
Get pre-processing time.
TTimeMeasure iterationTime
Iteration time of the simulation.
static void ErrorAndTerminate(const std::string &errorMessage)
Log an error and terminate the execution.
Definition: Logger.cpp:94
TOutputMessage OUT_FMT_FINSIH_NO_DONE
Output message - finish line without done.
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.
TMatrixName Ny_NAME
Ny variable name.
Definition: MatrixNames.h:73
TRealMatrix & Get_dyudyn_sgy()
Get the dyudyn_sgy matrix from the container.
void WriteScalarValue(const hid_t parentGroup, TMatrixName &datasetName, const float value)
Write the scalar value under a specified group, float value.
Definition: HDF5_File.cpp:748
static void Create_FFT_Plan_1DZ_C2R(const TDimensionSizes &outMatrixDims)
Create static cuFFT plan for Complex-to-Real in the Z dimension.
virtual size_t GetElementCount() const
Get element count of the matrix.
TCUFFTComplexMatrix & Get_cufft_x_temp()
Get the CUFFT_X_temp from the container.
TRealMatrix & Get_p_source_input()
Get the p_source_input matrix from the container.
virtual void WriteDataToHDF5File(THDF5_File &file, TMatrixName &matrixName, const size_t compressionLevel)
Write data into the HDF5 file.
Definition: RealMatrix.cpp:102
TOutputMessage OUT_FMT_POST_PROCESSING
Output message.
static void Create_FFT_Plan_3D_C2R(const TDimensionSizes &outMatrixDims)
Create static cuFFT plan for Complex-to-Real.
bool Get_c0_scalar_flag() const
Get c0_scalar_flag value.
Definition: Parameters.h:185
TOutputMessage OUT_FMT_SSE41
Print version output message.
TTimeMeasure totalTime
Total time of the simulation.
void Close()
Close file.
Definition: HDF5_File.cpp:175
void FreeStreams()
Free all streams - destroy them.
TOutputMessage OUT_FMT_DATA_LOADING
Output message.
TRealMatrix & Get_temp_3_real_3D()
Get the Temp_3_RS3D matrix from the container.
TOutputMessage OUT_FMT_CUDA_GRID_SHAPE_FORMAT
Output message.
void Compute_c2()
Calculate square of velocity.
void PrintFullNameCodeAndLicense() const
Print the code name and license.
TOutputMessage OUT_FMT_CUDA_CAPABILITY
Print version output message.
TRealMatrix & Get_transducer_source_input()
Get the transducer_source_input matrix from the container.
TOutputMessage OUT_FMT_CUDA_CODE_ARCH
Print version output message.
void SetCodeName(const std::string &codeName)
Set code name.
Definition: HDF5_File.h:808
bool Get_BonA_scalar_flag() const
Get BonA_scalar_flag value.
Definition: Parameters.h:195
bool IsOpen() const
Is the file opened?
Definition: HDF5_File.h:545
virtual float * GetHostData()
Get raw CPU data out of the class (for direct CPU kernel access).
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 SumPressureNonlinearLossless()
Sum sub-terms for new p, linear lossless case.
TComplexMatrix & Get_z_shift_neg_r()
Get the y_shift_neg_r from the container.
The header file containing the main class of the project responsible for the entire 3D fluid simulati...
const std::string GetCodeName() const
Get code name - release code version.
void AddStreams(TMatrixContainer &matrixContainer)
Add all streams into the container.
void CheckCheckpointFile()
Check the checkpoint file has the correct format and version.
void RecomputeIndicesToMatlab()
Recompute indices C++ -> MATLAB.
static void Create_FFT_Plan_1DX_C2R(const TDimensionSizes &outMatrixDims)
Create static cuFFT plan for Complex-to-Real in the X dimension.
TRealMatrix & Get_p0_source_input()
Get the p0_source_input from the container.
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 ReadHeaderFromCheckpointFile(THDF5_File &checkpointFile)
Read Header from checkpoint file (necessary for checkpoint-restart).
Definition: HDF5_File.cpp:1336
void Compute_FFT_1DX_R2C(TRealMatrix &inMatrix)
Compute 1D out-of-place Real-to-Complex FFT in the X dimension.
TRealMatrix & Get_p()
Get the p matrix from the container.
TComplexMatrix & Get_ddy_k_shift_pos()
Get the ddy_k_shift_pos matrix from the container.
void WriteHeaderToCheckpointFile(THDF5_File &checkpointFile)
Write header to the output file.
Definition: HDF5_File.cpp:1386
TOutputMessage OUT_FMT_FFT_PLANS
Output message.
TRealMatrix & Get_BonA()
Get the BonA matrix from the container.
TRealMatrix & Get_dxudxn()
Get the dxudxn matrix from the container.
Full level of verbosity.
Definition: Logger.h:67
TRealMatrix & Get_uz_sgz()
Get the uz_sgz matrix from the container.
static std::string GetCurrentHDF5_MajorVersion()
Get string version of current Major version.
Definition: HDF5_File.h:821
TMatrixName sensor_mask_corners_NAME
sensor_mask_corners variable name
Definition: MatrixNames.h:170
double GetCumulatedElapsedTimeOverPreviousLegs() const
Get time spent in previous legs.
Definition: TimeMeasure.h:159
TMatrixContainer matrixContainer
Matrix container with all the matrix classes.
virtual TDimensionSizes GetDimensionSizes() const
Get dimension sizes of the matrix.
void WriteHeaderToOutputFile(THDF5_File &outputFile)
Write header to the output file.
Definition: HDF5_File.cpp:1369
TMatrixName uy_final_NAME
uy_final variable name
Definition: MatrixNames.h:317
bool Get_alpha_coeff_scalar_flag() const
Get alpha_coeff_scalar_flag value.
Definition: Parameters.h:180
static void Create_FFT_Plan_1DY_R2C(const TDimensionSizes &inMatrixDims)
Create static cuFFT plan for Real-to-Complex in the Y dimension.
Class responsible for CUDA runtime setup.
std::string GetDeviceName() const
Get the name of the device used.
TOutputMessage OUT_FMT_SIMULATION_HEADER
Output message.
void CreateMatrices()
Create all matrices in the container.
TErrorMessage ERR_FMT_OUTPUT_DIMENSIONS_NOT_MATCH
KSpaceFirstOrder3DSolver error message.
TRealMatrix & Get_dzudzn()
Get the dzudzn matrix from the container.
TOutputMessage OUT_FMT_SIMULATOIN_END_SEPARATOR
Output message.
size_t Get_t_index() const
Get simulation time step.
Definition: Parameters.h:102
TOutputMessage OUT_FMT_SSE3
Print version output message.
void ReadHeaderFromOutputFile(THDF5_File &outputFile)
Read Header from output file (necessary for checkpoint-restart).
Definition: HDF5_File.cpp:1290
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).
double GetSimulationTime() const
Get simulation time (time loop).
void Calculate_p0_source()
Calculate p0_ ource.
TRealMatrix & Get_absorb_nabla2()
Get the absorb_nabla2 matrix from the container.
TDimensionSizes GetFullDimensionSizes() const
Full dimension sizes of the simulation (real classes).
Definition: Parameters.h:95
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 SetHostName()
Set host name.
Definition: HDF5_File.cpp:1503
The header file containing a class responsible for printing out info and error messages (stdout...
TOutputMessage OUT_FMT_AVX2
Print version output message.
void GenerateKappaAndNablas()
Generate kappa matrix, absorb_nabla1, absorb_nabla2 for absorbing media.
void Compute_FFT_3D_R2C(TRealMatrix &inMatrix)
Compute 3D out-of-place Real-to-Complex FFT.
void StoreSensorData()
Store sensor data.
TRealMatrix & Get_rho0()
Get the rho0 matrix from the container.
void WriteOutputDataInfo()
Write statistics and header into the output file.
TOutputMessage OUT_FMT_CUDA_DEVICE
Print version output message.
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.
TRealMatrix & Get_uz_source_input()
Get the uz_source_input matrix from the container.
void Compute_FFT_1DX_C2R(TRealMatrix &outMatrix)
Compute 1D out-of-place Complex-to-Real FFT in the X dimension.
static void Flush(const TLogLevel queryLevel)
Flush output messages.
Definition: Logger.cpp:109
void SumPressureTermsLinear(TRealMatrix &absorb_tau_temp, TRealMatrix &absorb_eta_temp, TRealMatrix &rhoxyz_sum)
Sum sub-terms to calculate new pressure, linear case.
TIndexMatrix & Get_sensor_mask_index()
Get the sensor_mask_index matrix from the container.
TOutputMessage OUT_FMT_ELAPSED_TIME
Output message.
void Start()
Get start timestamp.
Definition: TimeMeasure.h:100
void ReadScalarValue(const hid_t parentGroup, TMatrixName &datasetName, float &value)
Read the scalar value under a specified group, float value.
Definition: HDF5_File.cpp:855
bool IsStore_u_non_staggered_raw() const
Is –u_non_staggered_raw set?
Definition: Parameters.h:252
Class storing all parameters of the simulation.
Definition: Parameters.h:49
TOutputMessage OUT_FMT_INTEL_COMPILER
Print version output message.
TIndexMatrix & Get_delay_mask()
Get the delay_mask matrix from the container.
TRealMatrix & Get_temp_2_real_3D()
Get the Temp_2_RS3D matrix from the container.
size_t Get_uy_source_flag() const
Get uy_source_flag value.
Definition: Parameters.h:139
const TCUDAParameters & GetCudaParameters() const
Get class with CUDA Parameters (runtime setup), cosnt version.
Definition: Parameters.h:76
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.
THDF5_FileHeader & GetFileHeader()
Get File header handle.
Definition: Parameters.h:87
void PostProcessing()
Post processing, and closing the output streams.
TRealMatrix & Get_dt_rho0_sgz()
Get the dt.*rho0_sgz matrix from the container.
static bool IsHDF5(const std::string &fileName)
Does the file exist?
Definition: HDF5_File.cpp:158
hid_t GetRootGroup() const
Get handle to the root group.
Definition: HDF5_File.h:573
TOutputMessage OUT_FMT_LINUX_BUILD
Print version output message.
TIndexMatrix & Get_p_source_index()
Get the p_source_index matrix from the container.
THDF5_File & GetCheckpointFile()
Gest checkpoint file handle.
Definition: Parameters.h:85
void SampleStreams()
Sample all streams (only sample, no disk operations).
TOutputMessage OUT_FMT_VERSION_GIT_HASH
Print version output message.
bool IsStore_p_final() const
Is –p_final specified at the command line?
Definition: Parameters.h:247
void ComputeMainLoop()
Compute the main time loop of the kspaceFirstOrder3D.
TRealMatrix & Get_dyudyn()
Get the dyudyn matrix from the container.
TRealMatrix & Get_pml_y_sgy()
Get the pml_y_sgy matrix from the container.
TRealMatrix & Get_pml_y()
Get the pml_y matrix from the container.
TIndexMatrix & Get_u_source_index()
Get the u_source_index matrix from the container.
void ComputePressureLinear()
Compute acoustic pressure for linear case.
int GetSamplerBlockSize1D() const
Get number of threads for the 1D data sampling kernels.
double GetCumulatedTotalTime() const
Get total simulation time cumulated over all legs.
The header file containing routines for error messages and error messages common for both linux and w...
TIndexMatrix & Get_sensor_mask_corners()
Get the sensor_mask_corners matrix from the container.
size_t Get_p0_source_flag() const
Get p0_source_flag value.
Definition: Parameters.h:150
void ReopenStreams()
Reopen streams after checkpoint file (datasets).
TRealMatrix & Get_duxdx()
Get the duxdx matrix from the container.
size_t ny
number of elements in the y direction
static void Create_FFT_Plan_3D_R2C(const TDimensionSizes &inMatrixDims)
Create static cuFFT plan for Real-to-Complex.
float Get_dt() const
Get dt value.
Definition: Parameters.h:109
TCUFFTComplexMatrix & Get_cufft_shift_temp()
Get the FFT_shift_temp the container.
float Get_dx() const
Get dx value.
Definition: Parameters.h:111
TOutputMessage OUT_FMT_CUDA_SAMPLER_GRID_SHAPE
Output message.
void SumPressureTermsNonlinear(TRealMatrix &absorb_tau_temp, TRealMatrix &absorb_eta_temp, TRealMatrix &BonA_temp)
Sum sub-terms to calculate new pressure, non-linear case.
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).
TOutputMessage OUT_FMT_READING_OUTPUT_FILE
Output message.
int GetSolverGridSize1D() const
Get number of block for 1D grid used by kSpaceSolver.
virtual void LoadInputData()
Load simulation data from the input file.
TMatrixName Nz_NAME
Nz variable name.
Definition: MatrixNames.h:75
float Get_dy() const
Get dy value.
Definition: Parameters.h:113
int GetDeviceIdx() const
Get Idx of the device being used.
void FlushRawStreams()
Flush streams to disk - only raw streams.
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.
size_t Get_nonlinear_flag() const
Get nonlinear_flag value.
Definition: Parameters.h:161
Basic (default) level of verbosity.
Definition: Logger.h:63
The class for real matrices.
Definition: RealMatrix.h:45
TOutputMessage OUT_FMT_CUDA_RUNTIME
Print version output message.
void RecomputeIndicesToCPP()
Recompute indices MATALAB->C++.
virtual ~TKSpaceFirstOrder3DSolver()
Destructor.
float & Get_absorb_eta_scalar()
Get absorb_eta_scalar value.
Definition: Parameters.h:190
void Compute_FFT_1DZ_C2R(TRealMatrix &outMatrix)
Compute 1D out-of-place Complex-to-Real FFT in the Z dimension.
void Create(const std::string &fileName, unsigned int flags=H5F_ACC_TRUNC)
Create a file.
Definition: HDF5_File.cpp:97
TRealMatrix & Get_dt_rho0_sgx()
Get the dt.*rho0_sgx matrix from the container.
TMatrixName Nx_NAME
Nx variable name.
Definition: MatrixNames.h:71
float & Get_rho0_sgy_scalar()
Get rho0_sgy_scalar value.
Definition: Parameters.h:206
float & Get_rho0_scalar()
Get rho0_scalar value.
Definition: Parameters.h:202
size_t Get_uz_source_flag() const
Get uz_source_flag value.
Definition: Parameters.h:141
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 ComputePressurePartsNonLinear(TRealMatrix &rho_part, TRealMatrix &BonA_part, TRealMatrix &du_part)
Calculate three temporary sums in the new pressure formula, non-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 SaveScalarsToFile(THDF5_File &outputFile)
Save scalar values into the output HDF5 file.
Definition: Parameters.cpp:432
void SetActualCreationTime()
Set creation time.
Definition: HDF5_File.cpp:1429
TOutputMessage OUT_FMT_CUDA_DEVICE_NAME
Print version output message.
double GetCumulatedPostProcessingTime() const
Get post-processing time cumulated over all legs.
TOutputMessage OUT_FMT_LAST_SEPARATOR
Output message -last separator.
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.
void CheckpointStreams()
Checkpoint streams.
bool CheckMajorFileVersion()
Check major file version.
Definition: HDF5_File.h:862
TErrorMessage ERR_FMT_BAD_MINOR_FILE_VERSION
Command line parameters error message.
bool CheckMinorFileVersion()
Check minor file version.
Definition: HDF5_File.h:872
THDF5_FileHeader::TFileType GetFileType()
Get File type.
Definition: HDF5_File.cpp:1472
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.
TTimeMeasure postProcessingTime
Post-processing time of the simulation.
void CheckOutputFile()
Check the output file has the correct format and version.
void GenerateTauAndEta()
Generate absorb_tau, absorb_eta for heterogenous media.
TRealMatrix & Get_uy_sgy()
Get the uy_sgy matrix from the container.
void ComputeAbsorbtionTerm(TCUFFTComplexMatrix &fft1, TCUFFTComplexMatrix &fft2, const TRealMatrix &absorb_nabla1, const TRealMatrix &absorb_nabla2)
Compute absorbing term with abosrb_nabla1 and absorb_nabla2.
TOutputMessage OUT_FMT_CUDA_SOLVER_GRID_SHAPE
Output message.
size_t Get_p_source_flag() const
Get p_source_flag value.
Definition: Parameters.h:148
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 LoadElapsedTimeFromOutputFile(THDF5_File &o1utputFile)
Reads the header of the output file and sets the cumulative elapsed time from the first log...
TParameters & parameters
Global parameters of the simulation.
size_t Get_transducer_source_flag() const
Get transducer_source_flag value.
Definition: Parameters.h:163
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 ComputePressureNonlinear()
Compute acoustic pressure for nonlinear case.
Name space for all CUDA kernels used in the 3D solver.
size_t GetCompressionLevel() const
Get compression level.
Definition: Parameters.h:218
TTimeMeasure preProcessingTime
Pre-processing time of the simulation.
static std::string FormatMessage(const std::string &format, Args...args)
C++-11 replacement for sprintf that works with std::string instead of char *.
Definition: Logger.h:126
The header file containing the matrix container and the related matrix record class.
int GetSolverBlockSize1D() const
Get number of threads for 1D block used by kSpaceSolver.
double GetCumulatedSimulationTime() const
Get simulation time (time loop) cumulated over all legs.
TOutputMessage OUT_FMT_WINDOWS_BUILD
Print version output message.
TOutputMessage OUT_FMT_SSE42
Print version output message.
void Stop()
Get stop timestamp.
Definition: TimeMeasure.h:118
void CopyMatricesToDevice()
Copy all matrices from host to device (CPU -> GPU).
void SetNumberOfCores()
Set number of cores.
Definition: HDF5_File.cpp:1593
size_t actPercent
Percentage of the simulation done.
TRealMatrix & Get_duzdz()
Get the duzdz matrix from the container.
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.
bool IsStore_u_final() const
Is –u_final specified at the command line.
Definition: Parameters.h:264
TOutputMessage OUT_FMT_CUDA_DRIVER
Print version output message.
virtual void Compute()
Compute the k-space simulation.
TErrorMessage ERR_FMT_BAD_MAJOR_File_Version
Command line parameters error message.
TRealMatrix & Get_rhoz()
Get the rhoz matrix from the container.
TComplexMatrix & Get_ddx_k_shift_neg()
Get the ddx_k_shift_neg matrix from the container.
double GetCumulatedDataLoadTime() const
Get data load time cumulated over all legs.
double GetPostProcessingTime() const
Get post-processing time.
virtual float * GetDeviceData()
Get raw GPU data out of the class (for direct GPU kernel access).
void ComputeVelocityShiftInZ(TCUFFTComplexMatrix &cufft_shift_temp, const TComplexMatrix &z_shift_neg_r)
Compute the velocity shift in Fourier space over the Z axis.
float & Get_absorb_tau_scalar()
Get absorb_tau_scalar value.
Definition: Parameters.h:192
TTimeMeasure simulationTime
Simulation time of the simulation.
void CreateStreams()
Create all streams - opens the datasets.
TComplexMatrix & Get_ddy_k_shift_neg()
Get the ddy_k_shift_neg matrix from the container.
bool Get_rho0_scalar_flag() const
Get rho0_scalar_flag value.
Definition: Parameters.h:200
virtual void AllocateMemory()
Memory allocation.
void SaveCheckpointData()
Save checkpoint data.
void CalculateShiftedVelocity()
Calculate ux_shifted, uy_shifted and uz_shifted.
double GetCumulatedPreProcessingTime() const
Get pre-processing time cumulated over all legs.
TErrorMessage ERR_FMT_CHECKPOINT_DIMENSIONS_NOT_MATCH
KSpaceFirstOrder3DSolver error message.
TComplexMatrix & Get_ddx_k_shift_pos()
Get the ddx_k_shift_pos matrix from the container.
THDF5_File & GetInputFile()
Get input file handle.
Definition: Parameters.h:81
TOutputMessage OUT_FMT_BUILD_NO_DATE_TIME
Print version output message.
float & Get_alpha_coeff_scalar()
Get alpha_coeff_scalar value. Note: cannot be const because of other optimizations.
Definition: Parameters.h:182
float Get_alpha_power() const
Get alpha_power value.
Definition: Parameters.h:120
TOutputMessage OUT_FMT_CREATING_OUTPUT_FILE
Output message.
size_t Get_absorbing_flag() const
Get absorbing_flag value.
Definition: Parameters.h:159
void StoreDataIntoCheckpointFile(THDF5_File &checkpointFile)
Store selected matrices into the checkpoint file.
TRealMatrix & Get_pml_x_sgx()
Get the pml_x_sgx matrix from the container.
double GetDataLoadTime() const
Get data load time.
size_t Get_nonuniform_grid_flag() const
Get nonuniform_grid_flag value.
Definition: Parameters.h:157
TOutputMessage OUT_FMT_MAC_OS_BUILD
Print version output message.
void SetKernelConfiguration()
Set kernel configurations based on the simulation parameters.
TRealMatrix & Get_rhox()
Get the rhox matrix from the container.
const cudaDeviceProp & GetDeviceProperties() const
Return properties of currently used GPU.
size_t nz
number of elements in the z direction
static void Log(const TLogLevel queryLevel, const std::string &format, Args...args)
Log desired activity for a given log level, version with string format.
Definition: Logger.h:85
TOutputMessage OUT_FMT_SEPARATOR
Output message - separator.
void ComputeDensityNonliner()
Compute new values of acoustic density for non-linear case.
virtual void ScalarDividedBy(const float scalar)
Divide scalar/ matrix_element[i].
double GetElapsedTime() const
Get elapsed time.
Definition: TimeMeasure.h:140
void CloseStreams()
Close all streams.
void Increment_t_index()
Increment simulation time step.
Definition: Parameters.h:106
TOutputMessage OUT_FMT_AVX
Print version output message.
virtual void CopyFromDevice()
Copy data from GPU -> CPU (Device -> Host).
TOutputMessage OUT_FMT_PRE_PROCESSING
Output message.
size_t Get_nt() const
Get Nt value.
Definition: Parameters.h:100
TComplexMatrix & Get_y_shift_neg_r()
Get the y_shift_neg_r from the container.
TRealMatrix & Get_temp_1_real_3D()
Get the Temp_1_RS3D matrix from the container.
std::string GetOutputFileName() const
Get output file name.
Definition: Parameters.h:213
TOutputMessage OUT_FMT_CREATING_CHECKPOINT
Output message.
TOutputMessage OUT_FMT_FAILED
Output message - failed message.
std::string GetGitHash() const
Get Git hash of the code.
Definition: Parameters.cpp:513
TOutputMessage OUT_FMT_COMP_RESOURCES_HEADER
Output message.
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.
TRealMatrix & Get_rhoy()
Get the rhoy matrix from the container.
float & Get_rho0_sgx_scalar()
Get rho0_sgx_scalar value.
Definition: Parameters.h:204
TSensorMaskType Get_sensor_mask_type() const
Get sensor mask type (linear or corners).
Definition: Parameters.h:166
bool IsCopySensorMask() const
is –copy_mask set
Definition: Parameters.h:267
void GenerateInitialDenisty()
Calculate dt ./ rho0 for non-uniform grids.
TRealMatrix & Get_absorb_eta()
Get the absorb_eta matrix from the container.
bool IsCheckpointEnabled() const
Is checkpoint enabled?
Definition: Parameters.h:228
TMatrixName ux_final_NAME
ux_final variable name
Definition: MatrixNames.h:315
Class wrapping the HDF5 routines.
Definition: HDF5_File.h:500
void ComputeVelocity()
Compute new values of acoustic velocity.
void ComputeVelocityShiftInY(TCUFFTComplexMatrix &cufft_shift_temp, const TComplexMatrix &y_shift_neg_r)
Compute the velocity shift in Fourier space over the Y axis.
Structure with 4D dimension sizes (3 in space and 1 in time).
size_t GetCheckpointInterval() const
Get checkpoint interval.
Definition: Parameters.h:230
TOutputMessage OUT_FMT_DONE
Output message - Done with two spaces.
TRealMatrix & Get_uy_shifted()
Get the uy_shifted matrix from the container.
THDF5_File & GetOutputFile()
Get output file handle.
Definition: Parameters.h:83
TOutputMessage OUT_FMT_CUDA_DEVICE_NAME_PADDING
Print version output message.
size_t GetHostMemoryUsageInMB()
Get memory usage in MB on the CPU side.
void SetMemoryConsumption(const size_t totalMemory)
Set memory consumption.
Definition: HDF5_File.cpp:1531
float & Get_rho0_sgz_scalar()
Get rho0_sgz_scalar value.
Definition: Parameters.h:208
void PostProcessStreams()
Post-process all streams and flush them to the file.
TOutputMessage OUT_FMT_STORING_CHECKPOINT_DATA
Output message.
TMatrixName t_index_NAME
t_index name
Definition: MatrixNames.h:50
TRealMatrix & Get_uy_source_input()
Get the uy_source_input matrix from the container.
void PreProcessingPhase()
Compute pre-processing phase.
double GetCumulatedElapsedTimeOverAllLegs() const
Get cumulated elapsed time over all simulation legs.
Definition: TimeMeasure.h:150