diff --git a/src/main.cpp b/src/main.cpp index cd52860096e989d60bb046b5bced2e0ac7ba623f..d68d94569eb35c1bbd943ed35772f3d2c736db8a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -123,6 +123,87 @@ int main(int argc, char *argv[]) if (filename != "") { std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n"; GmshReader gmsh_reader(filename); + + typedef Mesh<Connectivity2D> MeshType; + typedef MeshData<MeshType> MeshDataType; + typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType; + + MeshType& mesh = gmsh_reader.mesh(); + + Kokkos::Timer timer; + timer.reset(); + MeshDataType mesh_data(mesh); + + std::vector<BoundaryConditionHandler> bc_list; + // { // quite dirty! + // SymmetryBoundaryCondition<MeshType::dimension>* sym_bc0 + // = new SymmetryBoundaryCondition<MeshType::dimension>(std::vector<unsigned int>({0u}), + // TinyVector<1>(-1)); + // std::shared_ptr<SymmetryBoundaryCondition<1>> bc0(sym_bc0); + // bc_list.push_back(BoundaryConditionHandler(bc0)); + + // PressureBoundaryCondition* pres_bc1 + // = new PressureBoundaryCondition(0.1, + // std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())})); + // std::shared_ptr<PressureBoundaryCondition> bc1(pres_bc1); + // bc_list.push_back(BoundaryConditionHandler(bc1)); + // } + + UnknownsType unknowns(mesh_data); + + unknowns.initializeSod(); + + AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list); + + typedef TinyVector<MeshType::dimension> Rd; + + 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; + } + + acoustic_solver.computeNextStep(t,dt, unknowns); + + block_eos.updatePandCFromRhoE(); + + t += dt; + ++iteration; + } + + std::ofstream gnuplot("sol.gnu"); + for (int j=0; j<mesh.numberOfCells(); ++j) { + for (int r=0; r<3; ++r) { + const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,r)]; + gnuplot << x[0] << ' ' << x[1] << '\n'; + } + const Rd& x = mesh.xr()[mesh.connectivity().cellNodes()(j,0)]; + gnuplot << x[0] << ' ' << x[1] << "\n\n"; + } + gnuplot.close(); + + 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(); + } else { // class for acoustic solver test Kokkos::Timer timer; @@ -207,8 +288,6 @@ int main(int argc, char *argv[]) Kokkos::finalize(); - std::cout << "----------------------\n"; - std::string::size_type size=0; for (const auto& method_cost : method_cost_map) { size = std::max(size, method_cost.first.size()); diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp index 55560c2f61aece668c5eb9419ace679576f6de44..640115039847a733a95ae80b2ff0d2e7926408b9 100644 --- a/src/mesh/Connectivity2D.hpp +++ b/src/mesh/Connectivity2D.hpp @@ -4,6 +4,10 @@ #include <Kokkos_Core.hpp> #include <TinyVector.hpp> +#include <vector> +#include <map> +#include <algorithm> + class Connectivity2D { public: @@ -14,20 +18,19 @@ private: size_t m_number_of_faces; size_t m_number_of_nodes; + const Kokkos::View<const unsigned short*> m_cell_nb_nodes; const Kokkos::View<const unsigned int**> m_cell_nodes; - Kokkos::View<unsigned int**> m_cell_faces; - - Kokkos::View<unsigned short*> m_cell_nb_nodes; Kokkos::View<double*> m_inv_cell_nb_nodes; + Kokkos::View<unsigned short*> m_cell_nb_faces; + Kokkos::View<unsigned int**> m_cell_faces; Kokkos::View<unsigned short*> m_node_nb_cells; - Kokkos::View<unsigned short*> m_face_nb_cells; - Kokkos::View<unsigned int**> m_node_cells; - Kokkos::View<unsigned int**> m_face_cells; - Kokkos::View<unsigned short**> m_node_cell_local_node; + + Kokkos::View<unsigned short*> m_face_nb_cells; + Kokkos::View<unsigned int**> m_face_cells; Kokkos::View<unsigned short**> m_face_cell_local_face; size_t m_max_nb_node_per_cell; @@ -110,13 +113,93 @@ public: Connectivity2D(const Connectivity2D&) = delete; - Connectivity2D(const Kokkos::View<const unsigned int**> cell_nodes) - : m_number_of_cells (cell_nodes.dimension_0()), + Connectivity2D(const Kokkos::View<const unsigned short*> cell_nb_nodes, + const Kokkos::View<const unsigned int**> cell_nodes) + : m_number_of_cells (cell_nodes.extent(0)), + m_cell_nb_nodes(cell_nb_nodes), m_cell_nodes (cell_nodes) { + { + Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", m_number_of_cells); + Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const int&j){ + inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j]; + }); + m_inv_cell_nb_nodes = inv_cell_nb_nodes; + } assert(m_number_of_cells>0); - std::cout << __PRETTY_FUNCTION__<< ": Number of cells = " << m_number_of_cells << '\n'; + // Computes inefficiently node->cells connectivity [Version 0] + + std::multimap<unsigned int, unsigned int> node_cells_map; + +#warning ONLY WORKS FOR TRIANGLES + // ------------------------------>>>> + m_max_nb_node_per_cell = 3; + for (unsigned int j=0; j<m_number_of_cells; ++j) { + for (unsigned int r=0; r<3; ++r) { + node_cells_map.insert(std::make_pair(cell_nodes(j,r),j)); + } + } + // <<<<------------------------------ + + std::vector<unsigned int> node_ids; + node_ids.reserve(node_cells_map.size()); + for (const auto& node_cell: node_cells_map) { + node_ids.push_back(node_cell.first); + } + std::vector<unsigned int>::iterator last_unique + = std::unique(node_ids.begin(), node_ids.end()); + node_ids.resize(std::distance(node_ids.begin(), last_unique)); + + m_number_of_nodes = node_ids.size(); + + std::cout << "node_ids.size()=" << node_ids.size() << '\n'; + if ((node_ids[0] != 0) or (node_ids[node_ids.size()-1] != node_ids.size()-1)) { + std::cerr << "sparse node numerotation NIY\n"; + std::exit(0); + } + + std::vector<std::vector<unsigned int>> node_cells_vector(node_ids.size()); + for (const auto& node_cell: node_cells_map) { + node_cells_vector[node_cell.first].push_back(node_cell.second); + } + + + Kokkos::View<unsigned short*> node_nb_cells("node_nb_cells", node_ids.size()); + int max_node_cells = 0; + for (int i=0; i<node_cells_vector.size(); ++i) { + const auto& cells_vector = node_cells_vector[i]; + const size_t nb_cells = cells_vector.size(); + node_nb_cells[i] = nb_cells; + if (nb_cells > max_node_cells) { + max_node_cells = nb_cells; + } + } + m_node_nb_cells = node_nb_cells; + Kokkos::View<unsigned int**> node_cells("node_cells", node_ids.size(), max_node_cells); + for (size_t i=0; i<node_cells_vector.size(); ++i) { + const auto& cells_vector = node_cells_vector[i]; + for (size_t j=0; j<cells_vector.size(); ++j) { + node_cells(i,j) = cells_vector[j]; + } + } + m_node_cells = node_cells; + + Kokkos::View<unsigned short**> node_cell_local_node("node_cell_local_node", + node_ids.size(), max_node_cells); + Kokkos::parallel_for(m_number_of_nodes, KOKKOS_LAMBDA(const int& r){ + for (int J=0; J<node_nb_cells[r]; ++J) { + const int j = node_cells(r,J); + for (int R=0; R<cell_nb_nodes[j]; ++R) { + if (cell_nodes(j,R) == r) { + node_cell_local_node(r,J)=R; + break; + } + } + } + }); + + m_node_cell_local_node = node_cell_local_node; } ~Connectivity2D() diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp index 140b9c6594caabb0154646799179157b50b27aaa..8090134e4b0263a0ccd11a04eb4db3fd92b0a0b7 100644 --- a/src/mesh/GmshReader.cpp +++ b/src/mesh/GmshReader.cpp @@ -761,30 +761,38 @@ GmshReader::__proceedData() std::exit(0); } else if ((dimension2_mask, elementNumber)>0) { if ((dimension1_mask, elementNumber)>0) { - throw ErrorHandler(__FILE__,__LINE__,"edges are not treated", ErrorHandler::normal); + // throw ErrorHandler(__FILE__,__LINE__,"edges are not treated", ErrorHandler::normal); } - if (__quadrangles.dimension_0()>0) { + if (__quadrangles.extent(0)>0) { std::cerr << "quadrangles NYI\n"; std::exit(0); } - Kokkos::View<unsigned int**> cell_nodes("cell_nodes", __triangles.dimension_0(), 3); - for (int j=0; j<__triangles.dimension_0(); ++j) { + + const size_t nb_cells = (dimension2_mask, elementNumber); + + std::cout << "nb_cells=" << nb_cells << '\n'; + const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, 3); + for (int j=0; j<nb_cells; ++j) { cell_nodes(j,0) = __triangles[j][0]; cell_nodes(j,1) = __triangles[j][1]; cell_nodes(j,2) = __triangles[j][2]; } - Connectivity2D connectivity(cell_nodes); - + const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells); + for (int j=0; j<nb_cells; ++j) { + cell_nb_nodes[j] = 3; + } + m_connectivity = new Connectivity2D(cell_nb_nodes, cell_nodes); + Connectivity2D& connectivity = *m_connectivity; typedef Mesh<Connectivity2D> MeshType; typedef TinyVector<2, double> Rd; - Kokkos::View<Rd*> xr("xr", __vertices.dimension_0()); - for (size_t i=0; i<__vertices.dimension_0(); ++i) { + Kokkos::View<Rd*> xr("xr", __vertices.extent(0)); + for (size_t i=0; i<__vertices.extent(0); ++i) { xr[i][0] = __vertices[i][0]; xr[i][1] = __vertices[i][1]; } - MeshType mesh(connectivity, xr); - + m_mesh = new MeshType(connectivity, xr); + MeshType& mesh = *m_mesh; std::ofstream gnuplot("mesh.gnu"); for (int j=0; j<mesh.numberOfCells(); ++j) { @@ -796,10 +804,10 @@ GmshReader::__proceedData() gnuplot << x[0] << ' ' << x[1] << "\n\n"; } gnuplot.close(); - std::exit(0); + } else if ((dimension1_mask, elementNumber)>0) { - std::cerr << "*** using a 1d mesh (NIY)\n"; - std::exit(0); + // std::cerr << "*** using a 1d mesh (NIY)\n"; + // std::exit(0); } else { std::cerr << "*** using a dimension 0 mesh is forbidden!\n"; std::exit(0); diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp index d916ad81ef2b611a47be57f2757050598c76efbe..5aec14118d6f7d08466e091f2c66201c375f480b 100644 --- a/src/mesh/GmshReader.hpp +++ b/src/mesh/GmshReader.hpp @@ -9,6 +9,9 @@ #include <Kokkos_Core.hpp> #include <TinyVector.hpp> +#include <Connectivity2D.hpp> +#include <Mesh.hpp> + class GmshReader { public: @@ -182,9 +185,26 @@ private: // void __writeReferences(const std::set<size_t>& references, // std::string objectName); + + + Connectivity2D* m_connectivity; + Mesh<Connectivity2D>* m_mesh; public: + + // Connectivity2D& connectivity() + // { + // return *m_connectivity; + // } + + Mesh<Connectivity2D>& mesh() + { + return *m_mesh; + } + GmshReader(const std::string& filename); ~GmshReader() = default; + + }; #endif // GMSH_READER_HPP diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp index 7fc8f8c806798747eae0190558c35b75beb6d232..02b5c0a39c3883cdca519cfcb7c771fcaff37172 100644 --- a/src/mesh/MeshData.hpp +++ b/src/mesh/MeshData.hpp @@ -58,7 +58,6 @@ private: { 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(); @@ -75,10 +74,43 @@ private: KOKKOS_INLINE_FUNCTION void _updateCjr() { - if(dimension == 1) { + if constexpr (dimension == 1) { // Cjr/njr/ljr are constant overtime + } + else if (dimension == 2) { + 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){ + for (int R=0; R<cell_nb_nodes[j]; ++R) { + int Rp1 = (R+1)%cell_nb_nodes[j]; + int Rm1 = (R+cell_nb_nodes[j]-1)%cell_nb_nodes[j]; + Rd half_xrp_xrm = 0.5*(xr(cell_nodes(j,Rp1))-xr(cell_nodes(j,Rm1))); + m_Cjr(j,R) = Rd{-half_xrp_xrm[1], half_xrp_xrm[0]}; + } + }); + + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + for (int R=0; R<cell_nb_nodes[j]; ++R) { + const Rd& Cjr = m_Cjr(j,R); + m_ljr(j,R) = sqrt((Cjr,Cjr)); + } + }); + + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + for (int R=0; R<cell_nb_nodes[j]; ++R) { + const Rd& Cjr = m_Cjr(j,R); + const double inv_ljr = 1./m_ljr(j,R); + m_njr(j,R) = inv_ljr*Cjr; + } + }); + } - static_assert(dimension==1, "only 1d is implemented"); + static_assert((dimension==1) or (dimension==2), "only 1d and 2d are implemented"); } public: @@ -114,9 +146,9 @@ public: void updateAllData() { + this->_updateCjr(); this->_updateCenter(); this->_updateVolume(); - this->_updateCjr(); } MeshData(const MeshType& mesh) @@ -127,7 +159,7 @@ public: m_xj("xj", mesh.numberOfCells()), m_Vj("Vj", mesh.numberOfCells()) { - if (dimension==1) { + if constexpr (dimension==1) { // in 1d Cjr are computed once for all Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { m_Cjr(j,0)=-1; diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index d5d64d77bd3ea00cd1e01dffd6304df5b9db3cef..3d2f79fe7391b8792b38f0a1fef7e77f52d611d4 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -228,9 +228,20 @@ private: void inverse(const Kokkos::View<const Rdd*>& A, Kokkos::View<Rdd*>& inv_A) const { - Kokkos::parallel_for(A.size(), KOKKOS_LAMBDA(const int& r) { - inv_A(r) = Rdd{1./(A(r)(0,0))}; - }); + if constexpr (dimension==1) { + Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) { + inv_A(r) = Rdd{1./(A(r)(0,0))}; + }); + } + else { + Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) { + const Rdd& M = A(r); + double det = M(0,0)*M(1,1)-M(1,0)*M(0,1); + double inv_det = 1./det; + inv_A(r) = inv_det*Rdd(M(1,1),-M(0,1), + -M(1,0), M(0,0)); + }); + } } void inverse(const Kokkos::View<const double*>& x,