Support analytical profile in Soliton test

[vivi235711]

Summary

Add Soliton_InitMode for the ELBDM Soliton test:

• 1 (default): Table-based profile (existing behavior)
• 2: Analytical soliton profile (no external table)

Changes

• Added Soliton_InitMode parameter and handling in initialization
• Implemented analytical density computation for mode 2

[boxianliu]

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

[vivi235711]

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

[hyschive]

@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.

[boxianliu]

I have made the following corrections, which you can refer to:

1. Use CoreRadiusRef = RadiusRef[b] + (DensCore - DensRef[b])*(RadiusRef[b+1]-RadiusRef[b])/(DensRef[b+1]-DensRef[b]); instead of CoreRadiusRef = 0.5*( RadiusRef[b] + RadiusRef[b+1] ); to improve interpolation accuracy.
2. Make 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 ); only print when Soliton_InitMode==1.
3. Aux_Message( stdout, " %7s %13s %13s %13s\n", "ID", "CoreRadius", "Center_X", "Center_Y", "Center_Z" ); is replaced with Aux_Message( stdout, " %7s %13s %13s %13s %13s\n", "ID", "CoreRadius", "Center_X", "Center_Y", "Center_Z" );
4. Adjust the overall arrangement of SetGridIC() so that const double *Table_Radius = Soliton_DensProf + 0*Soliton_DensProf_NBin; const double *Table_Density = Soliton_DensProf + 1*Soliton_DensProf_NBin; only executes when Soliton_InitMode==1 (when Soliton_InitMode==2, Soliton_DensProf == NULL).
5. Change const double m22 = ELBDM_MASS*UNIT_M/(Const_eV/SQR(Const_c))/1.0e-22; This should be moved outside the for loop.

In addition, I have the following suggestions:

1. In Description : Set the extenral boundary condition, external was mistakenly written as extenral.
2. In SolitonDensityProfile_Lambda0.0, Dens is actually |ψ|^2. In Init_TestProb_ELBDM_Soliton.cpp, it will be converted back to Dens using SCALE_D. This part might need to be noted, otherwise the user-inputted table may be incorrect.
3. In Soliton_InitMode==2 mode, code units and particle mass must be correctly set. Init_TestProb_ELBDM_Soliton.cpp will first convert the units to kpc and Msun/pc^3, so OPT__UNIT must not be set as 0. This can also be noted to remind users.

[boxianliu]

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

[vivi235711]

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.

[vivi235711]

@boxianliu Thanks for your detailed review! I've updated according to your comments.

In SolitonDensityProfile_Lambda0.0, Dens is actually |ψ|^2. In Init_TestProb_ELBDM_Soliton.cpp, it will be converted back to Dens using SCALE_D. This part might need to be noted, otherwise the user-inputted table may be incorrect.

Added a note in the README.

In Soliton_InitMode==2 mode, code units and particle mass must be correctly set. Init_TestProb_ELBDM_Soliton.cpp will first convert the units to kpc and Msun/pc^3, so OPT__UNIT must not be set as 0. This can also be noted to remind users.

I've added a check if the !OPT__UNIT and Soliton_InitMode == 2 the program would stop and output error.

Please take a look when you are able.

[boxianliu]

Thank you for your revisions. I have checked them and I think there are no other issues!