Skip to content

Fix a bug with rescaling the soliton density profile in the HaloMerger test problem #137

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
hsinhaoHHuang wants to merge 6 commits into
base: psidm
Choose a base branch
from
Open

Fix a bug with rescaling the soliton density profile in the HaloMerger test problem #137

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
3fdb2b2
Fix bug of rescaling soliton density profile
hsinhaoHHuang Oct 8, 2024
a587362
Use ScaleC to rescale soliton density profile
hsinhaoHHuang Oct 8, 2024
9368c5c
Merge branch 'psidm' into psidm-FixBugSolitonRescale
hsinhaoHHuang Oct 25, 2024
e161256
Update comments according to the review
hsinhaoHHuang Oct 25, 2024
9e1041f
Merge branch 'psidm' into psidm-FixBugSolitonRescale
hsinhaoHHuang Jan 27, 2025
3bbe099
Merge branch 'psidm' into psidm-FixBugSolitonRescale
hsinhaoHHuang Apr 7, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions example/test_problem/ELBDM/HaloMerger/Input__TestProb_Soliton
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
36 changes: 30 additions & 6 deletions src/TestProblem/ELBDM/HaloMerger/Init_TestProb_ELBDM_HaloMerger.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -326,6 +328,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 )
Expand Down Expand Up @@ -353,6 +357,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)

Expand All @@ -373,6 +378,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 )
{
Expand Down Expand Up @@ -401,6 +407,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 )
{
Expand Down Expand Up @@ -690,6 +697,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
Expand All @@ -699,7 +715,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 )
Expand All @@ -711,13 +727,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;
Expand Down Expand Up @@ -1098,14 +1119,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<HaloMerger_Soliton_Num; index_soliton++)
Aux_Message( stdout, " %7d %35s %16d %16d %16.6e %16.6e\n",
Aux_Message( stdout, " %7d %35s %16d %16d %16.6e %16.6e %16.6e %16.6e\n",
index_soliton+1, HaloMerger_Soliton_DensProf_Filename[index_soliton], HaloMerger_Soliton_DensProf_NBin[index_soliton],
HaloMerger_Soliton_DensProf_Rescale[index_soliton], HaloMerger_Soliton_DensProf_ScaleL[index_soliton], HaloMerger_Soliton_DensProf_ScaleD[index_soliton] );
HaloMerger_Soliton_DensProf_Rescale[index_soliton], HaloMerger_Soliton_DensProf_PhyConst[index_soliton],
HaloMerger_Soliton_DensProf_ScaleC[index_soliton], HaloMerger_Soliton_DensProf_ScaleL[index_soliton], HaloMerger_Soliton_DensProf_ScaleD[index_soliton] );

} // if ( HaloMerger_Soliton_InitMode == 1 )
else if ( HaloMerger_Soliton_InitMode == 2 )
Expand Down Expand Up @@ -1410,6 +1432,8 @@ void End_HaloMerger()
delete [] HaloMerger_Soliton_DensProf;
delete [] HaloMerger_Soliton_DensProf_NBin;
delete [] HaloMerger_Soliton_DensProf_Rescale;
delete [] HaloMerger_Soliton_DensProf_PhyConst;
delete [] HaloMerger_Soliton_DensProf_ScaleC;
delete [] HaloMerger_Soliton_DensProf_ScaleL;
delete [] HaloMerger_Soliton_DensProf_ScaleD;

Expand Down