Support feedback from resolved Type II supernovae

doc/wiki/Runtime-Parameters-related/[Runtime-Parameters]-All.md

@@ -114,7 +114,11 @@ For variables with `Default/Min/Max` labeled as `Depend`, click the parameter na
 
 | Name                                                                                                 |         Default |             Min |             Max | Short description |
 | :---                                                                                                 |            :--- |            :--- |            :--- | :--- |
-| [FB_LEVEL](%5BRuntime-Parameters%5D-Feedback#FB_LEVEL)                                               |              -1 |            None |       TOP_LEVEL | AMR level to apply feedback (must be MAX_LEVEL for now; <0=auto -> MAX_LEVEL) [-1] |
+| [FB_MIN_LEVEL](%5BRuntime-Parameters%5D-Feedback#FB_MIN_LEVEL)                                       |               0 |            None |       TOP_LEVEL | minimum AMR level allowed to apply feedback (<0=auto -> MAX_LEVEL) [0] |
+| [FB_RESOLVED_SNEII](%5BRuntime-Parameters%5D-Feedback#FB_RESOLVED_SNEII)                             |               0 |            None |            None | Type II supernovae feedback for high enough grid and particle resolution (must add one extra PAR_NATT_FLT_USER) |
+| [FB_RESOLVED_SNEII_MIN_M_GAS](%5BRuntime-Parameters%5D-Feedback#FB_RESOLVED_SNEII_MIN_M_GAS)         |             0.0 |             0.0 |            None | minimum mass of gas in the environment of SNeII to apply feedback (in Msun) (for FB_RESOLVED_SNEII only) [0.0] |
+| [FB_RESOLVED_SNEII_N_PER_MASS](%5BRuntime-Parameters%5D-Feedback#FB_RESOLVED_SNEII_N_PER_MASS)       |          1.0e-2 |             0.0 |            None | number of SNeII per stellar mass (in 1/Msun) (for FB_RESOLVED_SNEII only) [0.01] |
+| [FB_RESOLVED_SNEII_RECORD](%5BRuntime-Parameters%5D-Feedback#FB_RESOLVED_SNEII_RECORD)               |               0 |            None |            None | record the SNeII explosion info (for FB_RESOLVED_SNEII only) [0] |
 | [FB_RSEED](%5BRuntime-Parameters%5D-Feedback#FB_RSEED)                                               |             456 |               0 |            None | random seed [456] |
 | [FB_SNE](%5BRuntime-Parameters%5D-Feedback#FB_SNE)                                                   |               0 |            None |            None | supernova explosion feedback [0] |
 | [FB_USER](%5BRuntime-Parameters%5D-Feedback#FB_USER)                                                 |               0 |            None |            None | user-defined feedback [0] |

doc/wiki/Runtime-Parameters-related/[Runtime-Parameters]-Feedback.md

@@ -1,19 +1,22 @@
 Parameters described on this page:
-[FB_LEVEL](#FB_LEVEL),  
+[FB_MIN_LEVEL](#FB_MIN_LEVEL),  
 [FB_RSEED](#FB_RSEED),  
 [FB_SNE](#FB_SNE),  
+[FB_RESOLVED_SNEII](#FB_RESOLVED_SNEII),  
 [FB_USER](#FB_USER)  
+[FB_RESOLVED_SNEII_N_PER_MASS](#FB_RESOLVED_SNEII_N_PER_MASS),  
+[FB_RESOLVED_SNEII_MIN_M_GAS](#FB_RESOLVED_SNEII_MIN_M_GAS),  
+[FB_RESOLVED_SNEII_RECORD](#FB_RESOLVED_SNEII_RECORD),  
 
 
 Parameters below are shown in the format:   **`Name`   (Valid Values)   [Default Value]**
 
 
-<a name="FB_LEVEL"></a>
-* #### `FB_LEVEL` &ensp; (0 &#8804; input < [[--nlevel | [Installation]-Option-List#--nlevel]]; <0 &#8594; set to [[MAX_LEVEL | [Runtime-Parameters]-Refinement#MAX_LEVEL ]]) &ensp; [-1]
+<a name="FB_MIN_LEVEL"></a>
+* #### `FB_MIN_LEVEL` &ensp; (0 &#8804; input < [[--nlevel | [Installation]-Option-List#--nlevel]]; <0 &#8594; set to [[MAX_LEVEL | [Runtime-Parameters]-Refinement#MAX_LEVEL ]]) &ensp; [0]
     * **Description:**
-AMR level to apply feedback.
+Minimum AMR level to apply feedback.
     * **Restriction:**
-Must be [[MAX_LEVEL | [Runtime-Parameters]-Refinement#MAX_LEVEL ]] for now.
 
 <a name="FB_RSEED"></a>
 * #### `FB_RSEED` &ensp; (&#8805;0) &ensp; [456]
@@ -28,13 +31,71 @@ Supernova explosion feedback.
     * **Restriction:**
 Not supported yet.
 
+<a name="FB_RESOLVED_SNEII"></a>
+* #### `FB_RESOLVED_SNEII` &ensp; (0=off, 1=on) &ensp; [0]
+    * **Description:**
+Resolved Type II supernovae feedback.
+When a star particle forms, it would be sampled to contain a SNII progenitor stochastically
+with a probability of star particle mass times
+number of SNII per stellar masses ([[FB_RESOLVED_SNEII_N_PER_MASS | Runtime Parameters:-Feedback#FB_RESOLVED_SNEII_N_PER_MASS ]]).
+The selected particle for supernova will then explode when the age of the star particle reaches
+the given explosion time.
+When the supernova explodes, it will uniformly inject
+the given amount of thermal energy, mass, and metal
+into fluid of the surrounding cells where the particle is located.
+The explosion delay time of SNeII after star formation, the ejected internal energy of each SNeII explosion,
+the ejected total mass of each SNeII explosion, and the ejected metal mass of each SNeII explosion should be set in
+`FeedbackYieldTable_Resolved_SNeII_*`.
+Note that the input value should always be in units of Myr for explosion time,
+in units of erg for energy, and in nits of Msun for mass.
+The actual amount of the ejected metal will not be greater than the ejected total mass.
+The actual amount of the ejected mass will not be greater than the mass of the particle.
+The number of cells to apply feedback will increase from one until the enclosed mass exceeds the given
+minimum gas mass ([[FB_RESOLVED_SNEII_MIN_M_GAS | Runtime Parameters:-Feedback#FB_RESOLVED_SNEII_MIN_M_GAS ]]), or
+the diameter equals `FB_GHOST_SIZE+1`, or the region crosses the coarse-fine boundary.
+See sec. 2.6 in [Chia-Yu Hu et al. 2023](https://doi.org/10.3847/1538-4357/accf9e) for reference.
+    * **Restriction:**
+Must set one extra particle attribute with [[ --par_attribute_flt | Installation:-Option-List#--par_attribute_flt ]].
+The star particle mass resolution must be high enough to have at most one supernova explosion per particle.
+The grid resolution must be high enough to resolve the Sedov phase of supernova explosion blast wave,
+so the kinetic (outward momentum) feedback is not needed.
+
 <a name="FB_USER"></a>
 * #### `FB_USER` &ensp; (0=off, 1=on) &ensp; [0]
     * **Description:**
 User-defined feedback.
 See [Add User-defined Feedback](#add-user-defined-feedback) for details.
     * **Restriction:**
 
+<a name="FB_RESOLVED_SNEII_N_PER_MASS"></a>
+* #### `FB_RESOLVED_SNEII_N_PER_MASS` &ensp; (&#8805;0.0) &ensp; [1.0e-2]
+    * **Description:**
+Number of SNeII per stellar mass.
+Note that the input value should always be in units of 1/Msun.
+    * **Restriction:**
+Only for [[FB_RESOLVED_SNEII | Runtime Parameters:-Feedback#FB_RESOLVED_SNEII ]].
+Its value times the maximum star particle mass cannot be greater than one.
+
+<a name="FB_RESOLVED_SNEII_MIN_M_GAS"></a>
+* #### `FB_RESOLVED_SNEII_MIN_M_GAS` &ensp; (&#8805;0.0) &ensp; [0.0]
+    * **Description:**
+Minimum mass of gas in the environment of SNeII to apply feedback.
+If the mass of the cell where the particle is in is sufficient, then the feedback is applied to a single cell.
+Otherwise, the diameter of the feedback region will increase from `1` to `FB_GHOST_SIZE+1`
+until the enclosed mass is higher than this given threshold.
+When the particle is too close to the AMR coarse-fine boundary,
+the feedback region will be limited such that no feedback will be applied in the coarse level.
+Note that the input value should always be in units of Msun.
+    * **Restriction:**
+Only for [[FB_RESOLVED_SNEII | Runtime Parameters:-Feedback#FB_RESOLVED_SNEII ]].
+
+<a name="FB_RESOLVED_SNEII_RECORD"></a>
+* #### `FB_RESOLVED_SNEII_RECORD` &ensp; (0=off, 1=on) &ensp; [0]
+    * **Description:**
+Record the information of the resolved Type II supernovae feedback.
+    * **Restriction:**
+Only for [[FB_RESOLVED_SNEII | Runtime Parameters:-Feedback#FB_RESOLVED_SNEII ]].
+
 
 ## Remarks
 

example/script/clean.sh

@@ -3,4 +3,4 @@ rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Du
       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__Center \
-      GRACKLE_INFO Record__DivB Record__Hybrid
+      GRACKLE_INFO Record__DivB Record__Hybrid Record__FB_Resolved_SNeII

example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input_Options/Input__Parameter.high-res

@@ -157,6 +157,17 @@ SF_CREATE_STAR_MIN_STAR_MASS  2.0e3       # minimum star particle mass for the s
 SF_CREATE_STAR_MAX_STAR_MFRAC 0.5         # maximum gas mass fraction allowed to convert to stars per substep [0.5]
 
 
+# feedback (FEEDBACK only)
+FB_MIN_LEVEL                 -1           # minimum AMR level allowed to apply feedback (<0=auto -> MAX_LEVEL) [0]
+FB_RSEED                      456         # random seed [456]
+FB_SNE                        0           # supernova explosion feedback [0]
+FB_RESOLVED_SNEII             1           # Type II supernovae feedback for high enough grid and particle resolution (must add one extra PAR_NATT_FLT_USER)
+FB_USER                       0           # user-defined feedback [0]
+FB_RESOLVED_SNEII_N_PER_MASS  1.0e-5      # number of SNeII per stellar mass (in 1/Msun) (for FB_RESOLVED_SNEII only) [0.01]
+FB_RESOLVED_SNEII_MIN_M_GAS   0.0         # minimum mass of gas in the environment of SNeII to apply feedback (in Msun) (for FB_RESOLVED_SNEII only) [0.0]
+FB_RESOLVED_SNEII_RECORD      0           # record the SNeII explosion info (for FB_RESOLVED_SNEII only) [0]
+
+
 # fluid solver in HYDRO (MODEL==HYDRO only)
 GAMMA                         1.666667    # ratio of specific heats (i.e., adiabatic index) [5.0/3.0]
 MOLECULAR_WEIGHT              0.6         # mean molecular weight -> currently only for post-processing [0.6]

example/test_problem/Hydro/AGORA_IsolatedGalaxy/Input_Options/Input__Parameter.low-res

@@ -157,6 +157,17 @@ SF_CREATE_STAR_MIN_STAR_MASS  3.2e4       # minimum star particle mass for the s
 SF_CREATE_STAR_MAX_STAR_MFRAC 0.5         # maximum gas mass fraction allowed to convert to stars per substep [0.5]
 
 
+# feedback (FEEDBACK only)
+FB_MIN_LEVEL                 -1           # minimum AMR level allowed to apply feedback (<0=auto -> MAX_LEVEL) [0]
+FB_RSEED                      456         # random seed [456]
+FB_SNE                        0           # supernova explosion feedback [0]
+FB_RESOLVED_SNEII             1           # Type II supernovae feedback for high enough grid and particle resolution (must add one extra PAR_NATT_FLT_USER)
+FB_USER                       0           # user-defined feedback [0]
+FB_RESOLVED_SNEII_N_PER_MASS  1.0e-5      # number of SNeII per stellar mass (in 1/Msun) (for FB_RESOLVED_SNEII only) [0.01]
+FB_RESOLVED_SNEII_MIN_M_GAS   0.0         # minimum mass of gas in the environment of SNeII to apply feedback (in Msun) (for FB_RESOLVED_SNEII only) [0.0]
+FB_RESOLVED_SNEII_RECORD      0           # record the SNeII explosion info (for FB_RESOLVED_SNEII only) [0]
+
+
 # fluid solver in HYDRO (MODEL==HYDRO only)
 GAMMA                         1.666667    # ratio of specific heats (i.e., adiabatic index) [5.0/3.0]
 MOLECULAR_WEIGHT              0.6         # mean molecular weight -> currently only for post-processing [0.6]

example/test_problem/Hydro/AGORA_IsolatedGalaxy/README

@@ -42,6 +42,10 @@ Note:
 3. The low-resolution initial condition files from "sh download_ic_low_res.sh" are the same initial conditions used
    in the "AgoraGalaxy" test problem of Enzo
 4. Some handy yt analysis scripts are put at "yt_script"
+5. Supernovae feedback (with compilation option `FEEDBACK` enabled) is not officially supported for now
+   - The default star particle masses are too large to use `FB_RESOLVED_SNEII` physically
+   - Vaslues of `FB_RESOLVED_SNEII_N_PER_MASS` and `FB_RESOLVED_SNEII_DELAY_TIME` are unphysical
+     and are only for testing the implementation of SN event rate (with `plot_SN_feedback_rate.py`)
 
 
 First-time GRACKLE installation guide (on NTU clusters as an example):

example/test_problem/Hydro/AGORA_IsolatedGalaxy/clean.sh

@@ -3,4 +3,4 @@ rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Du
       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 \
-      GRACKLE_INFO
+      GRACKLE_INFO Record__FB_Resolved_SNeII

example/test_problem/Hydro/AGORA_IsolatedGalaxy/yt_script/plot_SN_feedback_rate.py

@@ -0,0 +1,92 @@
+# ref: https://yt-project.org/docs/dev/cookbook/calculating_information.html#using-particle-filters-to-calculate-star-formation-rates
+
+import yt
+import numpy as np
+from yt.data_objects.particle_filters import add_particle_filter
+import matplotlib
+matplotlib.use('Agg')
+from matplotlib import pyplot as plt
+
+
+filein  = "../Data_000050"
+fileout = "fig__SN_feedback_rate"
+nbin    = 50
+dpi     = 150
+
+
+# load data
+ds = yt.load( filein )
+
+
+# define the particle filter for the newly formed stars
+def new_star( pfilter, data ):
+   filter = data[ "all", "ParCreTime" ] > 0
+   return filter
+
+add_particle_filter( "new_star", function=new_star, filtered_type="all", requires=["ParCreTime"] )
+ds.add_particle_filter( "new_star" )
+
+
+# define the particle filter for the exploded SNe
+def exploded_SNe( pfilter, data ):
+   filter = data[ "all", "ParSNIITime" ] <= 0 # ParSNIITime is set as negative value of explosion time for exploded particle
+   return filter
+
+add_particle_filter( "exploded_SNe", function=exploded_SNe, filtered_type="all", requires=["ParSNIITime"] )
+ds.add_particle_filter( "exploded_SNe" )
+
+
+# get the creation time of the new stars and explosion time of teh SNe
+ad            = ds.all_data()
+
+star_mass     = ad[ "new_star", "ParMass" ].in_units( "Msun" )
+creation_time = ad[ "new_star", "ParCreTime" ].in_units( "Myr" )
+
+SNe_mass      = ad[ "exploded_SNe", "ParMass" ].in_units( "Msun" )
+SNe_expl_time = -1.0*( ad[ "exploded_SNe", "ParSNIITime" ]*ds.units.code_time ).in_units( "Myr" ) # ParSNIITime is set as negative value of explosion time for exploded particle
+
+print( 'Total number of stars = %d'%len(creation_time) )
+print( 'Total number of SNe   = %d'%len(SNe_expl_time) )
+
+
+# bin the data
+t_start        = 0.0
+t_end          = ds.current_time.in_units( "Myr" )
+t_bin          = np.linspace( start=t_start, stop=t_end, num=nbin+1 )
+time           = 0.5*( t_bin[:-1] + t_bin[1:] )
+star_upper_idx = np.digitize( creation_time, bins=t_bin, right=True )
+SNe_upper_idx  = np.digitize( SNe_expl_time, bins=t_bin, right=True )
+
+
+assert np.all( star_upper_idx > 0 ) and np.all( star_upper_idx < len(t_bin) ), "incorrect star_upper_idx !!"
+assert np.all( SNe_upper_idx > 0 )  and np.all( SNe_upper_idx < len(t_bin) ),  "incorrect SNe_upper_idx !!"
+
+
+# calculate the star formation and SNe mass rate
+Myr2yr = 1.0e6
+sfr    = np.array(  [ star_mass[star_upper_idx == j+1].sum() / ( (t_bin[j+1] - t_bin[j])*Myr2yr ) for j in range(len(time)) ]  )
+sfr[sfr == 0] = np.nan
+
+SNr     = np.array(  [ SNe_mass[SNe_upper_idx == j+1].sum()  / ( (t_bin[j+1] - t_bin[j])*Myr2yr ) for j in range(len(time)) ]  )
+SNr[SNr == 0] = np.nan
+
+
+# calulate the conversion factor
+StarsPerSN  = 1.0/(ds.parameters['FB_ResolvedSNeII_NPerMass']*ds.parameters['SF_CreateStar_MinStarMass'])
+SNDelayTime = ds.quan( ds.parameters['FB_ResolvedSNeII_DelayTime'], 'code_time' ).in_units('Myr').d
+
+
+# plot
+plt.plot( time,               sfr,                  label='Stars' )
+plt.plot( time,               SNr,                  label='SNe'   )
+plt.plot( time,               SNr*StarsPerSN, '--', label=r'SNe, $\times$ %.2f'%(StarsPerSN) )
+plt.plot( time.d-SNDelayTime, SNr*StarsPerSN, '--', label=r'SNe, $\times$ %.2f, shifted %.1f Myr'%(StarsPerSN, SNDelayTime) )
+plt.ylim( 0.0, 1.0e1 )
+plt.legend()
+plt.xlabel( "$\mathrm{t\ [Myr]}$",        fontsize="large" )
+plt.ylabel( "$\mathrm{[M_\odot yr^{-1}]}$", fontsize="large" )
+
+
+# show/save figure
+plt.savefig( fileout+".png", bbox_inches="tight", pad_inches=0.05, dpi=dpi )
+#plt.show()

example/test_problem/Hydro/Plummer/Input__Parameter

@@ -114,7 +114,7 @@ OPT__MINIMIZE_MPI_BARRIER     0           # minimize MPI barriers to improve loa
 
 
 # feedback (FEEDBACK only)
-FB_LEVEL                     -1           # AMR level to apply feedback (must be MAX_LEVEL for now; <0=auto -> MAX_LEVEL) [-1]
+FB_MIN_LEVEL                 -1           # minimum AMR level allowed to apply feedback (<0=auto -> MAX_LEVEL) [0]
 FB_RSEED                      456         # random seed [456]
 FB_USER                       0           # user-defined feedback [0]
 

example/test_problem/Hydro/SNFeedbackBlastWave/FeedbackYieldTable_Resolved_SNeII_N000001

@@ -0,0 +1,2 @@
+# explosion_time  progenitor_mass     ejected_mass   ejected_metals   ejected_energy
+      1.0000e+01       2.0000e+01       1.5000e+01       3.0000e+00       1.0000e+51

example/test_problem/Hydro/SNFeedbackBlastWave/ForBlastWave/Input__Flag_Lohner

@@ -0,0 +1,4 @@
+# Level    Threshold_Refine  Threshold_Derefine              Filter              Soften          MinDensity
+      0                0.80                0.80                0.01                0.00                0.00
+      1                0.80                0.80                0.01                0.00                0.00
+      2                0.80                0.80                0.01                0.00                0.00

example/test_problem/Hydro/SNFeedbackBlastWave/ForBlastWave/Input__Flag_User

@@ -0,0 +1,4 @@
+# Level           User-defined_criteria
+      0                          4.0e-2
+      1                          2.0e-2
+      2                          1.0e-2