diff --git a/example/test_problem/ELBDM/HaloMerger/Input__TestProb_Soliton b/example/test_problem/ELBDM/HaloMerger/Input__TestProb_Soliton index 869f4e061..5ae5e55c9 100644 --- a/example/test_problem/ELBDM/HaloMerger/Input__TestProb_Soliton +++ b/example/test_problem/ELBDM/HaloMerger/Input__TestProb_Soliton @@ -10,6 +10,7 @@ HaloMerger_Soliton_1_VelocityY 0.5 HaloMerger_Soliton_1_VelocityZ 0.0 HaloMerger_Soliton_1_DensProf_Filename SolitonDensityProfile_Lambda0.0 # filename of the density profile table for the 1st soliton (HaloMerger_Soliton_InitMode == 1 only) HaloMerger_Soliton_1_DensProf_Rescale 1 # whether to scale the density profile table for the 1st soliton (HaloMerger_Soliton_InitMode == 1 only) [1] +HaloMerger_Soliton_1_DensProf_PhyConst 1.0 # value of the dimensional constant 4*pi*G*(ELBDM_MASS/hbar)^2 in the units of density profile (HaloMerger_Soliton_InitMode == 1 only) (<=0.0=same as in the simulation units) [-1] HaloMerger_Soliton_1_OuterSlope -8.0 # outer slope of the analytical density profile of the 1st soliton (HaloMerger_Soliton_InitMode == 2 only) [-8.0] # 2nd soliton HaloMerger_Soliton_2_CoreRadius -1.0 @@ -22,4 +23,5 @@ HaloMerger_Soliton_2_VelocityY -0.5 HaloMerger_Soliton_2_VelocityZ 0.0 HaloMerger_Soliton_2_DensProf_Filename SolitonDensityProfile HaloMerger_Soliton_2_DensProf_Rescale 0 +HaloMerger_Soliton_2_DensProf_PhyConst -1.0 HaloMerger_Soliton_2_OuterSlope -8.0 diff --git a/src/TestProblem/ELBDM/HaloMerger/Init_TestProb_ELBDM_HaloMerger.cpp b/src/TestProblem/ELBDM/HaloMerger/Init_TestProb_ELBDM_HaloMerger.cpp index 91c4cd7f9..0e0d06d30 100644 --- a/src/TestProblem/ELBDM/HaloMerger/Init_TestProb_ELBDM_HaloMerger.cpp +++ b/src/TestProblem/ELBDM/HaloMerger/Init_TestProb_ELBDM_HaloMerger.cpp @@ -42,6 +42,8 @@ static double *HaloMerger_Soliton_OuterSlope = NULL; // ou static char (*HaloMerger_Soliton_DensProf_Filename)[MAX_STRING] = NULL; // filename of the density profile table of each soliton static int *HaloMerger_Soliton_DensProf_NBin = NULL; // number of bins of the density profile table static bool *HaloMerger_Soliton_DensProf_Rescale = NULL; // whether to scale the density profile table of each soliton +static double *HaloMerger_Soliton_DensProf_PhyConst = NULL; // value of the dimensional constant 4*pi*G*(ELBDM_MASS/hbar)^2 in the units of density profile +static double *HaloMerger_Soliton_DensProf_ScaleC = NULL; // ratio between the values of 4*pi*G*(ELBDM_MASS/hbar)^2 in the units of the simulation and input table static double *HaloMerger_Soliton_DensProf_ScaleL = NULL; // L/D: length/density scale factors of each soliton static double *HaloMerger_Soliton_DensProf_ScaleD = NULL; // (defined as the ratio between the core radii/peak // density of the target and reference soliton profiles) @@ -322,6 +324,8 @@ void SetParameter() HaloMerger_Soliton_DensProf = new double* [HaloMerger_Soliton_Num]; HaloMerger_Soliton_DensProf_NBin = new int [HaloMerger_Soliton_Num]; HaloMerger_Soliton_DensProf_Rescale = new bool [HaloMerger_Soliton_Num]; + HaloMerger_Soliton_DensProf_PhyConst = new double [HaloMerger_Soliton_Num]; + HaloMerger_Soliton_DensProf_ScaleC = new double [HaloMerger_Soliton_Num]; HaloMerger_Soliton_DensProf_ScaleL = new double [HaloMerger_Soliton_Num]; HaloMerger_Soliton_DensProf_ScaleD = new double [HaloMerger_Soliton_Num]; } // if ( HaloMerger_Soliton_InitMode == 1 ) @@ -349,6 +353,7 @@ void SetParameter() char HaloMerger_Soliton_i_DensProf_Filename[MAX_STRING]; // filename of the density profile table for the i-th soliton (HaloMerger_Soliton_InitMode == 1 only) char HaloMerger_Soliton_i_DensProf_Rescale[MAX_STRING]; // whether to scale the density profile table for the i-th soliton (HaloMerger_Soliton_InitMode == 1 only) + char HaloMerger_Soliton_i_DensProf_PhyConst[MAX_STRING]; // value of the dimensional constant 4*pi*G*(ELBDM_MASS/hbar)^2 in the units of density profile for the i-th soliton (HaloMerger_Soliton_InitMode == 1 only) (<=0.0=same as in the simulation units) char HaloMerger_Soliton_i_OuterSlope[MAX_STRING]; // outer slope of the analytical density profile of the i-th soliton (HaloMerger_Soliton_InitMode == 2 only) @@ -369,6 +374,7 @@ void SetParameter() { sprintf( HaloMerger_Soliton_i_DensProf_Filename, "HaloMerger_Soliton_%d_DensProf_Filename", index_soliton_input ); sprintf( HaloMerger_Soliton_i_DensProf_Rescale, "HaloMerger_Soliton_%d_DensProf_Rescale", index_soliton_input ); + sprintf( HaloMerger_Soliton_i_DensProf_PhyConst, "HaloMerger_Soliton_%d_DensProf_PhyConst", index_soliton_input ); } // if ( HaloMerger_Soliton_InitMode == 1 ) else if ( HaloMerger_Soliton_InitMode == 2 ) { @@ -397,6 +403,7 @@ void SetParameter() { ReadPara_Soliton->Add( HaloMerger_Soliton_i_DensProf_Filename, HaloMerger_Soliton_DensProf_Filename[index_soliton], NoDef_str, Useless_str, Useless_str ); ReadPara_Soliton->Add( HaloMerger_Soliton_i_DensProf_Rescale, &HaloMerger_Soliton_DensProf_Rescale[index_soliton], true, Useless_bool, Useless_bool ); + ReadPara_Soliton->Add( HaloMerger_Soliton_i_DensProf_PhyConst, &HaloMerger_Soliton_DensProf_PhyConst[index_soliton], -1.0, NoMin_double, NoMax_double ); } // if ( HaloMerger_Soliton_InitMode == 1 ) else if ( HaloMerger_Soliton_InitMode == 2 ) { @@ -686,6 +693,15 @@ void SetParameter() if ( CoreRadiusRef == NULL_REAL ) Aux_Error( ERROR_INFO, "cannot determine the reference core radius !!\n" ); + // set the default physical constant + if ( HaloMerger_Soliton_DensProf_PhyConst[index_soliton] <= 0.0 ) + HaloMerger_Soliton_DensProf_PhyConst[index_soliton] = 4.0*M_PI*NEWTON_G*SQR(ELBDM_ETA); + + // set the ratio of the physical constant between simulation and input table + HaloMerger_Soliton_DensProf_ScaleC[index_soliton] = (4.0*M_PI*NEWTON_G*SQR(ELBDM_ETA))/HaloMerger_Soliton_DensProf_PhyConst[index_soliton]; + + // rescale the soliton density profile by following the scaling relation: + // (4*pi*G*ELBDM_ETA^2)*(CoreRho)*(CoreRadius^4) = a dimensionless constnat -> (ScaleC)*(ScaleD)*(ScaleL^4) = 1 if ( HaloMerger_Soliton_DensProf_Rescale[index_soliton] ) { // evaluate the scale factors of each soliton @@ -695,7 +711,7 @@ void SetParameter() { // overwrite the core radius by the value calculated from the peak density HaloMerger_Soliton_DensProf_ScaleD[index_soliton] = HaloMerger_Soliton_CoreRho[index_soliton] / DensRef[0]; - HaloMerger_Soliton_DensProf_ScaleL[index_soliton] = sqrt( sqrt( 1.0 / (4.0*M_PI*NEWTON_G*SQR(ELBDM_ETA)*HaloMerger_Soliton_DensProf_ScaleD[index_soliton]) ) ); + HaloMerger_Soliton_DensProf_ScaleL[index_soliton] = sqrt( sqrt( 1.0 / (HaloMerger_Soliton_DensProf_ScaleC[index_soliton]*HaloMerger_Soliton_DensProf_ScaleD[index_soliton]) ) ); HaloMerger_Soliton_CoreRadius [index_soliton] = CoreRadiusRef*HaloMerger_Soliton_DensProf_ScaleL[index_soliton]; } else // if ( HaloMerger_Soliton_CoreRho[index_soliton] > 0.0 ) @@ -707,13 +723,18 @@ void SetParameter() { // overwrite the peak density by the value calculated from the core radius HaloMerger_Soliton_DensProf_ScaleL[index_soliton] = HaloMerger_Soliton_CoreRadius[index_soliton] / CoreRadiusRef; - HaloMerger_Soliton_DensProf_ScaleD[index_soliton] = 1.0 / ( 4.0*M_PI*NEWTON_G*SQR(ELBDM_ETA)*POW4(HaloMerger_Soliton_DensProf_ScaleL[index_soliton]) ); + HaloMerger_Soliton_DensProf_ScaleD[index_soliton] = 1.0 / ( HaloMerger_Soliton_DensProf_ScaleC[index_soliton]*POW4(HaloMerger_Soliton_DensProf_ScaleL[index_soliton]) ); HaloMerger_Soliton_CoreRho [index_soliton] = HaloMerger_Soliton_DensProf_ScaleD[index_soliton]*DensRef[0]; } // if ( HaloMerger_Soliton_CoreRadius[index_soliton] <= 0.0 ) ... else } // if ( HaloMerger_Soliton_DensProf_Rescale[index_soliton] ) else { + // check the units of table are consistent with the simulation + if ( ! Mis_CompareRealValue( HaloMerger_Soliton_DensProf_PhyConst[index_soliton], 4.0*M_PI*NEWTON_G*SQR(ELBDM_ETA), NULL, false ) ) + Aux_Error( ERROR_INFO, "HaloMerger_Soliton_%d_DensProf_PhyConst (%13.6e) in input table is not consistent with the value of 4.0*M_PI*NEWTON_G*SQR(ELBDM_ETA) = %13.6e in simulation units !!\n", + index_soliton_input, HaloMerger_Soliton_DensProf_PhyConst[index_soliton], 4.0*M_PI*NEWTON_G*SQR(ELBDM_ETA) ); + // overwrite the peak density and core radius from the table HaloMerger_Soliton_DensProf_ScaleL[index_soliton] = 1.0; HaloMerger_Soliton_DensProf_ScaleD[index_soliton] = 1.0; @@ -1094,14 +1115,15 @@ void SetParameter() if ( HaloMerger_Soliton_InitMode == 1 ) { Aux_Message( stdout, "\n soliton density profile information:\n" ); - Aux_Message( stdout, " %7s %35s %16s %16s %16s %16s\n", + Aux_Message( stdout, " %7s %35s %16s %16s %16s %16s %16s %16s\n", "ID", "DensProf_Filename", "NBin", - "Rescale", "ScaleL", "ScaleD" ); + "Rescale", "PhyConst", "ScaleC", "ScaleL", "ScaleD" ); for (int index_soliton=0; index_soliton