diff --git a/CMakeLists.txt b/CMakeLists.txt
index 341828ab792657d2c2d45a5da6bc1b4160f06a36..e369d66949e771f78d44b82d0ff7cc0cd23e400e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -136,8 +136,6 @@ include_directories(src/output)
include_directories(src/utils)
include_directories(src/scheme)
-include_directories(src/experimental)
-
# Pastis tests
set(CATCH_MODULE_PATH "${PASTIS_SOURCE_DIR}/packages/Catch2")
@@ -167,7 +165,6 @@ if("${CMAKE_BUILD_TYPE}" STREQUAL "Coverage")
endif()
set (GCOVR_EXCLUDE
- -e "${PASTIS_SOURCE_DIR}/src/experimental"
-e "${PASTIS_SOURCE_DIR}/src/main.cpp"
-e "${PASTIS_SOURCE_DIR}/src/utils/BacktraceManager.cpp"
-e "${PASTIS_SOURCE_DIR}/src/utils/BacktraceManager.hpp"
@@ -216,5 +213,4 @@ target_link_libraries(
pastis
kokkos
PastisUtils
- PastisMesh
- PastisExperimental)
+ PastisMesh)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index fcc3e886e01edf90d0aca9447099d64056d2b91c..ebf2a69d275051fcd7e0bdfd6a0a0718ae8e37ea 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -21,7 +21,3 @@ include_directories(scheme)
# Pastis output
include_directories(output)
-
-# Pastis experimental
-add_subdirectory(experimental)
-include_directories(experimental)
diff --git a/src/experimental/AcousticSolverClass.cpp b/src/experimental/AcousticSolverClass.cpp
deleted file mode 100644
index 051e148bb21babf5eb65c22b2bcf820df37632b6..0000000000000000000000000000000000000000
--- a/src/experimental/AcousticSolverClass.cpp
+++ /dev/null
@@ -1,381 +0,0 @@
-#include <AcousticSolverClass.hpp>
-#include <rang.hpp>
-
-#include <BlockPerfectGas.hpp>
-
-typedef const double my_double;
-
-struct AcousticSolverClass::ReduceMin
-{
-private:
- const Kokkos::View<my_double*> x_;
-
-public:
- typedef Kokkos::View<my_double*>::non_const_value_type value_type;
-
- ReduceMin(const Kokkos::View<my_double*>& x) : x_ (x) {}
-
- typedef Kokkos::View<my_double*>::size_type size_type;
-
- KOKKOS_INLINE_FUNCTION void
- operator() (const size_type i, value_type& update) const
- {
- if (x_(i) < update) {
- update = x_(i);
- }
- }
-
- KOKKOS_INLINE_FUNCTION void
- join (volatile value_type& dst,
- const volatile value_type& src) const
- {
- if (src < dst) {
- dst = src;
- }
- }
-
- KOKKOS_INLINE_FUNCTION void
- init (value_type& dst) const
- { // The identity under max is -Inf.
- dst= Kokkos::reduction_identity<value_type>::min();
- }
-};
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*>
-AcousticSolverClass::computeRhoCj(const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& cj)
-{
- Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
- m_rhocj[j] = rhoj[j]*cj[j];
- });
- return m_rhocj;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*[2]>
-AcousticSolverClass::computeAjr(const Kokkos::View<const double*>& rhocj,
- const Kokkos::View<const double*[2]>& Cjr)
-{
- Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
- for (int r=0; r<2; ++r) {
- m_Ajr(j,r) = rhocj(j)*Cjr(j,r)*Cjr(j,r);
- }
- });
-
- return m_Ajr;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*>
-AcousticSolverClass::computeAr(const Kokkos::View<const double*[2]>& Ajr,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- const Kokkos::View<const int*>& node_nb_cells)
-{
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
- double sum = 0;
- for (int j=0; j<node_nb_cells(r); ++j) {
- const int J = node_cells(r,j);
- const int R = node_cell_local_node(r,j);
- sum += Ajr(J,R);
- }
- m_Ar(r) = sum;
- });
-
- return m_Ar;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*>
-AcousticSolverClass::computeBr(const Kokkos::View<const double*[2]>& Ajr,
- const Kokkos::View<const double*[2]>& Cjr,
- const Kokkos::View<const double*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- const Kokkos::View<const int*>& node_nb_cells)
-{
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
- double sum = 0;
- for (int j=0; j<node_nb_cells(r); ++j) {
- const int J = node_cells(r,j);
- const int R = node_cell_local_node(r,j);
- sum += Ajr(J,R)*uj(J) + Cjr(J,R)*pj(J);
- }
- m_br(r) = sum;
- });
-
- return m_br;
-}
-
-KOKKOS_INLINE_FUNCTION
-Kokkos::View<double*>
-AcousticSolverClass::computeUr(const Kokkos::View<const double*>& Ar,
- const Kokkos::View<const double*>& br)
-{
- inverse(Ar, m_inv_Ar);
- const Kokkos::View<const double*> invAr = m_inv_Ar;
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
- m_ur[r]=br(r)*invAr(r);
- });
- m_ur[0]=0;
- m_ur[m_nr-1]=0;
-
- return m_ur;
-}
-
-KOKKOS_INLINE_FUNCTION
-Kokkos::View<double*[2]>
-AcousticSolverClass::computeFjr(const Kokkos::View<const double*[2]>& Ajr,
- const Kokkos::View<const double*>& ur,
- const Kokkos::View<const double*[2]>& Cjr,
- const Kokkos::View<const double*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const int*[2]>& cell_nodes)
-{
- Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
- for (int r=0; r<2; ++r) {
- m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+Cjr(j,r)*pj(j);
- }
- });
-
- return m_Fjr;
-}
-
-KOKKOS_INLINE_FUNCTION
-double AcousticSolverClass::
-acoustic_dt(const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*>& cj) const
-{
- Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){
- m_Vj_over_cj[j] = Vj[j]/cj[j];
- });
-
- double dt = std::numeric_limits<double>::max();
- Kokkos::parallel_reduce(m_nj, ReduceMin(m_Vj_over_cj), dt);
-
- return dt;
-}
-
-KOKKOS_INLINE_FUNCTION
-void
-AcousticSolverClass::inverse(const Kokkos::View<const double*>& x,
- Kokkos::View<double*>& inv_x) const
-{
- Kokkos::parallel_for(x.size(), KOKKOS_LAMBDA(const int& r) {
- inv_x(r) = 1./x(r);
- });
-}
-
-KOKKOS_INLINE_FUNCTION
-void AcousticSolverClass::computeExplicitFluxes(const Kokkos::View<const double*>& xr,
- const Kokkos::View<const double*>& xj,
- const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const double*>& cj,
- const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*[2]>& Cjr,
- const Kokkos::View<const int*[2]>& cell_nodes,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*>& node_nb_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- Kokkos::View<double*>& ur,
- Kokkos::View<double*[2]>& Fjr)
-{
- const Kokkos::View<const double*> rhocj = computeRhoCj(rhoj, cj);
- const Kokkos::View<const double*[2]> Ajr = computeAjr(rhocj, Cjr);
-
- const Kokkos::View<const double*> Ar = computeAr(Ajr, node_cells, node_cell_local_node, node_nb_cells);
- const Kokkos::View<const double*> br = computeBr(Ajr, Cjr, uj, pj,
- node_cells, node_cell_local_node, node_nb_cells);
-
- ur = computeUr(Ar, br);
- Fjr = computeFjr(Ajr, ur, Cjr, uj, pj, cell_nodes);
-}
-
-AcousticSolverClass::AcousticSolverClass(const long int& nj)
- : m_nj(nj),
- m_nr(nj+1),
- m_br("br", m_nr),
- m_Ajr("Ajr", m_nj),
- m_Ar("Ar", m_nr),
- m_inv_Ar("inv_Ar", m_nr),
- m_Fjr("Fjr", m_nj),
- m_ur("ur", m_nr),
- m_rhocj("rho_c", m_nj),
- m_Vj_over_cj("Vj_over_cj", m_nj)
-{
- Kokkos::View<double*> xj("xj",m_nj);
- Kokkos::View<double*> rhoj("rhoj",m_nj);
-
- Kokkos::View<double*> uj("uj",m_nj);
-
- 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<double*> xr("xr", m_nr);
-
- const double delta_x = 1./m_nj;
-
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
- xr[r] = r*delta_x;
- });
-
- Kokkos::View<int*[2]> cell_nodes("cell_nodes",m_nj,2);
- Kokkos::View<int*[2]> node_cells("node_cells",m_nr,2);
- Kokkos::View<int*[2]> node_cell_local_node("node_cells",m_nr,2);
- Kokkos::View<int*> node_nb_cells("node_cells",m_nr);
-
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
- node_nb_cells(r) = 2;
- });
- node_nb_cells(0) = 1;
- node_nb_cells(m_nr-1) = 1;
-
- node_cells(0,0) = 0;
- Kokkos::parallel_for(m_nr-2, KOKKOS_LAMBDA(const int& r){
- node_cells(r+1,0) = r;
- node_cells(r+1,1) = r+1;
- });
- node_cells(m_nr-1,0) =m_nj-1;
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- cell_nodes(j,0) = j;
- cell_nodes(j,1) = j+1;
- });
-
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
- for (int J=0; J<node_nb_cells(r); ++J) {
- int j = node_cells(r,J);
- for (int R=0; R<2; ++R) {
- if (cell_nodes(j,R) == r) {
- node_cell_local_node(r,J) = R;
- }
- }
- }
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]);
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Vj[j] = xr[cell_nodes(j,1)]-xr[cell_nodes(j,0)];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- if (xj[j]<0.5) {
- rhoj[j]=1;
- } else {
- rhoj[j]=0.125;
- }
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- if (xj[j]<0.5) {
- pj[j]=1;
- } else {
- pj[j]=0.1;
- }
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- gammaj[j] = 1.4;
- });
-
- BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-
- block_eos.updateEandCFromRhoP();
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Ej[j] = ej[j]+0.5*uj[j]*uj[j];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- mj[j] = rhoj[j] * Vj[j];
- });
-
- Kokkos::View<double*> inv_mj("inv_mj",m_nj);
- inverse(mj, inv_mj);
-
- const double tmax=0.2;
- double t=0;
-
- int itermax=std::numeric_limits<int>::max();
- int iteration=0;
-
- Kokkos::View<double*[2]> Cjr("Cjr",m_nj);
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
- Cjr(j,0)=-1;
- Cjr(j,1)= 1;
- });
-
- while((t<tmax) and (iteration<itermax)) {
- double dt = 0.4*acoustic_dt(Vj, cj);
- if (t+dt<tmax) {
- t+=dt;
- } else {
- dt=tmax-t;
- t=tmax;
- }
-
- computeExplicitFluxes(xr, xj,
- rhoj, uj, pj, cj, Vj, Cjr,
- cell_nodes, node_cells, node_nb_cells, node_cell_local_node,
- m_ur, m_Fjr);
-
- const Kokkos::View<const double*[2]> Fjr = m_Fjr;
- const Kokkos::View<const double*> ur = m_ur;
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- double momentum_fluxes = 0;
- double energy_fluxes = 0;
- for (int R=0; R<2; ++R) {
- const int r=cell_nodes(j,R);
- momentum_fluxes += Fjr(j,R);
- energy_fluxes += Fjr(j,R)*ur[r];
- }
- uj[j] -= dt*inv_mj[j]*(momentum_fluxes);
- Ej[j] -= dt*inv_mj[j]*(energy_fluxes);
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
- ej[j] = Ej[j] - 0.5 * uj[j]*uj[j];
- });
-
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
- xr[r] += dt*ur[r];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- xj[j] = 0.5*(xr[j]+xr[j+1]);
- Vj[j] = xr[j+1]-xr[j];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- rhoj[j] = mj[j]/Vj[j];
- });
-
- block_eos.updatePandCFromRhoE();
-
- ++iteration;
- }
-
- // {
- // std::ofstream fout("rho");
- // for (int j=0; j<nj; ++j) {
- // fout << xj[j] << ' ' << rhoj[j] << '\n';
- // }
- // }
-
- std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
- << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
-}
diff --git a/src/experimental/AcousticSolverClass.hpp b/src/experimental/AcousticSolverClass.hpp
deleted file mode 100644
index 2c2f8c0eccad97536b2e7f840aa01d8ccca058bf..0000000000000000000000000000000000000000
--- a/src/experimental/AcousticSolverClass.hpp
+++ /dev/null
@@ -1,84 +0,0 @@
-#ifndef ACOUSTIC_SOLVER_CLASS_HPP
-#define ACOUSTIC_SOLVER_CLASS_HPP
-
-#include <Kokkos_Core.hpp>
-
-class AcousticSolverClass
-{
-private:
- inline const Kokkos::View<const double*>
- computeRhoCj(const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& cj);
-
- inline const Kokkos::View<const double*[2]>
- computeAjr(const Kokkos::View<const double*>& rhocj,
- const Kokkos::View<const double*[2]>& Cjr);
-
- inline const Kokkos::View<const double*>
- computeAr(const Kokkos::View<const double*[2]>& Ajr,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- const Kokkos::View<const int*>& node_nb_cells);
-
- inline const Kokkos::View<const double*>
- computeBr(const Kokkos::View<const double*[2]>& Ajr,
- const Kokkos::View<const double*[2]>& Cjr,
- const Kokkos::View<const double*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- const Kokkos::View<const int*>& node_nb_cells);
-
- Kokkos::View<double*>
- computeUr(const Kokkos::View<const double*>& Ar,
- const Kokkos::View<const double*>& br);
-
- Kokkos::View<double*[2]>
- computeFjr(const Kokkos::View<const double*[2]>& Ajr,
- const Kokkos::View<const double*>& ur,
- const Kokkos::View<const double*[2]>& Cjr,
- const Kokkos::View<const double*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const int*[2]>& cell_nodes);
-
- void inverse(const Kokkos::View<const double*>& x,
- Kokkos::View<double*>& inv_x) const;
-
- inline double acoustic_dt(const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*>& cj) const;
-
- inline void computeExplicitFluxes(const Kokkos::View<const double*>& xr,
- const Kokkos::View<const double*>& xj,
- const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const double*>& cj,
- const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*[2]>& Cjr,
- const Kokkos::View<const int*[2]>& cell_nodes,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*>& node_nb_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- Kokkos::View<double*>& ur,
- Kokkos::View<double*[2]>& Fjr);
-
- struct ReduceMin;
-
- const long int m_nj;
- const long int m_nr;
-
- Kokkos::View<double*> m_br;
- Kokkos::View<double*[2]> m_Ajr;
- Kokkos::View<double*> m_Ar;
- Kokkos::View<double*> m_inv_Ar;
- Kokkos::View<double*[2]> m_Fjr;
- Kokkos::View<double*> m_ur;
- Kokkos::View<double*> m_rhocj;
- Kokkos::View<double*> m_Vj_over_cj;
-
-
-public:
- AcousticSolverClass(const long int& nj);
-};
-
-#endif // ACOUSTIC_SOLVER_CLASS_HPP
diff --git a/src/experimental/AcousticSolverTest.cpp b/src/experimental/AcousticSolverTest.cpp
deleted file mode 100644
index aaa9c0dbfefdbb603a95d657b81e23813ad5e69c..0000000000000000000000000000000000000000
--- a/src/experimental/AcousticSolverTest.cpp
+++ /dev/null
@@ -1,402 +0,0 @@
-#include <AcousticSolverTest.hpp>
-#include <rang.hpp>
-
-#include <BlockPerfectGas.hpp>
-#include <fstream>
-
-typedef const double my_double;
-
-struct AcousticSolverTest::ReduceMin
-{
-private:
- const Kokkos::View<my_double*> x_;
-
-public:
- typedef Kokkos::View<my_double*>::non_const_value_type value_type;
-
- ReduceMin(const Kokkos::View<my_double*>& x) : x_ (x) {}
-
- typedef Kokkos::View<my_double*>::size_type size_type;
-
- KOKKOS_INLINE_FUNCTION void
- operator() (const size_type i, value_type& update) const
- {
- if (x_(i) < update) {
- update = x_(i);
- }
- }
-
- KOKKOS_INLINE_FUNCTION void
- join (volatile value_type& dst,
- const volatile value_type& src) const
- {
- if (src < dst) {
- dst = src;
- }
- }
-
- KOKKOS_INLINE_FUNCTION void
- init (value_type& dst) const
- { // The identity under max is -Inf.
- dst= Kokkos::reduction_identity<value_type>::min();
- }
-};
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const double*>
-AcousticSolverTest::computeRhoCj(const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& cj)
-{
- Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
- m_rhocj[j] = rhoj[j]*cj[j];
- });
- return m_rhocj;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const Rdd*[2]>
-AcousticSolverTest::computeAjr(const Kokkos::View<const double*>& rhocj,
- const Kokkos::View<const Rd*[2]>& Cjr)
-{
- Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
- for (int r=0; r<2; ++r) {
- m_Ajr(j,r) = tensorProduct(rhocj(j)*Cjr(j,r), Cjr(j,r));
- }
- });
-
- return m_Ajr;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const Rdd*>
-AcousticSolverTest::computeAr(const Kokkos::View<const Rdd*[2]>& Ajr,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- const Kokkos::View<const int*>& node_nb_cells)
-{
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
- Rdd sum = zero;
- for (int j=0; j<node_nb_cells(r); ++j) {
- const int J = node_cells(r,j);
- const int R = node_cell_local_node(r,j);
- sum += Ajr(J,R);
- }
- m_Ar(r) = sum;
- });
-
- return m_Ar;
-}
-
-KOKKOS_INLINE_FUNCTION
-const Kokkos::View<const Rd*>
-AcousticSolverTest::computeBr(const Kokkos::View<const Rdd*[2]>& Ajr,
- const Kokkos::View<const Rd*[2]>& Cjr,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- const Kokkos::View<const int*>& node_nb_cells)
-{
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
- Rd& br = m_br(r);
- br = zero;
- for (int j=0; j<node_nb_cells(r); ++j) {
- const int J = node_cells(r,j);
- const int R = node_cell_local_node(r,j);
- br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R);
- }
- });
-
- return m_br;
-}
-
-KOKKOS_INLINE_FUNCTION
-Kokkos::View<Rd*>
-AcousticSolverTest::computeUr(const Kokkos::View<const Rdd*>& Ar,
- const Kokkos::View<const Rd*>& br)
-{
- inverse(Ar, m_inv_Ar);
- const Kokkos::View<const Rdd*> invAr = m_inv_Ar;
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r) {
- m_ur[r]=invAr(r)*br(r);
- });
- m_ur[0]=zero;
- m_ur[m_nr-1]=zero;
-
- return m_ur;
-}
-
-KOKKOS_INLINE_FUNCTION
-Kokkos::View<Rd*[2]>
-AcousticSolverTest::computeFjr(const Kokkos::View<const Rdd*[2]>& Ajr,
- const Kokkos::View<const Rd*>& ur,
- const Kokkos::View<const Rd*[2]>& Cjr,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const int*[2]>& cell_nodes)
-{
- Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j) {
- for (int r=0; r<2; ++r) {
- m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+pj(j)*Cjr(j,r);
- }
- });
-
- return m_Fjr;
-}
-
-KOKKOS_INLINE_FUNCTION
-double AcousticSolverTest::
-acoustic_dt(const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*>& cj) const
-{
- Kokkos::parallel_for(m_nj, KOKKOS_LAMBDA(const int& j){
- m_Vj_over_cj[j] = Vj[j]/cj[j];
- });
-
- double dt = std::numeric_limits<double>::max();
- Kokkos::parallel_reduce(m_nj, ReduceMin(m_Vj_over_cj), dt);
-
- return dt;
-}
-
-KOKKOS_INLINE_FUNCTION
-void
-AcousticSolverTest::inverse(const Kokkos::View<const double*>& x,
- Kokkos::View<double*>& inv_x) const
-{
- Kokkos::parallel_for(x.size(), KOKKOS_LAMBDA(const int& r) {
- inv_x(r) = 1./x(r);
- });
-}
-
-KOKKOS_INLINE_FUNCTION
-void
-AcousticSolverTest::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))};
- });
-}
-
-
-KOKKOS_INLINE_FUNCTION
-void AcousticSolverTest::computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
- const Kokkos::View<const Rd*>& xj,
- const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const double*>& cj,
- const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const Rd*[2]>& Cjr,
- const Kokkos::View<const int*[2]>& cell_nodes,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*>& node_nb_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- Kokkos::View<Rd*>& ur,
- Kokkos::View<Rd*[2]>& Fjr)
-{
- const Kokkos::View<const double*> rhocj = computeRhoCj(rhoj, cj);
- const Kokkos::View<const Rdd*[2]> Ajr = computeAjr(rhocj, Cjr);
-
- const Kokkos::View<const Rdd*> Ar = computeAr(Ajr, node_cells, node_cell_local_node, node_nb_cells);
- const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj,
- node_cells, node_cell_local_node, node_nb_cells);
-
- ur = computeUr(Ar, br);
- Fjr = computeFjr(Ajr, ur, Cjr, uj, pj, cell_nodes);
-}
-
-AcousticSolverTest::AcousticSolverTest(const long int& nj)
- : m_nj(nj),
- m_nr(nj+1),
- m_br("br", m_nr),
- m_Ajr("Ajr", m_nj),
- m_Ar("Ar", m_nr),
- m_inv_Ar("inv_Ar", m_nr),
- m_Fjr("Fjr", m_nj),
- m_ur("ur", m_nr),
- 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);
-
- 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("xr", m_nr);
-
- const double delta_x = 1./m_nj;
-
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
- xr[r][0] = r*delta_x;
- });
-
- Kokkos::View<int*[2]> cell_nodes("cell_nodes",m_nj,2);
- Kokkos::View<int*[2]> node_cells("node_cells",m_nr,2);
- Kokkos::View<int*[2]> node_cell_local_node("node_cells",m_nr,2);
- Kokkos::View<int*> node_nb_cells("node_cells",m_nr);
-
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
- node_nb_cells(r) = 2;
- });
- node_nb_cells(0) = 1;
- node_nb_cells(m_nr-1) = 1;
-
- node_cells(0,0) = 0;
- Kokkos::parallel_for(m_nr-2, KOKKOS_LAMBDA(const int& r){
- node_cells(r+1,0) = r;
- node_cells(r+1,1) = r+1;
- });
- node_cells(m_nr-1,0) =m_nj-1;
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- cell_nodes(j,0) = j;
- cell_nodes(j,1) = j+1;
- });
-
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
- for (int J=0; J<node_nb_cells(r); ++J) {
- int j = node_cells(r,J);
- for (int R=0; R<2; ++R) {
- if (cell_nodes(j,R) == r) {
- node_cell_local_node(r,J) = R;
- }
- }
- }
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]);
- });
-
- Kokkos::View<Rd*[2]> Cjr("Cjr",m_nj);
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
- Cjr(j,0)=-1;
- Cjr(j,1)= 1;
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Vj[j] = (xr[cell_nodes(j,1)], Cjr(j,1))
- + (xr[cell_nodes(j,0)], Cjr(j,0));
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- if (xj[j][0]<0.5) {
- rhoj[j]=1;
- } else {
- rhoj[j]=0.125;
- }
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- if (xj[j][0]<0.5) {
- pj[j]=1;
- } else {
- pj[j]=0.1;
- }
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- uj[j] = zero;
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- gammaj[j] = 1.4;
- });
-
- BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-
- block_eos.updateEandCFromRhoP();
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Ej[j] = ej[j]+0.5*(uj[j],uj[j]);
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- mj[j] = rhoj[j] * Vj[j];
- });
-
- Kokkos::View<double*> inv_mj("inv_mj",m_nj);
- inverse(mj, inv_mj);
-
- const double tmax=0.2;
- double t=0;
-
- int itermax=std::numeric_limits<int>::max();
- int iteration=0;
-
- while((t<tmax) and (iteration<itermax)) {
- double dt = 0.4*acoustic_dt(Vj, cj);
- if (t+dt<tmax) {
- t+=dt;
- } else {
- dt=tmax-t;
- t=tmax;
- }
-
- computeExplicitFluxes(xr, xj,
- rhoj, uj, pj, cj, Vj, Cjr,
- cell_nodes, node_cells, node_nb_cells, node_cell_local_node,
- m_ur, m_Fjr);
-
- const Kokkos::View<const Rd*[2]> Fjr = m_Fjr;
- const Kokkos::View<const Rd*> ur = m_ur;
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
- Rd momentum_fluxes = zero;
- double energy_fluxes = 0;
- for (int R=0; R<2; ++R) {
- const int r=cell_nodes(j,R);
- momentum_fluxes += Fjr(j,R);
- energy_fluxes += ((Fjr(j,R), ur[r]));
- }
- uj[j] -= (dt*inv_mj[j]) * momentum_fluxes;
- Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
- ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
- });
-
- Kokkos::parallel_for(m_nr, KOKKOS_LAMBDA(const int& r){
- xr[r] += dt*ur[r];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- xj[j] = 0.5*(xr[j]+xr[j+1]);
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Vj[j] = (xr[cell_nodes(j,1)], Cjr(j,1))
- + (xr[cell_nodes(j,0)], Cjr(j,0));
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- rhoj[j] = mj[j]/Vj[j];
- });
-
- block_eos.updatePandCFromRhoE();
-
- ++iteration;
- }
-
- {
- std::ofstream fout("rho");
- for (int j=0; j<nj; ++j) {
- fout << xj[j][0] << ' ' << rhoj[j] << '\n';
- }
- }
-
- std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
- << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
-}
diff --git a/src/experimental/AcousticSolverTest.hpp b/src/experimental/AcousticSolverTest.hpp
deleted file mode 100644
index 209a18747bb42eb467e6a992faaf8f64e706f856..0000000000000000000000000000000000000000
--- a/src/experimental/AcousticSolverTest.hpp
+++ /dev/null
@@ -1,96 +0,0 @@
-#ifndef ACOUSTIC_SOLVER_TEST_HPP
-#define ACOUSTIC_SOLVER_TEST_HPP
-
-#include <Kokkos_Core.hpp>
-
-#include <TinyVector.hpp>
-#include <TinyMatrix.hpp>
-
-constexpr static size_t dimension = 1;
-
-typedef TinyVector<dimension> Rd;
-typedef TinyMatrix<dimension> Rdd;
-
-class AcousticSolverTest
-{
-private:
-
- inline const Kokkos::View<const double*>
- computeRhoCj(const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& cj);
-
- inline const Kokkos::View<const Rdd*[2]>
- computeAjr(const Kokkos::View<const double*>& rhocj,
- const Kokkos::View<const Rd*[2]>& Cjr);
-
- inline const Kokkos::View<const Rdd*>
- computeAr(const Kokkos::View<const Rdd*[2]>& Ajr,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- const Kokkos::View<const int*>& node_nb_cells);
-
- inline const Kokkos::View<const Rd*>
- computeBr(const Kokkos::View<const Rdd*[2]>& Ajr,
- const Kokkos::View<const Rd*[2]>& Cjr,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- const Kokkos::View<const int*>& node_nb_cells);
-
- Kokkos::View<Rd*>
- computeUr(const Kokkos::View<const Rdd*>& Ar,
- const Kokkos::View<const Rd*>& br);
-
- Kokkos::View<Rd*[2]>
- computeFjr(const Kokkos::View<const Rdd*[2]>& Ajr,
- const Kokkos::View<const Rd*>& ur,
- const Kokkos::View<const Rd*[2]>& Cjr,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const int*[2]>& cell_nodes);
-
- void inverse(const Kokkos::View<const Rdd*>& A,
- Kokkos::View<Rdd*>& inv_A) const;
-
- void inverse(const Kokkos::View<const double*>& x,
- Kokkos::View<double*>& inv_x) const;
-
- inline double acoustic_dt(const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*>& cj) const;
-
- inline void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
- const Kokkos::View<const Rd*>& xj,
- const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const double*>& cj,
- const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const Rd*[2]>& Cjr,
- const Kokkos::View<const int*[2]>& cell_nodes,
- const Kokkos::View<const int*[2]>& node_cells,
- const Kokkos::View<const int*>& node_nb_cells,
- const Kokkos::View<const int*[2]>& node_cell_local_node,
- Kokkos::View<Rd*>& ur,
- Kokkos::View<Rd*[2]>& Fjr);
-
- struct ReduceMin;
-
- const size_t m_nj;
- const size_t m_nr;
-
- Kokkos::View<Rd*> m_br;
- Kokkos::View<Rdd*[2]> m_Ajr;
- Kokkos::View<Rdd*> m_Ar;
- Kokkos::View<Rdd*> m_inv_Ar;
- Kokkos::View<Rd*[2]> m_Fjr;
- Kokkos::View<Rd*> m_ur;
- Kokkos::View<double*> m_rhocj;
- Kokkos::View<double*> m_Vj_over_cj;
-
-
-public:
- AcousticSolverTest(const long int& nj);
-};
-
-#endif // ACOUSTIC_SOLVER_TEST_HPP
diff --git a/src/experimental/CMakeLists.txt b/src/experimental/CMakeLists.txt
deleted file mode 100644
index 1c6043787ec8d6809e233214fc57a465618eddc8..0000000000000000000000000000000000000000
--- a/src/experimental/CMakeLists.txt
+++ /dev/null
@@ -1,17 +0,0 @@
-include_directories(${CMAKE_CURRENT_SOURCE_DIR})
-include_directories(${CMAKE_CURRENT_BINARY_DIR})
-
-include_directories(${Kokkos_INCLUDE_DIRS_RET})
-include_directories(${PASTIS_SOURCE_DIR}/packages/rang/include)
-
-# ------------------- Source files --------------------
-
-add_library(
- PastisExperimental
- AcousticSolverClass.cpp
- AcousticSolverTest.cpp
- RawKokkosAcousticSolver.cpp
- MeshLessAcousticSolver.cpp
- )
-
-
diff --git a/src/experimental/MeshLessAcousticSolver.cpp b/src/experimental/MeshLessAcousticSolver.cpp
deleted file mode 100644
index bf501448130f9e2291a2692d956f2757f5ccfc05..0000000000000000000000000000000000000000
--- a/src/experimental/MeshLessAcousticSolver.cpp
+++ /dev/null
@@ -1,244 +0,0 @@
-#include <MeshLessAcousticSolver.hpp>
-#include <rang.hpp>
-
-#include <memory>
-
-#include <BlockPerfectGas.hpp>
-
-typedef const double my_double;
-
-struct MeshLessAcousticSolver::ReduceMin {
-private:
- const Kokkos::View<my_double*> x_;
-
-public:
- typedef Kokkos::View<my_double*>::non_const_value_type value_type;
-
- ReduceMin(const Kokkos::View<my_double*>& x) : x_ (x) {}
-
- typedef Kokkos::View<my_double*>::size_type size_type;
-
- KOKKOS_INLINE_FUNCTION void
- operator() (const size_type i, value_type& update) const
- {
- if (x_(i) < update) {
- update = x_(i);
- }
- }
-
- KOKKOS_INLINE_FUNCTION void
- join (volatile value_type& dst,
- const volatile value_type& src) const
- {
- if (src < dst) {
- dst = src;
- }
- }
-
- KOKKOS_INLINE_FUNCTION void
- init (value_type& dst) const
- { // The identity under max is -Inf.
- dst= Kokkos::reduction_identity<value_type>::min();
- }
-};
-
-
-KOKKOS_INLINE_FUNCTION
-double MeshLessAcousticSolver::
-acoustic_dt(const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*>& cj) const
-{
- const size_t nj = Vj.size();
- double dt = std::numeric_limits<double>::max();
-
- Kokkos::View<double*> Vj_cj("Vj_cj", nj);
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Vj_cj[j] = Vj[j]/cj[j];
- });
-
- Kokkos::parallel_reduce(nj, ReduceMin(Vj_cj), dt);
-
- return dt;
-}
-
-
-KOKKOS_INLINE_FUNCTION
-void MeshLessAcousticSolver::computeExplicitFluxes(const Kokkos::View<const double*>& xr,
- const Kokkos::View<const double*>& xj,
- const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const double*>& cj,
- const Kokkos::View<const double*>& Vj,
- Kokkos::View<double*>& ur,
- Kokkos::View<double*>& pr) const
-{
- // calcul de ur
- ur[0]=0;
- const size_t nr = ur.size();
- const size_t nj = uj.size();
-
- Kokkos::parallel_for(nj-1, KOKKOS_LAMBDA(const int& j) {
- const int r = j+1;
- const int k = r;
-
- const double ujr = uj[j];
- const double ukr = uj[k];
- const double pjr = pj[j];
- const double pkr = pj[k];
-
- ur[r]=(rhoj[j]*cj[j]*ujr + rhoj[k]*cj[k]*ukr + pjr-pkr)/(rhoj[j]*cj[j]+rhoj[k]*cj[k]);
- });
- ur[nr-1]=0;
-
- // calcul de pr
- pr[0] = pj[0] + rhoj[0]*cj[0]*(ur[0] - uj[0]);
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
- const int r = j+1;
-
- const double ujr = uj[j];
- const double pjr = pj[j];
-
- pr[r]=pjr+rhoj[j]*cj[j]*(ujr-ur[r]);
- });
-}
-
-
-MeshLessAcousticSolver::MeshLessAcousticSolver(const long int& nj)
-{
- Kokkos::View<double*> xj("xj", nj);
- Kokkos::View<double*> rhoj("rhoj", nj);
-
- Kokkos::View<double*> uj("uj", nj);
-
- Kokkos::View<double*> Ej("Ej", nj);
- Kokkos::View<double*> ej("ej", nj);
- Kokkos::View<double*> pj("pj", nj);
- Kokkos::View<double*> Vj("Vj", nj);
- Kokkos::View<double*> gammaj("gammaj", nj);
- Kokkos::View<double*> cj("cj", nj);
- Kokkos::View<double*> mj("mj", nj);
- Kokkos::View<double*> inv_mj("inv_mj", nj);
-
- const int nr=nj+1;
-
- Kokkos::View<double*> xr("xr", nr);
-
- const double delta_x = 1./nj;
-
- Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){
- xr[r] = r*delta_x;
- });
-
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- xj[j] = 0.5*(xr[j]+xr[j+1]);
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Vj[j] = xr[j+1]-xr[j];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- if (xj[j]<0.5) {
- rhoj[j]=1;
- } else {
- rhoj[j]=0.125;
- }
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- if (xj[j]<0.5) {
- pj[j]=1;
- } else {
- pj[j]=0.1;
- }
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- gammaj[j] = 1.4;
- });
-
- BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-
- block_eos.updateEandCFromRhoP();
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Ej[j] = ej[j]+0.5*uj[j]*uj[j];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- mj[j] = rhoj[j] * Vj[j];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- inv_mj[j] = 1./mj[j];
- });
-
- const double tmax=0.2;
- double t=0;
-
- int itermax=std::numeric_limits<int>::max();
- int iteration=0;
-
- while((t<tmax) and (iteration<itermax)) {
- double dt = 0.4*acoustic_dt(Vj, cj);
- if (t+dt<tmax) {
- t+=dt;
- } else {
- dt=tmax-t;
- t=tmax;
- }
-
- Kokkos::View<double*> ur("ur", nr);
- Kokkos::View<double*> pr("pr", nr);
-
- computeExplicitFluxes(xr, xj,
- rhoj, uj, pj, cj, Vj,
- ur, pr);
-
- Kokkos::View<double*> new_uj("new_uj", nj);
- Kokkos::View<double*> new_Ej("new_Ej", nj);
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- int rm=j;
- int rp=j+1;
-
- new_uj[j] = uj[j] + dt*inv_mj[j]*(pr[rm]-pr[rp]);
- new_Ej[j] = Ej[j] + dt*inv_mj[j]*(pr[rm]*ur[rm]-pr[rp]*ur[rp]);
- });
-
- uj=new_uj;
- Ej=new_Ej;
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- ej[j] = Ej[j] - 0.5 * uj[j]*uj[j];
- });
-
- Kokkos::View<double*> xr_new("new_xr", nr);
-
- Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){
- xr_new[r] = xr[r] + dt*ur[r];
- });
-
- xr = xr_new;
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- xj[j] = 0.5*(xr[j]+xr[j+1]);
- Vj[j] = xr[j+1]-xr[j];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- rhoj[j] = mj[j]/Vj[j];
- });
-
- block_eos.updatePandCFromRhoE();
-
- ++iteration;
- }
-
- std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
- << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
-}
diff --git a/src/experimental/MeshLessAcousticSolver.hpp b/src/experimental/MeshLessAcousticSolver.hpp
deleted file mode 100644
index ccdba58383fa691ed713887a3b0fd09e37b92446..0000000000000000000000000000000000000000
--- a/src/experimental/MeshLessAcousticSolver.hpp
+++ /dev/null
@@ -1,31 +0,0 @@
-#ifndef MESH_LESS_ACOUSTIC_SOLVER_HPP
-#define MESH_LESS_ACOUSTIC_SOLVER_HPP
-
-#include <Kokkos_Core.hpp>
-
-class MeshLessAcousticSolver
-{
-private:
-
- inline double e(double rho, double p, double gamma) const;
- inline double p(double rho, double e, double gamma) const;
-
- inline double acoustic_dt(const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*>& cj) const;
-
- inline void computeExplicitFluxes(const Kokkos::View<const double*>& xr,
- const Kokkos::View<const double*>& xj,
- const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const double*>& cj,
- const Kokkos::View<const double*>& Vj,
- Kokkos::View<double*>& ur,
- Kokkos::View<double*>& pr) const;
-
- struct ReduceMin;
-public:
- MeshLessAcousticSolver(const long int& nj);
-};
-
-#endif // MESH_LESS_ACOUSTIC_SOLVER_HPP
diff --git a/src/experimental/RawKokkosAcousticSolver.cpp b/src/experimental/RawKokkosAcousticSolver.cpp
deleted file mode 100644
index d27549d349183d6eb714be477750ca4fee6381b5..0000000000000000000000000000000000000000
--- a/src/experimental/RawKokkosAcousticSolver.cpp
+++ /dev/null
@@ -1,248 +0,0 @@
-#include <RawKokkosAcousticSolver.hpp>
-#include <Kokkos_Core.hpp>
-#include <rang.hpp>
-
-namespace RawKokkos
-{
-
-inline double e(double rho, double p, double gamma)
-{
- return p/(rho*(gamma-1));
-}
-
-inline double p(double rho, double e, double gamma)
-{
- return (gamma-1)*rho*e;
-}
-
-typedef const double my_double;
-
-struct ReduceMin {
-private:
- const Kokkos::View<my_double*> x_;
-
-public:
- typedef Kokkos::View<my_double*>::non_const_value_type value_type;
-
- ReduceMin(const Kokkos::View<my_double*>& x) : x_ (x) {}
-
- typedef Kokkos::View<my_double*>::size_type size_type;
-
- KOKKOS_INLINE_FUNCTION void
- operator() (const size_type i, value_type& update) const
- {
- if (x_(i) < update) {
- update = x_(i);
- }
- }
-
- KOKKOS_INLINE_FUNCTION void
- join (volatile value_type& dst,
- const volatile value_type& src) const
- {
- if (src < dst) {
- dst = src;
- }
- }
-
- KOKKOS_INLINE_FUNCTION void
- init (value_type& dst) const
- { // The identity under max is -Inf.
- dst= Kokkos::reduction_identity<value_type>::min();
- }
-};
-
-
-double acoustic_dt(const Kokkos::View<double*>& Vj,
- const Kokkos::View<double*>& cj)
-{
- const size_t nj = Vj.size();
- double dt = std::numeric_limits<double>::max();
-
- Kokkos::View<double*> Vj_cj("Vj_cj", nj);
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Vj_cj[j] = Vj[j]/cj[j];
- });
-
- Kokkos::parallel_reduce(nj, ReduceMin(Vj_cj), dt);
-
- // Kokkos::parallel_reduce(n, KOKKOS_LAMBDA(const long i, long& lcount) {
- // lcount += (i % 2) == 0;
- // }, count);
- return dt;
-}
-
-
-void computeExplicitFluxes(const Kokkos::View<double*>& xr,
- const Kokkos::View<double*>& xj,
- const Kokkos::View<double*>& rhoj,
- const Kokkos::View<double*>& uj,
- const Kokkos::View<double*>& pj,
- const Kokkos::View<double*>& cj,
- const Kokkos::View<double*>& Vj,
- Kokkos::View<double*>& ur,
- Kokkos::View<double*>& pr)
-{
- // calcul de ur
- ur[0]=0;
- const size_t nr = ur.size();
- const size_t nj = uj.size();
-
- Kokkos::parallel_for(nj-1, KOKKOS_LAMBDA(const int& j) {
- const int r = j+1;
- const int k = r;
- const double ujr = uj[j];
- const double ukr = uj[k];
- const double pjr = pj[j];
- const double pkr = pj[k];
-
- ur[r]=(rhoj[j]*cj[j]*ujr + rhoj[k]*cj[k]*ukr + pjr-pkr)/(rhoj[j]*cj[j]+rhoj[k]*cj[k]);
- });
- ur[nr-1]=0;
-
- // calcul de pr
- pr[0] = pj[0] + rhoj[0]*cj[0]*(ur[0] - uj[0]);
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j) {
- const int r = j+1;
-
- const double ujr = uj[j];
- const double pjr = pj[j];
-
- pr[r]=pjr+rhoj[j]*cj[j]*(ujr-ur[r]);
- });
-}
-
-
-void AcousticSolver(const long int& nj)
-{
- Kokkos::View<double*> xj("xj", nj);
- Kokkos::View<double*> rhoj("rhoj", nj);
-
- Kokkos::View<double*> uj("uj", nj);
-
- Kokkos::View<double*> Ej("Ej", nj);
- Kokkos::View<double*> ej("ej", nj);
- Kokkos::View<double*> pj("pj", nj);
- Kokkos::View<double*> Vj("Vj", nj);
- Kokkos::View<double*> gammaj("gammaj", nj);
- Kokkos::View<double*> cj("cj", nj);
- Kokkos::View<double*> mj("mj", nj);
- Kokkos::View<double*> inv_mj("inv_mj", nj);
-
- const int nr=nj+1;
-
- Kokkos::View<double*> xr("xr", nr);
-
- const double delta_x = 1./nj;
-
- Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){
- xr[r] = r*delta_x;
- });
-
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- xj[j] = 0.5*(xr[j]+xr[j+1]);
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Vj[j] = xr[j+1]-xr[j];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- if (xj[j]<0.5) {
- rhoj[j]=1;
- } else {
- rhoj[j]=0.125;
- }
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- if (xj[j]<0.5) {
- pj[j]=1;
- } else {
- pj[j]=0.1;
- }
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- gammaj[j] = 1.4;
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- ej[j] = e(rhoj[j],pj[j],gammaj[j]);
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- Ej[j] = ej[j]+0.5*uj[j]*uj[j];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- cj[j] = std::sqrt(gammaj[j]*pj[j]/rhoj[j]);
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- mj[j] = rhoj[j] * Vj[j];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- inv_mj[j] = 1./mj[j];
- });
-
- const double tmax=0.2;
- double t=0;
-
- int itermax=std::numeric_limits<int>::max();
- int iteration=0;
-
- while((t<tmax) and (iteration<itermax)) {
- double dt = 0.4*acoustic_dt(Vj, cj);
- if (t+dt<tmax) {
- t+=dt;
- } else {
- dt=tmax-t;
- t=tmax;
- }
-
- Kokkos::View<double*> ur("ur", nr);
- Kokkos::View<double*> pr("pr", nr);
-
- computeExplicitFluxes(xr, xj,
- rhoj, uj, pj, cj, Vj,
- ur, pr);
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- int rm=j;
- int rp=j+1;
-
- uj[j] += dt*inv_mj[j]*(pr[rm]-pr[rp]);
- Ej[j] += dt*inv_mj[j]*(pr[rm]*ur[rm]-pr[rp]*ur[rp]);
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- ej[j] = Ej[j] - 0.5 * uj[j]*uj[j];
- });
-
- Kokkos::parallel_for(nr, KOKKOS_LAMBDA(const int& r){
- xr[r] += dt*ur[r];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- xj[j] = 0.5*(xr[j]+xr[j+1]);
- Vj[j] = xr[j+1]-xr[j];
- });
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- rhoj[j] = mj[j]/Vj[j];
- pj[j] = p(rhoj[j], ej[j], gammaj[j]);
- cj[j] = std::sqrt(gammaj[j]*pj[j]/rhoj[j]); // inv_mj*vj
- });
-
- ++iteration;
- }
-
- std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
- << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-}
-
-}
diff --git a/src/experimental/RawKokkosAcousticSolver.hpp b/src/experimental/RawKokkosAcousticSolver.hpp
deleted file mode 100644
index e9c36ff437a7ef0f0c28c72e93b0a367d8142185..0000000000000000000000000000000000000000
--- a/src/experimental/RawKokkosAcousticSolver.hpp
+++ /dev/null
@@ -1,9 +0,0 @@
-#ifndef RAW_KOKKOS_ACOUSTIC_SOLVER_HPP
-#define RAW_KOKKOS_ACOUSTIC_SOLVER_HPP
-
-namespace RawKokkos
-{
- void AcousticSolver(const long int& nj);
-}
-
-#endif // RAW_KOKKOS_ACOUSTIC_SOLVER_HPP
diff --git a/src/main.cpp b/src/main.cpp
index b1a3decd13b36a4b738c49c9815a0c4bfb327d5f..ba6f6c504ae02f7859d722de46295e192f92fffe 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -147,10 +147,10 @@ int main(int argc, char *argv[])
bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
}
- typedef Connectivity1D ConnectivityType;
- typedef Mesh<ConnectivityType> MeshType;
- typedef MeshData<MeshType> MeshDataType;
- typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;
+ using ConnectivityType = Connectivity1D;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<MeshType>;
+ using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
@@ -192,9 +192,9 @@ int main(int argc, char *argv[])
AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
- typedef TinyVector<MeshType::dimension> Rd;
+ using Rd = TinyVector<MeshType::dimension>;
- const Kokkos::View<const double*> Vj = mesh_data.Vj();
+ const CellValue<const double>& Vj = mesh_data.Vj();
const double tmax=0.2;
double t=0;
@@ -202,11 +202,11 @@ int main(int argc, char *argv[])
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();
+ CellValue<double>& rhoj = unknowns.rhoj();
+ CellValue<double>& ej = unknowns.ej();
+ CellValue<double>& pj = unknowns.pj();
+ CellValue<double>& gammaj = unknowns.gammaj();
+ CellValue<double>& cj = unknowns.cj();
BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
@@ -233,10 +233,10 @@ int main(int argc, char *argv[])
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();
+ const CellValue<const Rd>& xj = mesh_data.xj();
+ const CellValue<const double>& rhoj = unknowns.rhoj();
std::ofstream fout("rho");
- for (size_t j=0; j<mesh.numberOfCells(); ++j) {
+ for (CellId j=0; j<mesh.numberOfCells(); ++j) {
fout << xj[j][0] << ' ' << rhoj[j] << '\n';
}
}
@@ -256,10 +256,10 @@ int main(int argc, char *argv[])
bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
}
- typedef Connectivity2D ConnectivityType;
- typedef Mesh<ConnectivityType> MeshType;
- typedef MeshData<MeshType> MeshDataType;
- typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;
+ using ConnectivityType = Connectivity2D;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<MeshType>;
+ using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
@@ -301,7 +301,7 @@ int main(int argc, char *argv[])
AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
- const Kokkos::View<const double*> Vj = mesh_data.Vj();
+ const CellValue<const double>& Vj = mesh_data.Vj();
const double tmax=0.2;
double t=0;
@@ -309,11 +309,11 @@ int main(int argc, char *argv[])
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();
+ CellValue<double>& rhoj = unknowns.rhoj();
+ CellValue<double>& ej = unknowns.ej();
+ CellValue<double>& pj = unknowns.pj();
+ CellValue<double>& gammaj = unknowns.gammaj();
+ CellValue<double>& cj = unknowns.cj();
BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
@@ -352,10 +352,10 @@ int main(int argc, char *argv[])
bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
}
- typedef Connectivity3D ConnectivityType;
- typedef Mesh<ConnectivityType> MeshType;
- typedef MeshData<MeshType> MeshDataType;
- typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;
+ using ConnectivityType = Connectivity3D;
+ using MeshType = Mesh<ConnectivityType>;
+ using MeshDataType = MeshData<MeshType>;
+ using UnknownsType = FiniteVolumesEulerUnknowns<MeshDataType>;
const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
@@ -397,7 +397,7 @@ int main(int argc, char *argv[])
AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
- const Kokkos::View<const double*> Vj = mesh_data.Vj();
+ const CellValue<const double>& Vj = mesh_data.Vj();
const double tmax=0.2;
double t=0;
@@ -405,11 +405,11 @@ int main(int argc, char *argv[])
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();
+ CellValue<double>& rhoj = unknowns.rhoj();
+ CellValue<double>& ej = unknowns.ej();
+ CellValue<double>& pj = unknowns.pj();
+ CellValue<double>& gammaj = unknowns.gammaj();
+ CellValue<double>& cj = unknowns.cj();
BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index eb5fb87b36947a5e4778f2e38161e6f35cb3f92d..f84f3835abe00acb25561ada3d071a29c5adc167 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -4,14 +4,14 @@
template<>
void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
{
- typedef std::tuple<unsigned int, unsigned short, bool> CellFaceInfo;
+ using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
const auto& cell_to_node_matrix
- = this->getMatrix(ItemType::cell, ItemType::node);
+ = this->_getMatrix(ItemType::cell, ItemType::node);
- Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", this->numberOfCells());
+ CellValue<unsigned short> cell_nb_faces(*this);
std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
- for (unsigned int j=0; j<this->numberOfCells(); ++j) {
+ for (CellId j=0; j<this->numberOfCells(); ++j) {
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
switch (m_cell_type[j]) {
@@ -97,12 +97,12 @@ void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
{
auto& cell_to_face_matrix
- = m_item_to_item_matrix[itemId(ItemType::cell)][itemId(ItemType::face)];
+ = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)];
std::vector<std::vector<unsigned int>> cell_to_face_vector(this->numberOfCells());
- for (size_t j=0; j<cell_to_face_vector.size(); ++j) {
+ for (CellId j=0; j<cell_to_face_vector.size(); ++j) {
cell_to_face_vector[j].resize(cell_nb_faces[j]);
}
- int l=0;
+ FaceId l=0;
for (const auto& face_cells_vector : face_cells_map) {
const auto& cells_vector = face_cells_vector.second;
for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
@@ -129,7 +129,7 @@ void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
{
auto& face_to_node_matrix
- = m_item_to_item_matrix[itemId(ItemType::face)][itemId(ItemType::node)];
+ = m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)];
std::vector<std::vector<unsigned int>> face_to_node_vector(face_cells_map.size());
int l=0;
@@ -157,16 +157,16 @@ template<>
void Connectivity<2>::_computeCellFaceAndFaceNodeConnectivities()
{
const auto& cell_to_node_matrix
- = this->getMatrix(ItemType::cell, ItemType::node);
+ = this->_getMatrix(ItemType::cell, ItemType::node);
// In 2D faces are simply define
- typedef std::pair<unsigned int, unsigned short> CellFaceId;
+ using CellFaceId = std::pair<CellId, unsigned short>;
std::map<Face, std::vector<CellFaceId>> face_cells_map;
- for (unsigned int j=0; j<this->numberOfCells(); ++j) {
+ for (CellId j=0; j<this->numberOfCells(); ++j) {
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
for (unsigned short r=0; r<cell_nodes.length; ++r) {
- unsigned int node0_id = cell_nodes(r);
- unsigned int node1_id = cell_nodes((r+1)%cell_nodes.length);
+ NodeId node0_id = cell_nodes(r);
+ NodeId node1_id = cell_nodes((r+1)%cell_nodes.length);
if (node1_id<node0_id) {
std::swap(node0_id, node1_id);
}
@@ -175,7 +175,7 @@ void Connectivity<2>::_computeCellFaceAndFaceNodeConnectivities()
}
{
- int l=0;
+ FaceId l=0;
for (const auto& face_cells_vector : face_cells_map) {
const Face& face = face_cells_vector.first;
m_face_number_map[face] = l;
@@ -192,7 +192,7 @@ void Connectivity<2>::_computeCellFaceAndFaceNodeConnectivities()
++l;
}
auto& face_to_node_matrix
- = m_item_to_item_matrix[itemId(ItemType::face)][itemId(ItemType::node)];
+ = m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)];
face_to_node_matrix = face_to_node_vector;
}
@@ -207,7 +207,7 @@ void Connectivity<2>::_computeCellFaceAndFaceNodeConnectivities()
++l;
}
auto& face_to_cell_matrix
- = m_item_to_item_matrix[itemId(ItemType::face)][itemId(ItemType::cell)];
+ = m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::cell)];
face_to_cell_matrix = face_to_cell_vector;
}
}
@@ -221,21 +221,21 @@ Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
Assert(cell_by_node_vector.size() == cell_type_vector.size());
auto& cell_to_node_matrix
- = m_item_to_item_matrix[itemId(ItemType::cell)][itemId(ItemType::node)];
+ = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];
cell_to_node_matrix = cell_by_node_vector;
Assert(this->numberOfCells()>0);
{
- Kokkos::View<CellType*> cell_type("cell_type", this->numberOfCells());
- Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ CellValue<CellType> cell_type(*this);
+ Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
cell_type[j] = cell_type_vector[j];
});
m_cell_type = cell_type;
}
{
- Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", this->numberOfCells());
- Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ CellValue<double> inv_cell_nb_nodes(*this);
+ Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
inv_cell_nb_nodes[j] = 1./cell_nodes.length;
});
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 0994fc9c84e4105c3f9de66b55231f3fd61ff8fe..cfcb7f67909e2c7668bb11d842eb5437e57734c2 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -5,10 +5,13 @@
#include <TinyVector.hpp>
#include <Kokkos_Core.hpp>
+#include <ItemValue.hpp>
#include <IConnectivity.hpp>
#include <ConnectivityMatrix.hpp>
+#include <ItemToItemMatrix.hpp>
+
#include <SubItemValuePerItem.hpp>
#include <ConnectivityComputer.hpp>
@@ -126,7 +129,7 @@ class ConnectivityFace<3>
};
bool m_reversed;
- std::vector<unsigned int> m_node_id_list;
+ std::vector<NodeId::base_type> m_node_id_list;
friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
{
@@ -222,19 +225,20 @@ class ConnectivityFace<3>
~ConnectivityFace() = default;
};
+
template <size_t Dimension>
class Connectivity final
: public IConnectivity
{
private:
- constexpr static auto& itemId = ItemId<Dimension>::itemId;
+ constexpr static auto& itemTId = ItemTypeId<Dimension>::itemTId;
public:
static constexpr size_t dimension = Dimension;
private:
ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1];
- Kokkos::View<const CellType*> m_cell_type;
+ CellValue<const CellType> m_cell_type;
FaceValuePerCell<const bool> m_cell_face_is_reversed;
@@ -259,11 +263,11 @@ class Connectivity final
std::vector<RefFaceList> m_ref_face_list;
std::vector<RefNodeList> m_ref_node_list;
- Kokkos::View<double*> m_inv_cell_nb_nodes;
+ CellValue<const double> m_inv_cell_nb_nodes;
using Face = ConnectivityFace<Dimension>;
- std::unordered_map<Face, unsigned int, typename Face::Hash> m_face_number_map;
+ std::unordered_map<Face, FaceId, typename Face::Hash> m_face_number_map;
void _computeCellFaceAndFaceNodeConnectivities();
@@ -281,6 +285,23 @@ class Connectivity final
return sub_item_value_per_item;
}
+ friend class ConnectivityComputer;
+
+ KOKKOS_INLINE_FUNCTION
+ const ConnectivityMatrix& _getMatrix(const ItemType& item_type_0,
+ const ItemType& item_type_1) const
+ {
+ const ConnectivityMatrix& connectivity_matrix
+ = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
+ if (not connectivity_matrix.isBuilt()) {
+ const_cast<ConnectivityMatrix&>(connectivity_matrix)
+ = m_connectivity_computer
+ . computeConnectivityMatrix(*this, item_type_0, item_type_1);
+
+ }
+ return connectivity_matrix;
+ }
+
public:
KOKKOS_INLINE_FUNCTION
@@ -288,25 +309,92 @@ class Connectivity final
const ItemType& item_type_1) const
{
const ConnectivityMatrix& connectivity_matrix
- = m_item_to_item_matrix[itemId(item_type_0)][itemId(item_type_1)];
+ = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
return connectivity_matrix.isBuilt();
}
+ template <ItemType source_item_type,
+ ItemType target_item_type>
+ KOKKOS_INLINE_FUNCTION
+ const auto getItemToItemMatrix() const
+ {
+ return ItemToItemMatrix<source_item_type,
+ target_item_type>(_getMatrix(source_item_type,
+ target_item_type));
+ }
KOKKOS_INLINE_FUNCTION
- const ConnectivityMatrix& getMatrix(const ItemType& item_type_0,
- const ItemType& item_type_1) const
+ const auto cellToFaceMatrix() const
{
- const ConnectivityMatrix& connectivity_matrix
- = m_item_to_item_matrix[itemId(item_type_0)][itemId(item_type_1)];
- if (not connectivity_matrix.isBuilt()) {
- const_cast<ConnectivityMatrix&>(connectivity_matrix)
- = m_connectivity_computer
- . computeConnectivityMatrix(*this, item_type_0, item_type_1);
+ return this->template getItemToItemMatrix<ItemType::cell, ItemType::face>();
+ }
- }
- return connectivity_matrix;
+ KOKKOS_INLINE_FUNCTION
+ const auto cellToEdgeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::cell, ItemType::edge>();
}
+ KOKKOS_INLINE_FUNCTION
+ const auto cellToNodeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::cell, ItemType::node>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto faceToCellMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::face, ItemType::cell>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto faceToEdgeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::face, ItemType::edge>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto faceToNodeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::face, ItemType::node>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto edgeToCellMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::edge, ItemType::cell>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto edgeToFaceMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::edge, ItemType::face>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto edgeToNodeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::edge, ItemType::node>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto nodeToCellMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::node, ItemType::cell>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto nodeToFaceMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::node, ItemType::face>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto nodeToEdgeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::node, ItemType::edge>();
+ }
+
+
KOKKOS_INLINE_FUNCTION
const auto& cellFaceIsReversed() const
{
@@ -439,30 +527,38 @@ class Connectivity final
}
KOKKOS_INLINE_FUNCTION
- size_t numberOfNodes() const
+ size_t numberOfNodes() const final
{
const auto& node_to_cell_matrix
- = this->getMatrix(ItemType::node,ItemType::cell);
+ = this->_getMatrix(ItemType::node,ItemType::cell);
return node_to_cell_matrix.numRows();
}
KOKKOS_INLINE_FUNCTION
- size_t numberOfFaces() const
+ size_t numberOfEdges() const final
+ {
+ const auto& edge_to_node_matrix
+ = this->_getMatrix(ItemType::edge,ItemType::node);
+ return edge_to_node_matrix.numRows();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfFaces() const final
{
const auto& face_to_node_matrix
- = this->getMatrix(ItemType::face,ItemType::cell);
+ = this->_getMatrix(ItemType::face,ItemType::cell);
return face_to_node_matrix.numRows();
}
KOKKOS_INLINE_FUNCTION
- size_t numberOfCells() const
+ size_t numberOfCells() const final
{
const auto& cell_to_node_matrix
- = this->getMatrix(ItemType::cell,ItemType::node);
+ = this->_getMatrix(ItemType::cell,ItemType::node);
return cell_to_node_matrix.numRows();
}
- const Kokkos::View<const double*> invCellNbNodes() const
+ const CellValue<const double>& invCellNbNodes() const
{
#warning add calculation on demand when variables will be defined
return m_inv_cell_nb_nodes;
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index 0c87ccb6d64f840cd450eb5ded67c6a337a5cf25..aa6f5c56ae612133afc3040d84fd6fe3c807b507 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -14,7 +14,7 @@ computeConnectivityMatrix(const ConnectivityType& connectivity,
ConnectivityMatrix item_to_child_item_matrix;
if (connectivity.isConnectivityMatrixBuilt(child_item_type, item_type)) {
const ConnectivityMatrix& child_to_item_matrix
- = connectivity.getMatrix(child_item_type, item_type);
+ = connectivity._getMatrix(child_item_type, item_type);
std::cout << "computing connectivity "
<< itemName(item_type) << " -> " << itemName(child_item_type) << '\n';
@@ -83,12 +83,12 @@ SubItemValuePerItem<const unsigned short, child_item_type, item_type>
ConnectivityComputer::computeLocalItemNumberInChildItem(const ConnectivityType& connectivity) const
{
const ConnectivityMatrix& child_item_to_items_matrix
- = connectivity.getMatrix(child_item_type, item_type);
+ = connectivity._getMatrix(child_item_type, item_type);
const ConnectivityMatrix& item_to_child_items_matrix
- = connectivity.getMatrix(item_type, child_item_type);
+ = connectivity._getMatrix(item_type, child_item_type);
SubItemValuePerItem<unsigned short, child_item_type, item_type> item_number_in_child_item(connectivity);
- for (unsigned int r=0; r<item_to_child_items_matrix.numRows(); ++r) {
+ for (ItemIdT<item_type> r=0; r<item_to_child_items_matrix.numRows(); ++r) {
const auto& item_to_child_items = item_to_child_items_matrix.rowConst(r);
for (unsigned short J=0; J<item_to_child_items.length; ++J) {
const unsigned int j = item_to_child_items(J);
diff --git a/src/mesh/ConnectivityMatrix.hpp b/src/mesh/ConnectivityMatrix.hpp
index 356ee6da56d474836180cc867318c2d4f912549a..4387070d0bc9ffc91e88f16c160b87a65e176925 100644
--- a/src/mesh/ConnectivityMatrix.hpp
+++ b/src/mesh/ConnectivityMatrix.hpp
@@ -7,13 +7,13 @@
class ConnectivityMatrix
{
private:
- typedef Kokkos::StaticCrsGraph<unsigned int, Kokkos::HostSpace> HostMatrix;
+ using HostMatrix = Kokkos::StaticCrsGraph<unsigned int, Kokkos::HostSpace>;
HostMatrix m_host_matrix;
bool m_is_built{false};
public:
- typedef HostMatrix::row_map_type HostRowType;
+ using HostRowType = HostMatrix::row_map_type;
const bool& isBuilt() const
{
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index f16ab2fa02e06aad3715eeb7c84028476c7565ae..ea003eb9b6c41d7f049c648ab62b01a01b0add26 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -310,7 +310,7 @@ void GmshReader::__readVertices()
}
__verticesNumbers.resize(numberOfVerices);
- __vertices = Kokkos::View<R3*>("vertices",numberOfVerices);
+ __vertices = Array<R3>(numberOfVerices);
if (not __binary) {
for (int i=0; i<numberOfVerices; ++i) {
@@ -552,38 +552,32 @@ GmshReader::__proceedData()
switch (i) {
// Supported entities
case 0: { // edges
- __edges = Kokkos::View<Edge*>("edges",
- elementNumber[i]);
+ __edges = Array<Edge>(elementNumber[i]);
__edges_ref.resize(elementNumber[i]);
break;
}
case 1: { // triangles
- __triangles = Kokkos::View<Triangle*>("triangles",
- elementNumber[i]);
+ __triangles = Array<Triangle>(elementNumber[i]);
__triangles_ref.resize(elementNumber[i]);
break;
}
case 2: { // quadrangles
- __quadrangles = Kokkos::View<Quadrangle*>("quadrangles",
- elementNumber[i]);
+ __quadrangles = Array<Quadrangle>(elementNumber[i]);
__quadrangles_ref.resize(elementNumber[i]);
break;
}
case 3: { // tetrahedra
- __tetrahedra = Kokkos::View<Tetrahedron*>("tetrahedra",
- elementNumber[i]);
+ __tetrahedra = Array<Tetrahedron>(elementNumber[i]);
__tetrahedra_ref.resize(elementNumber[i]);
break;
}
case 4: {// hexahedra
- __hexahedra = Kokkos::View<Hexahedron*>("hexahedra",
- elementNumber[i]);
+ __hexahedra = Array<Hexahedron>(elementNumber[i]);
__hexahedra_ref.resize(elementNumber[i]);
}
// Ignored entities
case 14: {// point
- __points = Kokkos::View<Point*>("points",
- elementNumber[i]);
+ __points = Array<Point>(elementNumber[i]);
__points_ref.resize(elementNumber[i]);
break;
}
@@ -798,7 +792,7 @@ GmshReader::__proceedData()
std::vector<CellType> cell_type_vector(nb_cells);
std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
- const size_t nb_tetrahedra = __tetrahedra.extent(0);
+ const size_t nb_tetrahedra = __tetrahedra.size();
for (size_t j=0; j<nb_tetrahedra; ++j) {
cell_by_node_vector[j].resize(4);
for (int r=0; r<4; ++r) {
@@ -806,7 +800,7 @@ GmshReader::__proceedData()
}
cell_type_vector[j] = CellType::Tetrahedron;
}
- const size_t nb_hexahedra = __hexahedra.extent(0);
+ const size_t nb_hexahedra = __hexahedra.size();
for (size_t j=0; j<nb_hexahedra; ++j) {
const size_t jh = nb_tetrahedra+j;
cell_by_node_vector[jh].resize(8);
@@ -821,13 +815,13 @@ GmshReader::__proceedData()
Connectivity3D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
- for (unsigned int f=0; f<__triangles.extent(0); ++f) {
+ for (unsigned int f=0; f<__triangles.size(); ++f) {
const unsigned int face_number
= connectivity.getFaceNumber({__triangles[f][0], __triangles[f][1], __triangles[f][2]});
const unsigned int& ref = __triangles_ref[f];
ref_faces_map[ref].push_back(face_number);
}
- for (unsigned int f=0; f<__quadrangles.extent(0); ++f) {
+ for (unsigned int f=0; f<__quadrangles.size(); ++f) {
const unsigned int face_number
= connectivity.getFaceNumber({__quadrangles[f][0], __quadrangles[f][1],
__quadrangles[f][2], __quadrangles[f][3]});
@@ -835,7 +829,7 @@ GmshReader::__proceedData()
ref_faces_map[ref].push_back(face_number);
}
for (const auto& ref_face_list : ref_faces_map) {
- Kokkos::View<unsigned int*> face_list("face_list", ref_face_list.second.size());
+ Array<FaceId> face_list(ref_face_list.second.size());
for (size_t j=0; j<ref_face_list.second.size(); ++j) {
face_list[j]=ref_face_list.second[j];
}
@@ -843,11 +837,11 @@ GmshReader::__proceedData()
connectivity.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
}
- typedef Mesh<Connectivity3D> MeshType;
- typedef TinyVector<3, double> Rd;
+ using MeshType = Mesh<Connectivity3D>;
+ using Rd = TinyVector<3, double>;
- Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
- for (size_t i=0; i<__vertices.extent(0); ++i) {
+ NodeValue<Rd> xr(connectivity);
+ for (NodeId i=0; i<__vertices.size(); ++i) {
xr[i][0] = __vertices[i][0];
xr[i][1] = __vertices[i][1];
xr[i][2] = __vertices[i][2];
@@ -860,7 +854,7 @@ GmshReader::__proceedData()
std::vector<CellType> cell_type_vector(nb_cells);
std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
- const size_t nb_triangles = __triangles.extent(0);
+ const size_t nb_triangles = __triangles.size();
for (size_t j=0; j<nb_triangles; ++j) {
cell_by_node_vector[j].resize(3);
for (int r=0; r<3; ++r) {
@@ -869,7 +863,7 @@ GmshReader::__proceedData()
cell_type_vector[j] = CellType::Triangle;
}
- const size_t nb_quadrangles = __quadrangles.extent(0);
+ const size_t nb_quadrangles = __quadrangles.size();
for (size_t j=0; j<nb_quadrangles; ++j) {
const size_t jq = j+nb_triangles;
cell_by_node_vector[jq].resize(4);
@@ -884,14 +878,14 @@ GmshReader::__proceedData()
Connectivity2D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
- for (unsigned int e=0; e<__edges.extent(0); ++e) {
+ for (unsigned int e=0; e<__edges.size(); ++e) {
const unsigned int edge_number = connectivity.getFaceNumber({__edges[e][0], __edges[e][1]});
const unsigned int& ref = __edges_ref[e];
ref_faces_map[ref].push_back(edge_number);
}
for (const auto& ref_face_list : ref_faces_map) {
- Kokkos::View<unsigned int*> face_list("face_list", ref_face_list.second.size());
+ Array<FaceId> face_list(ref_face_list.second.size());
for (size_t j=0; j<ref_face_list.second.size(); ++j) {
face_list[j]=ref_face_list.second[j];
}
@@ -900,14 +894,14 @@ GmshReader::__proceedData()
}
std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
- for (unsigned int r=0; r<__points.extent(0); ++r) {
+ for (unsigned int r=0; r<__points.size(); ++r) {
const unsigned int point_number = __points[r];
const unsigned int& ref = __points_ref[r];
ref_points_map[ref].push_back(point_number);
}
for (const auto& ref_point_list : ref_points_map) {
- Kokkos::View<unsigned int*> point_list("point_list", ref_point_list.second.size());
+ Array<NodeId> point_list(ref_point_list.second.size());
for (size_t j=0; j<ref_point_list.second.size(); ++j) {
point_list[j]=ref_point_list.second[j];
}
@@ -915,11 +909,11 @@ GmshReader::__proceedData()
connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_list));
}
- typedef Mesh<Connectivity2D> MeshType;
- typedef TinyVector<2, double> Rd;
+ using MeshType = Mesh<Connectivity2D>;
+ using Rd = TinyVector<2, double>;
- Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
- for (size_t i=0; i<__vertices.extent(0); ++i) {
+ NodeValue<Rd> xr(connectivity);
+ for (NodeId i=0; i<__vertices.size(); ++i) {
xr[i][0] = __vertices[i][0];
xr[i][1] = __vertices[i][1];
}
@@ -944,14 +938,14 @@ GmshReader::__proceedData()
Connectivity1D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
- for (unsigned int r=0; r<__points.extent(0); ++r) {
+ for (unsigned int r=0; r<__points.size(); ++r) {
const unsigned int point_number = __points[r];
const unsigned int& ref = __points_ref[r];
ref_points_map[ref].push_back(point_number);
}
for (const auto& ref_point_list : ref_points_map) {
- Kokkos::View<unsigned int*> point_list("point_list", ref_point_list.second.size());
+ Array<NodeId> point_list(ref_point_list.second.size());
for (size_t j=0; j<ref_point_list.second.size(); ++j) {
point_list[j]=ref_point_list.second[j];
}
@@ -959,11 +953,11 @@ GmshReader::__proceedData()
connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_list));
}
- typedef Mesh<Connectivity1D> MeshType;
- typedef TinyVector<1, double> Rd;
+ using MeshType = Mesh<Connectivity1D>;
+ using Rd = TinyVector<1, double>;
- Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
- for (size_t i=0; i<__vertices.extent(0); ++i) {
+ NodeValue<Rd> xr(connectivity);
+ for (NodeId i=0; i<__vertices.size(); ++i) {
xr[i][0] = __vertices[i][0];
}
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index 8825a67df3d6cf1593ff54bb6af8b7a96b079b5b..ba9ffe2502b4df46f59308f411747c7c2deaaac0 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -1,7 +1,7 @@
#ifndef GMSH_READER_HPP
#define GMSH_READER_HPP
-#include <Kokkos_Core.hpp>
+#include <Array.hpp>
#include <TinyVector.hpp>
#include <Mesh.hpp>
@@ -17,7 +17,7 @@
class GmshReader
{
public:
- typedef TinyVector<3, double> R3;
+ using R3 = TinyVector<3, double>;
class PhysicalRefId
{
@@ -78,25 +78,30 @@ private:
*/
std::vector<int> __verticesNumbers;
- Kokkos::View<R3*> __vertices;
+ Array<R3> __vertices;
- typedef unsigned int Point;
- Kokkos::View<Point*> __points;
+ using Point = unsigned int;
+ Array<Point> __points;
std::vector<int> __points_ref;
- typedef TinyVector<2,unsigned int> Edge;
- Kokkos::View<Edge*> __edges;
+
+ using Edge = TinyVector<2,unsigned int>;
+ Array<Edge> __edges;
std::vector<int> __edges_ref;
- typedef TinyVector<3,unsigned int> Triangle;
- Kokkos::View<Triangle*> __triangles;
+
+ using Triangle = TinyVector<3,unsigned int>;
+ Array<Triangle> __triangles;
std::vector<int> __triangles_ref;
- typedef TinyVector<4,unsigned int> Quadrangle;
- Kokkos::View<Quadrangle*> __quadrangles;
+
+ using Quadrangle = TinyVector<4,unsigned int>;
+ Array<Quadrangle> __quadrangles;
std::vector<int> __quadrangles_ref;
- typedef TinyVector<4,unsigned int> Tetrahedron;
- Kokkos::View<Tetrahedron*> __tetrahedra;
+
+ using Tetrahedron = TinyVector<4,unsigned int>;
+ Array<Tetrahedron> __tetrahedra;
std::vector<int> __tetrahedra_ref;
- typedef TinyVector<8,unsigned int> Hexahedron;
- Kokkos::View<Hexahedron*> __hexahedra;
+
+ using Hexahedron = TinyVector<8,unsigned int>;
+ Array<Hexahedron> __hexahedra;
std::vector<int> __hexahedra_ref;
/**
@@ -194,7 +199,7 @@ private:
* List of known keywords
*
*/
- typedef std::map<std::string, int> KeywordList;
+ using KeywordList = std::map<std::string, int>;
KeywordList __keywordList; /**< The keyword list */
@@ -202,7 +207,7 @@ private:
* Type for keyword
*
*/
- typedef std::pair<std::string, int> Keyword;
+ using Keyword = std::pair<std::string, int>;
/**
* Skips next comments if exists
diff --git a/src/mesh/IConnectivity.hpp b/src/mesh/IConnectivity.hpp
index 49432fdea276d770c993960a2ffb48da2f7407b0..e27c2ef08bf13e7e9946c095da1cb628027e3a44 100644
--- a/src/mesh/IConnectivity.hpp
+++ b/src/mesh/IConnectivity.hpp
@@ -4,15 +4,60 @@
#include <ItemType.hpp>
#include <ConnectivityMatrix.hpp>
-struct IConnectivity
+class IConnectivity
{
+ protected:
+ template <typename DataType,
+ ItemType sub_item_type,
+ ItemType item_type,
+ typename Allowed>
+ friend
+ class SubItemValuePerItem;
+
virtual const ConnectivityMatrix&
- getMatrix(const ItemType& item_type_0,
- const ItemType& item_type_1) const = 0;
+ _getMatrix(const ItemType& item_type_0,
+ const ItemType& item_type_1) const = 0;
+
+ public:
+ virtual size_t numberOfNodes() const = 0;
+ virtual size_t numberOfEdges() const = 0;
+ virtual size_t numberOfFaces() const = 0;
+ virtual size_t numberOfCells() const = 0;
+
+ template <ItemType item_type>
+ size_t numberOf() const = delete;
IConnectivity() = default;
IConnectivity(const IConnectivity&) = delete;
~IConnectivity() = default;
};
+template <>
+KOKKOS_INLINE_FUNCTION
+size_t IConnectivity::numberOf<ItemType::node>() const
+{
+ return this->numberOfNodes();
+}
+
+template <>
+KOKKOS_INLINE_FUNCTION
+size_t IConnectivity::numberOf<ItemType::edge>() const
+{
+ return this->numberOfEdges();
+}
+
+template <>
+KOKKOS_INLINE_FUNCTION
+size_t IConnectivity::numberOf<ItemType::face>() const
+{
+ return this->numberOfFaces();
+}
+
+template <>
+KOKKOS_INLINE_FUNCTION
+size_t IConnectivity::numberOf<ItemType::cell>() const
+{
+ return this->numberOfCells();
+}
+
#endif // ICONNECTIVITY_HPP
diff --git a/src/mesh/ItemId.hpp b/src/mesh/ItemId.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..11b87f7082bddbbbef76826bdc32b6272bc87d27
--- /dev/null
+++ b/src/mesh/ItemId.hpp
@@ -0,0 +1,81 @@
+#ifndef ITEM_ID_HPP
+#define ITEM_ID_HPP
+
+#include <ItemType.hpp>
+
+template <ItemType item_type>
+class ItemIdT
+{
+ public:
+ using base_type = unsigned int;
+
+ private:
+ base_type m_id;
+
+ public:
+ operator const base_type&() const
+ {
+ return m_id;
+ }
+
+ ItemIdT operator++(int)
+ {
+ ItemIdT item_id(m_id);
+ ++m_id;
+ return std::move(item_id);
+ }
+
+ ItemIdT& operator++()
+ {
+ ++m_id;
+ return *this;
+ }
+
+ ItemIdT& operator=(const base_type& id)
+ {
+ m_id = id;
+ return *this;
+ }
+
+ ItemIdT& operator=(const ItemIdT&) = default;
+ ItemIdT& operator=(ItemIdT&&) = default;
+
+ ItemIdT(const base_type& id) : m_id{id}{}
+
+ ItemIdT(const ItemIdT&) = default;
+
+ ItemIdT(ItemIdT&&) = default;
+ ItemIdT() = default;
+ ~ItemIdT() = default;
+
+ // forbidden rules
+ template <ItemType other_item_type>
+ ItemIdT(const ItemIdT<other_item_type>&)
+ {
+ static_assert(other_item_type == item_type,
+ "wrong type of item");
+ }
+
+ template <ItemType other_item_type>
+ ItemIdT& operator=(const ItemIdT<other_item_type>&)
+ {
+ static_assert(other_item_type == item_type,
+ "wrong type of item");
+ }
+
+
+ template <ItemType other_item_type>
+ ItemIdT& operator=(ItemIdT<other_item_type>&&)
+ {
+ static_assert(other_item_type == item_type,
+ "wrong type of item");
+ }
+
+};
+
+using NodeId = ItemIdT<ItemType::node>;
+using EdgeId = ItemIdT<ItemType::edge>;
+using FaceId = ItemIdT<ItemType::face>;
+using CellId = ItemIdT<ItemType::cell>;
+
+#endif // ITEM_ID_HPP
diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..510c421fc9c1a67d0565d4b5e5e3e08c91c4ccd9
--- /dev/null
+++ b/src/mesh/ItemToItemMatrix.hpp
@@ -0,0 +1,105 @@
+#ifndef ITEM_TO_ITEM_MATRIX_HPP
+#define ITEM_TO_ITEM_MATRIX_HPP
+
+#include <ItemType.hpp>
+#include <ItemId.hpp>
+
+#include <Kokkos_Core.hpp>
+#include <ConnectivityMatrix.hpp>
+
+template <ItemType SourceItemType,
+ ItemType TargetItemType>
+class ItemToItemMatrix
+{
+ public:
+ using SourceItemId = ItemIdT<SourceItemType>;
+ using TargetItemId = ItemIdT<TargetItemType>;
+
+ template <typename RowType>
+ class SubItemList
+ {
+ private:
+ const RowType m_row;
+
+ public:
+
+ KOKKOS_INLINE_FUNCTION
+ size_t size() const
+ {
+ return m_row.length;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TargetItemId operator[](const size_t& j) const
+ {
+ return m_row(j);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ SubItemList(const RowType& row)
+ : m_row{row}
+ {
+ ;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ SubItemList& operator=(const SubItemList&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ SubItemList& operator=(SubItemList&&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ SubItemList(const SubItemList&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ SubItemList(SubItemList&&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ~SubItemList() = default;
+ };
+
+ private:
+ const ConnectivityMatrix& m_connectivity_matrix;
+
+ public:
+ KOKKOS_INLINE_FUNCTION
+ const auto operator[](const SourceItemId& source_id) const
+ {
+ using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
+ return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
+ }
+
+ template <typename IndexType>
+ KOKKOS_INLINE_FUNCTION
+ const auto operator[](const IndexType& source_id) const
+ {
+ static_assert(std::is_same<IndexType, SourceItemId>(),
+ "ItemToItemMatrix must be indexed using correct ItemId");
+ using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
+ return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ ItemToItemMatrix(const ConnectivityMatrix& connectivity_matrix)
+ : m_connectivity_matrix{connectivity_matrix}
+ {
+ ;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ ItemToItemMatrix& operator=(const ItemToItemMatrix&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ItemToItemMatrix& operator=(ItemToItemMatrix&&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ItemToItemMatrix(ItemToItemMatrix&&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ItemToItemMatrix(const ItemToItemMatrix&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ~ItemToItemMatrix() = default;
+};
+
+#endif // ITEM_TO_ITEM_MATRIX_HPP
diff --git a/src/mesh/ItemType.hpp b/src/mesh/ItemType.hpp
index d37c02fae15120b872c008f939c5d4907debb085..965f888b1b7315f28f5132b0457b86ca1e85b210 100644
--- a/src/mesh/ItemType.hpp
+++ b/src/mesh/ItemType.hpp
@@ -39,12 +39,12 @@ std::string_view itemName(const ItemType& item_type)
}
template <size_t Dimension>
-struct ItemId {};
+struct ItemTypeId {};
template <>
-struct ItemId<1>
+struct ItemTypeId<1>
{
- inline static constexpr size_t itemId(const ItemType& item_type) {
+ inline static constexpr size_t itemTId(const ItemType& item_type) {
size_t i = std::numeric_limits<size_t>::max();
switch(item_type) {
case ItemType::cell: {
@@ -64,9 +64,9 @@ struct ItemId<1>
};
template <>
-struct ItemId<2>
+struct ItemTypeId<2>
{
- inline static constexpr size_t itemId(const ItemType& item_type) {
+ inline static constexpr size_t itemTId(const ItemType& item_type) {
size_t i = std::numeric_limits<size_t>::max();
switch(item_type) {
case ItemType::cell: {
@@ -89,9 +89,9 @@ struct ItemId<2>
};
template <>
-struct ItemId<3>
+struct ItemTypeId<3>
{
- inline static constexpr size_t itemId(const ItemType& item_type) {
+ inline static constexpr size_t itemTId(const ItemType& item_type) {
size_t i = std::numeric_limits<size_t>::max();
switch(item_type) {
case ItemType::cell: {
diff --git a/src/mesh/ItemValue.hpp b/src/mesh/ItemValue.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..77eb1cb37d9998fb1af6e4a27d05ff7fcdb60077
--- /dev/null
+++ b/src/mesh/ItemValue.hpp
@@ -0,0 +1,122 @@
+#ifndef ITEM_VALUE_HPP
+#define ITEM_VALUE_HPP
+
+#include <PastisAssert.hpp>
+
+#include <Array.hpp>
+
+#include <ItemType.hpp>
+#include <ItemId.hpp>
+
+#include <IConnectivity.hpp>
+
+template <typename DataType,
+ ItemType item_type>
+class ItemValue
+{
+ public:
+ static const ItemType item_t{item_type};
+ using data_type = DataType;
+
+ using ItemId = ItemIdT<item_type>;
+ using index_type = ItemId;
+
+ private:
+ bool m_is_built{false};
+
+ Array<DataType> m_values;
+
+ // Allows const version to access our data
+ friend ItemValue<std::add_const_t<DataType>,
+ item_type>;
+
+ public:
+ KOKKOS_FORCEINLINE_FUNCTION
+ const bool& isBuilt() const
+ {
+ return m_is_built;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ size_t size() const
+ {
+ Assert(m_is_built);
+ return m_values.size();
+ }
+
+ // Following Kokkos logic, these classes are view and const view does allow
+ // changes in data
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator[](const ItemId& i) const
+ {
+ Assert(m_is_built);
+ return m_values[i];
+ }
+
+ template <typename IndexType>
+ DataType& operator[](const IndexType& i) const
+ {
+ static_assert(std::is_same<IndexType,ItemId>(),
+ "ItemValue must be indexed by ItemId");
+ return m_values[i];
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfItems() const
+ {
+ Assert(m_is_built);
+ return m_values.size();
+ }
+
+ template <typename DataType2>
+ KOKKOS_INLINE_FUNCTION
+ ItemValue&
+ operator=(const ItemValue<DataType2, item_type>& value_per_item)
+ {
+ // ensures that DataType is the same as source DataType2
+ static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+ "Cannot assign ItemValue of different type");
+ // ensures that const is not lost through copy
+ static_assert(((std::is_const<DataType2>() and std::is_const<DataType>())
+ or not std::is_const<DataType2>()),
+ "Cannot assign ItemValue of const to ItemValue of non-const");
+
+ m_values = value_per_item.m_values;
+ m_is_built = value_per_item.m_is_built;
+
+ return *this;
+ }
+
+ template <typename DataType2>
+ KOKKOS_INLINE_FUNCTION
+ ItemValue(const ItemValue<DataType2, item_type>& value_per_item)
+ {
+ this->operator=(value_per_item);
+ }
+
+ ItemValue() = default;
+
+ ItemValue(const IConnectivity& connectivity)
+ : m_is_built{true},
+ m_values(connectivity.numberOf<item_type>())
+ {
+ static_assert(not std::is_const<DataType>(),
+ "Cannot allocate ItemValue of const data: only view is supported"); ;
+ }
+
+ ~ItemValue() = default;
+};
+
+template <typename DataType>
+using NodeValue = ItemValue<DataType, ItemType::node>;
+
+template <typename DataType>
+using EdgeValue = ItemValue<DataType, ItemType::edge>;
+
+template <typename DataType>
+using FaceValue = ItemValue<DataType, ItemType::face>;
+
+template <typename DataType>
+using CellValue = ItemValue<DataType, ItemType::cell>;
+
+#endif // ITEM_VALUE_HPP
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 3aee2c287d5cdd06c4a754115caf930973c4f2cb..0188d507cb2d463b899e3e318b5a75a002b152e4 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -1,7 +1,7 @@
#ifndef MESH_HPP
#define MESH_HPP
-#include <Kokkos_Core.hpp>
+#include <ItemValue.hpp>
#include <TinyVector.hpp>
#include <memory>
@@ -16,21 +16,24 @@ template <typename ConnectivityType>
class Mesh final : public IMesh
{
public:
- typedef ConnectivityType Connectivity;
+ using Connectivity = ConnectivityType;
+
static constexpr size_t dimension = ConnectivityType::dimension;
- typedef TinyVector<dimension> Rd;
+ using Rd = TinyVector<dimension>;
private:
const std::shared_ptr<Connectivity> m_connectivity;
- Kokkos::View<Rd*> m_xr;
+ NodeValue<const Rd> m_xr;
+ NodeValue<Rd> m_mutable_xr;
public:
-
+ KOKKOS_INLINE_FUNCTION
const size_t meshDimension() const
{
return dimension;
}
+ KOKKOS_INLINE_FUNCTION
const Connectivity& connectivity() const
{
return *m_connectivity;
@@ -54,28 +57,25 @@ public:
return m_connectivity->numberOfCells();
}
- Kokkos::View<Rd*> xr() const
+ KOKKOS_INLINE_FUNCTION
+ const NodeValue<const Rd>& xr() const
{
return m_xr;
}
+ [[deprecated("should rework this class: quite ugly")]]
KOKKOS_INLINE_FUNCTION
- Mesh(const std::shared_ptr<Connectivity>& connectivity)
- : m_connectivity(connectivity),
- m_xr("xr", connectivity->numberOfNodes())
+ NodeValue<Rd> mutableXr() const
{
- static_assert (Connectivity::dimension ==1,"no automatic calculation of vertices in dim>1");
- const double delta_x = 1./connectivity->numberOfCells();
- Kokkos::parallel_for(connectivity->numberOfNodes(), KOKKOS_LAMBDA(const int& r){
- m_xr[r][0] = r*delta_x;
- });
+ return m_mutable_xr;
}
KOKKOS_INLINE_FUNCTION
Mesh(const std::shared_ptr<Connectivity>& connectivity,
- Kokkos::View<Rd*>& xr)
- : m_connectivity(connectivity),
- m_xr(xr)
+ NodeValue<Rd>& xr)
+ : m_connectivity{connectivity},
+ m_xr{xr},
+ m_mutable_xr{xr}
{
;
}
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 0fe98256bb87d8d02b29899906e475737f4e835d..1a07f1bcc0e477a0141e0ab7d4f7065aad633336 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -4,6 +4,7 @@
#include <Kokkos_Core.hpp>
#include <TinyVector.hpp>
+#include <ItemValue.hpp>
#include <SubItemValuePerItem.hpp>
#include <map>
@@ -12,74 +13,77 @@ template <typename M>
class MeshData
{
public:
- typedef M MeshType;
+ using MeshType = M;
static constexpr size_t dimension = MeshType::dimension;
static_assert(dimension>0, "dimension must be strictly positive");
- typedef TinyVector<dimension> Rd;
+ using Rd = TinyVector<dimension>;
static constexpr double inv_dimension = 1./dimension;
private:
const MeshType& m_mesh;
- NodeValuePerCell<Rd> m_Cjr;
- NodeValuePerCell<double> m_ljr;
- NodeValuePerCell<Rd> m_njr;
- Kokkos::View<Rd*> m_xj;
- Kokkos::View<double*> m_Vj;
+ NodeValuePerCell<const Rd> m_Cjr;
+ NodeValuePerCell<const double> m_ljr;
+ NodeValuePerCell<const Rd> m_njr;
+ CellValue<const Rd> m_xj;
+ CellValue<const double> m_Vj;
KOKKOS_INLINE_FUNCTION
void _updateCenter()
{ // Computes vertices isobarycenter
if constexpr (dimension == 1) {
- const Kokkos::View<const Rd*> xr = m_mesh.xr();
+ const NodeValue<const Rd>& xr = m_mesh.xr();
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- m_xj[j] = 0.5*(xr[cell_nodes(0)]+xr[cell_nodes(1)]);
+ CellValue<Rd> xj(m_mesh.connectivity());
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ xj[j] = 0.5*(xr[cell_nodes[0]]+xr[cell_nodes[1]]);
});
+ m_xj = xj;
} else {
- const Kokkos::View<const Rd*> xr = m_mesh.xr();
- const Kokkos::View<const double*>& inv_cell_nb_nodes
+ const NodeValue<const Rd>& xr = m_mesh.xr();
+
+ const CellValue<const double>& inv_cell_nb_nodes
= m_mesh.connectivity().invCellNbNodes();
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
-
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ = m_mesh.connectivity().cellToNodeMatrix();
+ CellValue<Rd> xj(m_mesh.connectivity());
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
Rd X = zero;
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- for (size_t R=0; R<cell_nodes.length; ++R) {
- X += xr[cell_nodes(R)];
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ for (size_t R=0; R<cell_nodes.size(); ++R) {
+ X += xr[cell_nodes[R]];
}
- m_xj[j] = inv_cell_nb_nodes[j]*X;
+ xj[j] = inv_cell_nb_nodes[j]*X;
});
+ m_xj = xj;
}
}
KOKKOS_INLINE_FUNCTION
void _updateVolume()
{
- const Kokkos::View<const Rd*> xr = m_mesh.xr();
+ const NodeValue<const Rd>& xr = m_mesh.xr();
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ CellValue<double> Vj(m_mesh.connectivity());
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
double sum_cjr_xr = 0;
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix[j];
- for (size_t R=0; R<cell_nodes.length; ++R) {
- sum_cjr_xr += (xr[cell_nodes(R)], m_Cjr(j,R));
+ for (size_t R=0; R<cell_nodes.size(); ++R) {
+ sum_cjr_xr += (xr[cell_nodes[R]], m_Cjr(j,R));
}
- m_Vj[j] = inv_dimension * sum_cjr_xr;
+ Vj[j] = inv_dimension * sum_cjr_xr;
});
+ m_Vj = Vj;
}
KOKKOS_INLINE_FUNCTION
@@ -88,117 +92,136 @@ class MeshData
// Cjr/njr/ljr are constant overtime
}
else if constexpr (dimension == 2) {
- const Kokkos::View<const Rd*> xr = m_mesh.xr();
+ const NodeValue<const Rd>& xr = m_mesh.xr();
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
-
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- for (size_t R=0; R<cell_nodes.length; ++R) {
- int Rp1 = (R+1)%cell_nodes.length;
- int Rm1 = (R+cell_nodes.length-1)%cell_nodes.length;
- Rd half_xrp_xrm = 0.5*(xr(cell_nodes(Rp1))-xr(cell_nodes(Rm1)));
- m_Cjr(j,R) = Rd{-half_xrp_xrm[1], half_xrp_xrm[0]};
- }
- });
-
- const NodeValuePerCell<Rd>& Cjr = m_Cjr;
- Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){
- m_ljr[j] = l2Norm(Cjr[j]);
- });
-
- const NodeValuePerCell<double>& ljr = m_ljr;
- Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){
- m_njr[j] = (1./ljr[j])*Cjr[j];
- });
-
+ = m_mesh.connectivity().cellToNodeMatrix();
+
+ {
+ NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ for (size_t R=0; R<cell_nodes.size(); ++R) {
+ int Rp1 = (R+1)%cell_nodes.size();
+ int Rm1 = (R+cell_nodes.size()-1)%cell_nodes.size();
+ Rd half_xrp_xrm = 0.5*(xr[cell_nodes[Rp1]]-xr[cell_nodes[Rm1]]);
+ Cjr(j,R) = Rd{-half_xrp_xrm[1], half_xrp_xrm[0]};
+ }
+ });
+ m_Cjr = Cjr;
+ }
+
+ {
+ NodeValuePerCell<double> ljr(m_mesh.connectivity());
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const size_t& jr){
+ ljr[jr] = l2Norm(m_Cjr[jr]);
+ });
+ m_ljr = ljr;
+ }
+
+ {
+ NodeValuePerCell<Rd> njr(m_mesh.connectivity());
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const size_t& jr){
+ njr[jr] = (1./m_ljr[jr])*m_Cjr[jr];
+ });
+ m_njr = njr;
+ }
} else if (dimension ==3) {
- const Kokkos::View<const Rd*> xr = m_mesh.xr();
+ const NodeValue<const Rd>& xr = m_mesh.xr();
NodeValuePerFace<Rd> Nlr(m_mesh.connectivity());
const auto& face_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::face,
- ItemType::node);
+ = m_mesh.connectivity().faceToNodeMatrix();
- Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
- const auto& face_nodes = face_to_node_matrix.rowConst(l);
- const size_t nb_nodes = face_nodes.length;
+ Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const FaceId& l) {
+ const auto& face_nodes = face_to_node_matrix[l];
+ const size_t nb_nodes = face_nodes.size();
std::vector<Rd> dxr(nb_nodes);
for (size_t r=0; r<nb_nodes; ++r) {
- dxr[r] = xr[face_nodes((r+1)%nb_nodes)] - xr[face_nodes((r+nb_nodes-1)%nb_nodes)];
+ dxr[r]
+ = xr[face_nodes[(r+1)%nb_nodes]]
+ - xr[face_nodes[(r+nb_nodes-1)%nb_nodes]];
}
const double inv_12_nb_nodes = 1./(12.*nb_nodes);
for (size_t r=0; r<nb_nodes; ++r) {
Rd Nr = zero;
const Rd two_dxr = 2*dxr[r];
for (size_t s=0; s<nb_nodes; ++s) {
- Nr += crossProduct((two_dxr - dxr[s]), xr[face_nodes(s)]);
+ Nr += crossProduct((two_dxr - dxr[s]), xr[face_nodes[s]]);
}
Nr *= inv_12_nb_nodes;
- Nr -= (1./6.)*crossProduct(dxr[r], xr[face_nodes(r)]);
+ Nr -= (1./6.)*crossProduct(dxr[r], xr[face_nodes[r]]);
Nlr(l,r) = Nr;
}
});
- Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
- m_Cjr[jr] = zero;
- });
-
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
const auto& cell_to_face_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::face);
+ = m_mesh.connectivity().cellToFaceMatrix();
const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- const auto& cell_faces = cell_to_face_matrix.rowConst(j);
- const auto& face_is_reversed = cell_face_is_reversed.itemValues(j);
+ {
+ NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
+ Kokkos::parallel_for(Cjr.numberOfValues(), KOKKOS_LAMBDA(const size_t& jr){
+ Cjr[jr] = zero;
+ });
+
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ const auto& cell_faces = cell_to_face_matrix[j];
+ const auto& face_is_reversed = cell_face_is_reversed.itemValues(j);
- for (size_t L=0; L<cell_faces.length; ++L) {
- const size_t l = cell_faces(L);
- const auto& face_nodes = face_to_node_matrix.rowConst(l);
+ for (size_t L=0; L<cell_faces.size(); ++L) {
+ const FaceId& l = cell_faces[L];
+ const auto& face_nodes = face_to_node_matrix[l];
#warning should this lambda be replaced by a precomputed correspondance?
- std::function local_node_number_in_cell
- = [&](const size_t& node_number) {
- for (size_t i_node=0; i_node<cell_nodes.length; ++i_node) {
- if (node_number == cell_nodes(i_node)) {
- return i_node;
- break;
+ std::function local_node_number_in_cell
+ = [&](const NodeId& node_number) {
+ for (size_t i_node=0; i_node<cell_nodes.size(); ++i_node) {
+ if (node_number == cell_nodes[i_node]) {
+ return i_node;
+ break;
+ }
}
- }
- return std::numeric_limits<size_t>::max();
- };
-
- if (face_is_reversed[L]) {
- for (size_t rl = 0; rl<face_nodes.length; ++rl) {
- const size_t R = local_node_number_in_cell(face_nodes(rl));
- m_Cjr(j, R) -= Nlr(l,rl);
- }
- } else {
- for (size_t rl = 0; rl<face_nodes.length; ++rl) {
- const size_t R = local_node_number_in_cell(face_nodes(rl));
- m_Cjr(j, R) += Nlr(l,rl);
+ return std::numeric_limits<size_t>::max();
+ };
+
+ if (face_is_reversed[L]) {
+ for (size_t rl = 0; rl<face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
+ Cjr(j, R) -= Nlr(l,rl);
+ }
+ } else {
+ for (size_t rl = 0; rl<face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
+ Cjr(j, R) += Nlr(l,rl);
+ }
}
}
- }
- });
-
- const NodeValuePerCell<Rd>& Cjr = m_Cjr;
- Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
- m_ljr[jr] = l2Norm(Cjr[jr]);
- });
-
- const NodeValuePerCell<double>& ljr = m_ljr;
- Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
- m_njr[jr] = (1./ljr[jr])*Cjr[jr];
- });
+ });
+
+ m_Cjr = Cjr;
+ }
+
+ {
+ NodeValuePerCell<double> ljr(m_mesh.connectivity());
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const size_t& jr){
+ ljr[jr] = l2Norm(m_Cjr[jr]);
+ });
+ m_ljr = ljr;
+ }
+
+ {
+ NodeValuePerCell<Rd> njr(m_mesh.connectivity());
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const size_t& jr){
+ njr[jr] = (1./m_ljr[jr])*m_Cjr[jr];
+ });
+ m_njr = njr;
+ }
}
static_assert((dimension<=3), "only 1d, 2d and 3d are implemented");
}
@@ -209,27 +232,27 @@ class MeshData
return m_mesh;
}
- const NodeValuePerCell<Rd>& Cjr() const
+ const NodeValuePerCell<const Rd>& Cjr() const
{
return m_Cjr;
}
- const NodeValuePerCell<double>& ljr() const
+ const NodeValuePerCell<const double>& ljr() const
{
return m_ljr;
}
- const NodeValuePerCell<Rd>& njr() const
+ const NodeValuePerCell<const Rd>& njr() const
{
return m_njr;
}
- const Kokkos::View<const Rd*> xj() const
+ const CellValue<const Rd>& xj() const
{
return m_xj;
}
- const Kokkos::View<const double*> Vj() const
+ const CellValue<const double>& Vj() const
{
return m_Vj;
}
@@ -242,25 +265,27 @@ class MeshData
}
MeshData(const MeshType& mesh)
- : m_mesh(mesh),
- m_Cjr(mesh.connectivity()),
- m_ljr(mesh.connectivity()),
- m_njr(mesh.connectivity()),
- m_xj("xj", mesh.numberOfCells()),
- m_Vj("Vj", mesh.numberOfCells())
+ : m_mesh(mesh)
{
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;
- m_Cjr(j,1)= 1;
- });
+ {
+ NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
+ Cjr(j,0)=-1;
+ Cjr(j,1)= 1;
+ });
+ m_Cjr = Cjr;
+ }
// in 1d njr=Cjr (here performs a shallow copy)
- 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;
+ m_njr = m_Cjr;
+ {
+ NodeValuePerCell<double> ljr(m_mesh.connectivity());
+ Kokkos::parallel_for(ljr.numberOfValues(), KOKKOS_LAMBDA(const size_t& jr){
+ ljr[jr] = 1;
});
+ m_ljr = ljr;
+ }
}
this->updateAllData();
}
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 87c537439bf1a27e8443865677cd6b58c132378c..0a4b9ec63ff74b71e25342144871c34e3236c163 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -1,7 +1,9 @@
#ifndef MESH_NODE_BOUNDARY_HPP
#define MESH_NODE_BOUNDARY_HPP
-#include <Kokkos_Core.hpp>
+#include <Array.hpp>
+#include <ItemValue.hpp>
+
#include <Kokkos_Vector.hpp>
#include <TinyVector.hpp>
@@ -16,12 +18,12 @@ template <size_t dimension>
class MeshNodeBoundary
{
protected:
- Kokkos::View<const unsigned int*> m_node_list;
+ Array<const NodeId> m_node_list;
public:
MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;
- const Kokkos::View<const unsigned int*>& nodeList() const
+ const Array<const NodeId>& nodeList() const
{
return m_node_list;
}
@@ -32,12 +34,12 @@ class MeshNodeBoundary
{
static_assert(dimension == MeshType::dimension);
const auto& face_to_cell_matrix
- = mesh.connectivity().getMatrix(ItemType::face, ItemType::cell);
+ = mesh.connectivity().faceToCellMatrix();
- const Kokkos::View<const unsigned int*>& face_list = ref_face_list.faceList();
- Kokkos::parallel_for(face_list.extent(0), KOKKOS_LAMBDA(const int& l){
- const auto& face_cells = face_to_cell_matrix.rowConst(face_list[l]);
- if (face_cells.length>1) {
+ const Array<const FaceId>& face_list = ref_face_list.faceList();
+ Kokkos::parallel_for(face_list.size(), KOKKOS_LAMBDA(const int& l){
+ const auto& face_cells = face_to_cell_matrix[face_list[l]];
+ if (face_cells.size()>1) {
std::cerr << "internal faces cannot be used to define mesh boundaries\n";
std::exit(1);
}
@@ -45,24 +47,23 @@ class MeshNodeBoundary
Kokkos::vector<unsigned int> node_ids;
// not enough but should reduce significantly the number of resizing
- node_ids.reserve(dimension*face_list.extent(0));
+ node_ids.reserve(dimension*face_list.size());
const auto& face_to_node_matrix
- = mesh.connectivity().getMatrix(ItemType::face,
- ItemType::node);
+ = mesh.connectivity().faceToNodeMatrix();
- for (size_t l=0; l<face_list.extent(0); ++l) {
- const size_t face_number = face_list[l];
- const auto& face_nodes = face_to_node_matrix.rowConst(face_number);
+ for (size_t l=0; l<face_list.size(); ++l) {
+ const FaceId face_number = face_list[l];
+ const auto& face_nodes = face_to_node_matrix[face_number];
- for (size_t r=0; r<face_nodes.length; ++r) {
- node_ids.push_back(face_nodes(r));
+ for (size_t r=0; r<face_nodes.size(); ++r) {
+ node_ids.push_back(face_nodes[r]);
}
}
std::sort(node_ids.begin(), node_ids.end());
auto last = std::unique(node_ids.begin(), node_ids.end());
node_ids.resize(std::distance(node_ids.begin(), last));
- Kokkos::View<unsigned int*> node_list("node_list", node_ids.size());
+ Array<NodeId> node_list(node_ids.size());
Kokkos::parallel_for(node_ids.size(), KOKKOS_LAMBDA(const int& r){
node_list[r] = node_ids[r];
});
@@ -88,7 +89,9 @@ template <size_t dimension>
class MeshFlatNodeBoundary
: public MeshNodeBoundary<dimension>
{
- typedef TinyVector<dimension, double> Rd;
+ public:
+ using Rd = TinyVector<dimension, double>;
+
private:
const Rd m_outgoing_normal;
@@ -145,13 +148,14 @@ _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
const MeshType& mesh) const
{
static_assert(MeshType::dimension == 2);
- typedef TinyVector<2,double> R2;
+ using R2 = TinyVector<2,double>;
+
const R2 origin = 0.5*(xmin+xmax);
const double length = l2Norm(xmax-xmin);
- const Kokkos::View<const R2*> xr = mesh.xr();
+ const NodeValue<const R2>& xr = mesh.xr();
- Kokkos::parallel_for(m_node_list.extent(0), KOKKOS_LAMBDA(const size_t& r) {
+ Kokkos::parallel_for(m_node_list.size(), KOKKOS_LAMBDA(const size_t& r) {
const R2& x = xr[m_node_list[r]];
if ((x-origin,normal)>1E-13*length) {
std::cerr << "this FlatBoundary is not flat!\n";
@@ -167,9 +171,9 @@ MeshFlatNodeBoundary<1>::
_getNormal(const MeshType& mesh)
{
static_assert(MeshType::dimension == 1);
- typedef TinyVector<1,double> R;
+ using R = TinyVector<1,double>;
- if (m_node_list.extent(0) != 1) {
+ if (m_node_list.size() != 1) {
std::cerr << "Node boundaries in 1D require to have exactly 1 node\n";
std::exit(1);
}
@@ -184,9 +188,9 @@ MeshFlatNodeBoundary<2>::
_getNormal(const MeshType& mesh)
{
static_assert(MeshType::dimension == 2);
- typedef TinyVector<2,double> R2;
+ using R2 = TinyVector<2,double>;
- const Kokkos::View<const R2*> xr = mesh.xr();
+ const NodeValue<const R2>& xr = mesh.xr();
R2 xmin(std::numeric_limits<double>::max(),
std::numeric_limits<double>::max());
@@ -194,7 +198,7 @@ _getNormal(const MeshType& mesh)
R2 xmax(-std::numeric_limits<double>::max(),
-std::numeric_limits<double>::max());
- for (size_t r=0; r<m_node_list.extent(0); ++r) {
+ for (size_t r=0; r<m_node_list.size(); ++r) {
const R2& x = xr[m_node_list[r]];
if ((x[0]<xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
xmin = x;
@@ -226,7 +230,7 @@ MeshFlatNodeBoundary<3>::
_getNormal(const MeshType& mesh)
{
static_assert(MeshType::dimension == 3);
- typedef TinyVector<3,double> R3;
+ using R3 = TinyVector<3,double>;
R3 xmin(std::numeric_limits<double>::max(),
@@ -239,9 +243,9 @@ _getNormal(const MeshType& mesh)
R3 ymax = xmax;
R3 zmax = xmax;
- const Kokkos::View<const R3*> xr = mesh.xr();
+ const NodeValue<const R3>& xr = mesh.xr();
- for (size_t r=0; r<m_node_list.extent(0); ++r) {
+ for (size_t r=0; r<m_node_list.size(); ++r) {
const R3& x = xr[m_node_list[r]];
if (x[0]<xmin[0]) {
xmin = x;
@@ -310,25 +314,23 @@ MeshFlatNodeBoundary<1>::
_getOutgoingNormal(const MeshType& mesh)
{
static_assert(MeshType::dimension == 1);
- typedef TinyVector<1,double> R;
+ using R = TinyVector<1,double>;
const R normal = this->_getNormal(mesh);
- const Kokkos::View<const R*>& xr = mesh.xr();
+ const NodeValue<const R>& xr = mesh.xr();
const auto& cell_to_node_matrix
- = mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = mesh.connectivity().cellToNodeMatrix();
const auto& node_to_cell_matrix
- = mesh.connectivity().getMatrix(ItemType::node,
- ItemType::cell);
+ = mesh.connectivity().nodeToCellMatrix();
- const size_t r0 = m_node_list[0];
- const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
- const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
+ const NodeId r0 = m_node_list[0];
+ const CellId j0 = node_to_cell_matrix[r0][0];
+ const auto& j0_nodes = cell_to_node_matrix[j0];
double max_height = 0;
- for (size_t r=0; r<j0_nodes.length; ++r) {
- const double height = (xr(j0_nodes(r))-xr(r0), normal);
+ for (size_t r=0; r<j0_nodes.size(); ++r) {
+ const double height = (xr[j0_nodes[r]]-xr[r0], normal);
if (std::abs(height) > std::abs(max_height)) {
max_height = height;
}
@@ -347,25 +349,23 @@ MeshFlatNodeBoundary<2>::
_getOutgoingNormal(const MeshType& mesh)
{
static_assert(MeshType::dimension == 2);
- typedef TinyVector<2,double> R2;
+ using R2 = TinyVector<2,double>;
const R2 normal = this->_getNormal(mesh);
- const Kokkos::View<const R2*>& xr = mesh.xr();
+ const NodeValue<const R2>& xr = mesh.xr();
const auto& cell_to_node_matrix
- = mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = mesh.connectivity().cellToNodeMatrix();
const auto& node_to_cell_matrix
- = mesh.connectivity().getMatrix(ItemType::node,
- ItemType::cell);
+ = mesh.connectivity().nodeToCellMatrix();
- const size_t r0 = m_node_list[0];
- const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
- const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
+ const NodeId r0 = m_node_list[0];
+ const CellId j0 = node_to_cell_matrix[r0][0];
+ const auto& j0_nodes = cell_to_node_matrix[j0];
double max_height = 0;
- for (size_t r=0; r<j0_nodes.length; ++r) {
- const double height = (xr(j0_nodes(r))-xr(r0), normal);
+ for (size_t r=0; r<j0_nodes.size(); ++r) {
+ const double height = (xr[j0_nodes[r]]-xr[r0], normal);
if (std::abs(height) > std::abs(max_height)) {
max_height = height;
}
@@ -384,25 +384,23 @@ MeshFlatNodeBoundary<3>::
_getOutgoingNormal(const MeshType& mesh)
{
static_assert(MeshType::dimension == 3);
- typedef TinyVector<3,double> R3;
+ using R3 = TinyVector<3,double>;
const R3 normal = this->_getNormal(mesh);
- const Kokkos::View<const R3*>& xr = mesh.xr();
+ const NodeValue<const R3>& xr = mesh.xr();
const auto& cell_to_node_matrix
- = mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = mesh.connectivity().cellToNodeMatrix();
const auto& node_to_cell_matrix
- = mesh.connectivity().getMatrix(ItemType::node,
- ItemType::cell);
+ = mesh.connectivity().nodeToCellMatrix();
- const size_t r0 = m_node_list[0];
- const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
- const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
+ const NodeId r0 = m_node_list[0];
+ const CellId j0 = node_to_cell_matrix[r0][0];
+ const auto& j0_nodes = cell_to_node_matrix[j0];
double max_height = 0;
- for (size_t r=0; r<j0_nodes.length; ++r) {
- const double height = (xr(j0_nodes(r))-xr(r0), normal);
+ for (size_t r=0; r<j0_nodes.size(); ++r) {
+ const double height = (xr[j0_nodes[r]]-xr[r0], normal);
if (std::abs(height) > std::abs(max_height)) {
max_height = height;
}
diff --git a/src/mesh/RefFaceList.hpp b/src/mesh/RefFaceList.hpp
index e891ea2581bc1dbb19a0e82f4178a2a172112574..ab66943936c091735c404b287d587a15d81a1a9f 100644
--- a/src/mesh/RefFaceList.hpp
+++ b/src/mesh/RefFaceList.hpp
@@ -1,14 +1,14 @@
#ifndef REF_FACE_LIST_HPP
#define REF_FACE_LIST_HPP
-#include <Kokkos_Core.hpp>
+#include <Array.hpp>
#include <RefId.hpp>
class RefFaceList
{
private:
const RefId m_ref_id;
- const Kokkos::View<const unsigned int*> m_face_id_list;
+ const Array<const FaceId> m_face_id_list;
public:
const RefId& refId() const
@@ -16,13 +16,13 @@ class RefFaceList
return m_ref_id;
}
- const Kokkos::View<const unsigned int*>& faceList() const
+ const Array<const FaceId>& faceList() const
{
return m_face_id_list;
}
RefFaceList(const RefId& ref_id,
- const Kokkos::View<const unsigned int*>& face_id_list)
+ const Array<const FaceId>& face_id_list)
: m_ref_id(ref_id),
m_face_id_list(face_id_list)
{
diff --git a/src/mesh/RefNodeList.hpp b/src/mesh/RefNodeList.hpp
index 7bdc36fac84b1dd9cb713be5a5a64e9d992232c0..462d23bb3d35a6384fcd5176d0584419a586abc3 100644
--- a/src/mesh/RefNodeList.hpp
+++ b/src/mesh/RefNodeList.hpp
@@ -1,14 +1,14 @@
#ifndef REF_NODE_LIST_HPP
#define REF_NODE_LIST_HPP
-#include <Kokkos_Core.hpp>
+#include <Array.hpp>
#include <RefId.hpp>
class RefNodeList
{
private:
RefId m_ref_id;
- Kokkos::View<const unsigned int*> m_node_id_list;
+ Array<const NodeId> m_node_id_list;
public:
const RefId& refId() const
@@ -16,13 +16,13 @@ class RefNodeList
return m_ref_id;
}
- const Kokkos::View<const unsigned int*>& nodeList() const
+ const Array<const NodeId>& nodeList() const
{
return m_node_id_list;
}
RefNodeList(const RefId& ref_id,
- const Kokkos::View<const unsigned int*>& node_id_list)
+ const Array<const NodeId>& node_id_list)
: m_ref_id(ref_id),
m_node_id_list(node_id_list)
{
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index 5778fc73989ca547130288d1737bf9b1d767c2ef..3853038a1822ccd5bc6e4eb6d5c2d2696c8876ba 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -2,9 +2,14 @@
#define SUBITEM_VALUE_PER_ITEM_HPP
#include <Kokkos_StaticCrsGraph.hpp>
-#include <ItemType.hpp>
+
#include <PastisAssert.hpp>
+#include <Array.hpp>
+#include <ItemType.hpp>
+
+#include <ItemId.hpp>
+
#include <ConnectivityMatrix.hpp>
#include <IConnectivity.hpp>
@@ -25,11 +30,16 @@ class SubItemValuePerItem<DataType,
public:
static const ItemType item_t{item_type};
static const ItemType sub_item_t{sub_item_type};
+
+ using data_type = DataType;
+ using ItemId = ItemIdT<item_type>;
+ using index_type = ItemId;
+
private:
bool m_is_built{false};
- ConnectivityMatrix::HostRowType m_host_row_map;
- Kokkos::View<DataType*> m_values;
+ typename ConnectivityMatrix::HostRowType m_host_row_map;
+ Array<DataType> m_values;
// Allows const version to access our data
friend SubItemValuePerItem<std::add_const_t<DataType>,
@@ -39,9 +49,13 @@ class SubItemValuePerItem<DataType,
public:
class SubView
{
+ public:
+ using data_type = DataType;
+
private:
KOKKOS_RESTRICT DataType* const m_sub_values;
const size_t m_size;
+
public:
KOKKOS_INLINE_FUNCTION
const DataType& operator[](const size_t& i) const
@@ -69,14 +83,14 @@ class SubItemValuePerItem<DataType,
SubView(SubView&&) = default;
KOKKOS_INLINE_FUNCTION
- SubView(const Kokkos::View<DataType*>& values,
+ SubView(const Array<DataType>& values,
const size_t& begin,
const size_t& end)
: m_sub_values(&(values[begin])),
m_size(end-begin)
{
Assert(begin<=end);
- Assert(end<=values.extent(0));
+ Assert(end<=values.size());
}
};
@@ -86,42 +100,45 @@ class SubItemValuePerItem<DataType,
return m_is_built;
}
+ // Following Kokkos logic, these classes are view and const view does allow
+ // changes in data
KOKKOS_FORCEINLINE_FUNCTION
- DataType& operator()(const size_t& j, const size_t& r)
+ DataType& operator()(const ItemId& j, const size_t& r) const
{
Assert(m_is_built);
- return m_values[m_host_row_map(j)+r];
+ return m_values[m_host_row_map(size_t{j})+r];
}
- // Following Kokkos logic, these classes are view and const view does allow
- // changes in data
+ template <typename IndexType>
KOKKOS_FORCEINLINE_FUNCTION
- DataType& operator()(const size_t& j, const size_t& r) const
+ DataType& operator()(const IndexType& j, const size_t& r) const
{
- Assert(m_is_built);
- return m_values[m_host_row_map(j)+r];
+ static_assert(std::is_same<IndexType, size_t>(),
+ "SubItemValuePerItem indexed by ItemId");
+ return m_values[m_host_row_map(size_t{j})+r];
}
KOKKOS_INLINE_FUNCTION
size_t numberOfValues() const
{
Assert(m_is_built);
- return m_values.extent(0);
+ return m_values.size();
}
+ // Following Kokkos logic, these classes are view and const view does allow
+ // changes in data
KOKKOS_FORCEINLINE_FUNCTION
- DataType& operator[](const size_t& i)
+ DataType& operator[](const size_t& i) const
{
Assert(m_is_built);
return m_values[i];
}
- // Following Kokkos logic, these classes are view and const view does allow
- // changes in data
- KOKKOS_FORCEINLINE_FUNCTION
- DataType& operator[](const size_t & i) const
+ template <typename IndexType>
+ DataType& operator[](const IndexType& i) const
{
- Assert(m_is_built);
+ static_assert(std::is_same<IndexType, size_t>(),
+ "Access to SubItemValuePerItem's array must be indexed by size_t");
return m_values[i];
}
@@ -171,7 +188,7 @@ class SubItemValuePerItem<DataType,
// ensures that const is not lost through copy
static_assert(((std::is_const<DataType2>() and std::is_const<DataType>())
or not std::is_const<DataType2>()),
- "Cannot assign const SubItemValuePerItem to a non const SubItemValuePerItem");
+ "Cannot assign SubItemValuePerItem of const data to SubItemValuePerItem of non-const data");
m_host_row_map = sub_item_value_per_item.m_host_row_map;
m_values = sub_item_value_per_item.m_values;
@@ -192,11 +209,14 @@ class SubItemValuePerItem<DataType,
SubItemValuePerItem(const IConnectivity& connectivity)
: m_is_built{true}
{
+ static_assert(not std::is_const<DataType>(),
+ "Cannot allocate SubItemValuePerItem of const data: only view is supported"); ;
+
ConnectivityMatrix connectivity_matrix
- = connectivity.getMatrix(item_type, sub_item_type);
+ = connectivity._getMatrix(item_type, sub_item_type);
m_host_row_map = connectivity_matrix.rowsMap();
- m_values = Kokkos::View<DataType*>("values", connectivity_matrix.numEntries());
+ m_values = Array<std::remove_const_t<DataType>>(connectivity_matrix.numEntries());
}
~SubItemValuePerItem() = default;
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index a33fef0bfbe2abf2c73b34e7c3b564fb2cd81461..276bad33d5624f018b9d512743f4f7695a6213e5 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -8,6 +8,8 @@
#include <TinyVector.hpp>
#include <IConnectivity.hpp>
+#include <ItemValue.hpp>
+
class VTKWriter
{
private:
@@ -39,25 +41,24 @@ class VTKWriter
fout << "<Points>\n";
fout << "<DataArray Name=\"Positions\" NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">\n";
- if constexpr(MeshType::dimension ==1) {
- const Kokkos::View<const TinyVector<1>*> xr = mesh.xr();
- for (unsigned int r=0; r<mesh.numberOfNodes(); ++r) {
+ using Rd = TinyVector<MeshType::dimension>;
+ const NodeValue<const Rd>& xr = mesh.xr();
+ if constexpr(MeshType::dimension == 1) {
+ for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
for (unsigned short i=0; i<1; ++i) {
fout << xr[r][i] << ' ';
}
fout << "0 0 "; // VTK requires 3 components
}
- } else if constexpr (MeshType::dimension ==2) {
- const Kokkos::View<const TinyVector<2>*> xr = mesh.xr();
- for (unsigned int r=0; r<mesh.numberOfNodes(); ++r) {
+ } else if constexpr (MeshType::dimension == 2) {
+ for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
for (unsigned short i=0; i<2; ++i) {
fout << xr[r][i] << ' ';
}
fout << "0 "; // VTK requires 3 components
}
} else {
- const Kokkos::View<const TinyVector<3>*> xr = mesh.xr();
- for (unsigned int r=0; r<mesh.numberOfNodes(); ++r) {
+ for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
for (unsigned short i=0; i<3; ++i) {
fout << xr[r][i] << ' ';
}
@@ -72,13 +73,12 @@ class VTKWriter
fout << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n";
const auto& cell_to_node_matrix
- = mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = mesh.connectivity().cellToNodeMatrix();
- for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- for (unsigned short r=0; r<cell_nodes.length; ++r) {
- fout << cell_nodes(r) << ' ';
+ for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ for (unsigned short r=0; r<cell_nodes.size(); ++r) {
+ fout << cell_nodes[r] << ' ';
}
}
fout << '\n';
@@ -87,9 +87,9 @@ class VTKWriter
fout << "<DataArray type=\"UInt32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">\n";
{
unsigned int offset=0;
- for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- offset += cell_nodes.length;
+ for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ offset += cell_nodes.size();
fout << offset << ' ';
}
}
@@ -97,9 +97,9 @@ class VTKWriter
fout << "</DataArray>\n";
fout << "<DataArray type=\"Int8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">\n";
- for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- switch (cell_nodes.length) {
+ for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ switch (cell_nodes.size()) {
case 2: {
fout << "3 ";
break;
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 0cc6ec0df7a97b91999376873ed7b35a0c5acf6c..375828f41a86eef02f1d7fe87e48bb6ff4388223 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -1,9 +1,10 @@
#ifndef ACOUSTIC_SOLVER_HPP
#define ACOUSTIC_SOLVER_HPP
-#include <Kokkos_Core.hpp>
#include <rang.hpp>
+#include <ArrayUtils.hpp>
+
#include <BlockPerfectGas.hpp>
#include <PastisAssert.hpp>
@@ -19,8 +20,8 @@
template<typename MeshData>
class AcousticSolver
{
- typedef typename MeshData::MeshType MeshType;
- typedef FiniteVolumesEulerUnknowns<MeshData> UnknownsType;
+ using MeshType = typename MeshData::MeshType;
+ using UnknownsType = FiniteVolumesEulerUnknowns<MeshData>;
MeshData& m_mesh_data;
const MeshType& m_mesh;
@@ -29,66 +30,30 @@ class AcousticSolver
constexpr static size_t dimension = MeshType::dimension;
- typedef TinyVector<dimension> Rd;
- typedef TinyMatrix<dimension> Rdd;
+ using Rd = TinyVector<dimension>;
+ using Rdd = TinyMatrix<dimension>;
private:
- struct ReduceMin
- {
- private:
- const Kokkos::View<const double*> x_;
-
- public:
- typedef Kokkos::View<const double*>::non_const_value_type value_type;
-
- ReduceMin(const Kokkos::View<const double*>& x) : x_ (x) {}
-
- typedef Kokkos::View<const double*>::size_type size_type;
-
- KOKKOS_INLINE_FUNCTION
- void operator() (const size_type i, value_type& update) const
- {
- if (x_(i) < update) {
- update = x_(i);
- }
- }
-
- KOKKOS_INLINE_FUNCTION
- void join (volatile value_type& dst,
- const volatile value_type& src) const
- {
- if (src < dst) {
- dst = src;
- }
- }
-
- KOKKOS_INLINE_FUNCTION void
- init (value_type& dst) const
- { // The identity under max is -Inf.
- dst= Kokkos::reduction_identity<value_type>::min();
- }
- };
-
KOKKOS_INLINE_FUNCTION
- const Kokkos::View<const double*>
- computeRhoCj(const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& cj)
+ const CellValue<const double>
+ computeRhoCj(const CellValue<const double>& rhoj,
+ const CellValue<const double>& cj)
{
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
m_rhocj[j] = rhoj[j]*cj[j];
});
return m_rhocj;
}
KOKKOS_INLINE_FUNCTION
- void computeAjr(const Kokkos::View<const double*>& rhocj,
- const NodeValuePerCell<Rd>& Cjr,
- const NodeValuePerCell<double>& ljr,
- const NodeValuePerCell<Rd>& njr)
+ void computeAjr(const CellValue<const double>& rhocj,
+ const NodeValuePerCell<const Rd>& Cjr,
+ const NodeValuePerCell<const double>& ljr,
+ const NodeValuePerCell<const Rd>& njr)
{
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
const size_t& nb_nodes =m_Ajr.numberOfSubValues(j);
- const double& rho_c = rhocj(j);
+ const double& rho_c = rhocj[j];
for (size_t r=0; r<nb_nodes; ++r) {
m_Ajr(j,r) = tensorProduct(rho_c*Cjr(j,r), njr(j,r));
}
@@ -96,53 +61,51 @@ class AcousticSolver
}
KOKKOS_INLINE_FUNCTION
- const Kokkos::View<const Rdd*>
- computeAr(const NodeValuePerCell<Rdd>& Ajr) {
+ const NodeValue<const Rdd>
+ computeAr(const NodeValuePerCell<const Rdd>& Ajr) {
const auto& node_to_cell_matrix
- = m_connectivity.getMatrix(ItemType::node,
- ItemType::cell);
+ = m_connectivity.nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells
= m_connectivity.nodeLocalNumbersInTheirCells();
- Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r) {
Rdd sum = zero;
- const auto& node_to_cell = node_to_cell_matrix.rowConst(r);
+ const auto& node_to_cell = node_to_cell_matrix[r];
const auto& node_local_number_in_its_cells
= node_local_numbers_in_their_cells.itemValues(r);
- for (size_t j=0; j<node_to_cell.length; ++j) {
- const unsigned int J = node_to_cell(j);
+ for (size_t j=0; j<node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
const unsigned int R = node_local_number_in_its_cells[j];
sum += Ajr(J,R);
}
- m_Ar(r) = sum;
+ m_Ar[r] = sum;
});
return m_Ar;
}
KOKKOS_INLINE_FUNCTION
- const Kokkos::View<const Rd*>
+ const NodeValue<const Rd>
computeBr(const NodeValuePerCell<Rdd>& Ajr,
- const NodeValuePerCell<Rd>& Cjr,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj) {
+ const NodeValuePerCell<const Rd>& Cjr,
+ const CellValue<const Rd>& uj,
+ const CellValue<const double>& pj) {
const auto& node_to_cell_matrix
- = m_connectivity.getMatrix(ItemType::node,
- ItemType::cell);
+ = m_connectivity.nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells
= m_connectivity.nodeLocalNumbersInTheirCells();
- Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
- Rd& br = m_br(r);
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r) {
+ Rd& br = m_br[r];
br = zero;
- const auto& node_to_cell = node_to_cell_matrix.rowConst(r);
+ const auto& node_to_cell = node_to_cell_matrix[r];
const auto& node_local_number_in_its_cells
= node_local_numbers_in_their_cells.itemValues(r);
- for (size_t j=0; j<node_to_cell.length; ++j) {
- const unsigned int J = node_to_cell(j);
+ for (size_t j=0; j<node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
const unsigned int R = node_local_number_in_its_cells[j];
- br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R);
+ br += Ajr(J,R)*uj[J] + pj[J]*Cjr(J,R);
}
});
@@ -179,11 +142,13 @@ class AcousticSolver
const Rdd nxn = tensorProduct(n,n);
const Rdd P = I-nxn;
+ const Array<const NodeId>& node_list
+ = symmetry_bc.nodeList();
Kokkos::parallel_for(symmetry_bc.numberOfNodes(), KOKKOS_LAMBDA(const int& r_number) {
- const int r = symmetry_bc.nodeList()[r_number];
+ const NodeId r = node_list[r_number];
- m_Ar(r) = P*m_Ar(r)*P + nxn;
- m_br(r) = P*m_br(r);
+ m_Ar[r] = P*m_Ar[r]*P + nxn;
+ m_br[r] = P*m_br[r];
});
break;
}
@@ -191,14 +156,14 @@ class AcousticSolver
}
}
- Kokkos::View<Rd*>
- computeUr(const Kokkos::View<const Rdd*>& Ar,
- const Kokkos::View<const Rd*>& br)
+ NodeValue<Rd>
+ computeUr(const NodeValue<const Rdd>& Ar,
+ const NodeValue<const Rd>& br)
{
inverse(Ar, m_inv_Ar);
- const Kokkos::View<const Rdd*> invAr = m_inv_Ar;
- Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
- m_ur[r]=invAr(r)*br(r);
+ const NodeValue<const Rdd> invAr = m_inv_Ar;
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r) {
+ m_ur[r]=invAr[r]*br[r];
});
return m_ur;
@@ -206,64 +171,64 @@ class AcousticSolver
void
computeFjr(const NodeValuePerCell<Rdd>& Ajr,
- const Kokkos::View<const Rd*>& ur,
- const NodeValuePerCell<Rd>& Cjr,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj)
+ const NodeValue<const Rd>& ur,
+ const NodeValuePerCell<const Rd>& Cjr,
+ const CellValue<const Rd>& uj,
+ const CellValue<const double>& pj)
{
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
- for (size_t r=0; r<cell_nodes.length; ++r) {
- m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(r)))+pj(j)*Cjr(j,r);
+ for (size_t r=0; r<cell_nodes.size(); ++r) {
+ m_Fjr(j,r) = Ajr(j,r)*(uj[j]-ur[cell_nodes[r]])+pj[j]*Cjr(j,r);
}
});
}
- void inverse(const Kokkos::View<const Rdd*>& A,
- Kokkos::View<Rdd*>& inv_A) const
+ void inverse(const NodeValue<const Rdd>& A,
+ NodeValue<Rdd>& inv_A) const
{
- Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
- inv_A(r) = ::inverse(A(r));
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r) {
+ inv_A[r] = ::inverse(A[r]);
});
}
KOKKOS_INLINE_FUNCTION
- void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
- const Kokkos::View<const Rd*>& xj,
- const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const double*>& cj,
- const Kokkos::View<const double*>& Vj,
- const NodeValuePerCell<Rd>& Cjr,
- const NodeValuePerCell<double>& ljr,
- const NodeValuePerCell<Rd>& njr) {
- const Kokkos::View<const double*> rhocj = computeRhoCj(rhoj, cj);
+ void computeExplicitFluxes(const NodeValue<const Rd>& xr,
+ const CellValue<const Rd>& xj,
+ const CellValue<const double>& rhoj,
+ const CellValue<const Rd>& uj,
+ const CellValue<const double>& pj,
+ const CellValue<const double>& cj,
+ const CellValue<const double>& Vj,
+ const NodeValuePerCell<const Rd>& Cjr,
+ const NodeValuePerCell<const double>& ljr,
+ const NodeValuePerCell<const Rd>& njr) {
+ const CellValue<const double> rhocj = computeRhoCj(rhoj, cj);
computeAjr(rhocj, Cjr, ljr, njr);
- const Kokkos::View<const Rdd*> Ar = computeAr(m_Ajr);
- const Kokkos::View<const Rd*> br = computeBr(m_Ajr, Cjr, uj, pj);
+ NodeValuePerCell<const Rdd> Ajr = m_Ajr;
+ const NodeValue<const Rdd> Ar = computeAr(Ajr);
+ const NodeValue<const Rd> br = computeBr(m_Ajr, Cjr, uj, pj);
this->applyBoundaryConditions();
- Kokkos::View<Rd*> ur = m_ur;
- ur = computeUr(Ar, br);
+ NodeValue<Rd>& ur = m_ur;
+ ur = computeUr(Ar, br);
computeFjr(m_Ajr, ur, Cjr, uj, pj);
}
- Kokkos::View<Rd*> m_br;
+ NodeValue<Rd> m_br;
NodeValuePerCell<Rdd> m_Ajr;
- Kokkos::View<Rdd*> m_Ar;
- Kokkos::View<Rdd*> m_inv_Ar;
+ NodeValue<Rdd> m_Ar;
+ NodeValue<Rdd> m_inv_Ar;
NodeValuePerCell<Rd> m_Fjr;
- Kokkos::View<Rd*> m_ur;
- Kokkos::View<double*> m_rhocj;
- Kokkos::View<double*> m_Vj_over_cj;
+ NodeValue<Rd> m_ur;
+ CellValue<double> m_rhocj;
+ CellValue<double> m_Vj_over_cj;
public:
AcousticSolver(MeshData& mesh_data,
@@ -273,79 +238,72 @@ class AcousticSolver
m_mesh(mesh_data.mesh()),
m_connectivity(m_mesh.connectivity()),
m_boundary_condition_list(bc_list),
- m_br("br", m_mesh.numberOfNodes()),
+ m_br(m_connectivity),
m_Ajr(m_connectivity),
- m_Ar("Ar", m_mesh.numberOfNodes()),
- m_inv_Ar("inv_Ar", m_mesh.numberOfNodes()),
+ m_Ar(m_connectivity),
+ m_inv_Ar(m_connectivity),
m_Fjr(m_connectivity),
- m_ur("ur", m_mesh.numberOfNodes()),
- m_rhocj("rho_c", m_mesh.numberOfCells()),
- m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells())
+ m_ur(m_connectivity),
+ m_rhocj(m_connectivity),
+ m_Vj_over_cj(m_connectivity)
{
;
}
KOKKOS_INLINE_FUNCTION
- double acoustic_dt(const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*>& cj) const
+ double acoustic_dt(const CellValue<const double>& Vj,
+ const CellValue<const double>& cj) const
{
- const NodeValuePerCell<double>& ljr = m_mesh_data.ljr();
+ const NodeValuePerCell<const double>& ljr = m_mesh_data.ljr();
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
+ const auto& cell_nodes = cell_to_node_matrix[j];
double S = 0;
- for (size_t r=0; r<cell_nodes.length; ++r) {
+ for (size_t r=0; r<cell_nodes.size(); ++r) {
S += ljr(j,r);
}
m_Vj_over_cj[j] = 2*Vj[j]/(S*cj[j]);
});
- double dt = std::numeric_limits<double>::max();
- Kokkos::parallel_reduce(m_mesh.numberOfCells(), ReduceMin(m_Vj_over_cj), dt);
-
- return dt;
+ return ReduceMin(m_Vj_over_cj);
}
-
void computeNextStep(const double& t, const double& dt,
UnknownsType& unknowns)
{
- Kokkos::View<double*> rhoj = unknowns.rhoj();
- Kokkos::View<Rd*> uj = unknowns.uj();
- Kokkos::View<double*> Ej = unknowns.Ej();
-
- Kokkos::View<double*> ej = unknowns.ej();
- Kokkos::View<double*> pj = unknowns.pj();
- Kokkos::View<double*> gammaj = unknowns.gammaj();
- Kokkos::View<double*> cj = unknowns.cj();
-
- const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
- const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
- const NodeValuePerCell<Rd>& Cjr = m_mesh_data.Cjr();
- const NodeValuePerCell<double>& ljr = m_mesh_data.ljr();
- const NodeValuePerCell<Rd>& njr = m_mesh_data.njr();
- Kokkos::View<Rd*> xr = m_mesh.xr();
+ CellValue<double>& rhoj = unknowns.rhoj();
+ CellValue<Rd>& uj = unknowns.uj();
+ CellValue<double>& Ej = unknowns.Ej();
+
+ CellValue<double>& ej = unknowns.ej();
+ CellValue<double>& pj = unknowns.pj();
+ CellValue<double>& cj = unknowns.cj();
+
+ const CellValue<const Rd>& xj = m_mesh_data.xj();
+ const CellValue<const double>& Vj = m_mesh_data.Vj();
+ const NodeValuePerCell<const Rd>& Cjr = m_mesh_data.Cjr();
+ const NodeValuePerCell<const double>& ljr = m_mesh_data.ljr();
+ const NodeValuePerCell<const Rd>& njr = m_mesh_data.njr();
+ const NodeValue<const Rd>& xr = m_mesh.xr();
computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr, ljr, njr);
const NodeValuePerCell<Rd>& Fjr = m_Fjr;
- const Kokkos::View<const Rd*> ur = m_ur;
+ const NodeValue<const Rd> ur = m_ur;
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
- const Kokkos::View<const double*> inv_mj = unknowns.invMj();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ const CellValue<const double> inv_mj = unknowns.invMj();
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
- for (size_t R=0; R<cell_nodes.length; ++R) {
- const unsigned int r=cell_nodes(R);
+ for (size_t R=0; R<cell_nodes.size(); ++R) {
+ const NodeId r=cell_nodes[R];
momentum_fluxes += Fjr(j,R);
energy_fluxes += (Fjr(j,R), ur[r]);
}
@@ -353,18 +311,18 @@ class AcousticSolver
Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
});
- Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
- xr[r] += dt*ur[r];
+ NodeValue<Rd> mutable_xr = m_mesh.mutableXr();
+ Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r){
+ mutable_xr[r] += dt*ur[r];
});
-
m_mesh_data.updateAllData();
- const Kokkos::View<const double*> mj = unknowns.mj();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ const CellValue<const double> mj = unknowns.mj();
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
rhoj[j] = mj[j]/Vj[j];
});
}
diff --git a/src/scheme/BlockPerfectGas.hpp b/src/scheme/BlockPerfectGas.hpp
index c018399e5d0644762a2474f3a99da6de20f82ffa..704787756bc97ac60b48699c8215c9444cb87bf2 100644
--- a/src/scheme/BlockPerfectGas.hpp
+++ b/src/scheme/BlockPerfectGas.hpp
@@ -1,58 +1,60 @@
#ifndef BLOCK_PERFECTGAS_HPP
#define BLOCK_PERFECTGAS_HPP
+#include <ItemValue.hpp>
+
struct BlockPerfectGas
{
private:
- Kokkos::View<double*> m_rhoj;
- Kokkos::View<double*> m_ej;
- Kokkos::View<double*> m_pj;
- Kokkos::View<double*> m_gammaj;
- Kokkos::View<double*> m_cj;
+ CellValue<double>& m_rhoj;
+ CellValue<double>& m_ej;
+ CellValue<double>& m_pj;
+ CellValue<double>& m_gammaj;
+ CellValue<double>& m_cj;
public:
- BlockPerfectGas(Kokkos::View<double*> rhoj,
- Kokkos::View<double*> ej,
- Kokkos::View<double*> pj,
- Kokkos::View<double*> gammaj,
- Kokkos::View<double*> cj)
- : m_rhoj(rhoj),
- m_ej(ej),
- m_pj(pj),
- m_gammaj(gammaj),
- m_cj(cj)
+ BlockPerfectGas(CellValue<double>& rhoj,
+ CellValue<double>& ej,
+ CellValue<double>& pj,
+ CellValue<double>& gammaj,
+ CellValue<double>& cj)
+ : m_rhoj(rhoj),
+ m_ej(ej),
+ m_pj(pj),
+ m_gammaj(gammaj),
+ m_cj(cj)
{
;
}
void updatePandCFromRhoE()
{
- const int nj = m_ej.size();
- const Kokkos::View<const double*> rho = m_rhoj;
- const Kokkos::View<const double*> e = m_ej;
- const Kokkos::View<const double*> gamma = m_gammaj;
-
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- const double gamma_minus_one = gamma[j]-1;
- m_pj[j] = gamma_minus_one*rho[j]*e[j];
- m_cj[j] = std::sqrt(gamma[j]*gamma_minus_one*e[j]);
+ const size_t nj = m_ej.size();
+ const CellValue<const double>& rho = m_rhoj;
+ const CellValue<const double>& e = m_ej;
+ const CellValue<const double>& gamma = m_gammaj;
+
+ Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const CellId& j){
+ const double gamma_minus_one = gamma[j]-1;
+ m_pj[j] = gamma_minus_one*rho[j]*e[j];
+ m_cj[j] = std::sqrt(gamma[j]*gamma_minus_one*e[j]);
});
}
void updateEandCFromRhoP()
{
- const int nj = m_ej.size();
- const Kokkos::View<const double*> rho = m_rhoj;
- const Kokkos::View<const double*> p = m_pj;
- const Kokkos::View<const double*> gamma = m_gammaj;
+ const size_t nj = m_ej.size();
+ const CellValue<const double>& rho = m_rhoj;
+ const CellValue<const double>& p = m_pj;
+ const CellValue<const double>& gamma = m_gammaj;
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- m_ej[j] = p[j]/(rho[j]*(gamma[j]-1));
+ Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const CellId& j){
+ m_ej[j] = p[j]/(rho[j]*(gamma[j]-1));
});
- const Kokkos::View<const double*> e = m_ej;
- Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const int& j){
- m_cj[j] = std::sqrt(gamma[j]*(gamma[j]-1)*e[j]);
+ const CellValue<const double>& e = m_ej;
+ Kokkos::parallel_for(nj, KOKKOS_LAMBDA(const CellId& j){
+ m_cj[j] = std::sqrt(gamma[j]*(gamma[j]-1)*e[j]);
});
}
};
diff --git a/src/scheme/BoundaryCondition.hpp b/src/scheme/BoundaryCondition.hpp
index c59eae3d347d6354a6772c3c2001bfbf412595d5..8747cb0c5cff48ca787e4186638cb03902966ca7 100644
--- a/src/scheme/BoundaryCondition.hpp
+++ b/src/scheme/BoundaryCondition.hpp
@@ -4,7 +4,7 @@
#include <vector>
#include <memory>
-#include <Kokkos_Core.hpp>
+#include <Array.hpp>
#include <RefNodeList.hpp>
#include <MeshNodeBoundary.hpp>
@@ -44,7 +44,7 @@ private:
const double m_value;
const size_t m_number_of_faces;
- Kokkos::View<unsigned int*> m_face_list;
+ Array<const unsigned int> m_face_list;
public:
double value() const
@@ -57,7 +57,7 @@ public:
return m_number_of_faces;
}
- const Kokkos::View<const unsigned int*> faceList() const
+ const Array<const unsigned int>& faceList() const
{
return m_face_list;
}
@@ -66,12 +66,13 @@ public:
const std::vector<unsigned int>& faces)
: BoundaryCondition(BoundaryCondition::pressure),
m_value(value),
- m_number_of_faces(faces.size()),
- m_face_list("faces_list", m_number_of_faces)
+ m_number_of_faces(faces.size())
{
+ Array<unsigned int> face_list(faces.size());
Kokkos::parallel_for(m_number_of_faces, KOKKOS_LAMBDA(const int& f){
- m_face_list[f]=faces[f];
+ face_list[f]=faces[f];
});
+ m_face_list = face_list;
}
~PressureBoundaryCondition() = default;
@@ -82,7 +83,7 @@ class SymmetryBoundaryCondition
: public BoundaryCondition
{
public:
- typedef TinyVector<dimension, double> Rd;
+ using Rd = TinyVector<dimension, double>;
private:
@@ -95,10 +96,10 @@ public:
size_t numberOfNodes() const
{
- return m_mesh_flat_node_boundary.nodeList().extent(0);
+ return m_mesh_flat_node_boundary.nodeList().size();
}
- const Kokkos::View<const unsigned int*> nodeList() const
+ const Array<const NodeId>& nodeList() const
{
return m_mesh_flat_node_boundary.nodeList();
}
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 617cb8b39759dcf4ae4ff6f4a41f09a7bc056584..2eaaf474dbbe520da30cc2c5b0ad5b4737e12538 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -1,179 +1,182 @@
#ifndef FINITE_VOLUMES_EULER_UNKNOWNS_HPP
#define FINITE_VOLUMES_EULER_UNKNOWNS_HPP
+#include <TinyVector.hpp>
+#include <ItemValue.hpp>
+
template <typename TMeshData>
class FiniteVolumesEulerUnknowns
{
public:
- typedef TMeshData MeshDataType;
- typedef typename MeshDataType::MeshType MeshType;
+ using MeshDataType = TMeshData;
+ using MeshType = typename MeshDataType::MeshType;
static constexpr size_t dimension = MeshType::dimension;
- typedef TinyVector<dimension> Rd;
+ using Rd = TinyVector<dimension>;
private:
const MeshDataType& m_mesh_data;
const MeshType& m_mesh;
- Kokkos::View<double*> m_rhoj;
- Kokkos::View<Rd*> m_uj;
- Kokkos::View<double*> m_Ej;
+ CellValue<double> m_rhoj;
+ CellValue<Rd> m_uj;
+ CellValue<double> m_Ej;
- Kokkos::View<double*> m_ej;
- Kokkos::View<double*> m_pj;
- Kokkos::View<double*> m_gammaj;
- Kokkos::View<double*> m_cj;
- Kokkos::View<double*> m_mj;
- Kokkos::View<double*> m_inv_mj;
+ CellValue<double> m_ej;
+ CellValue<double> m_pj;
+ CellValue<double> m_gammaj;
+ CellValue<double> m_cj;
+ CellValue<double> m_mj;
+ CellValue<double> m_inv_mj;
public:
- Kokkos::View<double*> rhoj()
+ CellValue<double>& rhoj()
{
return m_rhoj;
}
- const Kokkos::View<const double*> rhoj() const
+ const CellValue<const double> rhoj() const
{
return m_rhoj;
}
- Kokkos::View<Rd*> uj()
+ CellValue<Rd>& uj()
{
return m_uj;
}
- const Kokkos::View<const Rd*> uj() const
+ const CellValue<const Rd> uj() const
{
return m_uj;
}
- Kokkos::View<double*> Ej()
+ CellValue<double>& Ej()
{
return m_Ej;
}
- const Kokkos::View<const double*> Ej() const
+ const CellValue<const double> Ej() const
{
return m_Ej;
}
- Kokkos::View<double*> ej()
+ CellValue<double>& ej()
{
return m_ej;
}
- const Kokkos::View<const double*> ej() const
+ const CellValue<const double> ej() const
{
return m_ej;
}
- Kokkos::View<double*> pj()
+ CellValue<double>& pj()
{
return m_pj;
}
- const Kokkos::View<const double*> pj() const
+ const CellValue<const double> pj() const
{
return m_pj;
}
- Kokkos::View<double*> gammaj()
+ CellValue<double>& gammaj()
{
return m_gammaj;
}
- const Kokkos::View<const double*> gammaj() const
+ const CellValue<const double> gammaj() const
{
return m_gammaj;
}
- Kokkos::View<double*> cj()
+ CellValue<double>& cj()
{
return m_cj;
}
- const Kokkos::View<const double*> cj() const
+ const CellValue<const double> cj() const
{
return m_cj;
}
- Kokkos::View<double*> mj()
+ CellValue<double>& mj()
{
return m_mj;
}
- const Kokkos::View<const double*> mj() const
+ const CellValue<const double> mj() const
{
return m_mj;
}
- Kokkos::View<double*> invMj()
+ CellValue<double>& invMj()
{
return m_inv_mj;
}
- const Kokkos::View<const double*> invMj() const
+ const CellValue<const double> invMj() const
{
return m_inv_mj;
}
void initializeSod()
{
- const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+ const CellValue<const Rd>& xj = m_mesh_data.xj();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- if (xj[j][0]<0.5) {
- m_rhoj[j]=1;
- } else {
- m_rhoj[j]=0.125;
- }
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
+ if (xj[j][0]<0.5) {
+ m_rhoj[j]=1;
+ } else {
+ m_rhoj[j]=0.125;
+ }
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- if (xj[j][0]<0.5) {
- m_pj[j]=1;
- } else {
- m_pj[j]=0.1;
- }
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
+ if (xj[j][0]<0.5) {
+ m_pj[j]=1;
+ } else {
+ m_pj[j]=0.1;
+ }
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_uj[j] = zero;
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
+ m_uj[j] = zero;
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_gammaj[j] = 1.4;
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
+ m_gammaj[j] = 1.4;
});
BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
block_eos.updateEandCFromRhoP();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
+ m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
});
- const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_mj[j] = m_rhoj[j] * Vj[j];
+ const CellValue<const double>& Vj = m_mesh_data.Vj();
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
+ m_mj[j] = m_rhoj[j] * Vj[j];
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_inv_mj[j] = 1./m_mj[j];
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
+ m_inv_mj[j] = 1./m_mj[j];
});
}
FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)
: m_mesh_data(mesh_data),
m_mesh(m_mesh_data.mesh()),
- m_rhoj("rhoj",m_mesh.numberOfCells()),
- m_uj("uj",m_mesh.numberOfCells()),
- m_Ej("Ej",m_mesh.numberOfCells()),
- m_ej("ej",m_mesh.numberOfCells()),
- m_pj("pj",m_mesh.numberOfCells()),
- m_gammaj("gammaj",m_mesh.numberOfCells()),
- m_cj("cj",m_mesh.numberOfCells()),
- m_mj("mj",m_mesh.numberOfCells()),
- m_inv_mj("inv_mj",m_mesh.numberOfCells())
+ m_rhoj(m_mesh.connectivity()),
+ m_uj(m_mesh.connectivity()),
+ m_Ej(m_mesh.connectivity()),
+ m_ej(m_mesh.connectivity()),
+ m_pj(m_mesh.connectivity()),
+ m_gammaj(m_mesh.connectivity()),
+ m_cj(m_mesh.connectivity()),
+ m_mj(m_mesh.connectivity()),
+ m_inv_mj(m_mesh.connectivity())
{
;
}
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e72efb042c5c9006960a15dfd20c87ed29c29d82
--- /dev/null
+++ b/src/utils/Array.hpp
@@ -0,0 +1,83 @@
+#ifndef ARRAY_HPP
+#define ARRAY_HPP
+
+#include <PastisAssert.hpp>
+#include <Kokkos_Core.hpp>
+
+template <typename DataType>
+class Array
+{
+ public:
+ using data_type = DataType;
+ using index_type = size_t;
+
+ private:
+ Kokkos::View<DataType*> m_values;
+
+ // Allows const version to access our data
+ friend Array<std::add_const_t<DataType>>;
+
+ public:
+ KOKKOS_INLINE_FUNCTION
+ size_t size() const
+ {
+ return m_values.extent(0);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ DataType& operator[](const index_type& i) const
+ {
+ Assert(i<m_values.extent(0));
+ return m_values[i];
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ Array(const size_t& size)
+ : m_values("anonymous", size)
+ {
+ static_assert(not std::is_const<DataType>(),
+ "Cannot allocate Array of const data: only view is supported");
+ }
+
+ template <typename DataType2>
+ KOKKOS_INLINE_FUNCTION
+ Array& operator=(const Array<DataType2>& array)
+ {
+ // ensures that DataType is the same as source DataType2
+ static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+ "Cannot assign Array of different type");
+ // ensures that const is not lost through copy
+ static_assert(((std::is_const<DataType2>() and std::is_const<DataType>())
+ or not std::is_const<DataType2>()),
+ "Cannot assign Array of const to Array of non-const");
+ m_values = array.m_values;
+ return *this;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ Array& operator=(const Array&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ Array& operator=(Array&&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ Array() = default;
+
+ template <typename DataType2>
+ KOKKOS_INLINE_FUNCTION
+ Array(const Array<DataType2>& array)
+ {
+ this->operator=(array);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ Array(Array&&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ Array(const Array&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ~Array() = default;
+};
+
+#endif // ARRAY_HPP
diff --git a/src/utils/ArrayUtils.hpp b/src/utils/ArrayUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9a55abfebde94b4ab53adff473c52907ec5385f8
--- /dev/null
+++ b/src/utils/ArrayUtils.hpp
@@ -0,0 +1,109 @@
+#ifndef ARRAY_UTILS_HPP
+#define ARRAY_UTILS_HPP
+
+#include <Kokkos_Core.hpp>
+
+template <typename ArrayType>
+class ReduceMin
+{
+ private:
+ const ArrayType& m_array;
+ using data_type = std::remove_const_t<typename ArrayType::data_type>;
+ using index_type = typename ArrayType::index_type;
+
+ public:
+ KOKKOS_INLINE_FUNCTION
+ operator data_type()
+ {
+ data_type reduced_value;
+ Kokkos::parallel_reduce(m_array.size(), *this, reduced_value);
+ return reduced_value;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ void operator()(const index_type& i, data_type& data) const
+ {
+ if (m_array[i] < data) {
+ data = m_array[i];
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ void join(volatile data_type& dst,
+ const volatile data_type& src) const
+ {
+ if (src < dst) {
+ dst = src;
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION void
+ init(data_type& value) const
+ {
+ value = std::numeric_limits<data_type>::max();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ ReduceMin(const ArrayType& array)
+ : m_array(array)
+ {
+ ;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ ~ReduceMin() = default;
+};
+
+
+template <typename ArrayType>
+class ReduceMax
+{
+ private:
+ const ArrayType& m_array;
+ using data_type = std::remove_const_t<typename ArrayType::data_type>;
+ using index_type = typename ArrayType::index_type;
+
+ public:
+ KOKKOS_INLINE_FUNCTION
+ operator data_type()
+ {
+ data_type reduced_value;
+ Kokkos::parallel_reduce(m_array.size(), *this, reduced_value);
+ return reduced_value;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ void operator()(const index_type& i, data_type& data) const
+ {
+ if (m_array[i] > data) {
+ data = m_array[i];
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ void join(volatile data_type& dst,
+ const volatile data_type& src) const
+ {
+ if (src > dst) {
+ dst = src;
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION void
+ init(data_type& value) const
+ {
+ value = -std::numeric_limits<data_type>::max();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ ReduceMax(const ArrayType& array)
+ : m_array(array)
+ {
+ ;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ ~ReduceMax() = default;
+};
+
+#endif //ARRAY_UTILS_HPP
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 3b4823f1eaa5d56c7bb2a7e7201d623e8308d60a..c36e9a2d570b0a385e679904a61dd3b260c25cd8 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -4,6 +4,8 @@ set(EXECUTABLE_OUTPUT_PATH ${PASTIS_BINARY_DIR})
add_executable (unit_tests
test_main.cpp
+ test_Array.cpp
+ test_ArrayUtils.cpp
test_ItemType.cpp
test_PastisAssert.cpp
test_RevisionInfo.cpp
@@ -15,6 +17,7 @@ target_include_directories(Catch INTERFACE ${CATCH_INCLUDE_DIR})
target_link_libraries (unit_tests
PastisUtils
+ kokkos
Catch
)
diff --git a/tests/test_Array.cpp b/tests/test_Array.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6cd02146c8775dca7a061672b8efed1fbdcad236
--- /dev/null
+++ b/tests/test_Array.cpp
@@ -0,0 +1,87 @@
+#include <catch.hpp>
+
+#include <PastisAssert.hpp>
+#include <Array.hpp>
+#include <Types.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class Array<int>;
+
+TEST_CASE("Array", "[utils]") {
+
+ Array<int> a(10);
+ REQUIRE(a.size() == 10);
+
+ for (size_t i=0; i<a.size(); ++i) {
+ a[i] = 2*i;
+ }
+
+ REQUIRE(((a[0] == 0) and (a[1] == 2) and
+ (a[2] == 4) and (a[3] == 6) and
+ (a[4] == 8) and (a[5] ==10) and
+ (a[6] ==12) and (a[7] ==14) and
+ (a[8] ==16) and (a[9] ==18)));
+
+ SECTION("checking for copies") {
+ Array<const int> b{a};
+
+ REQUIRE(((b[0] == 0) and (b[1] == 2) and
+ (b[2] == 4) and (b[3] == 6) and
+ (b[4] == 8) and (b[5] ==10) and
+ (b[6] ==12) and (b[7] ==14) and
+ (b[8] ==16) and (b[9] ==18)));
+
+ Array<int> c{a};
+
+ REQUIRE(((c[0] == 0) and (c[1] == 2) and
+ (c[2] == 4) and (c[3] == 6) and
+ (c[4] == 8) and (c[5] ==10) and
+ (c[6] ==12) and (c[7] ==14) and
+ (c[8] ==16) and (c[9] ==18)));
+
+
+ Array<int> d = std::move(c);
+
+ REQUIRE(((d[0] == 0) and (d[1] == 2) and
+ (d[2] == 4) and (d[3] == 6) and
+ (d[4] == 8) and (d[5] ==10) and
+ (d[6] ==12) and (d[7] ==14) and
+ (d[8] ==16) and (d[9] ==18)));
+
+ }
+
+ SECTION("checking for affectations") {
+ Array<const int> b;
+ b = a;
+
+ REQUIRE(((b[0] == 0) and (b[1] == 2) and
+ (b[2] == 4) and (b[3] == 6) and
+ (b[4] == 8) and (b[5] ==10) and
+ (b[6] ==12) and (b[7] ==14) and
+ (b[8] ==16) and (b[9] ==18)));
+
+
+ Array<int> c;
+ c = a;
+
+ REQUIRE(((c[0] == 0) and (c[1] == 2) and
+ (c[2] == 4) and (c[3] == 6) and
+ (c[4] == 8) and (c[5] ==10) and
+ (c[6] ==12) and (c[7] ==14) and
+ (c[8] ==16) and (c[9] ==18)));
+
+ Array<int> d;
+ d = std::move(c);
+
+ REQUIRE(((d[0] == 0) and (d[1] == 2) and
+ (d[2] == 4) and (d[3] == 6) and
+ (d[4] == 8) and (d[5] ==10) and
+ (d[6] ==12) and (d[7] ==14) and
+ (d[8] ==16) and (d[9] ==18)));
+ }
+#ifndef NDEBUG
+ SECTION("checking for bounds violation") {
+ REQUIRE_THROWS_AS(a[10], AssertError);
+ }
+#endif // NDEBUG
+}
diff --git a/tests/test_ArrayUtils.cpp b/tests/test_ArrayUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8183e7157a2c39d3f27b3fbd703917fcdf567ff1
--- /dev/null
+++ b/tests/test_ArrayUtils.cpp
@@ -0,0 +1,96 @@
+#include <catch.hpp>
+
+#include <PastisAssert.hpp>
+#include <Array.hpp>
+#include <ArrayUtils.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class Array<int>;
+
+TEST_CASE("ArrayUtils", "[utils]") {
+
+ SECTION("checking for reductions") {
+ using ArrayType = Array<int>;
+ using data_type = ArrayType::data_type;
+
+ Array<int> a(10);
+ a[0] =13;
+ a[1] = 1;
+ a[2] = 8;
+ a[3] =-3;
+ a[4] =23;
+ a[5] =-1;
+ a[6] =13;
+ a[7] = 0;
+ a[8] =12;
+ a[9] = 9;
+
+ SECTION("ReduceMin") {
+ ReduceMin r_min(a);
+ data_type data;
+
+ r_min.init(data);
+ REQUIRE(data == std::numeric_limits<data_type>::max());
+
+ r_min(2,data);
+ REQUIRE(data == 8);
+
+ r_min(9,data);
+ REQUIRE(data == 8);
+
+ r_min(5,data);
+ REQUIRE(data == -1);
+
+ {
+ data = 0;
+ data_type r_value{3};
+ r_min.join(data, r_value);
+ REQUIRE(data == 0);
+ }
+
+ {
+ data = 0;
+ data_type r_value{-5};
+ r_min.join(data, r_value);
+ REQUIRE(data == -5);
+ }
+
+ REQUIRE((r_min == -3));
+
+ REQUIRE((ReduceMin(a) == -3));
+ }
+
+ SECTION("ReduceMax") {
+ ReduceMax r_max(a);
+ data_type data;
+
+ r_max.init(data);
+ REQUIRE(data == -std::numeric_limits<data_type>::max());
+
+ r_max(3,data);
+ REQUIRE(data == -3);
+ r_max(5,data);
+ REQUIRE(data == -1);
+ r_max(8,data);
+ REQUIRE(data == 12);
+
+ {
+ data = 0;
+ data_type r_value{3};
+ r_max.join(data, r_value);
+ REQUIRE(data == 3);
+ }
+
+ {
+ data = 0;
+ data_type r_value{-5};
+ r_max.join(data, r_value);
+ REQUIRE(data == 0);
+ }
+
+ REQUIRE((r_max == 23));
+
+ REQUIRE((ReduceMax(a) == 23));
+ }
+ }
+}
diff --git a/tests/test_ItemType.cpp b/tests/test_ItemType.cpp
index 7c39963613dd3bae8d60e5fb4c9a4b0915ce6343..096dbdaf5b8a85099c8ec892a15c7f6ceb1000ca 100644
--- a/tests/test_ItemType.cpp
+++ b/tests/test_ItemType.cpp
@@ -26,23 +26,23 @@ TEST_CASE("ItemType", "[connectivity]") {
}
SECTION("checking for item ids in 1d") {
- REQUIRE(ItemId<1>::itemId(cell_type)==0);
- REQUIRE(ItemId<1>::itemId(edge_type)==1);
- REQUIRE(ItemId<1>::itemId(face_type)==1);
- REQUIRE(ItemId<1>::itemId(node_type)==1);
+ REQUIRE(ItemTypeId<1>::itemTId(cell_type)==0);
+ REQUIRE(ItemTypeId<1>::itemTId(edge_type)==1);
+ REQUIRE(ItemTypeId<1>::itemTId(face_type)==1);
+ REQUIRE(ItemTypeId<1>::itemTId(node_type)==1);
}
SECTION("checking for item ids in 2d") {
- REQUIRE(ItemId<2>::itemId(cell_type)==0);
- REQUIRE(ItemId<2>::itemId(edge_type)==1);
- REQUIRE(ItemId<2>::itemId(face_type)==1);
- REQUIRE(ItemId<2>::itemId(node_type)==2);
+ REQUIRE(ItemTypeId<2>::itemTId(cell_type)==0);
+ REQUIRE(ItemTypeId<2>::itemTId(edge_type)==1);
+ REQUIRE(ItemTypeId<2>::itemTId(face_type)==1);
+ REQUIRE(ItemTypeId<2>::itemTId(node_type)==2);
}
SECTION("checking for item ids in 3d") {
- REQUIRE(ItemId<3>::itemId(cell_type)==0);
- REQUIRE(ItemId<3>::itemId(edge_type)==1);
- REQUIRE(ItemId<3>::itemId(face_type)==2);
- REQUIRE(ItemId<3>::itemId(node_type)==3);
+ REQUIRE(ItemTypeId<3>::itemTId(cell_type)==0);
+ REQUIRE(ItemTypeId<3>::itemTId(edge_type)==1);
+ REQUIRE(ItemTypeId<3>::itemTId(face_type)==2);
+ REQUIRE(ItemTypeId<3>::itemTId(node_type)==3);
}
}
diff --git a/tests/test_main.cpp b/tests/test_main.cpp
index b3143fbb1788b7daef8a6cc20d8a0e9bafb308a5..98ed156f241538f4f0f47e967576a190ad949cf5 100644
--- a/tests/test_main.cpp
+++ b/tests/test_main.cpp
@@ -1,2 +1,15 @@
-#define CATCH_CONFIG_MAIN
+#define CATCH_CONFIG_RUNNER
#include <catch.hpp>
+
+#include <Kokkos_Core.hpp>
+
+int main( int argc, char* argv[] )
+{
+ Kokkos::initialize({4,-1,-1,true});
+
+ int result = Catch::Session().run( argc, argv );
+
+ Kokkos::finalize();
+
+ return result;
+}