From ef042b2b9a17f1cc98b3e4e529f77f2573c6799d Mon Sep 17 00:00:00 2001
From: Darwin Roduit <darwin.roduit@epfl.ch>
Date: Mon, 25 Nov 2024 21:20:30 +0100
Subject: [PATCH] Add adaptive softening struct to gpart

---
 src/adaptive_softening_struct.h           | 48 +++++++++++++++++++++++
 src/gravity/MultiSoftening/gravity_part.h |  4 ++
 2 files changed, 52 insertions(+)

diff --git a/src/adaptive_softening_struct.h b/src/adaptive_softening_struct.h
index c8301af8e8..b4c2468154 100644
--- a/src/adaptive_softening_struct.h
+++ b/src/adaptive_softening_struct.h
@@ -33,6 +33,49 @@ struct adaptive_softening_part_data {
   float zeta;
 };
 
+/**
+ * @brief Gravity particle-carried fields for the adaptive softening scheme.
+ */
+
+/* 1) Compute zeta and omega during the tree walk with the tidal tensor of the
+  previous timestep. We can drift them to et an estimate of the current
+  timestep, but we'll do this later. At the same time, compute the new tidal
+  tensor. (store the old value in another variable.
+  ---> create a new adaptive_softening_add_correction_term() fct and call it
+       in runner_iact_grav_pp_full() and _truncated().
+
+       To finish the computations, we might need to add a task.
+
+  2) Compute the Gamma_ab terms during the force interaction, i.e. when we
+  compute the acceleration. Compute d_epsilon_dt. Pay attention to the
+  interaction with gas particles -> not the same softening rule.
+  --> create a new adaptive_softening_get_acc_term() and call it in
+
+  3) Update the softening with the new tidal tensor. We can do this probably in
+  grav_end_force or in kick2() ? See where the gas updates its value.
+
+
+  Notes:
+  We only need to compute zeta and omega over all particles in the softening
+  kernel. Outside the softening kernel (r>H), the d Phi/d_epsilon = 0.
+
+  When computing the tidal tensor, use do not use the same Phi. See Hopkins
+  paper.
+
+  At the end, add a minimal softening and maximal softening.
+*/
+struct adaptive_softening_gpart_data {
+
+  /*! Correction term for energy conservation eq 15 in Hopkins 2023.*/
+  float zeta;
+
+  /*! Correction term for energy conservation (second term) eq 35 */
+  float omega;
+
+  /*! Time derivative of the softening eq 37 */
+  float depsilon_dt;
+};
+
 #else
 
 /**
@@ -40,6 +83,11 @@ struct adaptive_softening_part_data {
  */
 struct adaptive_softening_part_data {};
 
+/**
+ * @brief Gravity particle-carried fields for the adaptive softening scheme.
+ */
+struct adaptive_softening_gpart_data {};
+
 #endif
 
 #endif /* SWIFT_ADAPTIVE_SOFTENING_STRUCT_H */
diff --git a/src/gravity/MultiSoftening/gravity_part.h b/src/gravity/MultiSoftening/gravity_part.h
index 0c0c8c8ca7..6276642db5 100644
--- a/src/gravity/MultiSoftening/gravity_part.h
+++ b/src/gravity/MultiSoftening/gravity_part.h
@@ -19,6 +19,7 @@
 #ifndef SWIFT_MULTI_SOFTENING_GRAVITY_PART_H
 #define SWIFT_MULTI_SOFTENING_GRAVITY_PART_H
 
+#include "adaptive_softening_struct.h"
 #include "csds.h"
 #include "fof_struct.h"
 
@@ -69,6 +70,9 @@ struct gpart {
   /*! Type of the #gpart (DM, gas, star, ...) */
   enum part_type type;
 
+  /*! Additional data used for adaptive softening */
+  struct adaptive_softening_gpart_data adaptive_softening_data;
+
 #ifdef WITH_CSDS
   /* Additional data for the particle csds */
   struct csds_part_data csds_data;
-- 
GitLab