diff --git a/example/test_problem/Hydro/Bondi/Input__TestProb b/example/test_problem/Hydro/Bondi/Input__TestProb index 60fce18c2..df4e5c205 100644 --- a/example/test_problem/Hydro/Bondi/Input__TestProb +++ b/example/test_problem/Hydro/Bondi/Input__TestProb @@ -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] @@ -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-SmoothRoff) [-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] diff --git a/example/test_problem/Hydro/Bondi/plot_diagonal.gpt b/example/test_problem/Hydro/Bondi/plot_diagonal.gpt index 54e94c852..a4872a825 100644 --- a/example/test_problem/Hydro/Bondi/plot_diagonal.gpt +++ b/example/test_problem/Hydro/Bondi/plot_diagonal.gpt @@ -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' @@ -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' diff --git a/src/TestProblem/Hydro/Bondi/ExtAcc_Bondi.cpp b/src/TestProblem/Hydro/Bondi/ExtAcc_Bondi.cpp index 44af678bb..913058f67 100644 --- a/src/TestProblem/Hydro/Bondi/ExtAcc_Bondi.cpp +++ b/src/TestProblem/Hydro/Bondi/ExtAcc_Bondi.cpp @@ -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 @@ -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__ @@ -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 ) diff --git a/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp b/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp index 9f8fffb30..c8a6a2bb3 100644 --- a/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp +++ b/src/TestProblem/Hydro/Bondi/Flu_ResetByUser_Bondi.cpp @@ -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; @@ -19,6 +18,9 @@ extern double Bondi_SinkEk; extern double Bondi_SinkEt; extern int Bondi_SinkNCell; +extern bool Bondi_void; +extern bool Bondi_dynBH; + @@ -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; @@ -113,6 +118,8 @@ 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; @@ -120,7 +127,7 @@ void Flu_ResetByUser_API_Bondi( const int lv, const int FluSg, const int MagSg, # 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; PIDNPatchComma[lv][1]; PID++) { x0 = amr->patch[0][lv][PID]->EdgeL[0] + 0.5*dh; @@ -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; PIDNPatchComma[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 diff --git a/src/TestProblem/Hydro/Bondi/Init_TestProb_Hydro_Bondi.cpp b/src/TestProblem/Hydro/Bondi/Init_TestProb_Hydro_Bondi.cpp index f57f4af15..ab63b7f30 100644 --- a/src/TestProblem/Hydro/Bondi/Init_TestProb_Hydro_Bondi.cpp +++ b/src/TestProblem/Hydro/Bondi/Init_TestProb_Hydro_Bondi.cpp @@ -4,76 +4,89 @@ // problem-specific global variables // ======================================================================================= - double Bondi_MassBH; // black hole mass -static double Bondi_Rho0; // background density -static double Bondi_T0; // background temperature - double Bondi_RefineRadius0; // refinement radius at the base level - // NOTE: refinement radius at Lv is set to Bondi_RefineRadius0*2^(-Lv) - // --> all refinement shells have roughly the same number of cells at its level - // (except Lv=MAX_LEVEL, which will have twice the number of cells along the radius - // unless Bondi_HalfMaxLvRefR is on) - bool Bondi_HalfMaxLvRefR; // halve the refinement radius at the maximum level - double Bondi_InBC_Rho; // density inside the void region -static double Bondi_InBC_T; // temperature inside the void region -static double Bondi_InBC_NCell; // number of finest cells for the radius of the void region -static double Bondi_Soften_NCell; // number of finest cells for the soften length (<=0.0 ==> disable) - - double Bondi_InBC_R; // radius of the void region (=Bondi_InBC_NCell*dh[MAX_LEVEL]) - double Bondi_InBC_P; // pressure inside the void region - double Bondi_InBC_E; // energy inside the void region - double Bondi_Soften_R; // soften length (=Bondi_Soften_NCell*dh[MAX_LEVEL]) -static double Bondi_P0; // background pressure -static double Bondi_Cs; // background sound speed -static double Bondi_RS; // Schwarzschild radius -static double Bondi_RB; // Bondi radius -static double Bondi_TimeB; // Bondi time - - double Bondi_SinkMass; // total mass in the void region removed in one global time-step - double Bondi_SinkMomX; // total x-momentum ... - double Bondi_SinkMomY; // total y-momentum ... - double Bondi_SinkMomZ; // total z-momentum ... - double Bondi_SinkMomXAbs; // total |x-momentum| ... - double Bondi_SinkMomYAbs; // total |y-momentum| ... - double Bondi_SinkMomZAbs; // total |z-momentum| ... - double Bondi_SinkEk; // total kinematic energy ... - double Bondi_SinkEt; // total thermal energy ... - int Bondi_SinkNCell; // total number of finest cells within the void region + double Bondi_MassBH; // black hole mass +static double Bondi_Rho0; // background density +static double Bondi_T0; // background temperature + double Bondi_RefineRadius0; // refinement radius at the base level + // NOTE: refinement radius at Lv is set to Bondi_RefineRadius0*2^(-Lv) + // --> all refinement shells have roughly the same number of cells at its level + // (except Lv=MAX_LEVEL, which will have twice the number of cells along the radius + // unless Bondi_HalfMaxLvRefR is on) + bool Bondi_HalfMaxLvRefR; // halve the refinement radius at the maximum level + double Bondi_InBC_Rho; // density inside the void region +static double Bondi_InBC_T; // temperature inside the void region +static double Bondi_InBC_NCell; // number of finest cells for the radius of the void region +static double Bondi_Soften_NCell; // number of finest cells for the soften length (<=0.0 ==> disable) + bool Bondi_void; // enable the void region + bool Bondi_dynBH; // dynamically increase BH mass + + double Bondi_InBC_R; // radius of the void region (=Bondi_InBC_NCell*dh[MAX_LEVEL]) + double Bondi_InBC_P; // pressure inside the void region + double Bondi_InBC_E; // energy inside the void region + double Bondi_Soften_R; // soften length (=Bondi_Soften_NCell*dh[MAX_LEVEL]) +static double Bondi_P0; // background pressure +static double Bondi_Cs; // background sound speed +static double Bondi_RS; // Schwarzschild radius +static double Bondi_RB; // Bondi radius +static double Bondi_TimeB; // Bondi time + + double Bondi_SinkMass; // total mass in the void region removed in one global time-step + double Bondi_SinkMomX; // total x-momentum ... + double Bondi_SinkMomY; // total y-momentum ... + double Bondi_SinkMomZ; // total z-momentum ... + double Bondi_SinkMomXAbs; // total |x-momentum| ... + double Bondi_SinkMomYAbs; // total |y-momentum| ... + double Bondi_SinkMomZAbs; // total |z-momentum| ... + double Bondi_SinkEk; // total kinematic energy ... + double Bondi_SinkEt; // total thermal energy ... + int Bondi_SinkNCell; // total number of finest cells within the void region // external units in cgs -const double UnitExt_L = Const_kpc; -const double UnitExt_D = 1.0; -const double UnitExt_M = Const_Msun; -const double UnitExt_E = Const_keV; +const double UnitExt_L = Const_kpc; +const double UnitExt_D = 1.0; +const double UnitExt_M = Const_Msun; +const double UnitExt_E = Const_keV; // hydrostatic equilibrium (HSE) -static bool Bondi_HSE; // enable HSE -static int Bondi_HSE_Mode; // initial configuration (1:T=Bondi_T0, 2:rho~1/r, 3:beta model) -static double Bondi_HSE_Dens_NormR; // normalize the density profile to density(r=NormR)=NormD -static double Bondi_HSE_Dens_NormD; // see Bondi_HSE_Dens_NormR +static bool Bondi_HSE; // enable HSE +static int Bondi_HSE_Mode; // initial configuration (1:T=Bondi_T0, 2:rho~1/r, 3:beta model) +static double Bondi_HSE_Dens_NormR; // normalize the density profile to density(r=NormR)=NormD +static double Bondi_HSE_Dens_NormD; // see Bondi_HSE_Dens_NormR // parameters for Bondi_HSE_Mode=1 -static int Bondi_HSE_Dens_NBin; // number of bins in the density profile table -static double Bondi_HSE_Dens_MinR; // minimum radius in the density profile -static double Bondi_HSE_Dens_MaxR; // maximum ... -static bool Bondi_HSE_Truncate; // truncate density within r adjust P2 in the pressure profile such that T(r=NormR)=Bondi_T0 -static double Bondi_HSE_Pres_NormP1; // P=P1*r^-2+P2 (P2=0 when Bondi_HSE_Pres_NormT=false) -static double Bondi_HSE_Pres_NormP2; +static bool Bondi_HSE_Pres_NormT; // true --> adjust P2 in the pressure profile such that T(r=NormR)=Bondi_T0 +static double Bondi_HSE_Pres_NormP1; // P=P1*r^-2+P2 (P2=0 when Bondi_HSE_Pres_NormT=false) +static double Bondi_HSE_Pres_NormP2; // parameters for Bondi_HSE_Mode=3 (beta model) -const double Bondi_HSE_Beta = 2.0/3.0; // beta (must be 2/3 for now) -static double Bondi_HSE_Beta_Rcore; // core radius (input parameter) -static double Bondi_HSE_Beta_Rho0; // peak density (set by Bondi_HSE_Dens_NormR/D) -static double Bondi_HSE_Beta_P1; // P(r) = P1*( 1/x + atan(x) ) + P2 assuming beta=2/3, where x=r/Rcore, -static double Bondi_HSE_Beta_P2; // P1=G*MassBH*Rho0/Rcore, and P2 currently fixed to -0.5*pi*P1 so that P(inf)=0 +const double Bondi_HSE_Beta = 2.0/3.0; // beta (must be 2/3 for now) +static double Bondi_HSE_Beta_Rcore; // core radius (input parameter) +static double Bondi_HSE_Beta_Rho0; // peak density (set by Bondi_HSE_Dens_NormR/D) +static double Bondi_HSE_Beta_P1; // P(r) = P1*( 1/x + atan(x) ) + P2 assuming beta=2/3, where x=r/Rcore, +static double Bondi_HSE_Beta_P2; // P1=G*MassBH*Rho0/Rcore, and P2 currently fixed to -0.5*pi*P1 so that P(inf)=0 + +// parameters for soliton + bool Bondi_Soliton; // add soliton external potential + double Bondi_Soliton_m22; // FDM particle mass in 1e-22 eV/c^2 for Bondi_Soliton + int Bondi_Soliton_type; // functional form for gradually introducing the soliton potential + // (0:unity, 1:arctan, 2:linear, 3:smooth step function, 4:sigmoid, 5:tanh) + double Bondi_Soliton_t; // characteristic time normalized to Bondi_TimeB for adding the soliton potential + double Bondi_Soliton_rc; // soliton radius for Bondi_Soliton + // (<0.0 --> compute from Bondi_Soliton_MassHalo/Redshift using the core-halo relation) +static double Bondi_Soliton_MassHalo; // halo mass for determining Bondi_Soliton_rc +static double Bondi_Soliton_Redshift; // redshift for determining Bondi_Soliton_rc // ======================================================================================= @@ -89,6 +102,35 @@ static void BondiBC( real Array[], const int ArraySize[], real fluid[], const in const int GhostSize, const int idx[], const double pos[], const double Time, const int lv, const int TFluVarIdxList[], double AuxArray[] ); +void SetExtAccAuxArray_Bondi( double [], const double ); + + + + +#ifdef GRAVITY +//------------------------------------------------------------------------------------------------------- +// Function : Poi_UserWorkBeforePoisson_Bondi +// Description : Call SetExtAccAuxArray_Bondi() to reset Bondi_MassBH and soliton parameters before +// invoking the Poisson solver +// +// Note : 1. Invoked by Gra_AdvanceDt() using the function pointer "Poi_UserWorkBeforePoisson_Ptr" +// +// Parameter : Time : Target physical time +// lv : Target refinement level +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void Poi_UserWorkBeforePoisson_Bondi( const double Time, const int lv ) +{ + + SetExtAccAuxArray_Bondi( ExtAcc_AuxArray, Time ); + +# ifdef GPU + CUAPI_SetConstMemory_ExtAccPot(); +# endif + +} // FUNCTION : Poi_UserWorkBeforePoisson_Bondi +#endif // #ifdef GRAVITY @@ -186,31 +228,41 @@ void SetParameter() // --> note that VARIABLE, DEFAULT, MIN, and MAX must have the same data type // --> some handy constants (e.g., NoMin_int, Eps_float, ...) are defined in "include/ReadPara.h" // ******************************************************************************************************************************** -// ReadPara->Add( "KEY_IN_THE_FILE", &VARIABLE, DEFAULT, MIN, MAX ); +// ReadPara->Add( "KEY_IN_THE_FILE", &VARIABLE, DEFAULT, MIN, MAX ); // ******************************************************************************************************************************** - ReadPara->Add( "Bondi_MassBH", &Bondi_MassBH, -1.0, Eps_double, NoMax_double ); - ReadPara->Add( "Bondi_Rho0", &Bondi_Rho0, -1.0, Eps_double, NoMax_double ); - ReadPara->Add( "Bondi_T0", &Bondi_T0, -1.0, Eps_double, NoMax_double ); - ReadPara->Add( "Bondi_RefineRadius0", &Bondi_RefineRadius0, -1.0, Eps_double, NoMax_double ); - ReadPara->Add( "Bondi_HalfMaxLvRefR", &Bondi_HalfMaxLvRefR, true, Useless_bool, Useless_bool ); - ReadPara->Add( "Bondi_InBC_Rho", &Bondi_InBC_Rho, -1.0, Eps_double, NoMax_double ); - ReadPara->Add( "Bondi_InBC_T", &Bondi_InBC_T, -1.0, Eps_double, NoMax_double ); - ReadPara->Add( "Bondi_InBC_NCell", &Bondi_InBC_NCell, -1.0, Eps_double, NoMax_double ); - ReadPara->Add( "Bondi_Soften_NCell", &Bondi_Soften_NCell, -1.0, NoMin_double, NoMax_double ); - - ReadPara->Add( "Bondi_HSE", &Bondi_HSE, false, Useless_bool, Useless_bool ); - ReadPara->Add( "Bondi_HSE_Mode", &Bondi_HSE_Mode, 1, 1, 3 ); - ReadPara->Add( "Bondi_HSE_Dens_NBin", &Bondi_HSE_Dens_NBin, 10000, 2, NoMax_int ); - ReadPara->Add( "Bondi_HSE_Dens_MinR", &Bondi_HSE_Dens_MinR, -1.0, NoMin_double, NoMax_double ); - ReadPara->Add( "Bondi_HSE_Dens_MaxR", &Bondi_HSE_Dens_MaxR, -1.0, NoMin_double, NoMax_double ); - ReadPara->Add( "Bondi_HSE_Dens_NormR", &Bondi_HSE_Dens_NormR, -1.0, NoMin_double, NoMax_double ); - ReadPara->Add( "Bondi_HSE_Dens_NormD", &Bondi_HSE_Dens_NormD, -1.0, Eps_double, NoMax_double ); - ReadPara->Add( "Bondi_HSE_Truncate", &Bondi_HSE_Truncate, true, Useless_bool, Useless_bool ); - ReadPara->Add( "Bondi_HSE_TrunR", &Bondi_HSE_TrunR, -1.0, NoMin_double, NoMax_double ); - ReadPara->Add( "Bondi_HSE_TrunD", &Bondi_HSE_TrunD, -1.0, Eps_double, NoMax_double ); - ReadPara->Add( "Bondi_HSE_TrunSmoothR",&Bondi_HSE_TrunSmoothR, -1.0, NoMin_double, NoMax_double ); - ReadPara->Add( "Bondi_HSE_Pres_NormT", &Bondi_HSE_Pres_NormT, false, Useless_bool, Useless_bool ); - ReadPara->Add( "Bondi_HSE_Beta_Rcore", &Bondi_HSE_Beta_Rcore, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_MassBH", &Bondi_MassBH, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_Rho0", &Bondi_Rho0, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_T0", &Bondi_T0, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_RefineRadius0", &Bondi_RefineRadius0, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_HalfMaxLvRefR", &Bondi_HalfMaxLvRefR, true, Useless_bool, Useless_bool ); + ReadPara->Add( "Bondi_InBC_Rho", &Bondi_InBC_Rho, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_InBC_T", &Bondi_InBC_T, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_InBC_NCell", &Bondi_InBC_NCell, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_Soften_NCell", &Bondi_Soften_NCell, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Bondi_void", &Bondi_void, true, Useless_bool, Useless_bool ); + ReadPara->Add( "Bondi_dynBH", &Bondi_dynBH, false, Useless_bool, Useless_bool ); + + ReadPara->Add( "Bondi_HSE", &Bondi_HSE, false, Useless_bool, Useless_bool ); + ReadPara->Add( "Bondi_HSE_Mode", &Bondi_HSE_Mode, 1, 1, 3 ); + ReadPara->Add( "Bondi_HSE_Dens_NBin", &Bondi_HSE_Dens_NBin, 10000, 2, NoMax_int ); + ReadPara->Add( "Bondi_HSE_Dens_MinR", &Bondi_HSE_Dens_MinR, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Bondi_HSE_Dens_MaxR", &Bondi_HSE_Dens_MaxR, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Bondi_HSE_Dens_NormR", &Bondi_HSE_Dens_NormR, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Bondi_HSE_Dens_NormD", &Bondi_HSE_Dens_NormD, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_HSE_Truncate", &Bondi_HSE_Truncate, true, Useless_bool, Useless_bool ); + ReadPara->Add( "Bondi_HSE_TrunR", &Bondi_HSE_TrunR, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Bondi_HSE_TrunD", &Bondi_HSE_TrunD, -1.0, Eps_double, NoMax_double ); + ReadPara->Add( "Bondi_HSE_TrunSmoothR", &Bondi_HSE_TrunSmoothR, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Bondi_HSE_Pres_NormT", &Bondi_HSE_Pres_NormT, false, Useless_bool, Useless_bool ); + ReadPara->Add( "Bondi_HSE_Beta_Rcore", &Bondi_HSE_Beta_Rcore, -1.0, Eps_double, NoMax_double ); + + ReadPara->Add( "Bondi_Soliton", &Bondi_Soliton, false, Useless_bool, Useless_bool ); + ReadPara->Add( "Bondi_Soliton_m22", &Bondi_Soliton_m22, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Bondi_Soliton_type", &Bondi_Soliton_type, 5, 0, 5 ); + ReadPara->Add( "Bondi_Soliton_t", &Bondi_Soliton_t, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Bondi_Soliton_rc", &Bondi_Soliton_rc, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Bondi_Soliton_MassHalo", &Bondi_Soliton_MassHalo, -1.0, NoMin_double, NoMax_double ); + ReadPara->Add( "Bondi_Soliton_Redshift", &Bondi_Soliton_Redshift, -1.0, NoMin_double, NoMax_double ); ReadPara->Read( FileName ); @@ -241,6 +293,14 @@ void SetParameter() Aux_Message( stderr, "WARNING : OPT__BC_FLU[%d] != BC_FLU_OUTFLOW for non-HSE setup !?\n", s ); } + if ( Bondi_Soliton ) + { + if ( Bondi_Soliton_m22 <= 0.0 ) Aux_Error( ERROR_INFO, "Bondi_Soliton_m22 (%14.7e) <= 0.0 !!\n", Bondi_Soliton_m22 ); + if ( Bondi_Soliton_t < 0.0 ) Aux_Error( ERROR_INFO, "Bondi_Soliton_t (%14.7e) < 0.0 !!\n", Bondi_Soliton_t ); + if ( Bondi_Soliton_MassHalo <= 0.0 ) Aux_Error( ERROR_INFO, "Bondi_Soliton_MassHalo (%14.7e) <= 0.0 !!\n", Bondi_Soliton_MassHalo ); + if ( Bondi_Soliton_Redshift < 0.0 ) Aux_Error( ERROR_INFO, "Bondi_Soliton_Redshift (%14.7e) < 0.0 !!\n", Bondi_Soliton_Redshift ); + } + // (2) set the problem-specific derived parameters // must initialize EoS first @@ -329,7 +389,33 @@ void SetParameter() } // if ( Bondi_HSE ) -// (4) reset other general-purpose parameters +// (4) initialize the soliton setup + if ( Bondi_Soliton ) + { +// compute the soliton radius using the core-halo relation from Eq. 7 in Schive et al., PRL 113, 261302 (2014) + if ( Bondi_Soliton_rc < 0.0 ) + { + const double z = Bondi_Soliton_Redshift; + const double Mh = Bondi_Soliton_MassHalo*UnitExt_M/Const_Msun; // convert to Msun + const double H0 = 67.66*Const_km/Const_Mpc; // hard-coded for now + const double Om0 = 0.3111; // hard-coded for now + const double a0 = Om0*SQR(H0)/( 2.47e-5*SQR( 1e7/(Const_kpc*1e3) ) ); + const double H_H0_z = 1.0/( 1.0-Om0*( 1.0-pow(1.0+z,3.0)*(1.0+z+a0)/a0 ) ); + const double Om_z = Om0*pow(1.0+z,3.0)*H_H0_z; + const double Om_0 = Om0/( 1.0-Om0*( 1.0-(1.0+a0)/a0 ) ); + const double zeta_z = ( 18.0*SQR(M_PI) + 82.0*(Om_z-1.0) - 39.0*SQR(Om_z-1.0) )/Om_z; + const double zeta_0 = ( 18.0*SQR(M_PI) + 82.0*(Om_0-1.0) - 39.0*SQR(Om_0-1.0) )/Om_0; + + Bondi_Soliton_rc = 1.6/Bondi_Soliton_m22*pow( 1.0+z, -1.0/2.0 )*pow( zeta_z/zeta_0, -1.0/6.0 )*pow( Mh/1e9, -1.0/3.0 ); // in kpc + Bondi_Soliton_rc *= Const_kpc/UnitExt_L; // convert to external units to match the manually input value + } + + Bondi_Soliton_rc *= UnitExt_L/UNIT_L; // convert to internal units + Bondi_Soliton_t *= Bondi_TimeB; // input Bondi_Soliton_t is normalized to the Bondi time + } // if ( Bondi_Soliton ) + + +// (5) reset other general-purpose parameters // --> a helper macro PRINT_RESET_PARA is defined in Macro.h const long End_Step_Default = __INT_MAX__; const double End_T_Default = 1.0e1*Bondi_TimeB; // 10 Bondi time @@ -349,40 +435,51 @@ void SetParameter() if ( MPI_Rank == 0 ) { Aux_Message( stdout, "=============================================================================\n" ); - Aux_Message( stdout, " test problem ID = %d\n", TESTPROB_ID ); - Aux_Message( stdout, " Bondi_MassBH = %13.7e (%13.7e Msun)\n", Bondi_MassBH, Bondi_MassBH*UNIT_M/Const_Msun ); - Aux_Message( stdout, " Bondi_Rho0 = %13.7e (%13.7e g/cm^3)\n", Bondi_Rho0, Bondi_Rho0*UNIT_D ); - Aux_Message( stdout, " Bondi_T0 = %13.7e (%13.7e keV)\n", Bondi_T0, Bondi_T0*UNIT_E/Const_keV ); - Aux_Message( stdout, " Bondi_RefineRadius0 = %13.7e (%13.7e kpc)\n", Bondi_RefineRadius0, Bondi_RefineRadius0*UNIT_L/Const_kpc ); - Aux_Message( stdout, " Bondi_HalfMaxLvRefR = %s\n", (Bondi_HalfMaxLvRefR)?"YES":"NO" ); - Aux_Message( stdout, " Bondi_InBC_Rho = %13.7e (%13.7e g/cm^3)\n", Bondi_InBC_Rho, Bondi_InBC_Rho*UNIT_D ); - Aux_Message( stdout, " Bondi_InBC_T = %13.7e (%13.7e keV)\n", Bondi_InBC_T, Bondi_InBC_T*UNIT_E/Const_keV ); - Aux_Message( stdout, " Bondi_InBC_NCell = %13.7e\n", Bondi_InBC_NCell ); - Aux_Message( stdout, " Bondi_InBC_R = %13.7e (%13.7e kpc)\n", Bondi_InBC_R, Bondi_InBC_R*UNIT_L/Const_kpc ); - Aux_Message( stdout, " Bondi_InBC_E = %13.7e\n", Bondi_InBC_E ); - Aux_Message( stdout, " Bondi_Soften_NCell = %13.7e\n", Bondi_Soften_NCell ); - Aux_Message( stdout, " Bondi_Soften_R = %13.7e (%13.7e kpc)\n", Bondi_Soften_R, Bondi_Soften_R*UNIT_L/Const_kpc ); - Aux_Message( stdout, " Bondi_Cs = %13.7e (%13.7e km/s)\n", Bondi_Cs, Bondi_Cs*UNIT_V/Const_km ); - Aux_Message( stdout, " Schwarzschild radius = %13.7e (%13.7e kpc)\n", Bondi_RS, Bondi_RS*UNIT_L/Const_kpc ); - Aux_Message( stdout, " Bondi radius = %13.7e (%13.7e kpc)\n", Bondi_RB, Bondi_RB*UNIT_L/Const_kpc ); - Aux_Message( stdout, " Bondi time = %13.7e (%13.7e Myr)\n", Bondi_TimeB, Bondi_TimeB*UNIT_T/Const_Myr ); - - Aux_Message( stdout, " Bondi_HSE = %s\n", (Bondi_HSE)?"YES":"NO" ); + Aux_Message( stdout, " test problem ID = %d\n", TESTPROB_ID ); + Aux_Message( stdout, " Bondi_MassBH = %13.7e (%13.7e Msun)\n", Bondi_MassBH, Bondi_MassBH*UNIT_M/Const_Msun ); + Aux_Message( stdout, " Bondi_Rho0 = %13.7e (%13.7e g/cm^3)\n", Bondi_Rho0, Bondi_Rho0*UNIT_D ); + Aux_Message( stdout, " Bondi_T0 = %13.7e (%13.7e keV)\n", Bondi_T0, Bondi_T0*UNIT_E/Const_keV ); + Aux_Message( stdout, " Bondi_RefineRadius0 = %13.7e (%13.7e kpc)\n", Bondi_RefineRadius0, Bondi_RefineRadius0*UNIT_L/Const_kpc ); + Aux_Message( stdout, " Bondi_HalfMaxLvRefR = %s\n", (Bondi_HalfMaxLvRefR)?"YES":"NO" ); + Aux_Message( stdout, " Bondi_InBC_Rho = %13.7e (%13.7e g/cm^3)\n", Bondi_InBC_Rho, Bondi_InBC_Rho*UNIT_D ); + Aux_Message( stdout, " Bondi_InBC_T = %13.7e (%13.7e keV)\n", Bondi_InBC_T, Bondi_InBC_T*UNIT_E/Const_keV ); + Aux_Message( stdout, " Bondi_InBC_NCell = %13.7e\n", Bondi_InBC_NCell ); + Aux_Message( stdout, " Bondi_InBC_R = %13.7e (%13.7e kpc)\n", Bondi_InBC_R, Bondi_InBC_R*UNIT_L/Const_kpc ); + Aux_Message( stdout, " Bondi_InBC_E = %13.7e\n", Bondi_InBC_E ); + Aux_Message( stdout, " Bondi_Soften_NCell = %13.7e\n", Bondi_Soften_NCell ); + Aux_Message( stdout, " Bondi_Soften_R = %13.7e (%13.7e kpc)\n", Bondi_Soften_R, Bondi_Soften_R*UNIT_L/Const_kpc ); + Aux_Message( stdout, " Bondi_Cs = %13.7e (%13.7e km/s)\n", Bondi_Cs, Bondi_Cs*UNIT_V/Const_km ); + Aux_Message( stdout, " Schwarzschild radius = %13.7e (%13.7e kpc)\n", Bondi_RS, Bondi_RS*UNIT_L/Const_kpc ); + Aux_Message( stdout, " Bondi radius = %13.7e (%13.7e kpc)\n", Bondi_RB, Bondi_RB*UNIT_L/Const_kpc ); + Aux_Message( stdout, " Bondi time = %13.7e (%13.7e Myr)\n", Bondi_TimeB, Bondi_TimeB*UNIT_T/Const_Myr ); + Aux_Message( stdout, " Bondi_void = %s\n", (Bondi_void)?"YES":"NO" ); + Aux_Message( stdout, " Bondi_dynBH = %s\n", (Bondi_dynBH)?"YES":"NO" ); + + Aux_Message( stdout, " Bondi_HSE = %s\n", (Bondi_HSE)?"YES":"NO" ); if ( Bondi_HSE ) { - Aux_Message( stdout, " Bondi_HSE_Mode = %d\n", Bondi_HSE_Mode ); - Aux_Message( stdout, " Bondi_HSE_Dens_NBin = %d\n", Bondi_HSE_Dens_NBin ); - Aux_Message( stdout, " Bondi_HSE_Dens_MinR = %13.7e (%13.7e kpc)\n", Bondi_HSE_Dens_MinR, Bondi_HSE_Dens_MinR*UNIT_L/Const_kpc ); - Aux_Message( stdout, " Bondi_HSE_Dens_MaxR = %13.7e (%13.7e kpc)\n", Bondi_HSE_Dens_MaxR, Bondi_HSE_Dens_MaxR*UNIT_L/Const_kpc ); - Aux_Message( stdout, " Bondi_HSE_Dens_NormR = %13.7e (%13.7e kpc)\n", Bondi_HSE_Dens_NormR, Bondi_HSE_Dens_NormR*UNIT_L/Const_kpc ); - Aux_Message( stdout, " Bondi_HSE_Dens_NormD = %13.7e (%13.7e g/cm^3)\n", Bondi_HSE_Dens_NormD, Bondi_HSE_Dens_NormD*UNIT_D ); - Aux_Message( stdout, " Bondi_HSE_Truncate = %s\n", (Bondi_HSE_Truncate)?"YES":"NO" ); - Aux_Message( stdout, " Bondi_HSE_TrunR = %13.7e (%13.7e kpc)\n", Bondi_HSE_TrunR, Bondi_HSE_TrunR*UNIT_L/Const_kpc ); - Aux_Message( stdout, " Bondi_HSE_TrunD = %13.7e (%13.7e g/cm^3)\n", Bondi_HSE_TrunD, Bondi_HSE_TrunD*UNIT_D ); - Aux_Message( stdout, " Bondi_HSE_TrunSmoothR = %13.7e (%13.7e kpc)\n", Bondi_HSE_TrunSmoothR, Bondi_HSE_TrunSmoothR*UNIT_L/Const_kpc ); - Aux_Message( stdout, " Bondi_HSE_Pres_NormT = %s\n", (Bondi_HSE_Pres_NormT)?"YES":"NO" ); - Aux_Message( stdout, " Bondi_HSE_Beta = %13.7e\n", Bondi_HSE_Beta ); - Aux_Message( stdout, " Bondi_HSE_Beta_Rho0 = %13.7e (%13.7e g/cm^3)\n", Bondi_HSE_Beta_Rho0, Bondi_HSE_Beta_Rho0*UNIT_D ); - Aux_Message( stdout, " Bondi_HSE_Beta_Rcore = %13.7e (%13.7e kpc)\n", Bondi_HSE_Beta_Rcore, Bondi_HSE_Beta_Rcore*UNIT_L/Const_kpc ); } + Aux_Message( stdout, " Bondi_HSE_Mode = %d\n", Bondi_HSE_Mode ); + Aux_Message( stdout, " Bondi_HSE_Dens_NBin = %d\n", Bondi_HSE_Dens_NBin ); + Aux_Message( stdout, " Bondi_HSE_Dens_MinR = %13.7e (%13.7e kpc)\n", Bondi_HSE_Dens_MinR, Bondi_HSE_Dens_MinR*UNIT_L/Const_kpc ); + Aux_Message( stdout, " Bondi_HSE_Dens_MaxR = %13.7e (%13.7e kpc)\n", Bondi_HSE_Dens_MaxR, Bondi_HSE_Dens_MaxR*UNIT_L/Const_kpc ); + Aux_Message( stdout, " Bondi_HSE_Dens_NormR = %13.7e (%13.7e kpc)\n", Bondi_HSE_Dens_NormR, Bondi_HSE_Dens_NormR*UNIT_L/Const_kpc ); + Aux_Message( stdout, " Bondi_HSE_Dens_NormD = %13.7e (%13.7e g/cm^3)\n", Bondi_HSE_Dens_NormD, Bondi_HSE_Dens_NormD*UNIT_D ); + Aux_Message( stdout, " Bondi_HSE_Truncate = %s\n", (Bondi_HSE_Truncate)?"YES":"NO" ); + Aux_Message( stdout, " Bondi_HSE_TrunR = %13.7e (%13.7e kpc)\n", Bondi_HSE_TrunR, Bondi_HSE_TrunR*UNIT_L/Const_kpc ); + Aux_Message( stdout, " Bondi_HSE_TrunD = %13.7e (%13.7e g/cm^3)\n", Bondi_HSE_TrunD, Bondi_HSE_TrunD*UNIT_D ); + Aux_Message( stdout, " Bondi_HSE_TrunSmoothR = %13.7e (%13.7e kpc)\n", Bondi_HSE_TrunSmoothR, Bondi_HSE_TrunSmoothR*UNIT_L/Const_kpc ); + Aux_Message( stdout, " Bondi_HSE_Pres_NormT = %s\n", (Bondi_HSE_Pres_NormT)?"YES":"NO" ); + Aux_Message( stdout, " Bondi_HSE_Beta = %13.7e\n", Bondi_HSE_Beta ); + Aux_Message( stdout, " Bondi_HSE_Beta_Rho0 = %13.7e (%13.7e g/cm^3)\n", Bondi_HSE_Beta_Rho0, Bondi_HSE_Beta_Rho0*UNIT_D ); + Aux_Message( stdout, " Bondi_HSE_Beta_Rcore = %13.7e (%13.7e kpc)\n", Bondi_HSE_Beta_Rcore, Bondi_HSE_Beta_Rcore*UNIT_L/Const_kpc ); } + + Aux_Message( stdout, " Bondi_Soliton = %s\n", (Bondi_Soliton)?"YES":"NO" ); + if( Bondi_Soliton ) { + Aux_Message( stdout, " Bondi_Soliton_m22 = %13.7e\n", Bondi_Soliton_m22 ); + Aux_Message( stdout, " Bondi_Soliton_type = %d\n", Bondi_Soliton_type ); + Aux_Message( stdout, " Bondi_Soliton_t = %13.7e (%13.7e Myr)\n", Bondi_Soliton_t, Bondi_Soliton_t*UNIT_T/Const_Myr ); + Aux_Message( stdout, " Bondi_Soliton_rc = %13.7e (%13.7e kpc)\n", Bondi_Soliton_rc, Bondi_Soliton_rc*UNIT_L/Const_kpc ); + Aux_Message( stdout, " Bondi_Soliton_MassHalo = %13.7e Msun\n", Bondi_Soliton_MassHalo*UnitExt_M/Const_Msun ); + Aux_Message( stdout, " Bondi_Soliton_Redshift = %13.7e\n", Bondi_Soliton_Redshift ); } Aux_Message( stdout, "=============================================================================\n" ); } // if ( MPI_Rank == 0 ) @@ -463,7 +560,6 @@ void SetGridIC( real fluid[], const double x, const double y, const double z, co Aux_Error( ERROR_INFO, "unsupported Bondi_HSE_Mode (%d) !!\n", Bondi_HSE_Mode ); } // if ( Bondi_HSE ) - // uniform background else { @@ -585,6 +681,7 @@ void HSE_SetDensProfileTable() } // FUNCTION : HSE_SetDensProfileTable + //------------------------------------------------------------------------------------------------------- // Function : End_Bondi // Description : Free memory before terminating the program @@ -665,15 +762,16 @@ void Init_TestProb_Hydro_Bondi() // set the function pointers of various problem-specific routines - Init_Function_User_Ptr = SetGridIC; - Flag_User_Ptr = Flag_Bondi; - Aux_Record_User_Ptr = Record_Bondi; - BC_User_Ptr = BondiBC; - Flu_ResetByUser_Func_Ptr = Flu_ResetByUser_Func_Bondi; - Flu_ResetByUser_API_Ptr = Flu_ResetByUser_API_Bondi; - End_User_Ptr = End_Bondi; + Init_Function_User_Ptr = SetGridIC; + Flag_User_Ptr = Flag_Bondi; + Aux_Record_User_Ptr = Record_Bondi; + BC_User_Ptr = BondiBC; + Flu_ResetByUser_Func_Ptr = Flu_ResetByUser_Func_Bondi; + Flu_ResetByUser_API_Ptr = Flu_ResetByUser_API_Bondi; + End_User_Ptr = End_Bondi; # ifdef GRAVITY - Init_ExtAcc_Ptr = Init_ExtAcc_Bondi; + Init_ExtAcc_Ptr = Init_ExtAcc_Bondi; + Poi_UserWorkBeforePoisson_Ptr = Poi_UserWorkBeforePoisson_Bondi; # endif # endif // #if ( MODEL == HYDRO && defined GRAVITY ) diff --git a/src/TestProblem/Hydro/Bondi/Record_Bondi.cpp b/src/TestProblem/Hydro/Bondi/Record_Bondi.cpp index e6173e281..e597ab698 100644 --- a/src/TestProblem/Hydro/Bondi/Record_Bondi.cpp +++ b/src/TestProblem/Hydro/Bondi/Record_Bondi.cpp @@ -6,6 +6,7 @@ // specific global variables declared in Init_TestProb_Hydro_Bondi.cpp // ======================================================================================= +extern double Bondi_MassBH; extern double Bondi_SinkMass; extern double Bondi_SinkMomX; extern double Bondi_SinkMomY; @@ -50,10 +51,10 @@ void Record_Bondi() if ( Aux_CheckFileExist(FileName) ) Aux_Message( stderr, "WARNING : file \"%s\" already exists !!\n", FileName ); FILE *File_User = fopen( FileName, "a" ); - fprintf( File_User, "#%9s%16s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s\n", + fprintf( File_User, "#%9s%16s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s%20s\n", "Step", "Time [yr]", "NVoidCell", "Mass [Msun]", "Time [yr]", "dM/dt [Msun/yr]", "MomX [g*cm/s]", "MomY [g*cm/s]", "MomZ [g*cm/s]", "MomXAbs [g*cm/s]", "MomYAbs [g*cm/s]", "MomZAbs [g*cm/s]", - "Ek [erg]", "Et [erg]" ); + "Ek [erg]", "Et [erg]", "BHMass [Msun]" ); fclose( File_User ); } @@ -70,7 +71,7 @@ void Record_Bondi() // get the total amount of sunk variables - double Mass_Sum, MomX_Sum, MomY_Sum, MomZ_Sum, MomXAbs_Sum, MomYAbs_Sum, MomZAbs_Sum, Ek_Sum, Et_Sum; + double Mass_Sum, MomX_Sum, MomY_Sum, MomZ_Sum, MomXAbs_Sum, MomYAbs_Sum, MomZAbs_Sum, Ek_Sum, Et_Sum, Mass_BH; MPI_Reduce( &Bondi_SinkMass, &Mass_Sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); MPI_Reduce( &Bondi_SinkMomX, &MomX_Sum, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD ); @@ -91,6 +92,7 @@ void Record_Bondi() MomZAbs_Sum *= UNIT_M*UNIT_L/UNIT_T; Ek_Sum *= UNIT_E; Et_Sum *= UNIT_E; + Mass_BH = Bondi_MassBH*UNIT_M/Const_Msun; if ( MPI_Rank == 0 ) { @@ -99,9 +101,9 @@ void Record_Bondi() dTime *= UNIT_T/Const_yr; FILE *File_User = fopen( FileName, "a" ); - fprintf( File_User, "%10ld%16.7e%20d%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e\n", + fprintf( File_User, "%10ld%16.7e%20d%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e%20.7e\n", Step, Time[0]*UNIT_T/Const_yr, SinkNCell_Sum, Mass_Sum, dTime, Mass_Sum/dTime, - MomX_Sum, MomY_Sum, MomZ_Sum, MomXAbs_Sum, MomYAbs_Sum, MomZAbs_Sum, Ek_Sum, Et_Sum ); + MomX_Sum, MomY_Sum, MomZ_Sum, MomXAbs_Sum, MomYAbs_Sum, MomZAbs_Sum, Ek_Sum, Et_Sum, Mass_BH ); fclose( File_User ); } } // if ( FirstTime ) ... else ...