Add AGORA subgrid model

README

@@ -66,6 +66,10 @@ Parameters:
                                       GEAR model. This is equivalent to --hydro
                                       --limiter --sync --self-gravity --stars
                                       --star-formation --cooling --feedback.
+    --agora                           Run with all the options needed for the
+                                      AGORA model. This is equivalent to --hydro
+                                      --limiter --sync --self-gravity --stars
+                                      --star-formation --cooling --feedback.
                                       
   Control options:
 

README.md

@@ -160,6 +160,10 @@ Parameters:
                                       GEAR model. This is equivalent to --hydro
                                       --limiter --sync --self-gravity --stars
                                       --star-formation --cooling --feedback.
+    --agora                           Run with all the options needed for the
+                                      GEAR model. This is equivalent to --hydro
+                                      --limiter --sync --self-gravity --stars
+                                      --star-formation --cooling --feedback.
                                       
   Control options:
 

configure.ac

@@ -1783,7 +1783,7 @@ AC_SUBST([SUNDIALS_INCS])
 # As an example for this, see the call to AC_ARG_WITH for cooling.
 AC_ARG_WITH([subgrid],
 	[AS_HELP_STRING([--with-subgrid=<subgrid>],
-		[Master switch for subgrid methods. Inexperienced user should start here. Options are: @<:@none, GEAR, QLA, QLA-EAGLE, EAGLE, EAGLE-XL, SPIN_JET_EAGLE default: none@:>@]
+		[Master switch for subgrid methods. Inexperienced user should start here. Options are: @<:@none, GEAR, AGORA, QLA, QLA-EAGLE, EAGLE, EAGLE-XL, SPIN_JET_EAGLE default: none@:>@]
 	)],
 	[with_subgrid="$withval"],
 	[with_subgrid=none]
@@ -1819,6 +1819,18 @@ case "$with_subgrid" in
 	with_subgrid_extra_io=none
 	enable_fof=no
    ;;
+   AGORA)
+	with_subgrid_cooling=grackle_0
+	with_subgrid_chemistry=AGORA
+	with_subgrid_pressure_floor=GEAR
+	with_subgrid_stars=GEAR
+	with_subgrid_star_formation=GEAR
+	with_subgrid_feedback=AGORA
+	with_subgrid_black_holes=none
+	with_subgrid_sink=none
+	with_subgrid_extra_io=none
+	enable_fof=no
+   ;;   
    QLA)
 	with_subgrid_cooling=QLA
 	with_subgrid_chemistry=QLA
@@ -2240,7 +2252,7 @@ fi
 #  chemistry function
 AC_ARG_WITH([chemistry],
    [AS_HELP_STRING([--with-chemistry=<function>],
-      [chemistry function @<:@none, GEAR_*, QLA, EAGLE default: none@:>@
+      [chemistry function @<:@none, GEAR_*, AGORA, QLA, EAGLE default: none@:>@
       For GEAR, you need to provide the number of elements (e.g. GEAR_10)]
    )],
    [with_chemistry="$withval"],
@@ -2267,6 +2279,10 @@ case "$with_chemistry" in
       with_chemistry_name="GEAR (with $number_element elements)"
       with_chemistry="GEAR"
    ;;
+   AGORA)
+      AC_DEFINE([CHEMISTRY_AGORA], [1], [Chemistry taken from the AGORA model])
+      with_chemistry_name="AGORA"
+   ;;   
    QLA)
       AC_DEFINE([CHEMISTRY_QLA], [1], [Chemistry taken from the Quick-Lyman-alpha model])
       with_chemistry_name="QLA"
@@ -2463,7 +2479,7 @@ esac
 # Feedback model
 AC_ARG_WITH([feedback],
    [AS_HELP_STRING([--with-feedback=<model>],
-      [Feedback model to use @<:@none, EAGLE, EAGLE-thermal, EAGLE-kinetic, GEAR default: none@:>@]
+      [Feedback model to use @<:@none, EAGLE, EAGLE-thermal, EAGLE-kinetic, GEAR, AGORA default: none@:>@]
    )],
    [with_feedback="$withval"],
    [with_feedback="none"]
@@ -2495,6 +2511,10 @@ case "$with_feedback" in
       AC_DEFINE([FEEDBACK_GEAR], [1], [GEAR stellar feedback and evolution model])
       with_feedback_name="GEAR"
    ;;
+   AGORA)
+      AC_DEFINE([FEEDBACK_AGORA], [1], [AGORA stellar feedback and evolution model])
+      with_feedback_name="AGORA"
+   ;;   
    none)
       AC_DEFINE([FEEDBACK_NONE], [1], [No feedback])
    ;;
@@ -2882,6 +2902,9 @@ AM_CONDITIONAL([HAVEGRACKLECOOLING], [test "$with_cooling" = "grackle"])
 # check if using gear feedback
 AM_CONDITIONAL([HAVEGEARFEEDBACK], [test "$with_feedback" = "GEAR"])
 
+# check if using gear feedback
+AM_CONDITIONAL([HAVEAGORAFEEDBACK], [test "$with_feedback" = "AGORA"])
+
 # check if using SPHENIX
 AM_CONDITIONAL([HAVE_SPHENIX], [test "$with_hydro" = "sphenix"])
 
@@ -2894,6 +2917,9 @@ AM_CONDITIONAL([HAVE_CHEMISTRY_NONE], [test "$with_chemistry" = "none"])
 # check if using GEAR chemistry
 AM_CONDITIONAL([HAVE_CHEMISTRY_GEAR], [test "$with_chemistry" = "GEAR" || test "$with_chemistry" = "GEAR_DIFFUSION"])
 
+# check if using AGORA chemistry
+AM_CONDITIONAL([HAVE_CHEMISTRY_AGORA], [test "$with_chemistry" = "AGORA" || test "$with_chemistry" = "GEAR_DIFFUSION"])
+
 # check if using default stars
 AM_CONDITIONAL([HAVE_STARS_BASIC], [test "$with_stars" = "basic"])
 

doc/RTD/source/CommandLineOptions/index.rst

@@ -63,6 +63,10 @@ can be found by typing ``./swift -h``:
                                       GEAR model. This is equivalent to --hydro
                                       --limiter --sync --self-gravity --stars
                                       --star-formation --cooling --feedback.
+    --agora                           Run with all the options needed for the
+                                      AGORA model. This is equivalent to --hydro
+                                      --limiter --sync --self-gravity --stars
+                                      --star-formation --cooling --feedback.
 
   Control options:
 

doc/RTD/source/SubgridModels/AGORA/index.rst

+.. AGORA sub-grid model
+   Yves Revaz 04.2023
+
+.. _agora:
+
+AGORA model
+===========
+
+AGORA's model corresponds to the subgrid model adopted by the 
+`AGORA Project <https://sites.google.com/site/santacruzcomparisonproject/>`_ (High-resolution Galaxy Simulations Comparison Project).
+Details of the model is described in 
+`Kim \& al. 2014 <https://ui.adsabs.harvard.edu/link_gateway/2014ApJS..210...14K/PUB_PDF>`_
+and
+`Kim \& al. 2016 <https://ui.adsabs.harvard.edu/link_gateway/2016ApJ...833..202K/PUB_PDF>`_.
+
+
+
+This model can be selected with the configuration option ``--with-subgrid=AGORA`` and run with the option ``--agora``. 
+A few examples exist and can be found in ``examples/AGORA``. 
+
+.. _agora_grackle_cooling:
+
+Gas cooling/heating: Grackle
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The AGORA model uses the Grackle cooling library which is also used by the :ref:`GEAR model <gear_grackle_cooling>`.
+
+
+
+
+.. _agora_star_formation:
+
+Star formation
+~~~~~~~~~~~~~~
+
+The AGORA model uses the :ref:`GEAR model <gear_star_formation>` scheme, however with the 
+``GEARStarFormation:star_formation_mode`` parameter set to ``agora``. Instead of requiring the gas
+density to reach the pressure floor, we simply require it to be denser than a density
+threshold defined by ``GEARStarFormation:density_threshold``.
+
+
+Recommended parameters for the AGORA model should be:
+
+
+.. code:: YAML
+
+  GEARStarFormation:
+    star_formation_mode: agora            
+    star_formation_efficiency: 0.01   
+    maximal_temperature:  1e10       
+    n_stars_per_particle: 1
+    min_mass_frac: 0.5
+    density_threshold:   1.67e-23   
+
+
+
+.. _agora_feedback_and_chemistry:
+
+Stellar Feedback and Chemistry
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+In the current implementation, only two `elements`, iron (Fe) and the sum of all elements heavier than helium (metallicity)
+are considered.
+Only stellar feeback from core collapse supernovae (CCSNe) is considered.
+The model assumes that a number ``AGORAFeedback:ccsne_per_solar_mass`` of CCSNe per solar mass will form out of 
+each stellar particle. Those supernovae will all explode after a time 
+``AGORAFeedback:supernovae_explosion_time_myr`` after the birth of the stellar particle.
+At this time, each stellar particle will expel:
+
+  - ``AGORAFeedback:ejected_mass_in_solar_mass_per_CCSN`` amount of gas (in solar mass), per CCSN formed
+  - ``AGORAFeedback:ejected_Fe_mass_in_solar_mass_per_CCSN`` amount of iron (in solar mass), per CCSN formed
+  - ``AGORAFeedback:ejected_metal_mass_in_solar_mass_per_CCSN`` amount of metals (in solar mass), per CCSN formed
+
+In addition stellar particles will release energy:
+
+  - ``AGORAFeedback:energy_in_erg_per_CCSN`` erg per CCSN formed
+
+The energy released effectively into surrounding gas particles can be mitigated with the 
+parameter ``AGORAFeedback:supernovae_efficiency``, used as a simple factor to ``AGORAFeedback:energy_in_erg_per_CCSN``.
+
+Both energy and mass are ejected into the surrounding gas according to the SPH kernel.
+
+The mass fraction of elements received by gas particles will be smoothed (smoothed metallicity scheme)
+by using the SPH kernel. Note that this scheme does not exchange any material between particles. 
+Snapshots can store both the smoothed metallicity (``SmoothedMetalMassFractions``)
+and/or the non-smoothed one (``MetalMassFractions``), i.e., the mass fraction of elements effectively received
+by the gas particles.
+
+
+The initial metallicity of the gas can be defined by the parameter ``AGORAChemistry:initial_metallicity``.
+A value less than 0 forces the code to take the gas metallicity from the initial condition file (snapshot).
+Instead, if ``AGORAChemistry:scale_initial_metallicity`` is different than 0, the initial mass fraction
+of elements will be set to:
+
+  - ``AGORAChemistry:initial_metallicity`` time ``AGORAChemistry:solar_abundance_Metals`` for the iron
+  - ``AGORAChemistry:initial_metallicity`` time ``AGORAChemistry:solar_abundance_Metals`` for the metals
+
+
+Instead, they will be set to:
+
+  - ``AGORAChemistry:initial_metallicity`` for the iron
+  - ``AGORAChemistry:initial_metallicity`` for the metals
+
+
+
+Recommended parameters for the AGORA model should be:
+
+.. code:: YAML
+
+   AGORAChemistry:
+    initial_metallicity: 1  
+    scale_initial_metallicity: 1            
+    solar_abundance_Fe: 0.001771
+    solar_abundance_Metals: 0.02  
+    
+
+.. code:: YAML
+   
+   AGORAFeedback:
+    energy_in_erg_per_CCSN: 1e51
+    supernovae_efficiency: 1
+    supernovae_explosion_time_myr: 5
+    ccsne_per_solar_mass : 0.010989
+    ejected_mass_in_solar_mass_per_CCSN : 14.8
+    ejected_Fe_mass_in_solar_mass_per_CCSN : 2.63
+    ejected_metal_mass_in_solar_mass_per_CCSN : 2.63
+
+
+
+
+
+
+.. _agora_pressure_floor:
+
+Pressure Floor
+~~~~~~~~~~~~~~
+
+The AGORA model uses precisely the same pressure floor than the :ref:`GEAR model <gear_pressure_floor>`.
+
+
+

doc/RTD/source/SubgridModels/GEAR/index.rst

@@ -8,6 +8,8 @@ GEAR model
 GEAR's model are mainly described in `Revaz \& Jablonka <https://ui.adsabs.harvard.edu/abs/2018A%26A...616A..96R/abstract>`_.
 This model can be selected with the configuration option ``--with-subgrid=GEAR`` and run with the option ``--gear``. A few examples exist and can be found in ``examples/GEAR``. 
 
+.. _gear_pressure_floor:
+
 Pressure Floor
 ~~~~~~~~~~~~~~
 
@@ -93,7 +95,7 @@ The self shielding method is defined by ``GrackleCooling:self_shielding_method``
     self_shielding_method: -1                    # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
     self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
 
-
+.. _gear_star_formation:
 
 Star formation
 ~~~~~~~~~~~~~~
@@ -103,8 +105,11 @@ The star formation is done in two steps: first we check if a particle is in the
 A particle is in the star forming regime if:
  - The velocity divergence is negative (:math:`\nabla\cdot v < 0`),
  - The temperature is lower than a threshold (:math:`T < T_t` where :math:`T_t` is defined with ``GEARStarFormation:maximal_temperature``),
+ - The gas density is higher than a threshold (:math:`\rho > \rho_t` where :math:`\rho_t` is defined with ``GEARStarFormation:density_threshold``)
  - The particle reaches the pressure floor (:math:`\rho > \frac{\pi}{4 G N_\textrm{Jeans}^{2/3} h^2}\frac{\gamma k_B T}{\mu m_p}` where :math:`N_\textrm{Jeans}` is defined in the pressure floor).
 
+If ``GEARStarFormation:star_formation_mode`` is set to ``agora``, the condition on the pressure floor is ignored. Its default value is ``default``.
+
 A star will be able to form if a randomly drawn number is below :math:`\frac{m_g}{m_\star}\left(1 - \exp\left(-c_\star \Delta t / t_\textrm{ff}\right)\right)` where :math:`t_\textrm{ff}` is the free fall time, :math:`\Delta t` is the time step of the particle and :math:`c_\star` is the star formation coefficient (``GEARStarFormation:star_formation_efficiency``), :math:`m_g` the mass of the gas particle and :math:`m_\star` the mass of the possible future star. The mass of the star is computed from the average gas mass in the initial conditions divided by the number of possible stars formed per gas particle (``GEARStarFormation:n_stars_per_particle``). When we cannot have enough mass to form a second star (defined with the fraction of mass ``GEARStarFormation:min_mass_frac``), we fully convert the gas particle into a stellar particle. Once the star is formed, we move it a bit in a random direction and fraction of the smoothing length in order to avoid any division by 0.
 
 Currently, only the following hydro schemes are compatible: SPHENIX and Gadget2.

doc/RTD/source/SubgridModels/index.rst

@@ -18,3 +18,4 @@ be use as an empty canvas to be copied to create additional models.
    GEAR/index
    Basic/index	     
    AGNSpinJets/index
+   AGORA/index

examples/AGORA/AgoraDisk/README

+To run this example, SWIFT must be configured with the following options:
+
+./configure --with-subgrid=AGORA
+
+or similarily 
+
+./configure --with-feedback=AGORA --with-chemistry=AGORA  --with-cooling=grackle_0  --with-pressure-floor=GEAR --with-stars=GEAR --with-star-formation=GEAR 
+
+
+An analysis script for this test case can be found in the AGORA project repository:
+https://bitbucket.org/mornkr/agora-analysis-script/
+
+A modified version of the script able to read in SWIFT data is available upon request.
+
+(Note to the SWIFT authors: The modified script is stored with the other reference solutions.)

examples/AGORA/AgoraDisk/agora_disk.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e21 # kpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
+
+Scheduler:
+  max_top_level_cells: 16
+  cell_extra_sparts:  5000       # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
+  cell_extra_gparts:  5000       # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
+  cell_max_size:             8000   # (Optional) Maximal number of interactions per task if we force the split (this is the default value).
+  cell_sub_size_pair_hydro:  2560 # (Optional) Maximal number of hydro-hydro interactions per sub-pair hydro/star task (this is the default value).
+  cell_sub_size_self_hydro:  3200     # (Optional) Maximal number of hydro-hydro interactions per sub-self hydro/star task (this is the default value).
+  cell_sub_size_pair_stars:  2560 # (Optional) Maximal number of hydro-star interactions per sub-pair hydro/star task (this is the default value).
+  cell_sub_size_self_stars:  3200     # (Optional) Maximal number of hydro-star interactions per sub-self hydro/star task (this is the default value).
+  cell_sub_size_pair_grav:   2560 # (Optional) Maximal number of interactions per sub-pair gravity task  (this is the default value).
+  cell_sub_size_self_grav:   3200     # (Optional) Maximal number of interactions per sub-self gravity task  (this is the default value).
+  cell_split_size:           100       # (Optional) Maximal number of particles per cell (this is the default value).
+  
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   0.25  # The end time of the simulation (in internal units).
+  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
+  dt_max:     0.1  # The maximal time-step size of the simulation (in internal units).
+  max_dt_RMS_factor: 0.25  # (Optional) Dimensionless factor for the maximal displacement allowed based on the RMS velocities.
+
+# Parameters governing the snapshots
+Snapshots:
+  basename:            snapshot # Common part of the name of output files
+  time_first:          0.    # Time of the first output (in internal units)
+  delta_time:          1e-3  # Time difference between consecutive outputs (in internal units)
+  compression:         4
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  delta_time:          1e-3 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                       0.05    # Constant dimensionless multiplier for time integration.
+  MAC:                       geometric
+  theta_cr:                  0.7     
+  use_tree_below_softening:  1
+  comoving_DM_softening:     0.08 # Comoving softening length (in internal units).
+  max_physical_DM_softening: 0.08    # Physical softening length (in internal units).
+  comoving_baryon_softening:     0.08 # Comoving softening length (in internal units).
+  max_physical_baryon_softening: 0.08    # Physical softening length (in internal units).
+  mesh_side_length:       32        # Number of cells along each axis for the periodic gravity mesh.
+  
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   10.      # Kelvin
+  h_min_ratio:           0.1095       # (Optional) Minimal allowed smoothing length in units of the softening. Defaults to 0 if unspecified.
+  h_max:                     10.      # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./agora_disk.hdf5     # The file to read
+  periodic:   0                     # Non-periodic BCs
+  cleanup_h_factors: 1              # Remove the h-factors inherited from Gadget
+  shift:    [674.1175, 674.1175, 674.1175]   # Centre the box
+
+
+# Cooling with Grackle 2.0
+GrackleCooling:
+  cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+  with_UV_background: 1 # Enable or not the UV background
+  redshift: 0 # Redshift to use (-1 means time based redshift)
+  with_metal_cooling: 1 # Enable or not the metal cooling
+  provide_volumetric_heating_rates: 0 # User provide volumetric heating rates
+  provide_specific_heating_rates: 0 # User provide specific heating rates
+  self_shielding_method: -1 # Grackle (<= 3) or Gear self shielding method
+  self_shielding_threshold_atom_per_cm3: 1e10  # Required only with GEAR's self shielding. Density threshold of the self shielding
+  max_steps: 1000
+  convergence_limit: 1e-2
+  thermal_time_myr: 0
+
+
+GEARStarFormation:
+  star_formation_mode: agora        # default or agora
+  star_formation_efficiency: 0.01   # star formation efficiency (c_*)
+  maximal_temperature:  1e10        # Upper limit to the temperature of a star forming particle
+  n_stars_per_particle: 1
+  min_mass_frac: 0.5
+  density_threshold:   1.67e-23     # Density threashold in g/cm3
+
+
+GEARPressureFloor:
+  jeans_factor: 10
+
+AGORAFeedback:
+  energy_in_erg_per_CCSN: 1e51
+  supernovae_efficiency: 1
+  supernovae_explosion_time_myr: 5 # In Myr
+  ccsne_per_solar_mass : 0.010989 # 1/91
+  ejected_mass_in_solar_mass_per_CCSN : 14.8
+  ejected_Fe_mass_in_solar_mass_per_CCSN : 2.63
+  ejected_metal_mass_in_solar_mass_per_CCSN : 2.63
+
+AGORAChemistry:
+  initial_metallicity: 1          # if less than 0, read the metallicity from the snapshot (1 means solar abundance)
+  solar_abundance_Metals: 0.02  
+  scale_initial_metallicity: 1    # scale the initial metallicity using solar_abundance_Metals ?
+
+
+
+
+
+Restarts:
+  delta_hours:   1.        # (Optional) decimal hours between dumps of restart files.
+

examples/AGORA/AgoraDisk/checkSolution.py

+import h5py
+import numpy as np
+
+expected_total_mass = 227413100.0
+
+f = h5py.File("snapshot_0025.hdf5", "r")
+
+# get the total stellar mass
+mass = f["StarsParticles/Masses"]
+total_mass = sum(mass) * 1e10
+
+error = np.fabs((expected_total_mass - total_mass) / expected_total_mass)
+
+if error < 0.01:
+    print(48 * "!")
+    print(r"Well done !")
+    print(r"The stellar mass is %g Msol," % (total_mass))
+    print(r"The expected stellar mass is %g Msol" % (expected_total_mass))
+    print(r"This represents an error of %d %% !" % (100 * error))
+    print(48 * "!")
+else:
+    print(48 * "!")
+    print(r"Too bad !")
+    print(r"The stellar mass is %g Msol," % (total_mass))
+    print(r"While the expected stellar mass is %g Msol" % (expected_total_mass))
+    print(r"This represents an error of %d %% !" % (100 * error))
+    print(r"This is beyond the requirements.")
+    print(48 * "!")

examples/AGORA/AgoraDisk/getGrackleCoolingTable.sh

+#!/bin/bash
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/CloudyData_UVB=HM2012.h5

examples/AGORA/AgoraDisk/getIC.sh

+#!/bin/bash
+
+if [ "$#" -ne 1 ]; then
+    echo "You need to provide the resolution (e.g. ./getIC.sh LOW)."
+    echo "The possible options are LOW and MED."
+    exit
+fi
+
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/AgoraDisk/snap.$1.hdf5
+mv snap.$1.hdf5 agora_disk.hdf5
+

examples/AGORA/AgoraDisk/run.sh

+#!/bin/bash
+
+# This example is based on the AGORA disk article (DOI: 10.3847/1538-4357/833/2/202)
+
+# set the resolution (LOW or MED)
+sim=LOW
+
+
+# make run.sh fail if a subcommand fails
+set -e
+
+# Generate the initial conditions if they are not present.
+if [ ! -e agora_disk.hdf5 ]
+then
+    echo "Fetching initial glass file for the Sedov blast example..."
+    ./getIC.sh $sim
+fi
+
+# Get the Grackle cooling table
+if [ ! -e CloudyData_UVB=HM2012.h5 ]
+then
+    echo "Fetching the Cloudy tables required by Grackle..."
+    ./getGrackleCoolingTable.sh
+fi
+
+
+
+
+# Run SWIFT
+../../../swift --sync --limiter --cooling --hydro --self-gravity --star-formation --feedback --stars --threads=8 agora_disk.yml 2>&1 | tee output.log
+
+
+
+

examples/GEAR/AgoraCosmo/agora_cosmo.yml

@@ -93,7 +93,7 @@ GEARStarFormation:
   maximal_temperature:  3e4         # Upper limit to the temperature of a star forming particle
   n_stars_per_particle: 1
   min_mass_frac: 0.5
-  density_threashold:   1e-30       # Density threashold (in addition to the pressure floor) in g/cm3
+  density_threshold:   1e-30       # Density threashold (in addition to the pressure floor) in g/cm3
 
 GEARPressureFloor:
   jeans_factor: 8.75

examples/GEAR/AgoraDisk/agora_disk.yml

@@ -90,7 +90,7 @@ GEARStarFormation:
   maximal_temperature:  1e10         # Upper limit to the temperature of a star forming particle
   n_stars_per_particle: 1
   min_mass_frac: 0.5
-  density_threashold:   1.67e-23   # Density threashold in g/cm3
+  density_threshold:   1.67e-23   # Density threashold in g/cm3
 
 
 GEARPressureFloor:

examples/GEAR/ZoomIn/zoom_in.yml

@@ -96,7 +96,7 @@ GrackleCooling:
 GEARStarFormation:
   star_formation_efficiency: 0.01   # star formation efficiency (c_*)
   maximal_temperature:  3e4         # Upper limit to the temperature of a star forming particle
-  density_threashold:   1.67e-25    # Density threashold in g/cm3
+  density_threshold:   1.67e-25    # Density threashold in g/cm3
   n_stars_per_particle: 4
   min_mass_frac: 0.5
 

examples/parameter_example.yml

@@ -526,6 +526,14 @@ GEARChemistry:
   scale_initial_metallicity: 1   # Should we scale the initial metallicity with the solar one?
   diffusion_coefficient: 1e-3    # Coefficient for the diffusion (see Shen et al. 2010; differs by gamma^2 due to h^2).
 
+# AGORA chemistry model (Kim et al. 2014)
+AGORAChemistry:
+ initial_metallicity: 1
+ scale_initial_metallicity: 1
+ solar_abundance_Fe: 0.001771
+ solar_abundance_Metals: 0.02
+
+
 # Parameters related to star formation models  -----------------------------------------------
 
 # EAGLE star formation model (Schaye and Dalla Vecchia 2008)
@@ -555,11 +563,12 @@ QLAStarFormation:
 
 # GEAR star formation model (Revaz and Jablonka 2018)
 GEARStarFormation:
+  star_formation_mode: default      # default (use the pressure floor limit) or agora (do not use the pressure floor limit)
   star_formation_efficiency: 0.01   # star formation efficiency (c_*)
   maximal_temperature:  3e4         # Upper limit to the temperature of a star forming particle
   n_stars_per_particle: 4           # Number of stars that an hydro particle can generate
   min_mass_frac: 0.5                # Minimal mass for a stellar particle as a fraction of the average mass for the stellar particles.
-
+  density_threshold: 1.67e-23       # Density threshold for star formation
 
 # Parameters related to feedback models  -----------------------------------------------
 
@@ -621,6 +630,16 @@ GEARFeedback:
   discrete_yields: 0                                       # Should we use discrete yields or the IMF integrated one?
   elements: [Fe, Mg, O, S, Zn, Sr, Y, Ba, Eu]              # Elements to read in the yields table. The number of element should be one less than the number of elements (N) requested during the configuration (--with-chemistry=GEAR_N).
 
+# AGORA feedback model
+AGORAFeedback:
+  energy_in_erg_per_CCSN: 1e51
+  supernovae_efficiency: 1
+  supernovae_explosion_time_myr: 5
+  ccsne_per_solar_mass : 0.010989
+  ejected_mass_in_solar_mass_per_CCSN : 14.8
+  ejected_Fe_mass_in_solar_mass_per_CCSN : 2.63
+  ejected_metal_mass_in_solar_mass_per_CCSN : 2.63
+ 
 
 # Parameters related to AGN models  -----------------------------------------------
 

src/Makefile.am

@@ -129,6 +129,12 @@ GEAR_FEEDBACK_SOURCES += feedback/GEAR/stellar_evolution.c feedback/GEAR/feedbac
 	feedback/GEAR/initial_mass_function.c feedback/GEAR/supernovae_ia.c feedback/GEAR/supernovae_ii.c
 endif
 
+# source files for AGORA feedback
+AGORA_FEEDBACK_SOURCES =
+if HAVEAGORAFEEDBACK
+AGORA_FEEDBACK_SOURCES += feedback/AGORA/feedback.c
+endif
+
 # source files for GEAR RT
 GEAR_RT_SOURCES =
 if HAVEGEARRT
@@ -185,6 +191,7 @@ AM_SOURCES += $(EAGLE_EXTRA_IO_SOURCES)
 AM_SOURCES += $(QLA_COOLING_SOURCES) $(QLA_EAGLE_COOLING_SOURCES) 
 AM_SOURCES += $(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES) 
 AM_SOURCES += $(GRACKLE_COOLING_SOURCES) $(GEAR_FEEDBACK_SOURCES) 
+AM_SOURCES += $(AGORA_FEEDBACK_SOURCES)
 AM_SOURCES += $(PS2020_COOLING_SOURCES)
 AM_SOURCES += $(SPHM1RT_RT_SOURCES)
 AM_SOURCES += $(GEAR_RT_SOURCES)
@@ -414,6 +421,12 @@ nobase_noinst_HEADERS += chemistry/GEAR_DIFFUSION/chemistry_io.h
 nobase_noinst_HEADERS += chemistry/GEAR_DIFFUSION/chemistry_struct.h
 nobase_noinst_HEADERS += chemistry/GEAR_DIFFUSION/chemistry_iact.h
 nobase_noinst_HEADERS += chemistry/GEAR_DIFFUSION/chemistry_debug.h
+nobase_noinst_HEADERS += chemistry/AGORA/chemistry.h 
+nobase_noinst_HEADERS += chemistry/AGORA/chemistry_io.h 
+nobase_noinst_HEADERS += chemistry/AGORA/chemistry_csds.h
+nobase_noinst_HEADERS += chemistry/AGORA/chemistry_struct.h 
+nobase_noinst_HEADERS += chemistry/AGORA/chemistry_iact.h 
+nobase_noinst_HEADERS += chemistry/AGORA/chemistry_debug.h 
 nobase_noinst_HEADERS += chemistry/EAGLE/chemistry.h 
 nobase_noinst_HEADERS += chemistry/EAGLE/chemistry_io.h 
 nobase_noinst_HEADERS += chemistry/EAGLE/chemistry_struct.h
@@ -435,6 +448,9 @@ nobase_noinst_HEADERS += extra_io/EAGLE/extra_io.h extra_io/EAGLE/extra.h
 nobase_noinst_HEADERS += feedback/none/feedback.h feedback/none/feedback_struct.h feedback/none/feedback_iact.h 
 nobase_noinst_HEADERS += feedback/none/feedback_properties.h feedback/none/feedback.h
 nobase_noinst_HEADERS += feedback/none/feedback_debug.h
+nobase_noinst_HEADERS += feedback/AGORA/feedback.h feedback/AGORA/feedback_struct.h feedback/AGORA/feedback_iact.h 
+nobase_noinst_HEADERS += feedback/AGORA/feedback_properties.h feedback/AGORA/feedback.h
+nobase_noinst_HEADERS += feedback/AGORA/feedback_debug.h
 nobase_noinst_HEADERS += feedback/EAGLE_kinetic/feedback.h feedback/EAGLE_kinetic/feedback_struct.h 
 nobase_noinst_HEADERS += feedback/EAGLE_kinetic/feedback_properties.h feedback/EAGLE_kinetic/feedback_iact.h 
 nobase_noinst_HEADERS += feedback/EAGLE_kinetic/feedback_debug.h

src/chemistry.h

@@ -40,6 +40,9 @@
 #elif defined(CHEMISTRY_GEAR_DIFFUSION)
 #include "./chemistry/GEAR_DIFFUSION/chemistry.h"
 #include "./chemistry/GEAR_DIFFUSION/chemistry_iact.h"
+#elif defined(CHEMISTRY_AGORA)
+#include "./chemistry/AGORA/chemistry.h"
+#include "./chemistry/AGORA/chemistry_iact.h"
 #elif defined(CHEMISTRY_QLA)
 #include "./chemistry/QLA/chemistry.h"
 #include "./chemistry/QLA/chemistry_iact.h"

src/chemistry/AGORA/chemistry.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2023 Yves Revaz (yves.revaz@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_CHEMISTRY_AGORA_H
+#define SWIFT_CHEMISTRY_AGORA_H
+
+/**
+ * @file src/chemistry/none/chemistry.h
+ * @brief Empty infrastructure for the cases without chemistry function
+ */
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+#include <string.h>
+
+/* Local includes. */
+#include "chemistry_struct.h"
+#include "cosmology.h"
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/**
+ * @brief Get the name of the element i.
+ *
+ * @param sm The #stellar_model.
+ * @param i The element indice.
+ */
+static INLINE const char* chemistry_get_element_name(
+    const struct chemistry_global_data* data, int i) {
+
+  return data->elements_name + i * AGORA_LABELS_SIZE;
+}
+
+/**
+ * @brief Copies the chemistry properties of the gas particle over to the
+ * star particle.
+ *
+ * @param p the gas particles.
+ * @param xp the additional properties of the gas particles.
+ * @param sp the new created star particle with its properties.
+ */
+INLINE static void chemistry_copy_star_formation_properties(
+    struct part* p, const struct xpart* xp, struct spart* sp) {
+
+  float mass = hydro_get_mass(p);
+
+  /* Store the chemistry struct in the star particle */
+  for (int i = 0; i < AGORA_CHEMISTRY_ELEMENT_COUNT; i++) {
+    sp->chemistry_data.metal_mass_fraction[i] =
+        p->chemistry_data.smoothed_metal_mass_fraction[i];
+
+    /* Remove the metals taken by the star. */
+    p->chemistry_data.metal_mass[i] *= mass / (mass + sp->mass);
+  }
+}
+
+/**
+ * @brief Initialises the chemistry properties.
+ *
+ * Nothing to do here.
+ *
+ * @param parameter_file The parsed parameter file.
+ * @param us The current internal system of units.
+ * @param phys_const The physical constants in internal units.
+ * @param data The global chemistry information (to be filled).
+ */
+static INLINE void chemistry_init_backend(struct swift_params* parameter_file,
+                                          const struct unit_system* us,
+                                          const struct phys_const* phys_const,
+                                          struct chemistry_global_data* data) {
+
+  /* read parameters */
+  const float initial_metallicity = parser_get_param_float(
+      parameter_file, "AGORAChemistry:initial_metallicity");
+
+  if (initial_metallicity < 0) {
+    message("Setting the initial metallicity from the snapshot.");
+  } else {
+    message("Setting the initial metallicity from the parameter file.");
+  }
+
+  /* Set the initial metallicities */
+  for (int i = 0; i < AGORA_CHEMISTRY_ELEMENT_COUNT; i++) {
+    data->initial_metallicities[i] = initial_metallicity;
+  }
+
+  /* Check if need to scale the initial metallicity */
+  const int scale_metallicity = parser_get_opt_param_int(
+      parameter_file, "AGORAChemistry:scale_initial_metallicity", 0);
+
+  /* Scale the metallicities if required */
+  if (scale_metallicity) {
+
+    /* Set element name (Metals goes to the end). */
+    strcpy(data->elements_name + (0) * AGORA_LABELS_SIZE, "Fe");
+    strcpy(data->elements_name +
+               (AGORA_CHEMISTRY_ELEMENT_COUNT - 1) * AGORA_LABELS_SIZE,
+           "Metals");
+
+    /* Read the solar abundances. */
+    data->solar_abundances[0] = parser_get_opt_param_double(
+        parameter_file, "AGORAChemistry:solar_abundance_Fe", 0.001771);
+    data->solar_abundances[AGORA_CHEMISTRY_ELEMENT_COUNT - 1] =
+        parser_get_opt_param_double(
+            parameter_file, "AGORAChemistry:solar_abundance_Metals", 0.02);
+
+    /* Scale the metallicities */
+    for (int i = 0; i < AGORA_CHEMISTRY_ELEMENT_COUNT; i++) {
+      data->initial_metallicities[i] *= data->solar_abundances[i];
+    }
+  }
+}
+
+/**
+ * @brief Prepares a particle for the smooth metal calculation.
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the various smooth metallicity tasks
+ *
+ * @param p The particle to act upon
+ * @param cd #chemistry_global_data containing chemistry informations.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_init_part(
+    struct part* restrict p, const struct chemistry_global_data* cd) {
+
+  struct chemistry_part_data* cpd = &p->chemistry_data;
+
+  for (int i = 0; i < AGORA_CHEMISTRY_ELEMENT_COUNT; i++) {
+    /* Reset the smoothed metallicity */
+    cpd->smoothed_metal_mass_fraction[i] = 0.f;
+  }
+}
+
+/**
+ * @brief Prints the properties of the chemistry model to stdout.
+ *
+ * @brief The #chemistry_global_data containing information about the current
+ * model.
+ */
+static INLINE void chemistry_print_backend(
+    const struct chemistry_global_data* data) {
+
+  message("Chemistry function is 'AGORA'.");
+}
+
+/**
+ * @brief Finishes the density calculation.
+ *
+ * @param p The particle to act upon
+ * @param cd The global chemistry information.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_density(
+    struct part* restrict p, const struct chemistry_global_data* cd,
+    const struct cosmology* cosmo) {
+
+  /* Some smoothing length multiples. */
+  const float h = p->h;
+  const float h_inv = 1.0f / h;                       /* 1/h */
+  const float factor = pow_dimension(h_inv) / p->rho; /* 1 / h^d * rho */
+
+  struct chemistry_part_data* cpd = &p->chemistry_data;
+
+  for (int i = 0; i < AGORA_CHEMISTRY_ELEMENT_COUNT; i++) {
+    /* Final operation on the density (add self-contribution). */
+    cpd->smoothed_metal_mass_fraction[i] += cpd->metal_mass[i] * kernel_root;
+
+    /* Finish the calculation by inserting the missing h-factors */
+    cpd->smoothed_metal_mass_fraction[i] *= factor;
+  }
+}
+
+/**
+ * @brief Updates to the chemistry data after the hydro force loop.
+ *
+ * Nothing to do here.
+ *
+ * @param p The particle to act upon.
+ * @param cosmo The current cosmological model.
+ * @param with_cosmology Are we running with the cosmology?
+ * @param time Current time of the simulation.
+ * @param dt Time step (in physical units).
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_force(
+    struct part* restrict p, const struct cosmology* cosmo,
+    const int with_cosmology, const double time, const double dt) {}
+
+/**
+ * @brief Computes the chemistry-related time-step constraint.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param cosmo The current cosmological model.
+ * @param us The internal system of units.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cd The global properties of the chemistry scheme.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float chemistry_timestep(
+    const struct phys_const* restrict phys_const,
+    const struct cosmology* restrict cosmo,
+    const struct unit_system* restrict us,
+    const struct hydro_props* hydro_props,
+    const struct chemistry_global_data* cd, const struct part* restrict p) {
+  return FLT_MAX;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param cd #chemistry_global_data containing chemistry informations.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void
+chemistry_part_has_no_neighbours(struct part* restrict p,
+                                 struct xpart* restrict xp,
+                                 const struct chemistry_global_data* cd,
+                                 const struct cosmology* cosmo) {
+
+  /* Set the smoothed fractions with the non smoothed fractions */
+  for (int i = 0; i < AGORA_CHEMISTRY_ELEMENT_COUNT; i++) {
+    p->chemistry_data.smoothed_metal_mass_fraction[i] =
+        p->chemistry_data.metal_mass[i] / hydro_get_mass(p);
+  }
+}
+
+/**
+ * @brief Sets the chemistry properties of the (x-)particles to a valid start
+ * state.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param cosmo The current cosmological model.
+ * @param data The global chemistry information used for this run.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_part(
+    const struct phys_const* restrict phys_const,
+    const struct unit_system* restrict us,
+    const struct cosmology* restrict cosmo,
+    const struct chemistry_global_data* data, struct part* restrict p,
+    struct xpart* restrict xp) {
+
+  for (int i = 0; i < AGORA_CHEMISTRY_ELEMENT_COUNT; i++) {
+    if (data->initial_metallicities[i] < 0) {
+      /* Use the value from the IC. We are reading the metal mass fraction. */
+      p->chemistry_data.metal_mass[i] *= hydro_get_mass(p);
+    } else {
+      /* Use the value from the parameter file */
+      p->chemistry_data.metal_mass[i] =
+          data->initial_metallicities[i] * hydro_get_mass(p);
+    }
+  }
+
+  chemistry_init_part(p, data);
+}
+
+/**
+ * @brief Sets the chemistry properties of the sparticles to a valid start
+ * state.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cosmo The current cosmological model.
+ * @param data The global chemistry information.
+ * @param sp Pointer to the sparticle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
+    const struct chemistry_global_data* data, struct spart* restrict sp) {
+
+  for (int i = 0; i < AGORA_CHEMISTRY_ELEMENT_COUNT; i++) {
+    sp->chemistry_data.metal_mass_fraction[i] = data->initial_metallicities[i];
+  }
+}
+
+/**
+ * @brief Initialise the chemistry properties of a black hole with
+ * the chemistry properties of the gas it is born from.
+ *
+ * Nothing to do here.
+ *
+ * @param bp_data The black hole data to initialise.
+ * @param p_data The gas data to use.
+ * @param gas_mass The mass of the gas particle.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_bpart_from_part(
+    struct chemistry_bpart_data* bp_data,
+    const struct chemistry_part_data* p_data, const double gas_mass) {
+  error("To be implemented.");
+}
+
+/**
+ * @brief Add the chemistry data of a gas particle to a black hole.
+ *
+ * Nothing to do here.
+ *
+ * @param bp_data The black hole data to add to.
+ * @param p_data The gas data to use.
+ * @param gas_mass The mass of the gas particle.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_add_part_to_bpart(
+    struct chemistry_bpart_data* bp_data,
+    const struct chemistry_part_data* p_data, const double gas_mass) {
+  error("To be implemented.");
+}
+
+/**
+ * @brief Transfer chemistry data of a gas particle to a black hole.
+ *
+ * Nothing to do here.
+ *
+ * @param bp_data The black hole data to add to.
+ * @param p_data The gas data to use.
+ * @param nibble_mass The mass to be removed from the gas particle.
+ * @param nibble_fraction The fraction of the (original) mass of the gas
+ *        particle that is removed.
+ */
+__attribute__((always_inline)) INLINE static void
+chemistry_transfer_part_to_bpart(struct chemistry_bpart_data* bp_data,
+                                 struct chemistry_part_data* p_data,
+                                 const double nibble_mass,
+                                 const double nibble_fraction) {
+  error("To be implemented.");
+}
+
+/**
+ * @brief Add the chemistry data of a black hole to another one.
+ *
+ * Nothing to do here.
+ *
+ * @param bp_data The black hole data to add to.
+ * @param swallowed_data The black hole data to use.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_add_bpart_to_bpart(
+    struct chemistry_bpart_data* bp_data,
+    const struct chemistry_bpart_data* swallowed_data) {
+  error("To be implemented.");
+}
+
+/**
+ * @brief Add the chemistry data of a sink particle to a sink.
+ *
+ * Nothing to do here.
+ *
+ * @param si_data The black hole data to add to.
+ * @param sj_data The gas data to use.
+ * @param gas_mass The mass of the gas particle.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_add_sink_to_sink(
+    struct chemistry_sink_data* si_data,
+    const struct chemistry_sink_data* sj_data) {
+  error("To be implemented.");
+}
+
+/**
+ * @brief Add the chemistry data of a gas particle to a sink.
+ *
+ * Nothing to do here.
+ *
+ * @param sp_data The sink data to add to.
+ * @param p_data The gas data to use.
+ * @param gas_mass The mass of the gas particle.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_add_part_to_sink(
+    struct chemistry_sink_data* sp_data,
+    const struct chemistry_part_data* p_data, const double gas_mass) {
+  error("To be implemented.");
+}
+
+/**
+ * @brief Split the metal content of a particle into n pieces
+ *
+ * Nothing to do here.
+ *
+ * @param p The #part.
+ * @param n The number of pieces to split into.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_split_part(
+    struct part* p, const double n) {
+  error("To be implemented.");
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in feedback/enrichment related routines.
+ *
+ * No metallicity treatment here -> return 0.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_fraction_for_feedback(
+    const struct part* restrict p) {
+  error("To be implemented.");
+  return 0.f;
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * gas particle to be used in feedback/enrichment related routines.
+ *
+ * No metallicity treatment here -> return NULL array.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float const*
+chemistry_get_metal_mass_fraction_for_feedback(const struct part* restrict p) {
+  error("To be implemented.");
+  return NULL;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * star particle to be used in feedback/enrichment related routines.
+ *
+ * No metallicity treatment here -> return 0.
+ *
+ * @param sp Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_star_total_metal_mass_fraction_for_feedback(
+    const struct spart* restrict sp) {
+
+  return sp->chemistry_data
+      .metal_mass_fraction[AGORA_CHEMISTRY_ELEMENT_COUNT - 1];
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * star particle to be used in feedback/enrichment related routines.
+ *
+ * No metallicity treatment here -> return NULL array.
+ *
+ * @param sp Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double const*
+chemistry_get_star_metal_mass_fraction_for_feedback(
+    const struct spart* restrict sp) {
+
+  return sp->chemistry_data.metal_mass_fraction;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in cooling related routines.
+ *
+ * No metallicity treatment here -> return 0.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_fraction_for_cooling(
+    const struct part* restrict p) {
+
+  return p->chemistry_data
+      .smoothed_metal_mass_fraction[AGORA_CHEMISTRY_ELEMENT_COUNT - 1];
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * gas particle to be used in cooling related routines.
+ *
+ * No metallicity treatment here -> return NULL array.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double const*
+chemistry_get_metal_mass_fraction_for_cooling(const struct part* restrict p) {
+
+  return p->chemistry_data.smoothed_metal_mass_fraction;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * gas particle to be used in star formation related routines.
+ *
+ * No metallicity treatment here -> return 0.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double
+chemistry_get_total_metal_mass_fraction_for_star_formation(
+    const struct part* restrict p) {
+
+  return p->chemistry_data
+      .smoothed_metal_mass_fraction[AGORA_CHEMISTRY_ELEMENT_COUNT - 1];
+}
+
+/**
+ * @brief Returns the abundance array (metal mass fractions) of the
+ * gas particle to be used in star formation related routines.
+ *
+ * No metallicity treatment here -> return NULL array.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double const*
+chemistry_get_metal_mass_fraction_for_star_formation(
+    const struct part* restrict p) {
+
+  return p->chemistry_data.smoothed_metal_mass_fraction;
+}
+
+/**
+ * @brief Returns the total metal mass of the
+ * gas particle to be used in the stats related routines.
+ *
+ * No metallicity treatment here -> return 0.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_total_metal_mass_for_stats(const struct part* restrict p) {
+
+  return p->chemistry_data.metal_mass[AGORA_CHEMISTRY_ELEMENT_COUNT - 1];
+}
+
+/**
+ * @brief Returns the total metal mass of the
+ * star particle to be used in the stats related routines.
+ *
+ * No metallicity treatment here -> return 0.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_star_total_metal_mass_for_stats(const struct spart* restrict sp) {
+
+  return sp->chemistry_data
+             .metal_mass_fraction[AGORA_CHEMISTRY_ELEMENT_COUNT - 1] *
+         sp->mass;
+}
+
+/**
+ * @brief Returns the total metal mass of the
+ * black hole particle to be used in the stats related routines.
+ *
+ * No metallicity treatment here -> return 0.
+ *
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_bh_total_metal_mass_for_stats(const struct bpart* restrict bp) {
+
+  error("To be implemented.");
+  return 0.f;
+}
+
+/**
+ * @brief Returns the total metallicity (metal mass fraction) of the
+ * star particle to be used in the luminosity calculations.
+ *
+ * No metallicity treatment here -> return 0.
+ *
+ * @param sp Pointer to the star particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+chemistry_get_star_total_metal_mass_fraction_for_luminosity(
+    const struct spart* restrict sp) {
+
+  error("To be implemented.");
+  return 0.f;
+}
+
+#endif /* SWIFT_CHEMISTRY_AGORA_H */