diff --git a/experimental/AcousticSolver.cpp b/experimental/AcousticSolver.cpp index 3f9984d2bb70e6f88a1d3dae379ed31b5cd60290..84f17ac608d6abed1b4724aa655a8b71a6a7bc5d 100644 --- a/experimental/AcousticSolver.cpp +++ b/experimental/AcousticSolver.cpp @@ -1,18 +1,9 @@ #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){ @@ -199,10 +186,6 @@ AcousticSolver::AcousticSolver(const long int& nj) dt=tmax-t; 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,9 +227,9 @@ 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; } diff --git a/experimental/BlockPerfectGas.hpp b/experimental/BlockPerfectGas.hpp new file mode 100644 index 0000000000000000000000000000000000000000..c018399e5d0644762a2474f3a99da6de20f82ffa --- /dev/null +++ b/experimental/BlockPerfectGas.hpp @@ -0,0 +1,61 @@ +#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