sound speed calculation; EOS with temperature and internal energy (#478)

* Add factor of gamma in sound speed calculation; eos with temperature or internal energy
* change polytropic eos
* avoid two kernels for u and temp
* EOS choice with enum; EOS choice also in std hydro

main/src/init/turbulence_init.hpp

@@ -60,7 +60,7 @@ InitSettings TurbulenceConstants()
             {"mTotal", 1.0},
             {"powerLawExp", 5. / 3},
             {"anglesExp", 2.0},
-            {"eosChoice", 0},
+            {"eosChoice", sph::EosType::idealGas},
             {"gamma", 1.001},
             {"muiConst", 0.62},
             {"soundSpeedConst", 1.0},

sph/include/sph/eos.hpp

@@ -1,14 +1,19 @@
 #pragma once
 
-#include <vector>
-
 #include "cstone/util/tuple.hpp"
 
 #include "kernels.hpp"
 
 namespace sph
 {
 
+enum EosType : int
+{
+    idealGas   = 0,
+    isothermal = 1,
+    polytropic = 2
+};
+
 //! @brief returns the heat capacity for given mean molecular weight
 template<class T1, class T2>
 HOST_DEVICE_FUN constexpr T1 idealGasCv(T1 mui, T2 gamma)
@@ -21,24 +26,29 @@ HOST_DEVICE_FUN constexpr T1 idealGasCv(T1 mui, T2 gamma)
  *
  * @param u     internal energy
  * @param rho   baryonic density
- * @param mui   mean molecular weight
  * @param gamma adiabatic index
  *
  * This EOS is used for simple cases where we don't need the temperature.
  * Returns pressure, speed of sound
  */
 template<class T1, class T2, class T3>
-HOST_DEVICE_FUN auto idealGasEOS(T1 temp, T2 rho, T3 mui, T1 gamma)
+HOST_DEVICE_FUN auto idealGasEOS_u(T1 u, T2 rho, T3 gamma)
 {
     using Tc = std::common_type_t<T1, T2, T3>;
 
-    Tc tmp = idealGasCv(mui, gamma) * temp * (gamma - Tc(1));
+    Tc tmp = u * (gamma - Tc(1));
     Tc p   = rho * tmp;
-    Tc c   = std::sqrt(tmp);
+    Tc c   = std::sqrt(gamma * tmp);
 
     return util::tuple<Tc, Tc>{p, c};
 }
 
+template<class T1, class T2, class T3>
+HOST_DEVICE_FUN auto idealGasEOS(T1 temp, T2 rho, T3 mui, T1 gamma)
+{
+    return idealGasEOS_u(idealGasCv(mui, gamma) * temp, rho, gamma);
+}
+
 /*! @brief Isothermal equation of state
  *
  * @param c     speed of sound
@@ -53,49 +63,26 @@ HOST_DEVICE_FUN auto isothermalEOS(T1 c, T2 rho)
     return p;
 }
 
-/*! @brief Polytropic EOS for a 1.4 M_sun and 12.8 km neutron star
+/*! @brief General polytropic equation of state.
+ * @param K_poly       polytropic constant
+ * @param gamma_poly   polytropic exponent
+ * @param rho          SPH density
  *
- * @param rho  baryonic density
+ * Returns pressure and sound speed
  *
- * Kpol is hardcoded for these NS characteristics and is not valid for
- * other NS masses and radius
- * Returns pressure, and speed of sound
+ * For a 1.4 M_sun and 12.8 km neutron star the values are
+ * K_poly = 2.246341237993810232e-10
+ * gammapol = 3.e0
  */
-template<class T>
-HOST_DEVICE_FUN auto polytropicEOS(T rho)
+template<typename T1, typename T2, typename T3>
+HOST_DEVICE_FUN auto polytropicEOS(T1 K_poly, T2 gamma_poly, T3 rho)
 {
-    constexpr T Kpol     = 2.246341237993810232e-10;
-    constexpr T gammapol = 3.e0;
-
-    T p = Kpol * std::pow(rho, gammapol);
-    T c = std::sqrt(gammapol * p / rho);
-
-    return util::tuple<T, T>{p, c};
-}
-
-/*! @brief Polytropic EOS interface for SPH where rho is computed on-the-fly
- *
- * @tparam Dataset
- * @param startIndex  index of first locally owned particle
- * @param endIndex    index of last locally owned particle
- * @param d           the dataset with the particle buffers
- */
-template<typename Dataset>
-void computeEOS_Polytropic(size_t startIndex, size_t endIndex, Dataset& d)
-{
-    const auto* kx = d.kx.data();
-    const auto* xm = d.xm.data();
-    const auto* m  = d.m.data();
+    using Tc = std::common_type_t<T1, T2, T3>;
 
-    auto* p = d.p.data();
-    auto* c = d.c.data();
+    Tc p = K_poly * std::pow(rho, gamma_poly);
+    Tc c = std::sqrt(gamma_poly * p / rho);
 
-#pragma omp parallel for schedule(static)
-    for (size_t i = startIndex; i < endIndex; ++i)
-    {
-        auto rho             = kx[i] * m[i] / xm[i];
-        std::tie(p[i], c[i]) = polytropicEOS(rho);
-    }
+    return util::tuple<Tc, Tc>{p, c};
 }
 
 } // namespace sph

sph/include/sph/hydro_std/eos.hpp

@@ -52,30 +52,103 @@ namespace sph
  * we could potentially avoid halo exchange of p and c in return for exchanging halos of u.
  */
 template<typename Dataset>
-void computeEOS_HydroStdImpl(size_t startIndex, size_t endIndex, Dataset& d)
+void computeIdealGasEOS_HydroStd_Impl(size_t startIndex, size_t endIndex, Dataset& d)
 {
+    const auto* u    = d.u.data();
     const auto* temp = d.temp.data();
-    auto*       rho  = d.rho.data();
+    const auto* rho  = d.rho.data();
 
     auto* p = d.p.data();
     auto* c = d.c.data();
 
+    if (d.u.size() == 0)
+    {
+#pragma omp parallel for schedule(static)
+        for (size_t i = startIndex; i < endIndex; ++i)
+        {
+            std::tie(p[i], c[i]) = idealGasEOS(temp[i], rho[i], d.muiConst, d.gamma);
+        }
+    }
+    else
+    {
+#pragma omp parallel for schedule(static)
+        for (size_t i = startIndex; i < endIndex; ++i)
+        {
+            std::tie(p[i], c[i]) = idealGasEOS_u(u[i], rho[i], d.gamma);
+        }
+    }
+}
+
+template<typename Dataset>
+void computeIsothermalEOS_HydroStd_Impl(size_t startIndex, size_t endIndex, Dataset& d)
+{
+    const auto* rho = d.rho.data();
+
+    auto* p      = d.p.data();
+    auto* temp   = d.temp.data();
+    auto  cConst = d.soundSpeedConst;
+
+#pragma omp parallel for schedule(static)
+    for (size_t i = startIndex; i < endIndex; ++i)
+    {
+        p[i] = isothermalEOS(cConst, rho[i]);
+        if (temp) { temp[i] = 0; }
+    }
+}
+
+template<typename Dataset>
+void computePolytropicEOS_HydroStd_Impl(size_t startIndex, size_t endIndex, Dataset& d)
+{
+    const auto* rho = d.rho.data();
+
+    auto* p    = d.p.data();
+    auto* temp = d.temp.data();
+    auto* c    = d.c.data();
+
 #pragma omp parallel for schedule(static)
     for (size_t i = startIndex; i < endIndex; ++i)
     {
-        std::tie(p[i], c[i]) = idealGasEOS(temp[i], rho[i], d.muiConst, d.gamma);
+        std::tie(p[i], c[i]) = polytropicEOS(d.polytropic_const, d.polytropic_index, rho[i]);
+        if (temp) { temp[i] = 0; }
     }
 }
 
 template<class Dataset>
-void computeEOS_HydroStd(size_t startIndex, size_t endIndex, Dataset& d)
+void computeIdealGasEOS_HydroStd(size_t startIndex, size_t endIndex, Dataset& d)
 {
     if constexpr (cstone::HaveGpu<typename Dataset::AcceleratorType>{})
     {
-        cuda::computeEOS_HydroStd(startIndex, endIndex, d.muiConst, d.gamma, rawPtr(d.devData.temp),
-                                  rawPtr(d.devData.m), rawPtr(d.devData.rho), rawPtr(d.devData.p), rawPtr(d.devData.c));
+        cuda::computeIdealGasEOS_HydroStd(startIndex, endIndex, d);
     }
-    else { computeEOS_HydroStdImpl(startIndex, endIndex, d); }
+    else { computeIdealGasEOS_HydroStd_Impl(startIndex, endIndex, d); }
+}
+
+template<class Dataset>
+void computeIsothermalEOS_HydroStd(size_t startIndex, size_t endIndex, Dataset& d)
+{
+    if constexpr (cstone::HaveGpu<typename Dataset::AcceleratorType>{})
+    {
+        cuda::computeIsothermalEOS_HydroStd(startIndex, endIndex, d);
+    }
+    else { computeIsothermalEOS_HydroStd_Impl(startIndex, endIndex, d); }
+}
+
+template<class Dataset>
+void computePolytropicEOS_HydroStd(size_t startIndex, size_t endIndex, Dataset& d)
+{
+    if constexpr (cstone::HaveGpu<typename Dataset::AcceleratorType>{})
+    {
+        cuda::computePolytropicEOS_HydroStd(startIndex, endIndex, d);
+    }
+    else { computePolytropicEOS_HydroStd_Impl(startIndex, endIndex, d); }
+}
+
+template<class Dataset>
+void computeEOS_HydroStd(size_t startIndex, size_t endIndex, Dataset& d)
+{
+    if (d.eosChoice == EosType::idealGas) { computeIdealGasEOS_HydroStd(startIndex, endIndex, d); }
+    else if (d.eosChoice == EosType::isothermal) { computeIsothermalEOS_HydroStd(startIndex, endIndex, d); }
+    else if (d.eosChoice == EosType::polytropic) { computePolytropicEOS_HydroStd(startIndex, endIndex, d); }
 }
 
 } // namespace sph

sph/include/sph/hydro_std/eos_gpu.cu

@@ -35,37 +35,94 @@
 
 #include "sph/sph_gpu.hpp"
 #include "sph/eos.hpp"
+#include "sph/particles_data.hpp"
 
 namespace sph
 {
 namespace cuda
 {
 
-template<class Tt, class Trho, class Tp, class Tc>
-__global__ void cudaEOS_HydroStd(size_t firstParticle, size_t lastParticle, Trho mui, Tt gamma, const Tt* temp,
-                                 const Trho* m, Trho* rho, Tp* p, Tc* c)
+template<class Tt, class Tm, class Thydro>
+__global__ void cudaComputeIdealGasEOS_HydroStd(size_t firstParticle, size_t lastParticle, Tm mui, Tt gamma,
+                                                const Tt* temp, const Tt* u, const Tm* m, Thydro* rho, Thydro* p,
+                                                Thydro* c)
 {
     unsigned i = firstParticle + blockDim.x * blockIdx.x + threadIdx.x;
     if (i >= lastParticle) return;
 
-    util::tie(p[i], c[i]) = idealGasEOS(temp[i], rho[i], mui, gamma);
+    if (u == nullptr) { util::tie(p[i], c[i]) = idealGasEOS(temp[i], rho[i], mui, gamma); }
+    else { util::tie(p[i], c[i]) = idealGasEOS_u(u[i], rho[i], gamma); }
 }
 
-template<class Tt, class Trho, class Tp, class Tc>
-void computeEOS_HydroStd(size_t firstParticle, size_t lastParticle, Trho mui, Tt gamma, const Tt* temp, const Trho* m,
-                         Trho* rho, Tp* p, Tc* c)
+template<class Dataset>
+void computeIdealGasEOS_HydroStd(size_t firstParticle, size_t lastParticle, Dataset& d)
 {
     if (firstParticle == lastParticle) { return; }
     unsigned numThreads = 256;
     unsigned numBlocks  = cstone::iceil(lastParticle - firstParticle, numThreads);
-    cudaEOS_HydroStd<<<numBlocks, numThreads>>>(firstParticle, lastParticle, mui, gamma, temp, m, rho, p, c);
+
+    cudaComputeIdealGasEOS_HydroStd<<<numBlocks, numThreads>>>(
+        firstParticle, lastParticle, d.muiConst, d.gamma, rawPtr(d.devData.temp), rawPtr(d.devData.u),
+        rawPtr(d.devData.m), rawPtr(d.devData.rho), rawPtr(d.devData.p), rawPtr(d.devData.c));
+
+    checkGpuErrors(cudaDeviceSynchronize());
+}
+
+template void computeIdealGasEOS_HydroStd(size_t, size_t, sphexa::ParticlesData<cstone::GpuTag>&);
+
+template<typename Th, typename Tu>
+__global__ void cudaComputeIsothermalEOS_HydroStd(size_t first, size_t last, Th cConst, Th* c, Th* rho, Th* p, Tu* temp)
+{
+    unsigned i = first + blockDim.x * blockIdx.x + threadIdx.x;
+    if (i >= last) return;
+
+    p[i] = isothermalEOS(cConst, rho[i]);
+    c[i] = cConst;
+    if (temp) { temp[i] = 0; }
+}
+
+template<typename Dataset>
+void computeIsothermalEOS_HydroStd(size_t first, size_t last, Dataset& d)
+{
+    if (first == last) { return; }
+    unsigned numThreads = 256;
+    unsigned numBlocks  = cstone::iceil(last - first, numThreads);
+
+    cudaComputeIsothermalEOS_HydroStd<<<numBlocks, numThreads>>>(first, last, d.soundSpeedConst, rawPtr(d.devData.c),
+                                                                 rawPtr(d.devData.rho), rawPtr(d.devData.p),
+                                                                 rawPtr(d.devData.temp));
+
+    checkGpuErrors(cudaDeviceSynchronize());
+}
+
+template void computeIsothermalEOS_HydroStd(size_t, size_t, sphexa::ParticlesData<cstone::GpuTag>& d);
+
+template<typename Th, typename Tt>
+__global__ void cudaComputePolytropicEOS_HydroStd(size_t first, size_t last, Tt polytropic_const, Tt polytropic_index,
+                                                  Th* rho, Th* p, Tt* temp, Th* c)
+{
+    unsigned i = first + blockDim.x * blockIdx.x + threadIdx.x;
+    if (i >= last) return;
+
+    util::tie(p[i], c[i]) = polytropicEOS(polytropic_const, polytropic_index, rho[i]);
+    if (temp) { temp[i] = 0; }
+}
+
+template<typename Dataset>
+void computePolytropicEOS_HydroStd(size_t first, size_t last, Dataset& d)
+{
+    if (first == last) { return; }
+    unsigned numThreads = 256;
+    unsigned numBlocks  = cstone::iceil(last - first, numThreads);
+
+    cudaComputePolytropicEOS_HydroStd<<<numBlocks, numThreads>>>(first, last, d.polytropic_const, d.polytropic_index,
+                                                                 rawPtr(d.devData.rho), rawPtr(d.devData.p),
+                                                                 rawPtr(d.devData.temp), rawPtr(d.devData.c));
+
     checkGpuErrors(cudaDeviceSynchronize());
 }
 
-template void computeEOS_HydroStd(size_t, size_t, double, double, const double*, const double*, double*, double*,
-                                  double*);
-template void computeEOS_HydroStd(size_t, size_t, float, double, const double*, const float*, float*, float*, float*);
-template void computeEOS_HydroStd(size_t, size_t, float, float, const float*, const float*, float*, float*, float*);
+template void computePolytropicEOS_HydroStd(size_t, size_t, sphexa::ParticlesData<cstone::GpuTag>&);
 
 } // namespace cuda
 } // namespace sph