Added a simple BlockPerfectGas class

experimental/AcousticSolver.cpp

 #include <AcousticSolver.hpp>
 #include <rang.hpp>
 
+#include <memory>
 
-KOKKOS_INLINE_FUNCTION
-double AcousticSolver::e(double rho, double p, double gamma) const
-{
-  return p/(rho*(gamma-1));
-}
-
-KOKKOS_INLINE_FUNCTION
-double AcousticSolver::p(double rho, double e, double gamma) const
-{
-  return (gamma-1)*rho*e;
-}
+#include <BlockPerfectGas.hpp>
 
 typedef const double my_double;
 
@@ -169,16 +160,12 @@ AcousticSolver::AcousticSolver(const long int& nj)
       gammaj[j] = 1.4;
     });
 
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    ej[j] = e(rhoj[j],pj[j],gammaj[j]);
-  });
+  BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
-  Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    Ej[j] = ej[j]+0.5*uj[j]*uj[j];
-  });
+  block_eos.updateEandCFromRhoP();
 
   Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
-    cj[j] = std::sqrt(gammaj[j]*pj[j]/rhoj[j]);
+      Ej[j] = ej[j]+0.5*uj[j]*uj[j];
   });
 
   Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
@@ -200,10 +187,6 @@ AcousticSolver::AcousticSolver(const long int& nj)
       t=tmax;
     }
     
-    // if (iteration%100 == 0) {
-    //   std::cout << "dt=" << dt << " t=" << t << " i=" << iteration << std::endl;
-    // }
-    
     Kokkos::View<double*> ur("ur", nr);
     Kokkos::View<double*> pr("pr", nr);
 
@@ -229,10 +212,14 @@ AcousticSolver::AcousticSolver(const long int& nj)
 	ej[j] = Ej[j] - 0.5 * uj[j]*uj[j];
       });
 
+    Kokkos::View<double*> xr_new("new_xr", nr);
+
     Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){
-      xr[r] += dt*ur[r];
+      xr_new[r] = xr[r] + dt*ur[r];
       });
 
+    xr = xr_new;
+
     Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
       xj[j] = 0.5*(xr[j]+xr[j+1]);
       Vj[j] = xr[j+1]-xr[j];
@@ -240,10 +227,10 @@ AcousticSolver::AcousticSolver(const long int& nj)
 
     Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
 	rhoj[j] = mj[j]/Vj[j];
-	pj[j] = p(rhoj[j], ej[j], gammaj[j]);
-	cj[j] = std::sqrt(gammaj[j]*pj[j]/rhoj[j]); // inv_mj*vj
       });
 
+    block_eos.updatePandCFromRhoE();    
+    
     ++iteration;
   }
 

experimental/BlockPerfectGas.hpp

+#ifndef BLOCK_PERFECTGAS_HPP
+#define BLOCK_PERFECTGAS_HPP
+
+struct BlockPerfectGas
+{
+private:
+  Kokkos::View<double*> m_rhoj;
+  Kokkos::View<double*> m_ej;
+  Kokkos::View<double*> m_pj;
+  Kokkos::View<double*> m_gammaj;
+  Kokkos::View<double*> m_cj;
+
+public:
+  BlockPerfectGas(Kokkos::View<double*> rhoj,
+		 Kokkos::View<double*> ej,
+		 Kokkos::View<double*> pj,
+		 Kokkos::View<double*> gammaj,
+		 Kokkos::View<double*> cj)
+    : m_rhoj(rhoj),
+      m_ej(ej),
+      m_pj(pj),
+      m_gammaj(gammaj),
+      m_cj(cj)
+  {
+    ;
+  }
+
+  void updatePandCFromRhoE()
+  {
+    const int nj = m_ej.size();
+    const Kokkos::View<const double*> rho = m_rhoj;
+    const Kokkos::View<const double*> e = m_ej;
+    const Kokkos::View<const double*> gamma = m_gammaj;
+
+    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
+	const double gamma_minus_one = gamma[j]-1;
+	m_pj[j] = gamma_minus_one*rho[j]*e[j];
+	m_cj[j] = std::sqrt(gamma[j]*gamma_minus_one*e[j]);
+      });
+  }
+
+  void updateEandCFromRhoP()
+  {
+    const int nj = m_ej.size();
+    const Kokkos::View<const double*> rho = m_rhoj;
+    const Kokkos::View<const double*> p = m_pj;
+    const Kokkos::View<const double*> gamma = m_gammaj;
+
+    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
+	m_ej[j] = p[j]/(rho[j]*(gamma[j]-1));
+      });
+
+    const Kokkos::View<const double*> e = m_ej;
+    Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
+	m_cj[j] = std::sqrt(gamma[j]*(gamma[j]-1)*e[j]);
+      });
+  }
+};
+
+
+#endif // BLOCK_PERFECTGAS_HPP