diff --git a/include/Profile.h b/include/Profile.h index f912b5db7a..ca6ad14621 100644 --- a/include/Profile.h +++ b/include/Profile.h @@ -20,6 +20,7 @@ void Aux_Error( const char *File, const int Line, const char *Func, const char * // Allocated : Whether or not member arrays such as Radius[] have been allocated // Radius : Radial coordinate at each bin // Data : Profile data at each bin +// Data_Sigma : Standard deviation of profile data at each bin // Weight : Total weighting at each bin // NCell : Number of cells at each bin // @@ -42,6 +43,7 @@ struct Profile_t double *Radius; double *Data; + double *Data_Sigma; double *Weight; long *NCell; @@ -57,12 +59,13 @@ struct Profile_t Profile_t() { - NBin = -1; - Allocated = false; - Radius = NULL; - Data = NULL; - Weight = NULL; - NCell = NULL; + NBin = -1; + Allocated = false; + Radius = NULL; + Data = NULL; + Data_Sigma = NULL; + Weight = NULL; + NCell = NULL; } // METHOD : Profile_t @@ -92,7 +95,7 @@ struct Profile_t // // Parameter : None // - // Return : Radius[], Data[], Weight[], NCell[], Allocated + // Return : Radius[], Data[], Data_Sigma[], Weight[], NCell[], Allocated //=================================================================================== void AllocateMemory() { @@ -103,10 +106,11 @@ struct Profile_t // --> allows the same Profile_t object to be reused without the need to manually free memory first if ( Allocated ) FreeMemory(); - Radius = new double [NBin]; - Data = new double [NBin]; - Weight = new double [NBin]; - NCell = new long [NBin]; + Radius = new double [NBin]; + Data = new double [NBin]; + Data_Sigma = new double [NBin]; + Weight = new double [NBin]; + NCell = new long [NBin]; Allocated = true; @@ -122,15 +126,16 @@ struct Profile_t // // Parameter : None // - // Return : Radius[], Data[], Weight[], NCell[] + // Return : Radius[], Data[], Data_Sigma[], Weight[], NCell[] //=================================================================================== void FreeMemory() { - delete [] Radius; Radius = NULL; - delete [] Data; Data = NULL; - delete [] Weight; Weight = NULL; - delete [] NCell; NCell = NULL; + delete [] Radius; Radius = NULL; + delete [] Data; Data = NULL; + delete [] Data_Sigma; Data_Sigma = NULL; + delete [] Weight; Weight = NULL; + delete [] NCell; NCell = NULL; Allocated = false; diff --git a/include/Prototype.h b/include/Prototype.h index af1ba95c9b..f45ee37035 100644 --- a/include/Prototype.h +++ b/include/Prototype.h @@ -41,7 +41,7 @@ int Aux_CountRow( const char *FileName ); void Aux_ComputeProfile( Profile_t *Prof[], const double Center[], const double r_max_input, const double dr_min, const bool LogBin, const double LogBinRatio, const bool RemoveEmpty, const long TVarBitIdx[], const int NProf, const int MinLv, const int MaxLv, const PatchType_t PatchType, - const double PrepTimeIn ); + const double PrepTimeIn, const bool GetSigma ); void Aux_FindExtrema( Extrema_t *Extrema, const ExtremaMode_t Mode, const int MinLv, const int MaxLv, const PatchType_t PatchType ); void Aux_FindWeightedAverageCenter( double WeightedAverageCenter[], const double Center_ref[], const double MaxR, const double MinWD, diff --git a/src/Auxiliary/Aux_ComputeProfile.cpp b/src/Auxiliary/Aux_ComputeProfile.cpp index b408e6a61e..b417172a9f 100644 --- a/src/Auxiliary/Aux_ComputeProfile.cpp +++ b/src/Auxiliary/Aux_ComputeProfile.cpp @@ -11,11 +11,12 @@ extern void SetTempIntPara( const int lv, const int Sg0, const double PrepTime, // Description : Compute the average radial profile of target field(s) // // Note : 1. Results will be stored in the input "Prof" object -// --> Prof->Radius[]: Radial coordinate at each bin -// Prof->Data []: Profile data at each bin -// Prof->Weight[]: Total weighting at each bin -// Prof->NCell []: Number of cells at each bin -// Prof->NBin : Total number of bins +// --> Prof->Radius []: Radial coordinate at each bin +// Prof->Data []: Profile data at each bin +// Prof->Data_Sigma[]: Standard deviation of profile data at each bin +// Prof->Weight []: Total weighting at each bin +// Prof->NCell []: Number of cells at each bin +// Prof->NBin : Total number of bins // --> See the "Profile_t" structure defined in "include/Profile.h" for details // --> These arrays will be free'd when deleting "Prof" // 2. The exact maximum radius adopted may be slightly larger than the input "r_max" @@ -44,7 +45,7 @@ extern void SetTempIntPara( const int lv, const int Sg0, const double PrepTime, // --> Right edge of log bin n = dr_min*LogBinRatio^n // RemoveEmpty : true --> remove empty bins from the data // false --> these empty bins will still be in the profile arrays with -// Data[empty_bin]=Weight[empty_bin]=NCell[empty_bin]=0 +// Data[empty_bin]=Data_Sigma[empty_bin]=Weight[empty_bin]=NCell[empty_bin]=0 // TVarBitIdx : Bitwise indices of target variables for computing the profiles // --> Supported indices (defined in Macro.h): // HYDRO : _DENS, _MOMX, _MOMY, _MOMZ, _ENGY, _VELX, _VELY, _VELZ, _VELR, @@ -63,6 +64,8 @@ extern void SetTempIntPara( const int lv, const int Sg0, const double PrepTime, // (i.e., MinLv ~ MaxLv) and non-leaf patches only on MaxLv // PrepTimeIn : Target physical time to prepare data // --> If PrepTimeIn<0, turn off temporal interpolation and always use the most recent data +// GetSigma : true --> compute the standard deviation of profile data +// false --> set Data_Sigma=0 // // Example : const double Center[3] = { amr->BoxCenter[0], amr->BoxCenter[1], amr->BoxCenter[2] }; // const double MaxRadius = 0.5*amr->BoxSize[0]; @@ -76,12 +79,13 @@ extern void SetTempIntPara( const int lv, const int Sg0, const double PrepTime, // const int MaxLv = MAX_LEVEL; // const PatchType_t PatchType = PATCH_LEAF_PLUS_MAXNONLEAF; // const double PrepTime = -1.0; +// const bool GetSigma = false; // // Profile_t Prof_Dens, Prof_Pres; // Profile_t *Prof[] = { &Prof_Dens, &Prof_Pres }; // // Aux_ComputeProfile( Prof, Center, MaxRadius, MinBinSize, LogBin, LogBinRatio, RemoveEmptyBin, -// TVar, NProf, MinLv, MaxLv, PatchType, PrepTime ); +// TVar, NProf, MinLv, MaxLv, PatchType, PrepTime, GetSigma ); // // if ( MPI_Rank == 0 ) // { @@ -90,10 +94,10 @@ extern void SetTempIntPara( const int lv, const int Sg0, const double PrepTime, // char Filename[MAX_STRING]; // sprintf( Filename, "Profile%d.txt", p ); // FILE *File = fopen( Filename, "w" ); -// fprintf( File, "#%19s %21s %21s %10s\n", "Radius", "Data", "Weight", "Cells" ); +// fprintf( File, "#%19s %21s %21s %21s %10s\n", "Radius", "Data", "Data_Sigma", "Weight", "Cells" ); // for (int b=0; bNBin; b++) -// fprintf( File, "%20.14e %21.14e %21.14e %10ld\n", -// Prof[p]->Radius[b], Prof[p]->Data[b], Prof[p]->Weight[b], Prof[p]->NCell[b] ); +// fprintf( File, "%20.14e %21.14e %21.14e %21.14e %10ld\n", +// Prof[p]->Radius[b], Prof[p]->Data[b], Prof[p]->Data_Sigma[b], Prof[p]->Weight[b], Prof[p]->NCell[b] ); // fclose( File ); // } // } @@ -103,7 +107,7 @@ extern void SetTempIntPara( const int lv, const int Sg0, const double PrepTime, void Aux_ComputeProfile( Profile_t *Prof[], const double Center[], const double r_max_input, const double dr_min, const bool LogBin, const double LogBinRatio, const bool RemoveEmpty, const long TVarBitIdx[], const int NProf, const int MinLv, const int MaxLv, const PatchType_t PatchType, - const double PrepTimeIn ) + const double PrepTimeIn, const bool GetSigma ) { // check @@ -213,21 +217,23 @@ void Aux_ComputeProfile( Profile_t *Prof[], const double Center[], const double const int NT = 1; # endif - double ***OMP_Data=NULL, ***OMP_Weight=NULL; + double ***OMP_Data=NULL, ***OMP_Data_Sigma=NULL, ***OMP_Weight=NULL; long ***OMP_NCell=NULL; - Aux_AllocateArray3D( OMP_Data, NProf, NT, Prof[0]->NBin ); - Aux_AllocateArray3D( OMP_Weight, NProf, NT, Prof[0]->NBin ); - Aux_AllocateArray3D( OMP_NCell, NProf, NT, Prof[0]->NBin ); + Aux_AllocateArray3D( OMP_Data, NProf, NT, Prof[0]->NBin ); + Aux_AllocateArray3D( OMP_Data_Sigma, NProf, NT, Prof[0]->NBin ); + Aux_AllocateArray3D( OMP_Weight, NProf, NT, Prof[0]->NBin ); + Aux_AllocateArray3D( OMP_NCell, NProf, NT, Prof[0]->NBin ); // initialize profile arrays for (int p=0; pNBin; b++) { - OMP_Data [p][t][b] = 0.0; - OMP_Weight[p][t][b] = 0.0; - OMP_NCell [p][t][b] = 0; + OMP_Data [p][t][b] = 0.0; + OMP_Data_Sigma[p][t][b] = 0.0; + OMP_Weight [p][t][b] = 0.0; + OMP_NCell [p][t][b] = 0; } real (*Patch_Data)[8][PS1][PS1][PS1] = new real [NT][8][PS1][PS1][PS1]; // field data of each cell @@ -538,6 +544,9 @@ void Aux_ComputeProfile( Profile_t *Prof[], const double Center[], const double OMP_Data [p][TID][bin] += Patch_Data[TID][LocalID][k][j][i]*Weight; OMP_Weight[p][TID][bin] += Weight; OMP_NCell [p][TID][bin] ++; + if ( GetSigma ) + OMP_Data_Sigma[p][TID][bin] += SQR( Patch_Data[TID][LocalID][k][j][i] )*Weight; + } // i,j,k } // for (int LocalID=0; LocalID<8; LocalID++) } // for (int p=0; pNBin; b++) { - Prof[p]->Data [b] = OMP_Data [p][0][b]; - Prof[p]->Weight[b] = OMP_Weight[p][0][b]; - Prof[p]->NCell [b] = OMP_NCell [p][0][b]; + Prof[p]->Data [b] = OMP_Data [p][0][b]; + Prof[p]->Weight[b] = OMP_Weight[p][0][b]; + Prof[p]->NCell [b] = OMP_NCell [p][0][b]; + if ( GetSigma ) + Prof[p]->Data_Sigma[b] = OMP_Data_Sigma[p][0][b]; + } for (int t=1; tData [b] += OMP_Data [p][t][b]; Prof[p]->Weight[b] += OMP_Weight[p][t][b]; Prof[p]->NCell [b] += OMP_NCell [p][t][b]; + if ( GetSigma ) + Prof[p]->Data_Sigma[b] += OMP_Data_Sigma[p][t][b]; } } // free per-thread arrays Aux_DeallocateArray3D( OMP_Data ); + Aux_DeallocateArray3D( OMP_Data_Sigma ); Aux_DeallocateArray3D( OMP_Weight ); Aux_DeallocateArray3D( OMP_NCell ); @@ -593,16 +608,21 @@ void Aux_ComputeProfile( Profile_t *Prof[], const double Center[], const double { if ( MPI_Rank == 0 ) { - MPI_Reduce( MPI_IN_PLACE, Prof[p]->Data, Prof[p]->NBin, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); - MPI_Reduce( MPI_IN_PLACE, Prof[p]->Weight, Prof[p]->NBin, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); - MPI_Reduce( MPI_IN_PLACE, Prof[p]->NCell , Prof[p]->NBin, MPI_LONG, MPI_SUM, 0, MPI_COMM_WORLD ); + MPI_Reduce( MPI_IN_PLACE, Prof[p]->Data, Prof[p]->NBin, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); + MPI_Reduce( MPI_IN_PLACE, Prof[p]->Weight, Prof[p]->NBin, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); + MPI_Reduce( MPI_IN_PLACE, Prof[p]->NCell, Prof[p]->NBin, MPI_LONG, MPI_SUM, 0, MPI_COMM_WORLD ); + if ( GetSigma ) + MPI_Reduce( MPI_IN_PLACE, Prof[p]->Data_Sigma, Prof[p]->NBin, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); + } else { - MPI_Reduce( Prof[p]->Data, NULL, Prof[p]->NBin, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); - MPI_Reduce( Prof[p]->Weight, NULL, Prof[p]->NBin, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); - MPI_Reduce( Prof[p]->NCell, NULL, Prof[p]->NBin, MPI_LONG, MPI_SUM, 0, MPI_COMM_WORLD ); + MPI_Reduce( Prof[p]->Data, NULL, Prof[p]->NBin, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); + MPI_Reduce( Prof[p]->Weight, NULL, Prof[p]->NBin, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); + MPI_Reduce( Prof[p]->NCell, NULL, Prof[p]->NBin, MPI_LONG, MPI_SUM, 0, MPI_COMM_WORLD ); + if ( GetSigma ) + MPI_Reduce( Prof[p]->Data_Sigma, NULL, Prof[p]->NBin, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); } } # endif @@ -615,7 +635,18 @@ void Aux_ComputeProfile( Profile_t *Prof[], const double Center[], const double for (int b=0; bNBin; b++) { // skip empty bins since both their data and weight are zero - if ( Prof[p]->NCell[b] > 0L ) Prof[p]->Data[b] /= Prof[p]->Weight[b]; + if ( Prof[p]->NCell[b] > 0L ) + { + Prof[p]->Data [b] /= Prof[p]->Weight[b]; + Prof[p]->Data_Sigma[b] = ( GetSigma ) + ? Prof[p]->Data_Sigma[b]/Prof[p]->Weight[b] - SQR( Prof[p]->Data[b] ) + : 0.0; + + if ( Prof[p]->Data_Sigma[b] < 0.0 ) + Aux_Error( ERROR_INFO, "Prof[%d]->Data_Sigma[%d] = %14.7e < 0.0 !!\n", p, b, Prof[p]->Data_Sigma[b] ); + else + Prof[p]->Data_Sigma[b] = sqrt( Prof[p]->Data_Sigma[b] ); + } } } @@ -623,9 +654,10 @@ void Aux_ComputeProfile( Profile_t *Prof[], const double Center[], const double // broadcast data to all ranks for (int p=0; pData, Prof[p]->NBin, MPI_DOUBLE, 0, MPI_COMM_WORLD ); - MPI_Bcast( Prof[p]->Weight, Prof[p]->NBin, MPI_DOUBLE, 0, MPI_COMM_WORLD ); - MPI_Bcast( Prof[p]->NCell, Prof[p]->NBin, MPI_LONG, 0, MPI_COMM_WORLD ); + MPI_Bcast( Prof[p]->Data, Prof[p]->NBin, MPI_DOUBLE, 0, MPI_COMM_WORLD ); + MPI_Bcast( Prof[p]->Data_Sigma, Prof[p]->NBin, MPI_DOUBLE, 0, MPI_COMM_WORLD ); + MPI_Bcast( Prof[p]->Weight, Prof[p]->NBin, MPI_DOUBLE, 0, MPI_COMM_WORLD ); + MPI_Bcast( Prof[p]->NCell, Prof[p]->NBin, MPI_LONG, 0, MPI_COMM_WORLD ); } @@ -650,10 +682,11 @@ void Aux_ComputeProfile( Profile_t *Prof[], const double Center[], const double for (int p=0; pRadius[b_up_ms] = Prof[p]->Radius[b_up]; - Prof[p]->Data [b_up_ms] = Prof[p]->Data [b_up]; - Prof[p]->Weight[b_up_ms] = Prof[p]->Weight[b_up]; - Prof[p]->NCell [b_up_ms] = Prof[p]->NCell [b_up]; + Prof[p]->Radius [b_up_ms] = Prof[p]->Radius [b_up]; + Prof[p]->Data [b_up_ms] = Prof[p]->Data [b_up]; + Prof[p]->Data_Sigma[b_up_ms] = Prof[p]->Data_Sigma[b_up]; + Prof[p]->Weight [b_up_ms] = Prof[p]->Weight [b_up]; + Prof[p]->NCell [b_up_ms] = Prof[p]->NCell [b_up]; } } diff --git a/src/SourceTerms/Deleptonization/CPU_Src_Deleptonization.cpp b/src/SourceTerms/Deleptonization/CPU_Src_Deleptonization.cpp index 0f6289667d..a941dc34b9 100644 --- a/src/SourceTerms/Deleptonization/CPU_Src_Deleptonization.cpp +++ b/src/SourceTerms/Deleptonization/CPU_Src_Deleptonization.cpp @@ -188,7 +188,7 @@ void Src_WorkBeforeMajorFunc_Deleptonization( const int lv, const double TimeNew for (int v=0; v