Add soliton potential for Bondi accretion problem

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
@@ -43,6 +48,38 @@ void SetExtAccAuxArray_Bondi( double AuxArray[], const double Time )
    AuxArray[2] = amr->BoxCenter[2];
    AuxArray[3] = NEWTON_G*Bondi_MassBH;   // gravitational_constant*point_source_mass
    AuxArray[4] = Bondi_Soften_R;          // soften_length (<=0.0 --> disable)
+   AuxArray[5] = Bondi_Soliton;
+   if( Bondi_Soliton )
+   {
+      AuxArray[6] = Bondi_Soliton_m22;
+      AuxArray[7] = Bondi_Soliton_rc;
+      AuxArray[8] = NEWTON_G*Const_Msun/UNIT_M;
+      AuxArray[9] = UNIT_L/Const_kpc;
+      if( Bondi_Soliton_type > 0 )
+      {
+          switch(Bondi_Soliton_type)
+          {
+              case 1: // arctan function
+                  AuxArray[8] *= 2/(real)3.14159265*ATAN(Time/Bondi_Soliton_t);
+                  break;
+              case 2: // linear function
+                  if( Time < Bondi_Soliton_t )
+                      AuxArray[8] *= Time/Bondi_Soliton_t;
+                  break;
+              case 3: // smooth step function
+                  if( Time < Bondi_Soliton_t )
+                      AuxArray[8] *= 3*SQR(Time/Bondi_Soliton_t)-2*CUBE(Time/Bondi_Soliton_t);
+                  break;
+              case 4: // sigmoid
+                  AuxArray[8] *= 2/(1+exp(-Time*log(3)/Bondi_Soliton_t))-1;
+                  break;
+              case 5: // tanh
+                  AuxArray[8] *= tanh(Time/Bondi_Soliton_t);
+                  break;
+
+          }
+      } 
+   }
 
 } // FUNCTION : SetExtAccAuxArray_Bondi
 #endif // #ifndef __CUDACC__
@@ -73,13 +110,31 @@ static void ExtAcc_Bondi( real Acc[], const double x, const double y, const doub
 {
 
    const double Cen[3] = { UserArray[0], UserArray[1], UserArray[2] };
-   const real GM       = (real)UserArray[3];
+         real GM       = (real)UserArray[3];
    const real eps      = (real)UserArray[4];
+   const bool SOL      = (real)UserArray[5];
    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( SOL )
+   {
+      const real m22      = (real)UserArray[6];
+      const real rc       = (real)UserArray[7];  // In code unit or kpc
+      const real Coeff    = (real)UserArray[8];
+      const real UNIT_L   = (real)UserArray[9];
+
+#ifdef Plummer
+      double M = GM*CUBE(r)/pow(SQR(r)+SQR(rc),1.5);
+#else
+      real a = SQRT(POW(2.0,1.0/8.0)-1)*(r/rc);
+      real M = (real)4.17e9/(SQR(m22/1e-1)*(rc*UNIT_L*1e3)*POW(SQR(a)+1, 7.0))*((real)3465*POW(a,13.0)+(real)23100*POW(a,11.0)+(real)65373*POW(a,9.0)+(real)101376*POW(a,7.0)+(real)92323*POW(a,5.0)+(real)48580*POW(a,3.0)-(real)3465*a+(real)3465*POW(SQR(a)+1, 7.0)*ATAN(a));
+      M *= Coeff;
+#endif
+      GM += M;
+   }
+
 // Plummer
 #  if   ( defined SOFTEN_PLUMMER )
    const real _r3 = ( eps <= (real)0.0 ) ? (real)1.0/CUBE(r) : POW( SQR(r)+SQR(eps), (real)-1.5 );

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,8 +18,8 @@ extern double Bondi_SinkEk;
 extern double Bondi_SinkEt;
 extern int    Bondi_SinkNCell;
 
-
-
+extern bool   Bondi_void;
+extern bool   Bondi_dynBH;
 
 //-------------------------------------------------------------------------------------------------------
 // Function    :  Flu_ResetByUser_Func_Bondi
@@ -52,6 +51,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 +115,17 @@ 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;
+// Define varible to record sink mass at every time step
+   double SinkMass_OneSubStep_ThisRank = 0;
+   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 +198,23 @@ 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++)
 
+   MPI_Allreduce( &SinkMass_OneSubStep_ThisRank, &SinkMass_OneSubStep_AllRank, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
+
+   if( Bondi_dynBH ){ Bondi_MassBH += SinkMass_OneSubStep_AllRank; }
+
 } // FUNCTION : Flu_ResetByUser_API_Bondi