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

doc/wiki/Adding-New-Simulations.md

@@ -23,6 +23,7 @@ Mandatory steps are marked by 📌.
    *  [External Potential](#external-potential)
    *  [Equation of State](#equation-of-state)
    *  [Feedback](#feedback)
+   *  [HDF5 Output](#hdf5-output)
 7. [Add Problem-specific Validators](#vii-add-problem-specific-validators)
 8. [Store Problem-specific Input Files](#viii-store-problem-specific-input-files)
 
@@ -198,7 +199,34 @@ during the runtime.
     }
     ```
 
-4. Add these parameters to the input file `Input__TestProb`
+4. Edit the function `Output_HDF5_TestProb()` to store the variables in HDF5 files.
+   ```c++
+   void Output_HDF5_TestProb( HDF5_Output_t *HDF5_InputTest )
+   {
+
+      HDF5_InputTest->Add( "var_bool",   &var_bool   );
+      HDF5_InputTest->Add( "var_double", &var_double );
+      HDF5_InputTest->Add( "var_int",    &var_int    );
+      HDF5_InputTest->Add( "var_str",     var_str    );
+
+   } // FUNCTION : Output_HDF5_TestProb
+   ```
+> [!NOTE]
+> You should contain all the variables in step 2.
+
+> [!CAUTION]
+> There should be at least one variable to store. Otherwise, it should be like
+> ```c++
+> void Output_HDF5_TestProb( HDF5_Output_t *HDF5_InputTest )
+> {
+>
+>    HDF5_InputTest->Add( "CDM_LSS_TestProb_ID",  &TESTPROB_ID );
+>
+> } // FUNCTION : Output_HDF5_TestProb
+> ```
+> See `src/TestProblem/Hydro/CDM_LSS/Init_TestProb_Hydro_CDM_LSS.cpp`.
+
+5. Add these parameters to the input file `Input__TestProb`
 (see [[Input__TestProb | Runtime-Parameters:-Input__TestProb]]
 for the file format).
 This file must be put in the same directory as the executable `gamer`
@@ -617,6 +645,20 @@ Add a user-specified feedback. See [[FB_USER | Runtime-Parameters:-Feedback#FB_U
 * **Example:**
    * `src/TestProblem/Hydro/Plummer/FB_Plummer.cpp`
 
+### HDF5 Output
+* **Description:**
+Store user-specified variables. Similar usage of as [Store Problem-specific Variables](#iv-store-problem-specific-variables).
+* **Prototype:**
+   * `void Output_HDF5_User_NewProblem( HDF5_Output_t *HDF5_OutUser );`
+* **Function Pointer:**
+   * `Output_HDF5_User_Ptr`
+* **Compilation Option:**
+[[SUPPORT_HDF5 | Installation: Simulation-Options#SUPPORT_HDF5]]
+* **Runtime Option:**
+None
+* **Example:**
+   * `Output/Output_DumData_Total_HDF5.cpp` --> `Output_HDF5_User_Template()`
+
 
 ## VII. Add Problem-specific Validators
 
@@ -677,7 +719,6 @@ and other relevant files such as README and analysis scripts.
 
 4. Add analysis scripts (if any).
 
-
 <br>
 
 ## Remarks

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## -> NOT supported yet!
+SRC_EC_TEF_N                  6001        # number of points for lambda(T) sampling in LOG [1501]
+SRC_EC_SUBCYCLING             1           # perform subcycling when the cooling time step is small (0=off, 1=on) [0]
+SRC_EC_DTCOEF                 0.5         # the coefficient of the cooling time step [0.5]
+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,10 +157,13 @@ 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##
-MIN_ENTR                      0.0         # minimum entropy         (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##
 
 
 # gravity solvers in all models
@@ -231,6 +245,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,62 @@
 # problem-specific runtime parameters
-Merger_Coll_NumHalos    2 # number of halos
+Merger_Coll_NumHalos    2                        # number of halos
 
 # parameters for cluster 1
-Merger_File_Prof1       profile1_gamer.h5    # profile table of 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
-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
+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_Prof2       profile2_gamer.h5        # profile table of cluster 2
 Merger_File_Par2        1to3_b0.0_gamerp_2.h5    # particle file of cluster 2
-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
+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
-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
+Merger_File_Prof3       none                     # profile table of cluster 3
+Merger_File_Par3        none                     # particle file of cluster 3
+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)
+
+# parameters of AGN jet feedback
+AGN_feedback            0                        # turn on/off (1/0) AGN 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
+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

@@ -2,4 +2,5 @@ rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Du
       Record__Conservation Data* stderr stdout log XYslice* YZslice* XZslice* Xline* Yline* Zline* \
       Diag* Box* 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
+      Record__ParticleCount Record__User Patch_* Record__NCorrUnphy FailedPatchGroup* *.pyc Record__LoadBalance \
+      Record__ClusterCenter

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} )