Skip to content

Support analytical profile in Soliton test #484

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
vivi235711 wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Support analytical profile in Soliton test #484

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
98d52ee
Support analytical profile in Soliton test
vivi235711 Oct 7, 2025
855f6d9
Update code based on @boxianliu's comment
vivi235711 Dec 1, 2025
330e587
Typo
vivi235711 Dec 4, 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
1 change: 1 addition & 0 deletions example/test_problem/ELBDM/Soliton/Input__TestProb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,4 +5,5 @@ Soliton_CoreRadiusAll 1.0e1 # core radius for ALL solitons (in code uni
Soliton_EmptyRegion 0.0 # soliton-free region [0.0]
# -> solitons only reside in [Soliton_EmptyRegion : BOX_SIZE-Soliton_EmptyRegion]
# -> useful only when Soliton_RSeed >= 0
Soliton_InitMode 1 # soliton initialization mode (1=table of the density profile, 2=analytical function of the density profile) [1]
Soliton_DensProf_Filename SolitonDensityProfile_Lambda0.0 # filename of the input soliton density profile
1 change: 1 addition & 0 deletions example/test_problem/ELBDM/Soliton/README
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,3 +18,4 @@ Note:
Soliton_N > 1: soliton merger
--> Set Soliton_RSeed >= 0
Set Soliton_EmptyRegion > 0.0 (e.g., 1.0e2)
2. Input density profile tables should use $| \psi |^2$ units. The code automatically converts this to mass density.
180 changes: 116 additions & 64 deletions src/TestProblem/ELBDM/Soliton/Init_TestProb_ELBDM_Soliton.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ static double Soliton_CoreRadiusAll; // core radius for all
// (<=0.0 --> hard coding each soliton)
static double Soliton_EmptyRegion; // soliton-free region from the boundary
// (useful only when Soliton_RSeed>=0)
static int Soliton_InitMode; // soliton initialization mode (1=table of the density profile, 2=analytical function of the density profile)
static char Soliton_DensProf_Filename[MAX_STRING]; // filename of the reference soliton density profile

static int Soliton_DensProf_NBin; // number of radial bins of the soliton density profile
Expand All @@ -20,7 +21,7 @@ static double (*Soliton_Center)[3] = NULL; // center coordinates o
static double *Soliton_ScaleL = NULL; // L/D: length/density scale factors of each soliton
// (defined as the ratio between the core radii/peak
// density of the target and reference soliton profiles)
static double *Soliton_ScaleD = NULL;
static double *Soliton_ScaleD = NULL; // ( for Soliton_InitMode==1 only)
// =======================================================================================

static void BC( real Array[], const int ArraySize[], real fluid[], const int NVar_Flu,
Expand Down Expand Up @@ -124,6 +125,7 @@ void LoadInputTestProb( const LoadParaMode_t load_mode, ReadPara_t *ReadPara, HD
LOAD_PARA( load_mode, "Soliton_RSeed", &Soliton_RSeed, 0, NoMin_int, NoMax_int );
LOAD_PARA( load_mode, "Soliton_CoreRadiusAll", &Soliton_CoreRadiusAll, NoDef_double, NoMin_double, NoMax_double );
LOAD_PARA( load_mode, "Soliton_EmptyRegion", &Soliton_EmptyRegion, 0.0, NoMin_double, NoMax_double );
LOAD_PARA( load_mode, "Soliton_InitMode", &Soliton_InitMode, 1, 1, 2 );
LOAD_PARA( load_mode, "Soliton_DensProf_Filename", Soliton_DensProf_Filename, NoDef_str, Useless_str, Useless_str );

} // FUNCITON : LoadInputTestProb
Expand Down Expand Up @@ -171,13 +173,21 @@ void SetParameter()
if ( Soliton_CoreRadiusAll == NoDef_double )
Aux_Error( ERROR_INFO, "Runtime parameter \"Soliton_CoreRadiusAll\" is not set !!\n" );

if ( !OPT__UNIT && Soliton_InitMode == 2 )
Aux_Error( ERROR_INFO, "OPT__UNIT must be enabled for Soliton_InitMode == 2 !!\n" );

// (2) set the problem-specific derived parameters
// (2-1) allocate memory
Soliton_CoreRadius = new double [Soliton_N];
Soliton_Center = new double [Soliton_N][3];
Soliton_ScaleL = new double [Soliton_N];
Soliton_ScaleD = new double [Soliton_N];
if ( Soliton_InitMode == 1 )
{
Soliton_ScaleL = new double [Soliton_N];
Soliton_ScaleD = new double [Soliton_N];
}
else if ( Soliton_InitMode == 2 ) {}
else
Aux_Error( ERROR_INFO, "unsupported Soliton_InitMode (%d) !!\n", Soliton_InitMode );

// (2-2) soliton core radii
if ( Soliton_CoreRadiusAll > 0.0 )
Expand Down Expand Up @@ -231,41 +241,47 @@ void SetParameter()
// (3) load the reference soliton density profile and evaluate the scale factors
if ( OPT__INIT != INIT_BY_RESTART )
{
// load the reference profile
const bool RowMajor_No = false; // load data into the column-major order
const bool AllocMem_Yes = true; // allocate memory for Soliton_DensProf
const int NCol = 2; // total number of columns to load
const int Col[NCol] = {0, 1}; // target columns: (radius, density)
if ( Soliton_InitMode == 1 )
{
// load the reference profile
const bool RowMajor_No = false; // load data into the column-major order
const bool AllocMem_Yes = true; // allocate memory for Soliton_DensProf
const int NCol = 2; // total number of columns to load
const int Col[NCol] = {0, 1}; // target columns: (radius, density)

Soliton_DensProf_NBin = Aux_LoadTable( Soliton_DensProf, Soliton_DensProf_Filename, NCol, Col, RowMajor_No, AllocMem_Yes );
Soliton_DensProf_NBin = Aux_LoadTable( Soliton_DensProf, Soliton_DensProf_Filename, NCol, Col, RowMajor_No, AllocMem_Yes );


// get the core radius of the reference profile
const double *RadiusRef = Soliton_DensProf + 0*Soliton_DensProf_NBin;
const double *DensRef = Soliton_DensProf + 1*Soliton_DensProf_NBin;
const double DensCore = 0.5*DensRef[0]; // define core radius as the half-density radius
// get the core radius of the reference profile
const double *RadiusRef = Soliton_DensProf + 0*Soliton_DensProf_NBin;
const double *DensRef = Soliton_DensProf + 1*Soliton_DensProf_NBin;
const double DensCore = 0.5*DensRef[0]; // define core radius as the half-density radius

double CoreRadiusRef = NULL_REAL;
double CoreRadiusRef = NULL_REAL;

for (int b=1; b<Soliton_DensProf_NBin-1; b++)
{
if ( DensRef[b] >= DensCore && DensRef[b+1] <= DensCore )
for (int b=1; b<Soliton_DensProf_NBin-1; b++)
{
CoreRadiusRef = 0.5*( RadiusRef[b] + RadiusRef[b+1] );
break;
if ( DensRef[b] >= DensCore && DensRef[b+1] <= DensCore )
{
CoreRadiusRef = RadiusRef[b] + (DensCore - DensRef[b])*(RadiusRef[b+1]-RadiusRef[b])/(DensRef[b+1]-DensRef[b]);
break;
}
}
}

if ( CoreRadiusRef == NULL_REAL )
Aux_Error( ERROR_INFO, "cannot determine the reference core radius !!\n" );
if ( CoreRadiusRef == NULL_REAL )
Aux_Error( ERROR_INFO, "cannot determine the reference core radius !!\n" );


// evaluate the scale factors of each soliton
for (int t=0; t<Soliton_N; t++)
{
Soliton_ScaleL[t] = Soliton_CoreRadius[t] / CoreRadiusRef;
Soliton_ScaleD[t] = 1.0 / ( 4.0*M_PI*NEWTON_G*SQR(ELBDM_ETA)*POW4(Soliton_ScaleL[t]) );
}
// evaluate the scale factors of each soliton
for (int t=0; t<Soliton_N; t++)
{
Soliton_ScaleL[t] = Soliton_CoreRadius[t] / CoreRadiusRef;
Soliton_ScaleD[t] = 1.0 / ( 4.0*M_PI*NEWTON_G*SQR(ELBDM_ETA)*POW4(Soliton_ScaleL[t]) );
}
} // if ( Soliton_InitMode == 1 )
else if ( Soliton_InitMode == 2 ) {}
else
Aux_Error( ERROR_INFO, "unsupported Soliton_InitMode (%d) !!\n", Soliton_InitMode );
} // if ( OPT__INIT != INIT_BY_RESTART )


Expand Down Expand Up @@ -293,16 +309,34 @@ void SetParameter()
Aux_Message( stdout, " total number of solitons = %d\n", Soliton_N );
Aux_Message( stdout, " random seed for setting the center coord. = %d\n", Soliton_RSeed );
Aux_Message( stdout, " size of the soliton-free zone = %13.7e\n", Soliton_EmptyRegion );
Aux_Message( stdout, " density profile filename = %s\n", Soliton_DensProf_Filename );
Aux_Message( stdout, " number of bins of the density profile = %d\n", Soliton_DensProf_NBin );
if (Soliton_InitMode==1)
{
Aux_Message( stdout, " density profile filename = %s\n", Soliton_DensProf_Filename );
Aux_Message( stdout, " number of bins of the density profile = %d\n", Soliton_DensProf_NBin );
}
Aux_Message( stdout, "\n" );
Aux_Message( stdout, " Soliton info:\n" );
Aux_Message( stdout, " %7s %13s %13s %13s %13s %13s %13s\n",
"ID", "CoreRadius", "ScaleL", "ScaleD", "Center_X", "Center_Y", "Center_Z" );
for (int t=0; t<Soliton_N; t++)
Aux_Message( stdout, " %7d %13.6e %13.6e %13.6e %13.6e %13.6e %13.6e\n",
t, Soliton_CoreRadius[t], Soliton_ScaleL[t], Soliton_ScaleD[t],
Soliton_Center[t][0], Soliton_Center[t][1], Soliton_Center[t][2] );
if ( Soliton_InitMode == 1 )
{
Aux_Message( stdout, " (InitMode = 1 --> use the table of the density profile)\n" );
Aux_Message( stdout, " %7s %13s %13s %13s %13s %13s %13s\n",
"ID", "CoreRadius", "ScaleL", "ScaleD", "Center_X", "Center_Y", "Center_Z" );
for (int t=0; t<Soliton_N; t++)
Aux_Message( stdout, " %7d %13.6e %13.6e %13.6e %13.6e %13.6e %13.6e\n",
t, Soliton_CoreRadius[t], Soliton_ScaleL[t], Soliton_ScaleD[t],
Soliton_Center[t][0], Soliton_Center[t][1], Soliton_Center[t][2] );
} // if ( Soliton_InitMode == 1 )
else if ( Soliton_InitMode == 2 )
{
Aux_Message( stdout, " (InitMode = 2 --> use the analytical function of the density profile)\n" );
Aux_Message( stdout, " %7s %13s %13s %13s %13s\n", "ID", "CoreRadius", "Center_X", "Center_Y", "Center_Z" );
for (int t=0; t<Soliton_N; t++)
Aux_Message( stdout, " %7d %13.6e %13.6e %13.6e %13.6e\n",
t, Soliton_CoreRadius[t],
Soliton_Center[t][0], Soliton_Center[t][1], Soliton_Center[t][2] );
} // if ( Soliton_InitMode == 2 )
else
Aux_Error( ERROR_INFO, "unsupported Soliton_InitMode (%d) !!\n", Soliton_InitMode );
Aux_Message( stdout, "======================================================================================\n" );
}

Expand Down Expand Up @@ -333,43 +367,61 @@ void SetParameter()
void SetGridIC( real fluid[], const double x, const double y, const double z, const double Time,
const int lv, double AuxArray[] )
Copy link
Contributor

@hyschive hyschive Nov 13, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@boxianliu The format of this comment looks a bit strange and too long. You may consider:

  • Using the Add a suggestion button to highlight the differences.
  • Previewing your comment with the Preview button before submitting.

{
double r_tar, r_ref, dens_ref, dens;

const double *Table_Radius = Soliton_DensProf + 0*Soliton_DensProf_NBin; // radius
const double *Table_Density = Soliton_DensProf + 1*Soliton_DensProf_NBin; // density

double r_tar, r_ref, dens_ref;


// initialize density as zero since there may be multiple solitons
// initialize density as zero since there may be multiple solitons
fluid[DENS] = 0.0;

// loop over all solitons to get the total density
for (int t=0; t<Soliton_N; t++)
if ( Soliton_InitMode == 1 )
{
r_tar = sqrt( SQR(x-Soliton_Center[t][0]) + SQR(y-Soliton_Center[t][1]) + SQR(z-Soliton_Center[t][2]) );
const double *Table_Radius = Soliton_DensProf + 0*Soliton_DensProf_NBin; // radius
const double *Table_Density = Soliton_DensProf + 1*Soliton_DensProf_NBin; // density

// rescale radius (target radius --> reference radius)
r_ref = r_tar / Soliton_ScaleL[t];
// loop over all solitons to get the total density
for (int t=0; t<Soliton_N; t++)
{
r_tar = sqrt( SQR(x-Soliton_Center[t][0]) + SQR(y-Soliton_Center[t][1]) + SQR(z-Soliton_Center[t][2]) );

// linear interpolation
dens_ref = Mis_InterpolateFromTable( Soliton_DensProf_NBin, Table_Radius, Table_Density, r_ref );
// rescale radius (target radius --> reference radius)
r_ref = r_tar / Soliton_ScaleL[t];

if ( dens_ref == NULL_REAL )
{
if ( r_ref < Table_Radius[0] )
dens_ref = Table_Density[0];
// linear interpolation
dens_ref = Mis_InterpolateFromTable( Soliton_DensProf_NBin, Table_Radius, Table_Density, r_ref );

else if ( r_ref >= Table_Radius[Soliton_DensProf_NBin-1] )
dens_ref = Table_Density[Soliton_DensProf_NBin-1];
if ( dens_ref == NULL_REAL )
{
if ( r_ref < Table_Radius[0] )
dens_ref = Table_Density[0];

else
Aux_Error( ERROR_INFO, "interpolation failed at radius %13.7e (min/max radius = %13.7e/%13.7e) !!\n",
r_ref, Table_Radius[0], Table_Radius[Soliton_DensProf_NBin-1] );
}
else if ( r_ref >= Table_Radius[Soliton_DensProf_NBin-1] )
dens_ref = Table_Density[Soliton_DensProf_NBin-1];

// rescale density (reference density --> target density) and add to the fluid array
fluid[DENS] += dens_ref*Soliton_ScaleD[t];
} // for (int t=0; t<Soliton_N; t++)
else
Aux_Error( ERROR_INFO, "interpolation failed at radius %13.7e (min/max radius = %13.7e/%13.7e) !!\n",
r_ref, Table_Radius[0], Table_Radius[Soliton_DensProf_NBin-1] );
}

// rescale density (reference density --> target density) and add to the fluid array
fluid[DENS] += dens_ref*Soliton_ScaleD[t];
} // for (int t=0; t<Soliton_N; t++)
} // if ( Soliton_InitMode == 1 )
else if ( Soliton_InitMode == 2 )
{
const double m22 = ELBDM_MASS*UNIT_M/(Const_eV/SQR(Const_c))/1.0e-22;

// loop over all solitons to get the total density
for (int t=0; t<Soliton_N; t++)
{
r_tar = sqrt( SQR(x-Soliton_Center[t][0]) + SQR(y-Soliton_Center[t][1]) + SQR(z-Soliton_Center[t][2]) );

const double rc_kpc = Soliton_CoreRadius[t]*UNIT_L/Const_kpc;
dens = 1.945/SQR(m22*10)/POW4(rc_kpc) / POW(1 + (POW(2.0,1.0/8.0)-1.0)*SQR(r_tar/Soliton_CoreRadius[t]), 8); // in unit of Msun/pc^3
dens = dens* Const_Msun/CUBE(Const_pc) / UNIT_D ; // code unit
fluid[DENS] += dens;
} // for (int t=0; t<Soliton_N; t++)
} // if ( Soliton_InitMode == 2 )
else
Aux_Error( ERROR_INFO, "unsupported Soliton_InitMode (%d) !!\n", Soliton_InitMode );


# if ( ELBDM_SCHEME == ELBDM_HYBRID )
Expand Down Expand Up @@ -411,7 +463,7 @@ void End_Soliton()

//-------------------------------------------------------------------------------------------------------
// Function : BC
// Description : Set the extenral boundary condition
// Description : Set the extenal boundary condition
//
// Note : 1. Linked to the function pointer "BC_User_Ptr"
//
Expand Down