fix underflow in enerygUpdate

sphexa-org · Aug 7, 2024 · 1f29ba8 · 1f29ba8
1 parent cf3407c
commit 1f29ba8
Show file tree

Hide file tree

Showing 2 changed files with 6 additions and 8 deletions.
diff --git 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,19 +132,18 @@ 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];
     }
 }
 
 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
@@ -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];
 }