[WIP] Sink Particle

example/test_problem/Hydro/SinkParTest/Input__Flag_Jeans

@@ -0,0 +1,13 @@
+# Level        Jeans length / cell size
+      0                            8.0
+      1                            8.0
+      2                            8.0
+      3                            8.0
+      4                            8.0
+      5                            8.0
+      6                            8.0
+      7                            8.0
+      8                            8.0
+      9                            8.0
+     10                            8.0
+     11                            8.0

example/test_problem/Hydro/SinkParTest/Input__TestProb

@@ -0,0 +1,12 @@
+# problem-specific runtime parameters
+# The following parameters follow Hubber et al. 2018
+R0                           3.0856776e+16  # [cm]
+Omega0                       1.6e-12        # [s^-1]
+Core_Mass                    1.0            # [M_sun]
+Delta_Dens                   0.5            # m=2 density perturbation
+Dens_Contrast                100            # density contrast outside R0
+B0                           5e-6           # B field strenght [G]
+theta_B                      0.0            # angle between B and z-axis [degree]
+Mach_num                     0.0            # Mach number for turbulence
+rho_AD_SinkParTest           1e-14          # transition density for Barotropic EOS [g/cm^-3]
+Tur_Table                    Tur_Table      # the turbulence table

example/test_problem/Hydro/SinkParTest/VG_turb.py

@@ -0,0 +1,130 @@
+from pylab import *
+import sys
+
+'''
+Usage: python VG_turb.py [power law index/2] [kmin] [random seed]
+'''
+
+nmodes = int(128)
+
+grid = nmodes+1
+
+A0 = 1
+
+Ax = zeros((grid,grid,grid))
+Ay = zeros((grid,grid,grid))
+Az = zeros((grid,grid,grid))
+phasex = zeros((grid,grid,grid))
+phasey = zeros((grid,grid,grid))
+phasez = zeros((grid,grid,grid))
+
+n = float(sys.argv[1])
+kmin = int(sys.argv[2])
+sd = int(sys.argv[3])
+print(n,kmin,sd)
+seed(sd)
+
+half_modes = int(nmodes/2 + 1)
+
+# The Amplitudes of the modes are governed by a power law spectrum of the form A(k) = A0 * k**(-n)
+
+for ii in range(0,half_modes):
+	kz = ii
+	print(kz)
+	if(kz == 0):
+		for jj in range(0,half_modes):
+			ky = jj
+			if(ky==0):
+				for ll in range(0,half_modes):
+					kx = ll
+					k = sqrt(kx**2 + ky**2 + kz**2)
+					if(kx==0 or k<kmin):
+						Ax[ii,jj,ll] = 0
+						Ay[ii,jj,ll] = 0
+						Az[ii,jj,ll] = 0
+					else:
+						Ax[ii,jj,ll] = A0 * normal() * k**n
+						Ay[ii,jj,ll] = A0 * normal() * k**n
+						Az[ii,jj,ll] = A0 * normal() * k**n
+
+					phasex[ii,jj,ll]  = 2*pi*random()
+					phasey[ii,jj,ll]  = 2*pi*random()
+					phasez[ii,jj,ll]  = 2*pi*random()
+			if(ky>0):
+				for ll in range(0,nmodes + 1):
+					kx = ll
+					if(ll > nmodes/2):
+						kx = ll - nmodes - 1
+					k = sqrt(kx**2 + ky**2 + kz**2)
+
+					Ax[ii,jj,ll] = A0 * normal() * k**n
+					Ay[ii,jj,ll] = A0 * normal() * k**n
+					Az[ii,jj,ll] = A0 * normal() * k**n
+
+					if(k<kmin):
+						Ax[ii,jj,ll] = 0
+						Ay[ii,jj,ll] = 0
+						Az[ii,jj,ll] = 0
+
+					phasex[ii,jj,ll]  = 2*pi*random()
+					phasey[ii,jj,ll]  = 2*pi*random()
+					phasez[ii,jj,ll]  = 2*pi*random()	
+	if(kz>0):
+		for jj in range(0,nmodes + 1):
+			ky= jj
+			if(jj > nmodes/2):
+				ky = jj - nmodes - 1
+			for ll in range(0,nmodes + 1):
+				kx = ll
+				if(ll > nmodes/2):
+					kx = ll - nmodes - 1
+				k = sqrt(kx**2 + ky**2 + kz**2)
+
+				Ax[ii,jj,ll] = A0 * normal() * k**n
+				Ay[ii,jj,ll] = A0 * normal() * k**n
+				Az[ii,jj,ll] = A0 * normal() * k**n
+
+				if(k<kmin):
+					Ax[ii,jj,ll] = 0
+					Ay[ii,jj,ll] = 0
+					Az[ii,jj,ll] = 0
+
+				phasex[ii,jj,ll]  = 2*pi*random()
+				phasey[ii,jj,ll]  = 2*pi*random()
+				phasez[ii,jj,ll]  = 2*pi*random()
+
+
+
+Bx1 = Ax*cos(phasex)
+Bx2 = 1j*Ax*sin(phasex)
+
+vx = ifftn(Bx1 + Bx2)
+vx = vx.real
+
+By1 = Ay*cos(phasey)
+By2 = 1j*Ay*sin(phasey)
+
+vy = ifftn(By1 + By2)
+vy = vy.real
+
+Bz1 = Az*cos(phasez)
+Bz2 = 1j*Az*sin(phasez)
+
+vz = ifftn(Bz1 + Bz2)
+vz = vz.real
+
+x = linspace(-1,1,grid)
+y = linspace(-1,1,grid)
+z = linspace(-1,1,grid)
+
+f=open("v_field.dat","w")
+
+for ii in range(0,grid):
+	for jj in range(0,grid):
+		for ll in range(0,grid):
+			f.write("{0:.5g} {1:.5g} {2:.5g} {3:.5g} {4:.5g} {5:.5g}\n".format(x[ll],  y[jj],  z[ii],  vx[ii,jj,ll], vy[ii,jj,ll], vz[ii,jj,ll]))
+
+f.close()
+
+
+

example/test_problem/Hydro/SinkParTest/clean.sh

@@ -0,0 +1,6 @@
+rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Dump Record__MemInfo Record__L1Err \
+      Record__Conservation Data* stderr stdout log XYslice* YZslice* XZslice* Xline* Yline* Zline* \
+      Diag* BaseXYslice* BaseYZslice* BaseXZslice* BaseXline* BaseYline* BaseZline* BaseDiag* \
+      PowerSpec_* Particle_* nohup.out Record__Performance Record__TimingMPI_* \
+      Record__ParticleCount Record__User Patch_* Record__NCorrUnphy FailedPatchGroup* *.pyc Record__LoadBalance \
+      GRACKLE_INFO

include/Field.h

@@ -85,6 +85,7 @@ SET_GLOBAL( FieldIdx_t Idx_ParVelY,       Idx_Undefined );
 SET_GLOBAL( FieldIdx_t Idx_ParVelZ,       Idx_Undefined );
 SET_GLOBAL( FieldIdx_t Idx_ParTime,       Idx_Undefined );
 SET_GLOBAL( FieldIdx_t Idx_ParType,       Idx_Undefined );
+SET_GLOBAL( FieldIdx_t Idx_ParID,         Idx_Undefined );
 #ifdef STORE_PAR_ACC
 SET_GLOBAL( FieldIdx_t Idx_ParAccX,       Idx_Undefined );
 SET_GLOBAL( FieldIdx_t Idx_ParAccY,       Idx_Undefined );

include/Global.h

@@ -288,6 +288,9 @@ extern double                SF_CREATE_STAR_MIN_GAS_DENS;
 extern double                SF_CREATE_STAR_MASS_EFF;
 extern double                SF_CREATE_STAR_MIN_STAR_MASS;
 extern double                SF_CREATE_STAR_MAX_STAR_MFRAC;
+extern double                SF_CREATE_SINK_MIN_GAS_DENS;
+extern double                SF_CREATE_SINK_ACC_RADIUS;
+extern int                   SF_CREATE_SINK_MAX_NPAR_MPI;
 #endif
 
 
@@ -355,7 +358,7 @@ extern void (*Flu_DerivedField_User_Ptr)( real Out[], const real FluIn[], const
 // =======================================================================================================
 #ifdef FEEDBACK
 extern int  FB_LEVEL, FB_RSEED;
-extern bool FB_SNE, FB_USER;
+extern bool FB_SNE, FB_ACC, FB_USER;
 extern bool FB_Any;
 extern int  FB_ParaBuf;
 #endif

include/HDF5_Typedef.h

@@ -705,13 +705,17 @@ struct InputPara_t
    double SF_CreateStar_MassEff;
    double SF_CreateStar_MinStarMass;
    double SF_CreateStar_MaxStarMFrac;
+   double SF_CreateStar_SinkMinGasDens;
+   double SF_CreateStar_SinkAccRadius;
+   int    SF_CreateStar_SinkMaxNParMPI;
 #  endif
 
 // feedback
 #  ifdef FEEDBACK
    int   FB_Level;
    int   FB_RSeed;
    int   FB_SNe;
+   int   FB_Acc;
    int   FB_User;
 #  endif
 

include/Macro.h

@@ -828,7 +828,7 @@
 // --> can be changed manually
 // --> set to 0 if applicable to improve performance
 #ifdef FEEDBACK
-#        define FB_GHOST_SIZE       3
+#        define FB_GHOST_SIZE       8
 #endif
 
 

include/Prototype.h

@@ -829,6 +829,8 @@ void SF_FreeRNG();
 void SF_CreateStar_AGORA( const int lv, const real TimeNew, const real dt, RandomNumber_t *RNG,
                           const real GasDensThres, const real Efficiency, const real MinStarMass, const real MaxStarMFrac,
                           const bool DetRandom, const bool UseMetal );
+void SF_CreateStar_SinkParticle( const int lv, const real TimeNew, const real GasDensThres, const real AccCellNum, 
+                                 const int MaxNewPar );
 #endif
 
 
@@ -838,7 +840,11 @@ void FB_AdvanceDt( const int lv, const double TimeNew, const double TimeOld, con
                    const int SaveSg_Flu, const int SaveSg_Mag );
 void FB_Init();
 void FB_End();
-int FB_Aux_CellPatchRelPos( const int ijk[] );
+int  FB_Aux_CellPatchRelPos( const int ijk[] );
+int  FB_Accretion( const int lv, const real GasDensThres, const real AccCellNum, const int NPar, const long *ParSortID, real_par *ParAttFlt[PAR_NATT_FLT_TOTAL], 
+                   long_par *ParAttInt[PAR_NATT_INT_TOTAL], real (*Fluid)[FB_NXT][FB_NXT][FB_NXT], const double EdgeL[], const double dh, bool CoarseFine[] );
+void FB_Init_Accretion();
+void FB_End_Accretion();
 #endif
 
 

include/Typedef.h

@@ -97,6 +97,7 @@ const TestProbID_t
    TESTPROB_HYDRO_JET_ICM_WALL                 =   52,
    TESTPROB_HYDRO_CDM_LSS                      =  100,
    TESTPROB_HYDRO_ZELDOVICH                    =  101,
+   TESTPROB_HYDRO_SinkParTest                  =  150,
    TESTPROB_ELBDM_EXTPOT                       = 1000,
    TESTPROB_ELBDM_JEANS_INSTABILITY_COMOVING   = 1001,
    TESTPROB_ELBDM_JEANS_INSTABILITY_PHYSICAL   = 1002,
@@ -112,6 +113,7 @@ const TestProbID_t
    TESTPROB_ELBDM_HALO_MERGER                  = 1012,
    TESTPROB_ELBDM_DISK_HEATING                 = 1013;
 
+
 // program initialization options
 typedef int OptInit_t;
 const OptInit_t
@@ -524,7 +526,8 @@ const GracklePriChe_t
 typedef int SF_CreateStarScheme_t;
 const SF_CreateStarScheme_t
    SF_CREATE_STAR_SCHEME_NONE  = 0,
-   SF_CREATE_STAR_SCHEME_AGORA = 1;
+   SF_CREATE_STAR_SCHEME_AGORA = 1,
+   SF_CREATE_STAR_SCHEME_SINKPARTICLE = 2;
 #endif
 
 

src/Auxiliary/Aux_TakeNote.cpp

@@ -1144,17 +1144,21 @@ void Aux_TakeNote()
 #     ifdef STAR_FORMATION
       fprintf( Note, "Parameters of Star Formation\n" );
       fprintf( Note, "***********************************************************************************\n" );
-      fprintf( Note, "SF_CREATE_STAR_SCHEME          % d\n",           SF_CREATE_STAR_SCHEME                          );
+      fprintf( Note, "SF_CREATE_STAR_SCHEME          % d\n",           SF_CREATE_STAR_SCHEME                                           );
       if ( SF_CREATE_STAR_SCHEME != SF_CREATE_STAR_SCHEME_NONE ) {
-      fprintf( Note, "SF_CREATE_STAR_RSEED           % d\n",           SF_CREATE_STAR_RSEED                           );
-      fprintf( Note, "SF_CREATE_STAR_DET_RANDOM      % d\n",           SF_CREATE_STAR_DET_RANDOM                      );
-      fprintf( Note, "SF_CREATE_STAR_MIN_LEVEL       % d\n",           SF_CREATE_STAR_MIN_LEVEL                       );
-      fprintf( Note, "SF_CREATE_STAR_MIN_GAS_DENS    % 14.7e\n",       SF_CREATE_STAR_MIN_GAS_DENS                    );
-      fprintf( Note, "                              =% 14.7e cm^-3\n", SF_CREATE_STAR_MIN_GAS_DENS*UNIT_D/Const_mH    );
-      fprintf( Note, "SF_CREATE_STAR_MASS_EFF        % 14.7e\n",       SF_CREATE_STAR_MASS_EFF                        );
-      fprintf( Note, "SF_CREATE_STAR_MIN_STAR_MASS   % 14.7e\n",       SF_CREATE_STAR_MIN_STAR_MASS                   );
-      fprintf( Note, "                              =% 14.7e Msun\n",  SF_CREATE_STAR_MIN_STAR_MASS*UNIT_M/Const_Msun );
-      fprintf( Note, "SF_CREATE_STAR_MAX_STAR_MFRAC  % 14.7e\n",       SF_CREATE_STAR_MAX_STAR_MFRAC                  ); }
+      fprintf( Note, "SF_CREATE_STAR_RSEED           % d\n",           SF_CREATE_STAR_RSEED                                            );
+      fprintf( Note, "SF_CREATE_STAR_DET_RANDOM      % d\n",           SF_CREATE_STAR_DET_RANDOM                                       );
+      fprintf( Note, "SF_CREATE_STAR_MIN_LEVEL       % d\n",           SF_CREATE_STAR_MIN_LEVEL                                        );
+      fprintf( Note, "SF_CREATE_STAR_MIN_GAS_DENS    % 14.7e\n",       SF_CREATE_STAR_MIN_GAS_DENS                                     );
+      fprintf( Note, "                              =% 14.7e cm^-3\n", SF_CREATE_STAR_MIN_GAS_DENS*UNIT_D/Const_mH                     );
+      fprintf( Note, "SF_CREATE_STAR_MASS_EFF        % 14.7e\n",       SF_CREATE_STAR_MASS_EFF                                         );
+      fprintf( Note, "SF_CREATE_STAR_MIN_STAR_MASS   % 14.7e\n",       SF_CREATE_STAR_MIN_STAR_MASS                                    );
+      fprintf( Note, "                              =% 14.7e Msun\n",  SF_CREATE_STAR_MIN_STAR_MASS*UNIT_M/Const_Msun                  );
+      fprintf( Note, "SF_CREATE_STAR_MAX_STAR_MFRAC  % 14.7e\n",       SF_CREATE_STAR_MAX_STAR_MFRAC                                   ); 
+      fprintf( Note, "SF_CREATE_SINK_MIN_GAS_DENS    % 14.7e\n",       SF_CREATE_SINK_MIN_GAS_DENS                                     );
+      fprintf( Note, "                              =% 14.7e cm^-3\n", SF_CREATE_SINK_MIN_GAS_DENS*UNIT_D/Const_mH/MOLECULAR_WEIGHT    ); 
+      fprintf( Note, "SF_CREATE_SINK_ACC_RADIUS      % 14.7e\n",       SF_CREATE_SINK_ACC_RADIUS                                       );
+      fprintf( Note, "SF_CREATE_SINK_MAX_NPAR_MPI    % d\n",           SF_CREATE_SINK_MAX_NPAR_MPI                                     ); }
       fprintf( Note, "***********************************************************************************\n" );
       fprintf( Note, "\n\n" );
 #     endif // #ifdef STAR_FORMATION
@@ -1167,6 +1171,7 @@ void Aux_TakeNote()
       fprintf( Note, "FB_LEVEL                       % d\n",      FB_LEVEL                );
       fprintf( Note, "FB_RSEED                       % d\n",      FB_RSEED                );
       fprintf( Note, "FB_SNE                         % d\n",      FB_SNE                  );
+      fprintf( Note, "FB_ACC                         % d\n",      FB_ACC                  );
       fprintf( Note, "FB_USER                        % d\n",      FB_USER                 );
       fprintf( Note, "***********************************************************************************\n" );
       fprintf( Note, "\n\n" );