Full implementation of AGN feedback for the cool-core destruction project

example/test_problem/Hydro/ClusterMerger/Input__Flag_User

@@ -0,0 +1,22 @@
+# Level           User-defined Criteria
+      0                             6.0
+      1                             6.0
+      2                             6.0
+      3                             6.0
+      4                             6.0
+      5                             6.0
+      6                             6.0
+      7                             6.0
+      8                             6.0
+      9                             6.0
+     10                             6.0
+     11                             6.0
+     12                             6.0
+     13                             6.0
+     14                             6.0
+     15                             6.0
+     16                             6.0
+     17                             6.0
+     18                             6.0
+     19                             6.0
+     20                             6.0

example/test_problem/Hydro/ClusterMerger/Input__Parameter

@@ -124,9 +124,20 @@ OPT__MINIMIZE_MPI_BARRIER     0           # minimize MPI barriers to improve loa
                                           # (STORE_POT_GHOST, PAR_IMPROVE_ACC=1, OPT__TIMING_BARRIER=0 only; recommend AUTO_REDUCE_DT=0)
 
 
+# source terms                                                                                            
+SRC_DELEPTONIZATION           0           # deleptonization (for simulations of stellar core collapse) [0] ##HYDRO ONLY##
+SRC_EXACTCOOLING              0           # exact cooling scheme from Gaspari (2009) [0] ##HYDRO ONLY##  
+SRC_EC_TEF_N                  6001        # number of points for lambda(T) sampling in LOG
+SRC_EC_SUBCYCLING             1           # perform subcycling when the cooling time step is small [0=off, 1=on] 
+SRC_EC_DTCOEF                 0.5         # the coefficient of the cooling time step
+SRC_USER                      0           # user-defined source terms -> edit "Src_User.cpp" [0]       
+SRC_GPU_NPGROUP              -1           # number of patch groups sent into the CPU/GPU source-term solver (<=0=auto) [-1]
+
+
 # fluid solver in HYDRO (MODEL==HYDRO only)
 GAMMA                         1.666666666667       # ratio of specific heats (i.e., adiabatic index) [5.0/3.0]
 MOLECULAR_WEIGHT              0.59242761692650336  # mean molecular weight -> currently only for post-processing [0.6]
+MU_NORM                       1.6737352238051868e-24   # m_H value, fully consistent with yt 
 MINMOD_COEFF                  2.0         # coefficient of the generalized MinMod limiter (1.0~2.0) [1.5]
 MINMOD_MAX_ITER               0           # maximum number of iterations to reduce MINMOD_COEFF when data reconstruction fails (0=off) [0]
 OPT__LR_LIMITER               1           # slope limiter of data reconstruction in the MHM/MHM_RP/CTU schemes:
@@ -146,9 +157,12 @@ OPT__FIXUP_RESTRICT           1           # correct coarse grids by averaging th
 OPT__CORR_AFTER_ALL_SYNC     -1           # apply various corrections after all levels are synchronized (see "Flu_CorrAfterAllSync"):
                                           # (-1=auto, 0=off, 1=every step, 2=before dump) [-1]
 OPT__NORMALIZE_PASSIVE        1           # ensure "sum(passive_scalar_density) == gas_density" [1]
-MIN_DENS                      1.0e-07     # minimum mass density (must >= 0.0) [0.0] ##HYDRO, MHD, and ELBDM ONLY##
-MIN_PRES                      1.0e-09     # minimum pressure     (must >= 0.0) [0.0] ##HYDRO and MHD ONLY##
-MIN_EINT                      1.0e-09     # minimum internal energy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY##
+OPT__RESET_FLUID              0           # reset fluid variables after each update -> edit "Flu_ResetByUser.cpp" [0]
+OPT__RESET_FLUID_INIT         0           # reset fluid variables during initialization (<0=auto -> OPT__RESET_FLUID, 0=off, 1=on) [-1]
+MIN_DENS                      0.0         # minimum mass density (must >= 0.0) [0.0] ##HYDRO, MHD, and ELBDM ONLY##
+MIN_PRES                      0.0         # minimum pressure     (must >= 0.0) [0.0] ##HYDRO and MHD ONLY##
+MIN_EINT                      0.0         # minimum internal energy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY##
+MIN_TEMP                      1.0e4       # minimum temperature in K (must >= 0.0) [0.0] ##HYDRO and MHD ONLY##
 MIN_ENTR                      0.0         # minimum entropy         (must >= 0.0) [0.0] ##HYDRO and MHD ONLY##
 
 
@@ -230,6 +244,7 @@ OPT__MANUAL_CONTROL           1           # support manually dump data or stop r
                                           # (by generating the file DUMP_GAMER_DUMP or STOP_GAMER_STOP) [1]
 OPT__RECORD_USER              0           # record the user-specified info -> edit "Aux_Record_User.cpp" [0]
 OPT__OPTIMIZE_AGGRESSIVE      0           # apply aggressive optimizations (experimental) [0]
+OPT__SORT_PATCH_BY_LBIDX      1           # sort patches to improve bitwise reproducibility [SERIAL:0, LOAD_BALACNE:1]
 
 
 # checks

example/test_problem/Hydro/ClusterMerger/Input__TestProb

@@ -1,32 +1,68 @@
 # problem-specific runtime parameters
 Merger_Coll_NumHalos    2 # number of halos
+AGN_feedback            0 # turn on/off (1/0) AGN feedback
 
 # parameters for cluster 1
 Merger_File_Prof1       profile1_gamer.h5    # profile table of cluster 1
 Merger_File_Par1        1to3_b0.0_gamerp_1.h5    # particle file of cluster 1
+Unit_R1                 1.0 # the unit of length in cluster 1's IC file (in cgs)           
+Unit_D1                 1.0 # the unit of density in cluster 1's IC file (in cgs)          
+Unit_P1                 1.0 # the unit of pressure in cluster 1's IC file (in cgs) 
 Merger_Coll_PosX1       6000.0    # X-center of cluster 1 in kpc
 Merger_Coll_PosY1       7500.0    # Y-center of cluster 1 in kpc
 Merger_Coll_VelX1       375.0    # X-velocity of cluster 1 in km/s
 Merger_Coll_VelY1       0.0    # Y-velocity of cluster 1 in km/s
 Merger_Coll_IsGas1      1 # If cluster 1 should have gas and not have DM only
+Bondi_MassBH1           3.4e8 # black hole mass (in Msun)
+Mdot_BH1                0.0 # accretion rate of cluster 1 in Msun/yr
+Jet_HalfHeight1         0.70 # jet 1: half height of the cylinder-shape jet source (in kpc) 
+Jet_Radius1             0.35 # jet 1: radius of the cylinder-shape jet source (in kpc)
 
 # parameters for cluster 2
 Merger_File_Prof2       profile2_gamer.h5    # profile table of cluster 2
 Merger_File_Par2        1to3_b0.0_gamerp_2.h5    # particle file of cluster 2
+Unit_R2                 1.0 # the unit of length in cluster 2's IC file (in cgs)
+Unit_D2                 1.0 # the unit of density in cluster 2's IC file (in cgs)
+Unit_P2                 1.0 # the unit of pressure in cluster 2's IC file (in cgs)
 Merger_Coll_PosX2       9000.0    # X-center of cluster 2 in kpc
 Merger_Coll_PosY2       7500.0    # Y-center of cluster 2 in kpc
 Merger_Coll_VelX2       -1125.0000000000002    # X-velocity of cluster 2 in km/s
 Merger_Coll_VelY2       0.0    # Y-velocity of cluster 2 in km/s
 Merger_Coll_IsGas2      1 # If cluster 2 should have gas and not have DM only
+Bondi_MassBH2           3.4e8 # black hole mass (in Msun)
+Mdot_BH2                0.0 # accretion rate of cluster 2 in Msun/yr
+Jet_HalfHeight2         0.70 # jet 2: half height of the cylinder-shape jet source (in kpc) 
+Jet_Radius2             0.35 # jet 2: radius of the cylinder-shape jet source (in kpc) 
 
 # parameters for cluster 3 (disabled for this example)
 Merger_File_Prof3       none # profile table of cluster 3
 Merger_File_Par3        none # particle file of cluster 3
+Unit_R3                 1.0 # the unit of length in cluster 3's IC file (in cgs)           
+Unit_D3                 1.0 # the unit of density in cluster 3's IC file (in cgs)          
+Unit_P3                 1.0 # the unit of pressure in cluster 3's IC file (in cgs) 
 Merger_Coll_PosX3       0 # X-center of cluster 3 in kpc
 Merger_Coll_PosY3       0.0 # Y-center of cluster 3 in kpc
 Merger_Coll_VelX3       0.0 # X-velocity of cluster 3 in km/s
 Merger_Coll_VelY3       0.0 # Y-velocity of cluster 3 in km/s
 Merger_Coll_IsGas3      0 # If cluster 3 should have gas and not have DM only
+Bondi_MassBH3           3.4e8 # black hole mass (in Msun)
+Mdot_BH3                0.0 # accretion rate of cluster 3 in Msun/yr
+Jet_HalfHeight3         0.70 # jet 3: half height of the cylinder-shape jet source (in kpc) 
+Jet_Radius3             0.35 # jet 3: radius of the cylinder-shape jet source (in kpc) 
 
-# other parameters
-Merger_Coll_UseMetals   0 # look for a metal field in the profile files
+# parameters of AGN jet feedback
+Accretion_Mode          1 # 1: hot mode; 2: code mode; 3: combine (hot + cold)
+eta                     1.0 # mass loading factor in jet feedback
+eps_f                   0.001 # the radiative efficiency in jet feedback
+eps_m                   0.0 # the fraction of total energy that goes into the thermal energy in jet feedback
+R_acc                   4.0 # accretion radius: compute the accretion rate (in kpc)
+R_dep                   4.0 # radius to deplete the accreted gas (in kpc)   
+JetDirection_case       1 # Methods for choosing the jet direction: 1. Fixed at x-axis; 2. Import from table (generate JetDirection.txt); 3. Align with angular momentum
+
+# other	parameters
+Merger_Coll_UseMetals   0 # look for a metal field in the profile files [1]
+Merger_Coll_LabelCenter 1 # label the particle closest to the center of each cluster [1]
+AdjustBHPos             1 # (true/false) --> Adjust the BH position
+AdjustBHVel             0 # (true/false) --> Adjust the BH velocity
+AdjustPeriod            10 # the time interval of adjustment (in Myr)
+fixBH                   0 # fix the BH at the simulation box center and set its velocity to be zero (1 cluster only)

example/test_problem/Hydro/ClusterMerger/README

@@ -24,6 +24,12 @@ Default setup:
 2. Default resolution ~ 15 kpc (MAX_LEVEL = 3)
 
 
+AGN feedback:
+========================================
+1. Enable OPT__RESET_FLUID in Input__Parameter and "AGN_feedback" in Input__TestProb.
+2. To refine regions around the SMBHs, enable OPT__FLAG_USER and set MAX_LEVEL accordingly.
+
+
 Note:
 ========================================
 1. Mimic the merging cluster simulation setup of Flash provided by John ZuHone

example/test_problem/Hydro/ClusterMerger/clean.sh

@@ -3,3 +3,5 @@ rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Du
       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
+
+rm -f Record__Center BH_variable.bin

example/test_problem/Hydro/ClusterMerger/yt_script/plot_profile.py

@@ -40,13 +40,13 @@
 
 
 # load data
-ds_gamer = yt.load( "../Data_000008" )
+ds_gamer = yt.load( "../Data_000000" )
 #ds_gamer = yt.load( "gamer/Data_000100" )
 #ds_flash = yt.load( "flash/fiducial_1to1_b0_hdf5_plt_cnt_0100" )
 
 
 # create a sphere on the max density location
-sp_gamer = ds_gamer.sphere( "max", (2.5e3, "kpc") )
+sp_gamer = ds_gamer.sphere( "max", (1.0e4, "kpc") )
 #sp_flash = ds_flash.sphere( "max", (2.5e3, "kpc") )
 
 
@@ -144,4 +144,4 @@
 
 # show/save figure
 plt.savefig( FileOut, bbox_inches='tight', pad_inches=0.05 )
-plt.show()
+#plt.show()

example/test_problem/Hydro/ClusterMerger/yt_script/plot_slice-z.py

@@ -35,11 +35,30 @@
 center_mode   = 'c'
 dpi           = 150
 
+
+# define the particle filters for SMBH particles
+def smbh1( pfilter, data ):
+   filter = data[ "all", "ParType" ] == 4
+   return filter
+
+def smbh2( pfilter, data ):
+   filter = data[ "all", "ParType" ] == 5
+   return filter
+
+yt.add_particle_filter( "smbh1", function=smbh1, filtered_type="all", requires=["ParType"] )
+yt.add_particle_filter( "smbh2", function=smbh2, filtered_type="all", requires=["ParType"] )
+
+
 yt.enable_parallelism()
 ts = yt.DatasetSeries( [ prefix+'/Data_%06d'%idx for idx in range(idx_start, idx_end+1, didx) ] )
 
 for ds in ts.piter():
 
+#  apply the particle filters to the dataset
+   ds.add_particle_filter( 'smbh1' )
+   ds.add_particle_filter( 'smbh2' )
+
+
 #  density
    field   = 'density'
    sz_dens = yt.SlicePlot( ds, 'z', field, center=center_mode )
@@ -74,6 +93,8 @@
    cdm_mass.set_cmap( field, colormap_cdm )
    cdm_mass.set_colorbar_label( field, 'Dark matter mass [$M_{\odot}$]' )
    cdm_mass.annotate_timestamp( time_unit='Myr', corner='upper_right' )
+   cdm_mass.annotate_particles( width=(10.0*Mpc), p_size=10.0, col='b', ptype='smbh1' )
+   cdm_mass.annotate_particles( width=(10.0*Mpc), p_size=10.0, col='r', ptype='smbh2' )
 
    cdm_mass.save( mpl_kwargs={"dpi":dpi} )