![]() |
kspaceFirstOrder3D-OMP 1.0
The C++ implementation of the k-wave toolbox for the time-domain simulation of acoustic wave fields in 3D
|
00001 00030 #include <omp.h> 00031 #include <iostream> 00032 #include <string.h> 00033 #include <sstream> 00034 #include <exception> 00035 #include <stdexcept> 00036 00037 #include <Parameters/Parameters.h> 00038 00039 #include <Utils/MatrixNames.h> 00040 #include <Utils/ErrorMessages.h> 00041 00042 using namespace std; 00043 00044 //----------------------------------------------------------------------------// 00045 // Constants // 00046 //----------------------------------------------------------------------------// 00047 00048 00049 00050 00051 //----------------------------------------------------------------------------// 00052 // Definitions // 00053 //----------------------------------------------------------------------------// 00054 00055 00056 00057 00058 00059 bool TParameters::ParametersInstanceFlag = false; 00060 00061 TParameters* TParameters::ParametersSingleInstance = NULL; 00062 00063 00064 //----------------------------------------------------------------------------// 00065 // Implementation // 00066 // public methods // 00067 //----------------------------------------------------------------------------// 00068 00072 TParameters* TParameters::GetInstance(){ 00073 00074 00075 if(!ParametersInstanceFlag) 00076 { 00077 ParametersSingleInstance = new TParameters(); 00078 ParametersInstanceFlag = true; 00079 return ParametersSingleInstance; 00080 } 00081 else 00082 { 00083 return ParametersSingleInstance; 00084 } 00085 00086 }// end of Create 00087 //------------------------------------------------------------------------------ 00088 00094 void TParameters::ParseCommandLine(int argc, char** argv){ 00095 00096 CommandLinesParameters.ParseCommandLine(argc, argv); 00097 00098 if (CommandLinesParameters.IsVersion()){ 00099 return; 00100 } 00101 00102 ReadScalarsFromHDF5InputFile(HDF5_InputFile); 00103 00104 if (CommandLinesParameters.IsBenchmarkFlag()) 00105 Nt = CommandLinesParameters.GetBenchmarkTimeStepsCount(); 00106 00107 if ((Nt <= CommandLinesParameters.GetStartTimeIndex()) || 00108 ( 0 > CommandLinesParameters.GetStartTimeIndex()) ){ 00109 fprintf(stderr,Parameters_ERR_FMT_Illegal_StartTime_value, (long) 1, Nt); 00110 CommandLinesParameters.PrintUsageAndExit(); 00111 } 00112 00113 }// end of ParseCommandLine 00114 //------------------------------------------------------------------------------ 00115 00116 00122 void TParameters::ReadScalarsFromHDF5InputFile(THDF5_File & HDF5_InputFile){ 00123 00124 00125 TDimensionSizes ScalarSizes(1,1,1); 00126 00127 if (!HDF5_InputFile.IsOpened()) { 00128 00129 // Open file 00130 try{ 00131 HDF5_InputFile.Open(CommandLinesParameters.GetInputFileName().c_str()); 00132 } catch (ios::failure e){ 00133 fprintf(stderr,"%s",e.what()); 00134 PrintUsageAndExit(); 00135 } 00136 00137 } 00138 00139 HDF5_FileHeader.ReadHeaderFromInputFile(HDF5_InputFile); 00140 00141 if (HDF5_FileHeader.GetFileType() != THDF5_FileHeader::hdf5_ft_input) { 00142 char ErrorMessage[256] = ""; 00143 sprintf(ErrorMessage,Parameters_ERR_FMT_IncorrectInputFileFormat,GetInputFileName().c_str()); 00144 throw ios::failure(ErrorMessage); 00145 } 00146 00147 if (!HDF5_FileHeader.CheckMajorFileVersion()) { 00148 char ErrorMessage[256] = ""; 00149 sprintf(ErrorMessage,Parameters_ERR_FMT_IncorrectMajorHDF5FileVersion,GetInputFileName().c_str(), 00150 HDF5_FileHeader.GetSupportedHDF5_MajorVersion().c_str()); 00151 throw ios::failure(ErrorMessage); 00152 } 00153 00154 if (!HDF5_FileHeader.CheckMinorFileVersion()) { 00155 char ErrorMessage[256] = ""; 00156 sprintf(ErrorMessage,Parameters_ERR_FMT_IncorrectMinorHDF5FileVersion,GetInputFileName().c_str(), 00157 HDF5_FileHeader.GetSupportedHDF5_MinorVersion().c_str()); 00158 throw ios::failure(ErrorMessage); 00159 } 00160 00161 00162 00163 HDF5_InputFile.ReadCompleteDataset(Nt_Name, ScalarSizes, &Nt); 00164 00165 HDF5_InputFile.ReadCompleteDataset(dt_Name, ScalarSizes, &dt); 00166 HDF5_InputFile.ReadCompleteDataset(dx_Name, ScalarSizes, &dx); 00167 HDF5_InputFile.ReadCompleteDataset(dy_Name, ScalarSizes, &dy); 00168 HDF5_InputFile.ReadCompleteDataset(dz_Name, ScalarSizes, &dz); 00169 00170 HDF5_InputFile.ReadCompleteDataset(c_ref_Name, ScalarSizes, &c_ref); 00171 00172 HDF5_InputFile.ReadCompleteDataset(pml_x_size_Name, ScalarSizes, &pml_x_size); 00173 HDF5_InputFile.ReadCompleteDataset(pml_y_size_Name, ScalarSizes, &pml_y_size); 00174 HDF5_InputFile.ReadCompleteDataset(pml_z_size_Name, ScalarSizes, &pml_z_size); 00175 00176 HDF5_InputFile.ReadCompleteDataset(pml_x_alpha_Name, ScalarSizes, &pml_x_alpha); 00177 HDF5_InputFile.ReadCompleteDataset(pml_y_alpha_Name, ScalarSizes, &pml_y_alpha); 00178 HDF5_InputFile.ReadCompleteDataset(pml_z_alpha_Name, ScalarSizes, &pml_z_alpha); 00179 00180 00181 00182 00183 long X, Y, Z; 00184 HDF5_InputFile.ReadCompleteDataset(Nx_Name , ScalarSizes, &X); 00185 HDF5_InputFile.ReadCompleteDataset(Ny_Name , ScalarSizes, &Y); 00186 HDF5_InputFile.ReadCompleteDataset(Nz_Name , ScalarSizes, &Z); 00187 00188 FullDimensionSizes.X = X; 00189 FullDimensionSizes.Y = Y; 00190 FullDimensionSizes.Z = Z; 00191 00192 ReducedDimensionSizes.X = ((X/2) + 1); 00193 ReducedDimensionSizes.Y = Y; 00194 ReducedDimensionSizes.Z = Z; 00195 00196 00197 sensor_mask_ind_size = HDF5_InputFile.GetDatasetElementCount(sensor_mask_index_Name); 00198 00199 // flags 00200 HDF5_InputFile.ReadCompleteDataset(ux_source_flag_Name , ScalarSizes, &ux_source_flag); 00201 HDF5_InputFile.ReadCompleteDataset(uy_source_flag_Name , ScalarSizes, &uy_source_flag); 00202 HDF5_InputFile.ReadCompleteDataset(uz_source_flag_Name , ScalarSizes, &uz_source_flag); 00203 HDF5_InputFile.ReadCompleteDataset(transducer_source_flag_Name, ScalarSizes, &transducer_source_flag); 00204 00205 HDF5_InputFile.ReadCompleteDataset(p_source_flag_Name , ScalarSizes, &p_source_flag); 00206 HDF5_InputFile.ReadCompleteDataset(p0_source_flag_Name , ScalarSizes, &p0_source_flag); 00207 00208 HDF5_InputFile.ReadCompleteDataset(nonuniform_grid_flag_Name , ScalarSizes, &nonuniform_grid_flag); 00209 HDF5_InputFile.ReadCompleteDataset(absorbing_flag_Name , ScalarSizes, &absorbing_flag); 00210 HDF5_InputFile.ReadCompleteDataset(nonlinear_flag_Name , ScalarSizes, &nonlinear_flag); 00211 00212 00213 00214 //--- Vector sizes ---// 00215 if (transducer_source_flag == 0) transducer_source_input_size = 0; 00216 else { 00217 transducer_source_input_size = HDF5_InputFile.GetDatasetElementCount(transducer_source_input_Name); 00218 } 00219 00220 if ((transducer_source_flag > 0) || (ux_source_flag > 0) || (uy_source_flag > 0) || (uz_source_flag > 0)){ 00221 00222 u_source_index_size = HDF5_InputFile.GetDatasetElementCount(u_source_index_Name); 00223 } 00224 00225 00226 //-- uxyz_source_flags --// 00227 if ((ux_source_flag > 0) || (uy_source_flag > 0) || (uz_source_flag > 0)){ 00228 00229 HDF5_InputFile.ReadCompleteDataset(u_source_many_Name, ScalarSizes, &u_source_many); 00230 HDF5_InputFile.ReadCompleteDataset(u_source_mode_Name, ScalarSizes, &u_source_mode); 00231 } else{ 00232 u_source_many = 0; 00233 u_source_mode = 0; 00234 } 00235 00236 00237 //-- p_source_flag --// 00238 if (p_source_flag != 0) { 00239 HDF5_InputFile.ReadCompleteDataset(p_source_many_Name, ScalarSizes, &p_source_many); 00240 HDF5_InputFile.ReadCompleteDataset(p_source_mode_Name, ScalarSizes, &p_source_mode); 00241 00242 p_source_index_size = HDF5_InputFile.GetDatasetElementCount(p_source_index_Name); 00243 00244 } else{ 00245 p_source_mode = 0; 00246 p_source_many = 0; 00247 p_source_index_size = 0; 00248 } 00249 00250 00251 // absorb flag 00252 if (absorbing_flag != 0) { 00253 HDF5_InputFile.ReadCompleteDataset(alpha_power_Name, ScalarSizes, &alpha_power); 00254 if (alpha_power == 1.0f) { 00255 fprintf(stderr,"%s", Parameters_ERR_FMT_Illegal_alpha_power_value); 00256 PrintUsageAndExit(); 00257 } 00258 00259 alpha_coeff_scalar_flag = HDF5_InputFile.GetDatasetDimensionSizes(alpha_coeff_Name) == ScalarSizes; 00260 if (alpha_coeff_scalar_flag) HDF5_InputFile.ReadCompleteDataset(alpha_coeff_Name, ScalarSizes, &alpha_coeff_scalar); 00261 } 00262 00263 00264 c0_scalar_flag = HDF5_InputFile.GetDatasetDimensionSizes(c0_Name) == ScalarSizes; 00265 if (c0_scalar_flag) HDF5_InputFile.ReadCompleteDataset(c0_Name, ScalarSizes, &c0_scalar); 00266 00267 if (nonlinear_flag) { 00268 BonA_scalar_flag = HDF5_InputFile.GetDatasetDimensionSizes(BonA_Name) == ScalarSizes; 00269 if (BonA_scalar_flag) HDF5_InputFile.ReadCompleteDataset(BonA_Name, ScalarSizes, &BonA_scalar); 00270 } 00271 00272 rho0_scalar_flag = HDF5_InputFile.GetDatasetDimensionSizes(rho0_Name) == ScalarSizes; 00273 if ( rho0_scalar_flag) { 00274 HDF5_InputFile.ReadCompleteDataset(rho0_Name, ScalarSizes, &rho0_scalar); 00275 HDF5_InputFile.ReadCompleteDataset(rho0_sgx_Name, ScalarSizes, &rho0_sgx_scalar); 00276 HDF5_InputFile.ReadCompleteDataset(rho0_sgy_Name, ScalarSizes, &rho0_sgy_scalar); 00277 HDF5_InputFile.ReadCompleteDataset(rho0_sgz_Name, ScalarSizes, &rho0_sgz_scalar); 00278 } 00279 00280 00281 00282 }// end of ReadScalarsFromMatlabInputFile 00283 //------------------------------------------------------------------------------ 00284 00285 00290 void TParameters::SaveScalarsToHDF5File(THDF5_File & HDF5_OutputFile){ 00291 00292 00293 // Write dimension sizes 00294 HDF5_OutputFile.WriteScalarValue(Nx_Name ,(long) FullDimensionSizes.X); 00295 HDF5_OutputFile.WriteScalarValue(Ny_Name ,(long) FullDimensionSizes.Y); 00296 HDF5_OutputFile.WriteScalarValue(Nz_Name ,(long) FullDimensionSizes.Z); 00297 00298 HDF5_OutputFile.WriteScalarValue(Nt_Name ,(long) Nt); 00299 00300 HDF5_OutputFile.WriteScalarValue(dt_Name , dt); 00301 HDF5_OutputFile.WriteScalarValue(dx_Name , dx); 00302 HDF5_OutputFile.WriteScalarValue(dy_Name , dy); 00303 HDF5_OutputFile.WriteScalarValue(dz_Name , dz); 00304 00305 HDF5_OutputFile.WriteScalarValue(c_ref_Name, c_ref); 00306 00307 00308 HDF5_OutputFile.WriteScalarValue(pml_x_size_Name , pml_x_size); 00309 HDF5_OutputFile.WriteScalarValue(pml_y_size_Name , pml_y_size); 00310 HDF5_OutputFile.WriteScalarValue(pml_z_size_Name , pml_z_size); 00311 00312 HDF5_OutputFile.WriteScalarValue(pml_x_alpha_Name, pml_x_alpha); 00313 HDF5_OutputFile.WriteScalarValue(pml_y_alpha_Name, pml_y_alpha); 00314 HDF5_OutputFile.WriteScalarValue(pml_z_alpha_Name, pml_z_alpha); 00315 00316 00317 00318 HDF5_OutputFile.WriteScalarValue(ux_source_flag_Name , ux_source_flag); 00319 HDF5_OutputFile.WriteScalarValue(uy_source_flag_Name , uy_source_flag); 00320 HDF5_OutputFile.WriteScalarValue(uz_source_flag_Name , uz_source_flag); 00321 HDF5_OutputFile.WriteScalarValue(transducer_source_flag_Name , transducer_source_flag); 00322 00323 HDF5_OutputFile.WriteScalarValue(p_source_flag_Name , p_source_flag); 00324 HDF5_OutputFile.WriteScalarValue(p0_source_flag_Name , p0_source_flag); 00325 00326 HDF5_OutputFile.WriteScalarValue(nonuniform_grid_flag_Name, nonuniform_grid_flag); 00327 HDF5_OutputFile.WriteScalarValue(absorbing_flag_Name , absorbing_flag); 00328 HDF5_OutputFile.WriteScalarValue(nonlinear_flag_Name , nonlinear_flag); 00329 00330 00331 //-- uxyz_source_flags --// 00332 if ((ux_source_flag > 0) || (uy_source_flag > 0) || (uz_source_flag > 0)){ 00333 00334 HDF5_OutputFile.WriteScalarValue(u_source_many_Name , u_source_many); 00335 HDF5_OutputFile.WriteScalarValue(u_source_mode_Name , u_source_mode); 00336 } 00337 00338 //-- p_source_flag --// 00339 if (p_source_flag != 0) { 00340 00341 HDF5_OutputFile.WriteScalarValue(p_source_many_Name , p_source_many); 00342 HDF5_OutputFile.WriteScalarValue(p_source_mode_Name , p_source_mode); 00343 00344 } 00345 00346 00347 // absorb flag 00348 if (absorbing_flag != 0) { 00349 HDF5_OutputFile.WriteScalarValue(alpha_power_Name, alpha_power); 00350 00351 } 00352 00353 00354 00355 }// end of SaveScalarsToHDF5File 00356 00357 00358 00359 00360 //----------------------------------------------------------------------------// 00361 // Implementation // 00362 // protected methods // 00363 //----------------------------------------------------------------------------// 00364 00365 00369 TParameters::TParameters() : 00370 HDF5_InputFile(), HDF5_OutputFile(), HDF5_FileHeader(), 00371 CommandLinesParameters(), 00372 Nt(0),dt(0.0f), 00373 dx(0.0f), dy(0.0f), dz(0.0f), 00374 c_ref(0.0f), alpha_power(0.0f), 00375 FullDimensionSizes(0,0,0), ReducedDimensionSizes(0,0,0), 00376 sensor_mask_ind_size (0), u_source_index_size(0), p_source_index_size(0), transducer_source_input_size(0), 00377 ux_source_flag(0), uy_source_flag(0), uz_source_flag(0), 00378 p_source_flag(0), p0_source_flag(0), transducer_source_flag(0), 00379 u_source_many(0), u_source_mode(0), p_source_mode(0), p_source_many(0), 00380 nonuniform_grid_flag(0), absorbing_flag(0), nonlinear_flag(0), 00381 pml_x_size(0), pml_y_size(0), pml_z_size(0), 00382 alpha_coeff_scalar_flag(false), alpha_coeff_scalar(0.0f), 00383 c0_scalar_flag(false), c0_scalar(0.0f), 00384 absorb_eta_scalar(0.0f), absorb_tau_scalar (0.0f), 00385 BonA_scalar_flag(false), BonA_scalar (0.0f), 00386 rho0_scalar_flag(false), rho0_scalar(0.0f), rho0_sgx_scalar(0.0f), rho0_sgy_scalar(0.0f), rho0_sgz_scalar(0.0f) 00387 00388 { 00389 00390 }// end of TFFT1DParameters 00391 //------------------------------------------------------------------------------ 00392 00393 00394 00395 //----------------------------------------------------------------------------// 00396 // Implementation // 00397 // private methods // 00398 //----------------------------------------------------------------------------// 00399 00403 void TParameters::PrintUsageAndExit(){ 00404 00405 CommandLinesParameters.PrintUsageAndExit(); 00406 00407 00408 }// end of PrintUsage 00409 //------------------------------------------------------------------------------ 00410