diff --git a/experimental/AcousticSolverWithMesh.hpp b/experimental/AcousticSolverWithMesh.hpp index 4755e8f0c3002035372d312a7861f98cdf7f1244..9be4585db426f135b411d0af4d38c0cc413016f3 100644 --- a/experimental/AcousticSolverWithMesh.hpp +++ b/experimental/AcousticSolverWithMesh.hpp @@ -9,11 +9,14 @@ #include <TinyVector.hpp> #include <TinyMatrix.hpp> #include <Mesh.hpp> +#include <MeshData.hpp> - -template<typename MeshType> +template<typename MeshData> class AcousticSolverWithMesh { + typedef typename MeshData::MeshType MeshType; + + MeshData& m_mesh_data; const MeshType& m_mesh; const typename MeshType::Connectivity& m_connectivity; @@ -225,11 +228,12 @@ private: public: - AcousticSolverWithMesh(MeshType& mesh) - : m_mesh(mesh), + AcousticSolverWithMesh(MeshData& mesh_data) + : m_mesh_data(mesh_data), + m_mesh(mesh_data.mesh()), m_connectivity(m_mesh.connectivity()), - m_nj(mesh.numberOfCells()), - m_nr(mesh.numberOfNodes()), + m_nj(m_mesh.numberOfCells()), + m_nr(m_mesh.numberOfNodes()), m_br("br", m_nr), m_Ajr("Ajr", m_nj, m_connectivity.maxNbNodePerCell()), m_Ar("Ar", m_nr), @@ -239,7 +243,6 @@ public: m_rhocj("rho_c", m_nj), m_Vj_over_cj("Vj_over_cj", m_nj) { - Kokkos::View<Rd*> xj("xj",m_nj); Kokkos::View<double*> rhoj("rhoj",m_nj); Kokkos::View<Rd*> uj("uj",m_nj); @@ -247,28 +250,17 @@ public: Kokkos::View<double*> Ej("Ej",m_nj); Kokkos::View<double*> ej("ej",m_nj); Kokkos::View<double*> pj("pj",m_nj); - Kokkos::View<double*> Vj("Vj",m_nj); Kokkos::View<double*> gammaj("gammaj",m_nj); Kokkos::View<double*> cj("cj",m_nj); Kokkos::View<double*> mj("mj",m_nj); - Kokkos::View<Rd*> xr = mesh.xr(); - - const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes(); - Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){ - xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]); - }); + Kokkos::View<Rd*> xr = m_mesh.xr(); - Kokkos::View<Rd*[2]> Cjr("Cjr",m_nj); - Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) { - Cjr(j,0)=-1; - Cjr(j,1)= 1; - }); + const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); + const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); - Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){ - Vj[j] = (xr[cell_nodes(j,1)], Cjr(j,1)) - + (xr[cell_nodes(j,0)], Cjr(j,0)); - }); + const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes(); + const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr(); Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){ if (xj[j][0]<0.5) { @@ -351,15 +343,7 @@ public: xr[r] += dt*ur[r]; }); - Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){ - xj[j] = 0.5*(xr[j]+xr[j+1]); - }); - - Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){ - Vj[j] - = (xr[cell_nodes(j,1)], Cjr(j,1)) - + (xr[cell_nodes(j,0)], Cjr(j,0)); - }); + m_mesh_data.updateAllData(); Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){ rhoj[j] = mj[j]/Vj[j]; diff --git a/experimental/Mesh.hpp b/experimental/Mesh.hpp index 98cd0af12499eaa291d890f549a3a42b5eac716e..f199266b5ce8efef07d946a40134690233ae8952 100644 --- a/experimental/Mesh.hpp +++ b/experimental/Mesh.hpp @@ -39,7 +39,7 @@ public: } #warning PASSER CES NOUVEAUX VECTEURS EN CONST - Kokkos::View<Rd*> xr() + Kokkos::View<Rd*> xr() const { return m_xr; } diff --git a/experimental/MeshData.hpp b/experimental/MeshData.hpp new file mode 100644 index 0000000000000000000000000000000000000000..c2125f5f0b04a47b7b982102949f81cf94b27d90 --- /dev/null +++ b/experimental/MeshData.hpp @@ -0,0 +1,138 @@ +#ifndef MESH_DATA_HPP +#define MESH_DATA_HPP + +#include <Kokkos_Core.hpp> +#include <TinyVector.hpp> + +template <typename M> +class MeshData +{ +public: + typedef M MeshType; + + typedef TinyVector<dimension> Rd; + typedef double R; + + static constexpr size_t dimension = MeshType::dimension; + static_assert(dimension>0, "dimension must be strictly positive"); + + static constexpr R inv_dimension = 1./dimension; + +private: + const MeshType& m_mesh; + Kokkos::View<Rd**> m_Cjr; + Kokkos::View<Rd*> m_xj; + Kokkos::View<R*> m_Vj; + + template<size_t Dim> + KOKKOS_INLINE_FUNCTION + void _updateCenter(); + + template <> + KOKKOS_INLINE_FUNCTION + void _updateCenter<1>() + { + const Kokkos::View<const Rd*> xr = m_mesh.xr(); + + const Kokkos::View<const unsigned int**>& cell_nodes + = m_mesh.connectivity().cellNodes(); + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + m_xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]); + }); + } + + template<size_t Dim> + KOKKOS_INLINE_FUNCTION + void _updateVolume() + { + const Kokkos::View<const unsigned int**>& cell_nodes + = m_mesh.connectivity().cellNodes(); + + const Kokkos::View<const unsigned short*> cell_nb_nodes + = m_mesh.connectivity().cellNbNodes(); + + const Kokkos::View<const Rd*> xr = m_mesh.xr(); + + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + R sum_cjr_xr = 0; + for (int R=0; R<cell_nb_nodes[j]; ++R) { + sum_cjr_xr += (xr[cell_nodes(j,R)], m_Cjr(j,R)); + } + m_Vj[j] = inv_dimension * sum_cjr_xr; + }); + } + + template<> + KOKKOS_INLINE_FUNCTION + void _updateVolume<1>() + { + const Kokkos::View<const Rd*> xr = m_mesh.xr(); + + const Kokkos::View<const unsigned int**>& cell_nodes + = m_mesh.connectivity().cellNodes(); + + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + m_Vj[j] + = (xr[cell_nodes(j,1)], m_Cjr(j,1)) + + (xr[cell_nodes(j,0)], m_Cjr(j,0)); + }); + } + + template<size_t Dim> + KOKKOS_INLINE_FUNCTION + void _updateCjr(); + + template<> + KOKKOS_INLINE_FUNCTION + void _updateCjr<1>() {} + +public: + + const MeshType& mesh() const + { + return m_mesh; + } + + + const Kokkos::View<const Rd**> Cjr() const + { + return m_Cjr; + } + + const Kokkos::View<const Rd*> xj() const + { + return m_xj; + } + + const Kokkos::View<const R*> Vj() const + { + return m_Vj; + } + + void updateAllData() + { + this->_updateCenter<dimension>(); + this->_updateVolume<dimension>(); + this->_updateCjr<dimension>(); + } + + MeshData(const MeshType& mesh) + : m_mesh(mesh), + m_Cjr("Cjr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()), + m_xj("xj", mesh.numberOfCells()), + m_Vj("Vj", mesh.numberOfCells()) + { + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { + m_Cjr(j,0)=-1; + m_Cjr(j,1)= 1; + }); + this->updateAllData(); + } + + ~MeshData() + { + ; + } +}; + +#endif // MESH_DATA_HPP diff --git a/main.cpp b/main.cpp index 66dbbd7d062e8433fa19dcc5e2196f02ddeba05f..2d42c5e22fc3024fb5db3639dd1fb42fb131249a 100644 --- a/main.cpp +++ b/main.cpp @@ -118,8 +118,13 @@ int main(int argc, char *argv[]) Kokkos::Timer timer; timer.reset(); Connectivity1D connectivity(number); - Mesh<Connectivity1D> mesh(connectivity); - AcousticSolverWithMesh<Mesh<Connectivity1D>> acoustic_solver(mesh); + typedef Mesh<Connectivity1D> MeshType; + typedef MeshData<MeshType> MeshDataType; + + MeshType mesh(connectivity); + MeshDataType mesh_data(mesh); + AcousticSolverWithMesh<MeshDataType> acoustic_solver(mesh_data); + method_cost_map["AcousticSolverWithMesh"] = timer.seconds(); }