diff --git a/src/mesh/Connectivity1D.hpp b/src/mesh/Connectivity1D.hpp
index fdebf6efac58d521889b656763a4f8a723bb1cb4..9b9e11842ed2871d81add5fc2777ac3a52e01cb9 100644
--- a/src/mesh/Connectivity1D.hpp
+++ b/src/mesh/Connectivity1D.hpp
@@ -20,6 +20,12 @@ public:
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;
+
private:
std::vector<RefNodeList> m_ref_node_list;
@@ -97,6 +103,23 @@ public:
});
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_max_nb_node_per_cell,
diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp
index 2ca5ad40fe650279cabef934c7f3e0ff32824b5a..99d313fd15b0d4164a1bfc75b0bd919d23867b3c 100644
--- a/src/mesh/Connectivity2D.hpp
+++ b/src/mesh/Connectivity2D.hpp
@@ -27,6 +27,12 @@ class Connectivity2D
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;
+
private:
std::vector<RefFaceList> m_ref_face_list;
std::vector<RefNodeList> m_ref_node_list;
@@ -236,6 +242,22 @@ class Connectivity2D
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_max_nb_node_per_cell,
diff --git a/src/mesh/Connectivity3D.hpp b/src/mesh/Connectivity3D.hpp
index 6b70560d86d73cbec5a9bfdd13673d6766e80a83..9e6afd7b279de1c88b659e80493d02415344c5a9 100644
--- a/src/mesh/Connectivity3D.hpp
+++ b/src/mesh/Connectivity3D.hpp
@@ -34,6 +34,12 @@ public:
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;
+
private:
std::vector<RefFaceList> m_ref_face_list;
std::vector<RefNodeList> m_ref_node_list;
@@ -479,6 +485,22 @@ private:
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",
+ node_id_per_cell_vector);
+ }
+
ConnectivityUtils utils;
utils.computeNodeCellConnectivity(m_cell_to_node_matrix,
m_max_nb_node_per_cell,
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 0ce9beef810467a35466247c04a1abccb233aa1e..bf8a8a3c5a09f5fc208143bff151221d5bd9cd5b 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -14,6 +14,8 @@
#include <FiniteVolumesEulerUnknowns.hpp>
#include <BoundaryCondition.hpp>
+#include <NodeByCellData.hpp>
+
template<typename MeshData>
class AcousticSolver
{
@@ -84,14 +86,86 @@ private:
const Kokkos::View<const double**>& ljr,
const Kokkos::View<const Rd**>& njr)
{
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ int cpt0 = 0;
+ Kokkos::parallel_for("T new Ajr (raw vector)", m_Ajr_new.numberOfValues(), [&](const int& j) {
+ m_Ajr_new.value(j) = tensorProduct(Rd(zero), Rd(zero));
+ cpt0 ++;
+ });
+
+
+ int cpt1 = 0;
+ Kokkos::parallel_for("T new Ajr", m_mesh.numberOfCells(), [&](const int& j) {
+ auto Aj = m_Ajr_new.myCellValues(j);
+ const size_t nb_nodes = Aj.size();
+ for (size_t r=0; r<nb_nodes; ++r) {
+ Aj[r] = tensorProduct(Rd(zero), Rd(zero));
+ cpt1++;
+ }
+ });
+
+ int cpt2 = 0;
+ Kokkos::parallel_for("T old Ajr", m_mesh.numberOfCells(), [&](const int& j) {
const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
+ for (size_t r=0; r<cell_nodes.length; ++r) {
+ m_Ajr(j,r) = tensorProduct(Rd(zero), Rd(zero));
+ cpt2++;
+ }
+ });
+ int cpt3 = 0;
+ Kokkos::parallel_for("new Ajr", m_mesh.numberOfCells(), [&](const int& j) {
+ auto Aj = m_Ajr_new.myCellValues(j);
+ const size_t nb_nodes = Aj.size();
+ for (size_t r=0; r<nb_nodes; ++r) {
+ Aj[r] = tensorProduct((rhocj(j)*ljr(j,r))*njr(j,r), njr(j,r));
+ cpt3++;
+ }
+ });
+
+ int cpt4 = 0;
+ Kokkos::parallel_for("old Ajr", m_mesh.numberOfCells(), [&](const int& j) {
+ const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
for (size_t r=0; r<cell_nodes.length; ++r) {
m_Ajr(j,r) = tensorProduct((rhocj(j)*ljr(j,r))*njr(j,r), njr(j,r));
+ cpt4++;
}
});
+ int cpt5 = 0;
+ Kokkos::View<Rdd*> simple_view("simpleview", m_Ajr_new.numberOfValues());
+ Kokkos::parallel_for("T new Ajr (simple View)", simple_view.extent(0), [&](const int& j) {
+ simple_view(j) = tensorProduct(Rd(zero), Rd(zero));
+ cpt5 ++;
+ });
+
+ Kokkos::View<unsigned int*> shift ("shift", m_mesh.numberOfCells());
+ shift(0)=0;
+ for (unsigned int j=1; j<m_mesh.numberOfCells(); ++j) {
+ shift[j] = shift[j-1]+m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j-1).length;
+ }
+ int cpt6 = 0;
+ Kokkos::View<Rdd*> simple_view2("simpleview2", m_Ajr_new.numberOfValues());
+ Kokkos::parallel_for("T new Ajr (simple View2)", m_mesh.numberOfCells(),
+ [&] (const int& j) {
+
+ for (unsigned r=0; r<m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j).length; ++r) {
+ simple_view2(shift[j]+r) = tensorProduct(Rd(zero), Rd(zero));
+ cpt6 ++;
+ }
+ });
+
+
+ std::cout << cpt0 << ' ' << cpt1 << ' ' << cpt2 << ' ' << cpt3 << ' ' << cpt4 << ' ' << cpt5 << ' ' << cpt6 << '\n';
+ // for (unsigned int j=0; j<m_mesh.numberOfCells(); ++j) {
+ // const auto& Aj = m_Ajr_new.cellValues(j);
+ // const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
+ // for (size_t r=0; r<<cell_nodes.length; ++r) {
+ // if (Aj[r] != m_Ajr(j,r)) {
+ // std::cerr << "viens !\n";
+ // }
+ // }
+ // }
+
return m_Ajr;
}
@@ -242,6 +316,7 @@ private:
}
Kokkos::View<Rd*> m_br;
+ NodeByCellData<Rdd> m_Ajr_new;
Kokkos::View<Rdd**> m_Ajr;
Kokkos::View<Rdd*> m_Ar;
Kokkos::View<Rdd*> m_inv_Ar;
@@ -259,6 +334,7 @@ public:
m_connectivity(m_mesh.connectivity()),
m_boundary_condition_list(bc_list),
m_br("br", m_mesh.numberOfNodes()),
+ m_Ajr_new(m_connectivity.m_node_id_per_cell_matrix),
m_Ajr("Ajr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
m_Ar("Ar", m_mesh.numberOfNodes()),
m_inv_Ar("inv_Ar", m_mesh.numberOfNodes()),
diff --git a/src/scheme/NodeByCellData.hpp b/src/scheme/NodeByCellData.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..bd8718092f88f826f41ed2ded895c94d28919376
--- /dev/null
+++ b/src/scheme/NodeByCellData.hpp
@@ -0,0 +1,131 @@
+#ifndef NODE_BY_CELL_DATA_HPP
+#define NODE_BY_CELL_DATA_HPP
+
+#include <Kokkos_StaticCrsGraph.hpp>
+
+typedef Kokkos::StaticCrsGraph<unsigned int, Kokkos::HostSpace> ConnectivityMatrix;
+
+template <typename DataType>
+class NodeByCellData
+{
+#warning should eventually filter const from DataType
+ private:
+ ConnectivityMatrix m_node_id_per_cell_matrix;
+
+ Kokkos::View<unsigned int*> m_begin;
+ Kokkos::View<unsigned int*> m_end;
+
+ Kokkos::View<DataType*> m_values;
+
+ public:
+ class SubView
+ {
+ private:
+ DataType* const m_values;
+ const size_t m_size;
+ public:
+ KOKKOS_FORCEINLINE_FUNCTION
+ const DataType& operator[](const size_t& i) const
+ {
+ return m_values[i];
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator[](const size_t& i)
+ {
+ return m_values[i];
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ const size_t& size() const
+ {
+ return m_size;
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ SubView(const SubView&) = delete;
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ SubView(SubView&&) = default;
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ SubView(Kokkos::View<DataType*> values,
+ const unsigned int& begin,
+ const unsigned int& end)
+ : m_values(&(values[begin])),
+ m_size(end-begin)
+ {
+
+ }
+ };
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ size_t numberOfValues() const
+ {
+ return m_values.extent(0);
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& value(const size_t& i)
+ {
+ return m_values[i];
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ const DataType& value(const size_t & i) const
+ {
+ return m_values[i];
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ size_t numberOfEntities() const
+ {
+ return m_node_id_per_cell_matrix.numRows();
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ SubView&& myCellValues(const size_t& i_cell)
+ {
+ ConnectivityMatrix::size_type cell_begin = m_begin[i_cell];
+ ConnectivityMatrix::size_type cell_end = m_begin[i_cell+1];
+ return std::move(SubView(m_values, cell_begin, cell_end));
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ Kokkos::View<DataType*>
+ cellValues(const size_t& i_cell)
+ {
+ ConnectivityMatrix::size_type cell_begin = m_begin[i_cell];
+ ConnectivityMatrix::size_type cell_end = m_begin[i_cell+1];
+ // ConnectivityMatrix::size_type cell_begin = m_node_id_per_cell_matrix.row_map(i_cell);
+ // ConnectivityMatrix::size_type cell_end = m_node_id_per_cell_matrix.row_map(i_cell+1);
+ return Kokkos::subview(m_values, std::make_pair(cell_begin, cell_end));
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ Kokkos::View<const DataType*>
+ cellValues(const size_t& i_cell) const
+ {
+ ConnectivityMatrix::size_type cell_begin = m_begin[i_cell];
+ ConnectivityMatrix::size_type cell_end = m_end[i_cell+1];
+ // ConnectivityMatrix::size_type cell_begin = m_node_id_per_cell_matrix.row_map(i_cell);
+ // ConnectivityMatrix::size_type cell_end = m_node_id_per_cell_matrix.row_map(i_cell+1);
+ return Kokkos::subview(m_values, std::make_pair(cell_begin, cell_end));
+ }
+
+ NodeByCellData(const ConnectivityMatrix& node_id_per_cell_matrix)
+ : m_node_id_per_cell_matrix(node_id_per_cell_matrix),
+ m_begin("begin", m_node_id_per_cell_matrix.numRows()+1),
+ // m_end("end", m_node_id_per_cell_matrix.numRows()),
+ m_values("values", m_node_id_per_cell_matrix.entries.extent(0))
+ {
+ for (unsigned int j=0; j<m_node_id_per_cell_matrix.numRows()+1; ++j) {
+ m_begin[j] = m_node_id_per_cell_matrix.row_map(j);
+ // m_end[j] = m_node_id_per_cell_matrix.row_map(j+1);
+ }
+ }
+
+ ~NodeByCellData() = default;
+};
+
+#endif // NODE_BY_CELL_DATA_HPP