diff --git a/src/main.cpp b/src/main.cpp index 34ad22fea1716fa49262894add80484df6662583..7420998163d5d2b2923c47a8732e360d904aaeee 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -11,7 +11,6 @@ // #include <AcousticSolverClass.hpp> // #include <AcousticSolverTest.hpp> -#include <Connectivity1D.hpp> #include <Connectivity3D.hpp> #include <Mesh.hpp> @@ -441,89 +440,8 @@ int main(int argc, char *argv[]) std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n"; } else { - // class for acoustic solver test - Kokkos::Timer timer; - timer.reset(); - std::shared_ptr<Connectivity1D>connectivity( new Connectivity1D(number)); - typedef Mesh<Connectivity1D> MeshType; - typedef MeshData<MeshType> MeshDataType; - typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType; - - MeshType mesh(connectivity); - 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_bc0 - = new PressureBoundaryCondition(1, - std::vector<unsigned int>({static_cast<unsigned int>(0u)})); - std::shared_ptr<PressureBoundaryCondition> bc0(pres_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 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::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 - const Kokkos::View<const Rd*> xj = mesh_data.xj(); - const Kokkos::View<const double*> rhoj = unknowns.rhoj(); - std::ofstream fout("rho"); - for (size_t j=0; j<mesh.numberOfCells(); ++j) { - fout << xj[j][0] << ' ' << rhoj[j] << '\n'; - } - } - + std::cerr << "Connectivity1D defined by number of nodes no more implemented\n"; + std::exit(0); } } catch (const AssertError& error) { diff --git a/src/mesh/Connectivity1D.hpp b/src/mesh/Connectivity1D.hpp deleted file mode 100644 index 51ca51a760feb20cdd72d5aa02b8fbde541d46de..0000000000000000000000000000000000000000 --- a/src/mesh/Connectivity1D.hpp +++ /dev/null @@ -1,134 +0,0 @@ -#ifndef CONNECTIVITY_1D_HPP -#define CONNECTIVITY_1D_HPP - -#include <Kokkos_Core.hpp> -#include <PastisAssert.hpp> - -#include <TinyVector.hpp> -#include <ConnectivityUtils.hpp> - -#include <TypeOfItem.hpp> -#include <RefId.hpp> -#include <RefNodeList.hpp> - -class Connectivity1D -{ -public: - static constexpr size_t dimension = 1; - - ConnectivityMatrix m_cell_to_node_matrix; - - ConnectivityMatrix m_node_to_cell_matrix; - ConnectivityMatrixShort m_node_to_cell_local_node_matrix; - - // Stores numbering of nodes of each cell. - // - // This is different from m_cell_to_node_matrix which return the global id of - // a local node in a cell - ConnectivityMatrix m_node_id_per_cell_matrix; - - ConnectivityMatrix subItemIdPerItemMatrix() const - { - return m_node_id_per_cell_matrix; - } - - private: - std::vector<RefNodeList> m_ref_node_list; - - Kokkos::View<double*> m_inv_cell_nb_nodes; - - std::vector<std::vector<unsigned int>> - _buildConnectivity(const size_t& number_of_cells) - { - std::vector<std::vector<unsigned int>> cell_by_node_vector(number_of_cells); - for (unsigned int j=0; j<number_of_cells; ++j) { - cell_by_node_vector[j] = {j, j+1}; - } - return cell_by_node_vector; - } -public: - void addRefNodeList(const RefNodeList& ref_node_list) - { - m_ref_node_list.push_back(ref_node_list); - } - - size_t numberOfRefNodeList() const - { - return m_ref_node_list.size(); - } - - const RefNodeList& refNodeList(const size_t& i) const - { - return m_ref_node_list[i]; - } - - KOKKOS_INLINE_FUNCTION - size_t numberOfNodes() const - { - return m_node_to_cell_matrix.numRows(); - } - - KOKKOS_INLINE_FUNCTION - size_t numberOfCells() const - { - return m_cell_to_node_matrix.numRows(); - } - - const Kokkos::View<const double*> invCellNbNodes() const - { - return m_inv_cell_nb_nodes; - } - - Connectivity1D(const Connectivity1D&) = delete; - - Connectivity1D(const size_t& number_of_cells) - : Connectivity1D(_buildConnectivity(number_of_cells)) - { - ; - } - - Connectivity1D(const std::vector<std::vector<unsigned int>>& cell_by_node_vector) - { - m_cell_to_node_matrix - = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("cell_to_node_matrix", - cell_by_node_vector); - - Assert(this->numberOfCells()>0); - { - Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", this->numberOfCells()); - Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const size_t& j) { - const auto& cell_nodes = m_cell_to_node_matrix.rowConst(j); - inv_cell_nb_nodes[j] = 1./cell_nodes.length; - }); - m_inv_cell_nb_nodes = inv_cell_nb_nodes; - } - - { - std::vector<std::vector<unsigned int>> node_id_per_cell_vector(this->numberOfCells()); - unsigned int id=0; - for (unsigned int j=0; j<this->numberOfCells(); ++j) { - const auto& cell_to_node = m_cell_to_node_matrix.rowConst(j); - auto& node_id_per_cell = node_id_per_cell_vector[j]; - node_id_per_cell.resize(cell_to_node.length); - for (size_t r=0; r<cell_to_node.length; ++r) { - node_id_per_cell[r] = id++; - } - } - m_node_id_per_cell_matrix - = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("node_id_per_cell_matrix", - cell_by_node_vector); - } - - ConnectivityUtils utils; - utils.computeNodeCellConnectivity(m_cell_to_node_matrix, - m_node_to_cell_matrix, - m_node_to_cell_local_node_matrix); - } - - ~Connectivity1D() - { - ; - } -}; - -#endif // CONNECTIVITY_1D_HPP diff --git a/src/mesh/Connectivity3D.hpp b/src/mesh/Connectivity3D.hpp index 10f8f2d1b6b74fb284573bc280dd9176a4c462ec..810be687128aa3f2fc86d592668bba634d2760b7 100644 --- a/src/mesh/Connectivity3D.hpp +++ b/src/mesh/Connectivity3D.hpp @@ -24,6 +24,18 @@ class Connectivity; template <size_t Dimension> class ConnectivityFace; + +template<> +class ConnectivityFace<1> +{ + public: + friend struct Hash; + struct Hash + { + size_t operator()(const ConnectivityFace& f) const; + }; +}; + template<> class ConnectivityFace<2> { @@ -364,7 +376,9 @@ private: m_node_to_cell_matrix, m_node_to_cell_local_node_matrix); - this->_computeFaceCellConnectivities(); + if constexpr (Dimension>1) { + this->_computeFaceCellConnectivities(); + } } ~Connectivity() @@ -663,4 +677,15 @@ Connectivity<2>::subItemIdPerItemMatrix<TypeOfItem::node, return m_node_id_per_cell_matrix; } +using Connectivity1D = Connectivity<1>; + +template <> +template <> +inline const ConnectivityMatrix& +Connectivity<1>::subItemIdPerItemMatrix<TypeOfItem::node, + TypeOfItem::cell>() const +{ + return m_node_id_per_cell_matrix; +} + #endif // CONNECTIVITY_3D_HPP diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp index 6e422a5e8381336623a8fdc3a891c33dcb0d173e..573ae0227107cefe62e4ee1b35b66fb4153adce4 100644 --- a/src/mesh/GmshReader.cpp +++ b/src/mesh/GmshReader.cpp @@ -5,7 +5,6 @@ #include <set> #include <rang.hpp> -#include <Connectivity1D.hpp> #include <Connectivity3D.hpp> #include <Mesh.hpp> diff --git a/src/scheme/SubItemValuePerItem.hpp b/src/scheme/SubItemValuePerItem.hpp index 80acb50a8e17119d21f6618414984c0f826bd1bf..202ffdf2085b60f859f0dc11320569721ac44f65 100644 --- a/src/scheme/SubItemValuePerItem.hpp +++ b/src/scheme/SubItemValuePerItem.hpp @@ -164,11 +164,7 @@ class SubItemValuePerItem template <typename ConnectivityType> SubItemValuePerItem(const ConnectivityType& connectivity) { - if constexpr (ConnectivityType::dimension > 1) { - m_subitem_id_per_item_matrix = connectivity.template subItemIdPerItemMatrix<SubItemType,ItemType>(); - } else { - m_subitem_id_per_item_matrix = connectivity.subItemIdPerItemMatrix(); - } + m_subitem_id_per_item_matrix = connectivity.template subItemIdPerItemMatrix<SubItemType,ItemType>(); m_values = Kokkos::View<DataType*>("values", m_subitem_id_per_item_matrix.entries.extent(0)); }