Skip to content

Commit

Permalink
Merge pull request #371 from sandy0216/soliton
Browse files Browse the repository at this point in the history 
Add soliton potential for Bondi accretion problem
  • Loading branch information
@hyschive
hyschive authored Jan 26, 2025
2 parents 5d779cd + 8a5149c commit b9d3fd3
Show file tree
Hide file tree
Showing 6 changed files with 344 additions and 143 deletions.
12 changes: 12 additions & 0 deletions example/test_problem/Hydro/Bondi/Input__TestProb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ Bondi_InBC_Rho 1.0e-32 # density inside the void region (in
Bondi_InBC_T 1.0e-7 # temperature inside the void region (in keV)
Bondi_InBC_NCell 2.0 # number of finest cells (can be a fractional number) for the inner BC
Bondi_Soften_NCell 0.0 # number of finest cells (can be a fractional number) for the soften length (<=0.0 -> disable)
Bondi_void 1 # enable the void region [1]
Bondi_dynBH 0 # dynamically increase BH mass [0]

Bondi_HSE 0 # enable HSE [0]
Bondi_HSE_Mode 1 # initial configuration (1:T=Bondi_T0, 2:rho~1/r, 3:beta model) [1]
Expand All @@ -22,3 +24,13 @@ Bondi_HSE_TrunD 1.6735328e-24 # see Bondi_HSE_Truncate (in g/cm^3)
Bondi_HSE_TrunSmoothR 5.0e-2 # smooth out density within TrunR-SmoothR<r<TrunR+SmoothR (in kpc; <0.0->off) [-1.0]
Bondi_HSE_Pres_NormT 0 # normalize pressure profile such that T(r=Dens_NormR)=Bondi_T0 [0]
Bondi_HSE_Beta_Rcore 1.0e-1 # core radius in the beta model (in kpc)

Bondi_Soliton 0 # add soliton external potential [0]
Bondi_Soliton_m22 1.0 # FDM particle mass in 1e-22 eV/c^2 for Bondi_Soliton [-1.0]
Bondi_Soliton_type 5 # functional form for gradually introducing the soliton potential
# (0:unity, 1:arctan, 2:linear, 3:smooth step function, 4:sigmoid, 5:tanh) [5]
Bondi_Soliton_t 4.0e3 # characteristic time normalized to Bondi_TimeB for adding the soliton potential [-1.0]
Bondi_Soliton_rc -1.0 # soliton radius for Bondi_Soliton (in kpc)
# (<0.0 --> compute from Bondi_Soliton_MassHalo/Redshift using the core-halo relation) [-1.0]
Bondi_Soliton_MassHalo 1.0e11 # halo mass for determining Bondi_Soliton_rc (in Msun) [-1.0]
Bondi_Soliton_Redshift 7.0 # redshift for determining Bondi_Soliton_rc [-1.0]
4 changes: 2 additions & 2 deletions example/test_problem/Hydro/Bondi/plot_diagonal.gpt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ do for [ID=ID_START:ID_END:ID_DELTA] {
set yrange [1.0e-3:2.0e1]

plot sprintf( '%s/Diag_%06d', FILE_SIMU, ID ) \
u (abs($4-CENTER)*3**0.5):16 w p pt 6 lc 6 tit 'Simulation' \
u (abs($4-CENTER)*3**0.5):17 w p pt 6 lc 6 tit 'Simulation' \
,FILE_BONDI u 2:7 w l lc -1 tit 'Analytical'


Expand All @@ -79,7 +79,7 @@ do for [ID=ID_START:ID_END:ID_DELTA] {
set yrange [1.0e-2:1.0e3]

plot sprintf( '%s/Diag_%06d', FILE_SIMU, ID ) \
u (abs($4-CENTER)*3**0.5):13 w p pt 6 lc 6 tit 'Simulation' \
u (abs($4-CENTER)*3**0.5):14 w p pt 6 lc 6 tit 'Simulation' \
,FILE_BONDI u 2:( ($6/$7)**2.0*$4/GAMMA ) w l lc -1 tit 'Analytical'


Expand Down
84 changes: 76 additions & 8 deletions src/TestProblem/Hydro/Bondi/ExtAcc_Bondi.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,11 @@

extern double Bondi_MassBH;
extern double Bondi_Soften_R;
extern bool Bondi_Soliton;
extern double Bondi_Soliton_m22;
extern double Bondi_Soliton_rc;
extern int Bondi_Soliton_type;
extern double Bondi_Soliton_t;

//-------------------------------------------------------------------------------------------------------
// Function : SetExtAccAuxArray_Bondi
Expand All @@ -41,9 +46,55 @@ void SetExtAccAuxArray_Bondi( double AuxArray[], const double Time )
AuxArray[0] = amr->BoxCenter[0];
AuxArray[1] = amr->BoxCenter[1];
AuxArray[2] = amr->BoxCenter[2];
AuxArray[3] = NEWTON_G*Bondi_MassBH; // gravitational_constant*point_source_mass
AuxArray[3] = NEWTON_G*Bondi_MassBH; // gravitational_constant*black_hole_mass (in code units)
AuxArray[4] = Bondi_Soften_R; // soften_length (<=0.0 --> disable)

double Coeff_t;
switch ( Bondi_Soliton_type )
{
// unity
case 0: Coeff_t = 1.0;
break;

// arctan function
case 1: Coeff_t = 2.0/M_PI*atan( Time/Bondi_Soliton_t );
break;

// linear function
case 2: Coeff_t = ( Time < Bondi_Soliton_t ) ? Time/Bondi_Soliton_t
: 1.0;
break;

// smooth step function
case 3: Coeff_t = ( Time < Bondi_Soliton_t ) ? 3.0*SQR( Time/Bondi_Soliton_t ) - 2.0*CUBE( Time/Bondi_Soliton_t )
: 1.0;
break;

// sigmoid
case 4: Coeff_t = 2.0 / ( 1.0 + exp( -Time*log(3.0)/Bondi_Soliton_t ) ) - 1.0;
break;

// tanh
case 5: Coeff_t = tanh( Time/Bondi_Soliton_t );
break;

default:
Aux_Error( ERROR_INFO, "unsupported Bondi_Soliton_type (%d) !!\n", Bondi_Soliton_type );
} // switch ( Bondi_Soliton_type )

if ( Bondi_Soliton )
{
AuxArray[5] = Coeff_t*NEWTON_G*( 4.17e9*Const_Msun/UNIT_M )/
( SQR( Bondi_Soliton_m22*(real)10.0 )*( Bondi_Soliton_rc*UNIT_L/Const_pc ) );
AuxArray[6] = Bondi_Soliton_rc;
}

else
{
AuxArray[5] = -1.0;
AuxArray[6] = -1.0;
}

} // FUNCTION : SetExtAccAuxArray_Bondi
#endif // #ifndef __CUDACC__

Expand Down Expand Up @@ -72,13 +123,30 @@ static void ExtAcc_Bondi( real Acc[], const double x, const double y, const doub
const double UserArray[] )
{

const double Cen[3] = { UserArray[0], UserArray[1], UserArray[2] };
const real GM = (real)UserArray[3];
const real eps = (real)UserArray[4];
const real dx = (real)(x - Cen[0]);
const real dy = (real)(y - Cen[1]);
const real dz = (real)(z - Cen[2]);
const real r = SQRT( dx*dx + dy*dy + dz*dz );
const double Cen[3] = { UserArray[0], UserArray[1], UserArray[2] };
real GM = (real)UserArray[3];
const real eps = (real)UserArray[4];
const real GM0_sol = (real)UserArray[5];
const real rc = (real)UserArray[6];
const real dx = (real)(x - Cen[0]);
const real dy = (real)(y - Cen[1]);
const real dz = (real)(z - Cen[2]);
const real r = SQRT( dx*dx + dy*dy + dz*dz );

if ( GM0_sol > (real)0.0 && rc > (real)0.0 )
{
const real a = SQRT( POW( (real)2.0, (real)1.0/(real)8.0 ) - (real)1.0 )*(r/rc);
const real GM_sol = (real)GM0_sol/( POW( SQR(a)+(real)1.0, (real)7.0 ) )*
( (real) 3465*POW( a, (real)13.0 )
+(real) 23100*POW( a, (real)11.0 )
+(real) 65373*POW( a, (real) 9.0 )
+(real)101376*POW( a, (real) 7.0 )
+(real) 92323*POW( a, (real) 5.0 )
+(real) 48580*POW( a, (real) 3.0 )
-(real) 3465*a
+(real) 3465*POW( SQR(a)+(real)1.0, (real)7.0 )*ATAN(a) );
GM += GM_sol;
} // if ( GM0_sol > 0.0 && rc > 0.0 )

// Plummer
# if ( defined SOFTEN_PLUMMER )
Expand Down
29 changes: 25 additions & 4 deletions src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,7 @@

#if ( MODEL == HYDRO && defined GRAVITY )



extern double Bondi_MassBH;
extern double Bondi_InBC_Rho;
extern double Bondi_InBC_R;
extern double Bondi_InBC_E;
Expand All @@ -19,6 +18,9 @@ extern double Bondi_SinkEk;
extern double Bondi_SinkEt;
extern int Bondi_SinkNCell;

extern bool Bondi_void;
extern bool Bondi_dynBH;




Expand Down Expand Up @@ -52,6 +54,9 @@ int Flu_ResetByUser_Func_Bondi( real fluid[], const double Emag, const double x,
const double dt, const int lv, double AuxArray[] )
{

if ( !Bondi_void ) return false;


const double Pos[3] = { x, y, z };
const double InBC_R2 = SQR( Bondi_InBC_R );
double dr2[3], r2;
Expand Down Expand Up @@ -113,14 +118,16 @@ void Flu_ResetByUser_API_Bondi( const int lv, const int FluSg, const int MagSg,
int Reset;
real fluid[NCOMP_TOTAL], fluid_bk[NCOMP_TOTAL];
double x, y, z, x0, y0, z0;
double SinkMass_OneSubStep_ThisRank = 0.0; // variables to record sink mass at every time step
double SinkMass_OneSubStep_AllRank;

// reset to 0 since we only want to record the number of void cells **for one sub-step**
Bondi_SinkNCell = 0;


# pragma omp parallel for private( Reset, fluid, fluid_bk, x, y, z, x0, y0, z0 ) schedule( runtime ) \
reduction(+:Bondi_SinkMass, Bondi_SinkMomX, Bondi_SinkMomY, Bondi_SinkMomZ, Bondi_SinkMomXAbs, Bondi_SinkMomYAbs, Bondi_SinkMomZAbs, \
Bondi_SinkEk, Bondi_SinkEt, Bondi_SinkNCell)
Bondi_SinkEk, Bondi_SinkEt, Bondi_SinkNCell, SinkMass_OneSubStep_ThisRank )
for (int PID=0; PID<amr->NPatchComma[lv][1]; PID++)
{
x0 = amr->patch[0][lv][PID]->EdgeL[0] + 0.5*dh;
Expand Down Expand Up @@ -193,12 +200,26 @@ void Flu_ResetByUser_API_Bondi( const int lv, const int FluSg, const int MagSg,
Bondi_SinkEk += dv*Ek;
Bondi_SinkEt += dv*Et;
Bondi_SinkNCell ++;

SinkMass_OneSubStep_ThisRank += dv*fluid_bk[DENS];
}
} // if ( Reset )

else if ( amr->patch[0][lv][PID]->son == -1 )
{
// void region must be completely refined to the max level
Aux_Error( ERROR_INFO, "void region lies outside the max-level region !!\n" );
}
} // if ( Reset )
}}} // i,j,k
} // for (int PID=0; PID<amr->NPatchComma[lv][1]; PID++)

if ( Bondi_dynBH )
{
MPI_Allreduce( &SinkMass_OneSubStep_ThisRank, &SinkMass_OneSubStep_AllRank, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );

Bondi_MassBH += SinkMass_OneSubStep_AllRank;
}

} // FUNCTION : Flu_ResetByUser_API_Bondi


Expand Down
Loading

0 comments on commit b9d3fd3

Please sign in to comment.