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
+118 −64
Open

Support analytical profile in Soliton test #484

Show file tree
Hide file tree
Changes from 1 commit
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
154 changes: 96 additions & 58 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 @@ -176,8 +178,14 @@ void SetParameter()
// (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 +239,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 = 0.5*( RadiusRef[b] + RadiusRef[b+1] );
Copy link

@boxianliu boxianliu Nov 13, 2025

Choose a reason for hiding this comment

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

           CoreRadiusRef = RadiusRef[b] + (DensCore - DensRef[b])*(RadiusRef[b+1]-RadiusRef[b])/(DensRef[b+1]-DensRef[b]);

Copy link
Contributor Author

@vivi235711 vivi235711 Dec 1, 2025

Choose a reason for hiding this comment

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

I agree with you that it will improve interpolation accuracy.
But the result will be a little different compare to previous version ( about 0.1% )
What do you think? @hyschive

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 @@ -297,12 +311,27 @@ void SetParameter()
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\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,42 +362,51 @@ 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
fluid[DENS] = 0.0;

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

// loop over all solitons to get the total density
// 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]) );
if ( Soliton_InitMode == 1 )
{
// rescale radius (target radius --> reference radius)
r_ref = r_tar / Soliton_ScaleL[t];

// rescale radius (target radius --> reference radius)
r_ref = r_tar / Soliton_ScaleL[t];

// linear interpolation
dens_ref = Mis_InterpolateFromTable( Soliton_DensProf_NBin, Table_Radius, Table_Density, r_ref );
// linear interpolation
dens_ref = Mis_InterpolateFromTable( Soliton_DensProf_NBin, Table_Radius, Table_Density, r_ref );

if ( dens_ref == NULL_REAL )
{
if ( r_ref < Table_Radius[0] )
dens_ref = Table_Density[0];
if ( dens_ref == NULL_REAL )
{
if ( r_ref < Table_Radius[0] )
dens_ref = Table_Density[0];

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

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
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];
// rescale density (reference density --> target density) and add to the fluid array
fluid[DENS] += dens_ref*Soliton_ScaleD[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;
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;
} // if ( Soliton_InitMode == 2 )
else
Aux_Error( ERROR_INFO, "unsupported Soliton_InitMode (%d) !!\n", Soliton_InitMode );
} // for (int t=0; t<Soliton_N; t++)
Copy link

@boxianliu boxianliu Nov 13, 2025

Choose a reason for hiding this comment

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

double r_tar, r_ref, dens_ref, dens;

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

if ( Soliton_InitMode == 1 )
{
    if ( Soliton_DensProf == NULL || Soliton_DensProf_NBin < 2 )
        Aux_Error( ERROR_INFO, "invalid Soliton_DensProf (ptr=%p, NBin=%d)\n",
                  (void*)Soliton_DensProf, Soliton_DensProf_NBin );
   
    const double *Table_Radius  = Soliton_DensProf + 0*Soliton_DensProf_NBin;  // radius
    const double *Table_Density = Soliton_DensProf + 1*Soliton_DensProf_NBin;  // density

    // 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]) );

        // rescale radius (target radius --> reference radius)
        r_ref = r_tar / Soliton_ScaleL[t];

        // linear interpolation
        dens_ref = Mis_InterpolateFromTable( Soliton_DensProf_NBin, Table_Radius, Table_Density, r_ref );

        if ( dens_ref == NULL_REAL )
        {
            if      ( r_ref <  Table_Radius[0] )
                dens_ref = Table_Density[0];

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

            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 );

Copy link
Contributor Author

@vivi235711 vivi235711 Dec 1, 2025

Choose a reason for hiding this comment

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

Your version is way better than mine. Thanks!
By the way, I removed the following error check from the code:

    if ( Soliton_DensProf == NULL || Soliton_DensProf_NBin < 2 )
        Aux_Error( ERROR_INFO, "invalid Soliton_DensProf (ptr=%p, NBin=%d)\n",
                  (void*)Soliton_DensProf, Soliton_DensProf_NBin );

This check seemed redundant since the error should already be caught within the Aux_LoadTable function.



Expand Down