From 1f29ba87304b27a26ffdbc9feb58b57014f0c7e6 Mon Sep 17 00:00:00 2001
From: Sebastian Keller <sebastian.f.keller@gmail.com>
Date: Wed, 17 Jul 2024 16:29:46 +0200
Subject: [PATCH] fix underflow in enerygUpdate

---
 sph/include/sph/positions.hpp    | 10 ++++------
 sph/include/sph/positions_gpu.cu |  4 ++--
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/sph/include/sph/positions.hpp b/sph/include/sph/positions.hpp
index 3c19827ba..f0c1e9a28 100644
--- a/sph/include/sph/positions.hpp
+++ b/sph/include/sph/positions.hpp
@@ -51,8 +51,8 @@ HOST_DEVICE_FUN bool fbcCheck(Tc coord, Th h, Tc top, Tc bottom, bool fbc)
 }
 
 //! @brief update the energy according to Adams-Bashforth (2nd order)
-template<class TU, class TD>
-HOST_DEVICE_FUN TU energyUpdate(TU u_old, double dt, double dt_m1, TD du, TD du_m1)
+template<class TU>
+HOST_DEVICE_FUN TU energyUpdate(TU u_old, double dt, double dt_m1, double du, double du_m1)
 {
     TU u_new = u_old + du * dt + 0.5 * (du - du_m1) / dt_m1 * std::abs(dt) * dt;
     // To prevent u < 0 (when cooling with GRACKLE is active)
@@ -124,7 +124,6 @@ void updatePositionsHost(size_t startIndex, size_t endIndex, Dataset& d, const c
 template<class Dataset>
 void updateTempHost(size_t startIndex, size_t endIndex, Dataset& d)
 {
-    using Tdu    = decltype(d.du[0]);
     bool haveMui = !d.mui.empty();
     auto constCv = idealGasCv(d.muiConst, d.gamma);
 
@@ -133,7 +132,7 @@ void updateTempHost(size_t startIndex, size_t endIndex, Dataset& d)
     {
         auto cv    = haveMui ? idealGasCv(d.mui[i], d.gamma) : constCv;
         auto u_old = cv * d.temp[i];
-        d.temp[i]  = energyUpdate(u_old, d.minDt, d.minDt_m1, d.du[i], Tdu(d.du_m1[i])) / cv;
+        d.temp[i]  = energyUpdate(u_old, d.minDt, d.minDt_m1, d.du[i], d.du_m1[i]) / cv;
         d.du_m1[i] = d.du[i];
     }
 }
@@ -141,11 +140,10 @@ void updateTempHost(size_t startIndex, size_t endIndex, Dataset& d)
 template<class Dataset>
 void updateIntEnergyHost(size_t startIndex, size_t endIndex, Dataset& d)
 {
-    using Tdu = decltype(d.du[0]);
 #pragma omp parallel for schedule(static)
     for (size_t i = startIndex; i < endIndex; i++)
     {
-        d.u[i]     = energyUpdate(d.u[i], d.minDt, d.minDt_m1, d.du[i], Tdu(d.du_m1[i]));
+        d.u[i]     = energyUpdate(d.u[i], d.minDt, d.minDt_m1, d.du[i], d.du_m1[i]);
         d.du_m1[i] = d.du[i];
     }
 }
diff --git a/sph/include/sph/positions_gpu.cu b/sph/include/sph/positions_gpu.cu
index 190108d31..5310e572b 100644
--- a/sph/include/sph/positions_gpu.cu
+++ b/sph/include/sph/positions_gpu.cu
@@ -157,9 +157,9 @@ __global__ void computePositionsKernel(GroupView grp, float dt, util::array<floa
     {
         Thydro cv    = (constCv < 0) ? idealGasCv(mui[i], gamma) : constCv;
         auto   u_old = temp[i] * cv;
-        temp[i]      = energyUpdate(u_old, dt, dt_m1_rung, du[i], Tdu(du_m1[i])) / cv;
+        temp[i]      = energyUpdate(u_old, dt, dt_m1_rung, du[i], du_m1[i]) / cv;
     }
-    else if (u != nullptr) { u[i] = energyUpdate(u[i], dt, dt_m1_rung, du[i], Tdu(du_m1[i])); }
+    else if (u != nullptr) { u[i] = energyUpdate(u[i], dt, dt_m1_rung, du[i], du_m1[i]); }
     du_m1[i] = du[i];
 }