Add soliton potential for Bondi accretion problem

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-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]

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'
 
 

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 )

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,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;
@@ -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