From 421afab0139a55756fe348cb4d833ab94e91f19a Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Thu, 26 Nov 2020 10:01:14 +0000
Subject: [PATCH] Store a_grav back in the xpart to avoid an incorrect drift
 operation where the gpart->a_grav was sometimes accessed after the next step
 had already started accumulating accelerations. That was a bad mistake
 introduced in the change to the gravity mesh integration, which primarilty
 affected the Planetary calculations due to their choice of system of units.

---
 src/drift.h                                           | 6 +++---
 src/hydro/AnarchyPU/hydro_part.h                      | 3 +++
 src/hydro/Gadget2/hydro_part.h                        | 3 +++
 src/hydro/Gizmo/hydro_part.h                          | 3 +++
 src/hydro/Minimal/hydro_part.h                        | 3 +++
 src/hydro/Phantom/hydro_part.h                        | 3 +++
 src/hydro/Planetary/hydro_part.h                      | 3 +++
 src/hydro/PressureEnergy/hydro_part.h                 | 3 +++
 src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h | 3 +++
 src/hydro/PressureEntropy/hydro_part.h                | 3 +++
 src/hydro/SPHENIX/hydro_part.h                        | 3 +++
 src/hydro/Shadowswift/hydro_part.h                    | 3 +++
 src/runner_time_integration.c                         | 8 ++++++++
 13 files changed, 44 insertions(+), 3 deletions(-)

diff --git a/src/drift.h b/src/drift.h
index e994e4fde7..7af873a542 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -151,9 +151,9 @@ __attribute__((always_inline)) INLINE static void drift_part(
 
   /* Predict velocities (for gravity terms) */
   if (p->gpart != NULL) {
-    p->v[0] += (p->gpart->a_grav[0] + p->gpart->a_grav_mesh[0]) * dt_kick_grav;
-    p->v[1] += (p->gpart->a_grav[1] + p->gpart->a_grav_mesh[1]) * dt_kick_grav;
-    p->v[2] += (p->gpart->a_grav[2] + p->gpart->a_grav_mesh[2]) * dt_kick_grav;
+    p->v[0] += xp->a_grav[0] * dt_kick_grav;
+    p->v[1] += xp->a_grav[1] * dt_kick_grav;
+    p->v[2] += xp->a_grav[2] * dt_kick_grav;
   }
 
   /* Predict the values of the extra fields */
diff --git a/src/hydro/AnarchyPU/hydro_part.h b/src/hydro/AnarchyPU/hydro_part.h
index d55307bf9f..185416f317 100644
--- a/src/hydro/AnarchyPU/hydro_part.h
+++ b/src/hydro/AnarchyPU/hydro_part.h
@@ -55,6 +55,9 @@ struct xpart {
   /*! Velocity at the last full step. */
   float v_full[3];
 
+  /*! Gravitational acceleration at the end of the last step */
+  float a_grav[3];
+
   /*! Internal energy at the last full step. */
   float u_full;
 
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index bb614c8975..fc34fa8026 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -54,6 +54,9 @@ struct xpart {
   /* Velocity at the last full step. */
   float v_full[3];
 
+  /*! Gravitational acceleration at the end of the last step */
+  float a_grav[3];
+
   /* Entropy at the last full step. */
   float entropy_full;
 
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index 3244aa1b6c..17b0179916 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -39,6 +39,9 @@ struct xpart {
   /* Velocity at the last full step. */
   float v_full[3];
 
+  /*! Gravitational acceleration at the end of the last step */
+  float a_grav[3];
+
   /* Additional data used to record cooling information */
   struct cooling_xpart_data cooling_data;
 
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index 6db84cdc01..94aaf51282 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -58,6 +58,9 @@ struct xpart {
   /*! Velocity at the last full step. */
   float v_full[3];
 
+  /*! Gravitational acceleration at the end of the last step */
+  float a_grav[3];
+
   /*! Internal energy at the last full step. */
   float u_full;
 
diff --git a/src/hydro/Phantom/hydro_part.h b/src/hydro/Phantom/hydro_part.h
index 170c0d9e43..73a9c4bde7 100644
--- a/src/hydro/Phantom/hydro_part.h
+++ b/src/hydro/Phantom/hydro_part.h
@@ -56,6 +56,9 @@ struct xpart {
   /*! Velocity at the last full step. */
   float v_full[3];
 
+  /*! Gravitational acceleration at the end of the last step */
+  float a_grav[3];
+
   /*! Internal energy at the last full step. */
   float u_full;
 
diff --git a/src/hydro/Planetary/hydro_part.h b/src/hydro/Planetary/hydro_part.h
index bcd905acc4..40d2f77835 100644
--- a/src/hydro/Planetary/hydro_part.h
+++ b/src/hydro/Planetary/hydro_part.h
@@ -61,6 +61,9 @@ struct xpart {
   /*! Velocity at the last full step. */
   float v_full[3];
 
+  /*! Gravitational acceleration at the end of the last step */
+  float a_grav[3];
+
   /*! Internal energy at the last full step. */
   float u_full;
 
diff --git a/src/hydro/PressureEnergy/hydro_part.h b/src/hydro/PressureEnergy/hydro_part.h
index 66b44d56c8..8be543aadc 100644
--- a/src/hydro/PressureEnergy/hydro_part.h
+++ b/src/hydro/PressureEnergy/hydro_part.h
@@ -58,6 +58,9 @@ struct xpart {
   /*! Velocity at the last full step. */
   float v_full[3];
 
+  /*! Gravitational acceleration at the end of the last step */
+  float a_grav[3];
+
   /*! Internal energy at the last full step. */
   float u_full;
 
diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
index 30db47c9d2..605bdfb8fb 100644
--- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
+++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro_part.h
@@ -58,6 +58,9 @@ struct xpart {
   /*! Velocity at the last full step. */
   float v_full[3];
 
+  /*! Gravitational acceleration at the end of the last step */
+  float a_grav[3];
+
   /*! Internal energy at the last full step. */
   float u_full;
 
diff --git a/src/hydro/PressureEntropy/hydro_part.h b/src/hydro/PressureEntropy/hydro_part.h
index 9a03c8862a..3177475b1c 100644
--- a/src/hydro/PressureEntropy/hydro_part.h
+++ b/src/hydro/PressureEntropy/hydro_part.h
@@ -50,6 +50,9 @@ struct xpart {
   /*! Velocity at the last full step. */
   float v_full[3];
 
+  /*! Gravitational acceleration at the end of the last step */
+  float a_grav[3];
+
   /*! Entropy at the last full step. */
   float entropy_full;
 
diff --git a/src/hydro/SPHENIX/hydro_part.h b/src/hydro/SPHENIX/hydro_part.h
index 4a3c287790..f13e29c0e3 100644
--- a/src/hydro/SPHENIX/hydro_part.h
+++ b/src/hydro/SPHENIX/hydro_part.h
@@ -53,6 +53,9 @@ struct xpart {
   /*! Velocity at the last full step. */
   float v_full[3];
 
+  /*! Gravitational acceleration at the end of the last step */
+  float a_grav[3];
+
   /*! Internal energy at the last full step. */
   float u_full;
 
diff --git a/src/hydro/Shadowswift/hydro_part.h b/src/hydro/Shadowswift/hydro_part.h
index bc374e5500..2ed26b8568 100644
--- a/src/hydro/Shadowswift/hydro_part.h
+++ b/src/hydro/Shadowswift/hydro_part.h
@@ -41,6 +41,9 @@ struct xpart {
   /* Velocity at the last full step. */
   float v_full[3];
 
+  /*! Gravitational acceleration at the end of the last step */
+  float a_grav[3];
+
   /* Additional data used to record cooling information */
   struct cooling_xpart_data cooling_data;
 
diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c
index b7b014b731..bfb2611e5f 100644
--- a/src/runner_time_integration.c
+++ b/src/runner_time_integration.c
@@ -177,6 +177,14 @@ void runner_do_kick1(struct runner *r, struct cell *c, const int timer) {
         kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_mesh_grav,
                   dt_kick_therm, dt_kick_corr, cosmo, hydro_props,
                   entropy_floor, ti_begin, ti_end, ti_begin_mesh, ti_end_mesh);
+
+        /* Update the accelerations to be used in the drift for hydro */
+        if (p->gpart != NULL) {
+
+          xp->a_grav[0] = p->gpart->a_grav[0] + p->gpart->a_grav_mesh[0];
+          xp->a_grav[1] = p->gpart->a_grav[1] + p->gpart->a_grav_mesh[1];
+          xp->a_grav[2] = p->gpart->a_grav[2] + p->gpart->a_grav_mesh[2];
+        }
       }
     }
 
-- 
GitLab