diff --git a/experimental/AcousticSolverWithMesh.hpp b/experimental/AcousticSolverWithMesh.hpp index e5e7839029953e105d23117bcf811ee5a4bfdc0c..f035ced50369f19aeb9a0427781e68f43f0bd4ff 100644 --- a/experimental/AcousticSolverWithMesh.hpp +++ b/experimental/AcousticSolverWithMesh.hpp @@ -182,19 +182,6 @@ private: }); } - KOKKOS_INLINE_FUNCTION - double acoustic_dt(const Kokkos::View<const double*>& Vj, - const Kokkos::View<const double*>& cj) const { - Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - m_Vj_over_cj[j] = Vj[j]/cj[j]; - }); - - double dt = std::numeric_limits<double>::max(); - Kokkos::parallel_reduce(m_mesh.numberOfCells(), ReduceMin(m_Vj_over_cj), dt); - - return dt; - } - KOKKOS_INLINE_FUNCTION void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr, const Kokkos::View<const Rd*>& xj, @@ -225,7 +212,6 @@ private: Kokkos::View<double*> m_rhocj; Kokkos::View<double*> m_Vj_over_cj; - public: AcousticSolverWithMesh(MeshData& mesh_data, UnknownsType& unknowns) @@ -240,6 +226,26 @@ public: m_ur("ur", m_mesh.numberOfNodes()), m_rhocj("rho_c", m_mesh.numberOfCells()), m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells()) + { + ; + } + + KOKKOS_INLINE_FUNCTION + double acoustic_dt(const Kokkos::View<const double*>& Vj, + const Kokkos::View<const double*>& cj) const { + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + m_Vj_over_cj[j] = Vj[j]/cj[j]; + }); + + double dt = std::numeric_limits<double>::max(); + Kokkos::parallel_reduce(m_mesh.numberOfCells(), ReduceMin(m_Vj_over_cj), dt); + + return dt; + } + + + void computeNextStep(const double& t, const double& dt, + UnknownsType& unknowns) { Kokkos::View<double*> rhoj = unknowns.rhoj(); Kokkos::View<Rd*> uj = unknowns.uj(); @@ -249,77 +255,48 @@ public: Kokkos::View<double*> pj = unknowns.pj(); Kokkos::View<double*> gammaj = unknowns.gammaj(); Kokkos::View<double*> cj = unknowns.cj(); - const Kokkos::View<const double*> mj = unknowns.mj(); - - Kokkos::View<Rd*> xr = m_mesh.xr(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); - - const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes(); const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr(); + Kokkos::View<Rd*> xr = m_mesh.xr(); - Kokkos::View<double*> inv_mj("inv_mj",m_mesh.numberOfCells()); - inverse(mj, inv_mj); - - const double tmax=0.2; - double t=0; - - int itermax=std::numeric_limits<int>::max(); - int iteration=0; - - BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); - - while((t<tmax) and (iteration<itermax)) { - double dt = 0.4*acoustic_dt(Vj, cj); - if (t+dt<tmax) { - t+=dt; - } else { - dt=tmax-t; - t=tmax; - } - - computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr); - - const Kokkos::View<const Rd**> Fjr = m_Fjr; - const Kokkos::View<const Rd*> ur = m_ur; - const Kokkos::View<const unsigned short*> cell_nb_nodes - = m_connectivity.cellNbNodes(); - - Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { - Rd momentum_fluxes = zero; - double energy_fluxes = 0; - for (int R=0; R<cell_nb_nodes[j]; ++R) { - const int r=cell_nodes(j,R); - momentum_fluxes += Fjr(j,R); - energy_fluxes += (Fjr(j,R), ur[r]); - } - uj[j] -= (dt*inv_mj[j]) * momentum_fluxes; - Ej[j] -= (dt*inv_mj[j]) * energy_fluxes; - }); - - Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { - ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]); - }); + computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr); - Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r){ - xr[r] += dt*ur[r]; - }); + const Kokkos::View<const Rd**> Fjr = m_Fjr; + const Kokkos::View<const Rd*> ur = m_ur; + const Kokkos::View<const unsigned short*> cell_nb_nodes + = m_connectivity.cellNbNodes(); + const Kokkos::View<const unsigned int**>& cell_nodes + = m_connectivity.cellNodes(); - m_mesh_data.updateAllData(); + const Kokkos::View<const double*> inv_mj = unknowns.invMj(); + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { + Rd momentum_fluxes = zero; + double energy_fluxes = 0; + for (int R=0; R<cell_nb_nodes[j]; ++R) { + const int r=cell_nodes(j,R); + momentum_fluxes += Fjr(j,R); + energy_fluxes += (Fjr(j,R), ur[r]); + } + uj[j] -= (dt*inv_mj[j]) * momentum_fluxes; + Ej[j] -= (dt*inv_mj[j]) * energy_fluxes; + }); - Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - rhoj[j] = mj[j]/Vj[j]; - }); + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { + ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]); + }); - block_eos.updatePandCFromRhoE(); - - ++iteration; - } + Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r){ + xr[r] += dt*ur[r]; + }); - std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset - << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n"; + m_mesh_data.updateAllData(); + const Kokkos::View<const double*> mj = unknowns.mj(); + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + rhoj[j] = mj[j]/Vj[j]; + }); } }; diff --git a/experimental/FiniteVolumesEulerUnknowns.hpp b/experimental/FiniteVolumesEulerUnknowns.hpp index 8f4ac4e572f487fa2c8937063fac949484bf2bf3..617cb8b39759dcf4ae4ff6f4a41f09a7bc056584 100644 --- a/experimental/FiniteVolumesEulerUnknowns.hpp +++ b/experimental/FiniteVolumesEulerUnknowns.hpp @@ -24,6 +24,7 @@ private: Kokkos::View<double*> m_gammaj; Kokkos::View<double*> m_cj; Kokkos::View<double*> m_mj; + Kokkos::View<double*> m_inv_mj; public: Kokkos::View<double*> rhoj() @@ -106,6 +107,16 @@ public: return m_mj; } + Kokkos::View<double*> invMj() + { + return m_inv_mj; + } + + const Kokkos::View<const double*> invMj() const + { + return m_inv_mj; + } + void initializeSod() { const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); @@ -146,6 +157,9 @@ public: m_mj[j] = m_rhoj[j] * Vj[j]; }); + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + m_inv_mj[j] = 1./m_mj[j]; + }); } FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data) @@ -158,7 +172,8 @@ public: m_pj("pj",m_mesh.numberOfCells()), m_gammaj("gammaj",m_mesh.numberOfCells()), m_cj("cj",m_mesh.numberOfCells()), - m_mj("mj",m_mesh.numberOfCells()) + m_mj("mj",m_mesh.numberOfCells()), + m_inv_mj("inv_mj",m_mesh.numberOfCells()) { ; } diff --git a/main.cpp b/main.cpp index 56081e71993150e70fbfdba010e6b2907961a74d..60a288255f8b3d00c8f20eebadc1053a2693c30c 100644 --- a/main.cpp +++ b/main.cpp @@ -125,11 +125,46 @@ int main(int argc, char *argv[]) MeshType mesh(connectivity); MeshDataType mesh_data(mesh); UnknownsType unknowns(mesh_data); - + unknowns.initializeSod(); AcousticSolverWithMesh<MeshDataType> acoustic_solver(mesh_data, unknowns); + const Kokkos::View<const double*> Vj = mesh_data.Vj(); + const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr(); + + const double tmax=0.2; + double t=0; + + int itermax=std::numeric_limits<int>::max(); + int iteration=0; + + Kokkos::View<double*> rhoj = unknowns.rhoj(); + Kokkos::View<double*> ej = unknowns.ej(); + Kokkos::View<double*> pj = unknowns.pj(); + Kokkos::View<double*> gammaj = unknowns.gammaj(); + Kokkos::View<double*> cj = unknowns.cj(); + + BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); + + while((t<tmax) and (iteration<itermax)) { + double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj); + if (t+dt>tmax) { + dt=tmax-t; + t=tmax; + } + + acoustic_solver.computeNextStep(t,dt, unknowns); + + block_eos.updatePandCFromRhoE(); + + t += dt; + ++iteration; + } + + std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset + << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n"; + method_cost_map["AcousticSolverWithMesh"] = timer.seconds(); { // gnuplot output for density