From 5f251335d2133786d341114d4e7a3f681fc9d77a Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Thu, 15 Aug 2019 14:45:08 +0200
Subject: [PATCH] Final updates to the chemistry schemes backported from the
 COLIBRE fork.

---
 configure.ac                    |  2 +-
 src/chemistry/EAGLE/chemistry.h |  8 ++++++--
 src/chemistry/GEAR/chemistry.h  | 30 ++++++++++++++++++++++++++++++
 src/chemistry/none/chemistry.h  | 30 ++++++++++++++++++++++++++++++
 src/runner.c                    |  1 +
 src/timestep.h                  | 10 ++++++++--
 6 files changed, 76 insertions(+), 5 deletions(-)

diff --git a/configure.ac b/configure.ac
index 45c0ab5224..6bbb800db7 100644
--- a/configure.ac
+++ b/configure.ac
@@ -2074,8 +2074,8 @@ AC_MSG_RESULT([
    Make gravity glass  : $gravity_glass_making
    External potential  : $with_potential
 
-   Entropy floor        : $with_entropy_floor
    Pressure floor       : $with_pressure_floor
+   Entropy floor        : $with_entropy_floor
    Cooling function     : $with_cooling
    Chemistry            : $with_chemistry
    Tracers              : $with_tracers
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index c06b29f441..a470eef3fa 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -243,8 +243,9 @@ static INLINE void chemistry_print_backend(
 }
 
 /**
- * @brief Updates the metal mass fractions after diffusion at the end of the
- * force loop.
+ * @brief Updates to the chemistry data after the hydro force loop.
+ *
+ * Nothing to do here in EAGLE.
  *
  * @param p The particle to act upon.
  * @param cosmo The current cosmological model.
@@ -255,10 +256,13 @@ __attribute__((always_inline)) INLINE static void chemistry_end_force(
 /**
  * @brief Computes the chemistry-related time-step constraint.
  *
+ * No constraints in the EAGLE model (no diffusion etc.) --> FLT_MAX
+ *
  * @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(
diff --git a/src/chemistry/GEAR/chemistry.h b/src/chemistry/GEAR/chemistry.h
index 951d565337..34c4d10e5f 100644
--- a/src/chemistry/GEAR/chemistry.h
+++ b/src/chemistry/GEAR/chemistry.h
@@ -135,6 +135,15 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density(
   }
 }
 
+/**
+ * @brief Updates to the chemistry data after the hydro force loop.
+ *
+ * @param p The particle to act upon.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_force(
+    struct part* restrict p, const struct cosmology* cosmo) {}
+
 /**
  * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
  *
@@ -151,6 +160,27 @@ chemistry_part_has_no_neighbours(struct part* restrict p,
   error("Needs implementing!");
 }
 
+/**
+ * @brief Computes the chemistry-related time-step constraint.
+ *
+ * No constraints in the GEAR model (no diffusion) --> FLT_MAX
+ *
+ * @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 the chemistry properties of the (x-)particles to a valid start
  * state.
diff --git a/src/chemistry/none/chemistry.h b/src/chemistry/none/chemistry.h
index 543a5e77ee..bb35ea2535 100644
--- a/src/chemistry/none/chemistry.h
+++ b/src/chemistry/none/chemistry.h
@@ -87,6 +87,36 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density(
     struct part* restrict p, const struct chemistry_global_data* cd,
     const struct cosmology* cosmo) {}
 
+/**
+ * @brief Updates to the chemistry data after the hydro force loop.
+ *
+ * Nothing to do here in EAGLE.
+ *
+ * @param p The particle to act upon.
+ * @param cosmo The current cosmological model.
+ */
+__attribute__((always_inline)) INLINE static void chemistry_end_force(
+    struct part* restrict p, const struct cosmology* cosmo) {}
+
+/**
+ * @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.
  *
diff --git a/src/runner.c b/src/runner.c
index 05dc33e97e..ba11945f83 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -3618,6 +3618,7 @@ void runner_do_end_hydro_force(struct runner *r, struct cell *c, int timer) {
 
         /* Finish the force loop */
         hydro_end_force(p, cosmo);
+        chemistry_end_force(p, cosmo);
 
 #ifdef SWIFT_BOUNDARY_PARTICLES
 
diff --git a/src/timestep.h b/src/timestep.h
index c2b1a10fcb..cd9faaea61 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -146,8 +146,14 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
     new_dt_grav = min(new_dt_self_grav, new_dt_ext_grav);
   }
 
-  /* Final time-step is minimum of hydro and gravity */
-  float new_dt = min3(new_dt_hydro, new_dt_cooling, new_dt_grav);
+  /* Compute the next timestep (chemistry condition, e.g. diffusion) */
+  const float new_dt_chemistry =
+      chemistry_timestep(e->physical_constants, e->cosmology, e->internal_units,
+                         e->hydro_properties, e->chemistry, p);
+
+  /* Final time-step is minimum of hydro, gravity and subgrid */
+  float new_dt =
+      min4(new_dt_hydro, new_dt_cooling, new_dt_grav, new_dt_chemistry);
 
   /* Limit change in smoothing length */
   const float dt_h_change =
-- 
GitLab