diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp index bed5cf97e0c7ed3a0ee60a14b6eaa35c83907d2d..7fc8f8c806798747eae0190558c35b75beb6d232 100644 --- a/src/mesh/MeshData.hpp +++ b/src/mesh/MeshData.hpp @@ -20,6 +20,8 @@ public: private: const MeshType& m_mesh; Kokkos::View<Rd**> m_Cjr; + Kokkos::View<double**> m_ljr; + Kokkos::View<Rd**> m_njr; Kokkos::View<Rd*> m_xj; Kokkos::View<double*> m_Vj; @@ -74,7 +76,7 @@ private: KOKKOS_INLINE_FUNCTION void _updateCjr() { if(dimension == 1) { - // Cjr are constant overtime + // Cjr/njr/ljr are constant overtime } static_assert(dimension==1, "only 1d is implemented"); } @@ -90,6 +92,16 @@ public: return m_Cjr; } + const Kokkos::View<const double**> ljr() const + { + return m_ljr; + } + + const Kokkos::View<const Rd**> njr() const + { + return m_njr; + } + const Kokkos::View<const Rd*> xj() const { return m_xj; @@ -110,6 +122,8 @@ public: MeshData(const MeshType& mesh) : m_mesh(mesh), m_Cjr("Cjr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()), + m_ljr("ljr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()), + m_njr("njr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()), m_xj("xj", mesh.numberOfCells()), m_Vj("Vj", mesh.numberOfCells()) { @@ -119,6 +133,12 @@ public: m_Cjr(j,0)=-1; m_Cjr(j,1)= 1; }); + // in 1d njr=Cjr + m_njr=m_Cjr; + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { + m_ljr(j,0)= 1; + m_ljr(j,1)= 1; + }); } this->updateAllData(); } diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index 1889f9f8381f4967ab2528958b37db5975492c19..6f6d9a86492e11b2bb80d40f9c3f4108637099ff 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -79,13 +79,14 @@ private: KOKKOS_INLINE_FUNCTION const Kokkos::View<const Rdd**> computeAjr(const Kokkos::View<const double*>& rhocj, - const Kokkos::View<const Rd**>& Cjr) { + const Kokkos::View<const double**>& ljr, + const Kokkos::View<const Rd**>& njr) { const Kokkos::View<const unsigned short*> cell_nb_nodes = m_connectivity.cellNbNodes(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { for (int r=0; r<cell_nb_nodes[j]; ++r) { - m_Ajr(j,r) = tensorProduct(rhocj(j)*Cjr(j,r), Cjr(j,r)); + m_Ajr(j,r) = tensorProduct((rhocj(j)*ljr(j,r))*njr(j,r), njr(j,r)); } }); @@ -190,9 +191,11 @@ private: const Kokkos::View<const double*>& pj, const Kokkos::View<const double*>& cj, const Kokkos::View<const double*>& Vj, - const Kokkos::View<const Rd**>& Cjr) { + const Kokkos::View<const Rd**>& Cjr, + const Kokkos::View<const double**>& ljr, + const Kokkos::View<const Rd**>& njr) { const Kokkos::View<const double*> rhocj = computeRhoCj(rhoj, cj); - const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, Cjr); + const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, ljr, njr); const Kokkos::View<const Rdd*> Ar = computeAr(Ajr); const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj); @@ -259,9 +262,11 @@ public: const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr(); + const Kokkos::View<const double**> ljr = m_mesh_data.ljr(); + const Kokkos::View<const Rd**> njr = m_mesh_data.njr(); Kokkos::View<Rd*> xr = m_mesh.xr(); - computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr); + computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr, ljr, njr); const Kokkos::View<const Rd**> Fjr = m_Fjr; const Kokkos::View<const Rd*> ur = m_ur;