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 8e640d841a1e3d82b4fabe8c255f7f038801005d..06aec114228cbc8858c0b6c7769b916702ff5936 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -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,8 +233,8 @@ 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) {
fout << xj[j][0] << ' ' << rhoj[j] << '\n';
@@ -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);
@@ -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..c3a2c01a9751db986e64647467e92240a7f49611 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -9,7 +9,7 @@ void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
const auto& cell_to_node_matrix
= 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) {
const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
@@ -227,14 +227,14 @@ Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
Assert(this->numberOfCells()>0);
{
- Kokkos::View<CellType*> cell_type("cell_type", this->numberOfCells());
+ CellValue<CellType> cell_type(*this);
Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& 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());
+ CellValue<double> inv_cell_nb_nodes(*this);
Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& 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 b72a38abdadae65a8aed62d4ec9ebd0f09cbf628..964daeb262b287c60cf9c90ad1ff400f813e5af2 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -5,6 +5,7 @@
#include <TinyVector.hpp>
#include <Kokkos_Core.hpp>
+#include <ItemValue.hpp>
#include <IConnectivity.hpp>
@@ -234,7 +235,7 @@ class Connectivity final
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,7 +260,7 @@ 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>;
@@ -470,7 +471,7 @@ class Connectivity final
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/GmshReader.cpp b/src/mesh/GmshReader.cpp
index f16ab2fa02e06aad3715eeb7c84028476c7565ae..d483ee92421713aa1d9a38b7e99118c7235a6b1a 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -846,7 +846,7 @@ GmshReader::__proceedData()
typedef Mesh<Connectivity3D> MeshType;
typedef TinyVector<3, double> Rd;
- Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
+ NodeValue<Rd> xr(connectivity);
for (size_t i=0; i<__vertices.extent(0); ++i) {
xr[i][0] = __vertices[i][0];
xr[i][1] = __vertices[i][1];
@@ -918,7 +918,7 @@ GmshReader::__proceedData()
typedef Mesh<Connectivity2D> MeshType;
typedef TinyVector<2, double> Rd;
- Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
+ NodeValue<Rd> xr(connectivity);
for (size_t i=0; i<__vertices.extent(0); ++i) {
xr[i][0] = __vertices[i][0];
xr[i][1] = __vertices[i][1];
@@ -962,7 +962,7 @@ GmshReader::__proceedData()
typedef Mesh<Connectivity1D> MeshType;
typedef TinyVector<1, double> Rd;
- Kokkos::View<Rd*> xr("xr", __vertices.extent(0));
+ NodeValue<Rd> xr(connectivity);
for (size_t i=0; i<__vertices.extent(0); ++i) {
xr[i][0] = __vertices[i][0];
}
diff --git a/src/mesh/ItemValue.hpp b/src/mesh/ItemValue.hpp
index a00cb3bfd773d0dc0990be6436d39f97b4c70918..5534df4893230d4abfa9d0502a334a28b77eb71f 100644
--- a/src/mesh/ItemValue.hpp
+++ b/src/mesh/ItemValue.hpp
@@ -1,6 +1,8 @@
#ifndef ITEM_VALUE_HPP
#define ITEM_VALUE_HPP
+#include <Kokkos_Core.hpp>
+
#include <ItemType.hpp>
#include <PastisAssert.hpp>
@@ -12,6 +14,7 @@ class ItemValue
{
public:
static const ItemType item_t{item_type};
+ typedef DataType data_type;
private:
bool m_is_built{false};
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 3aee2c287d5cdd06c4a754115caf930973c4f2cb..383e92294d6064484e6eaf0b777caa9ce06fbb0e 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>
@@ -22,15 +22,17 @@ public:
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 +56,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 254ff42b1ac47f57c377762ed95ac0dfc9af7865..aa24493e89a9ede90f1d46cb2c3e935db3d4fa35 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -24,50 +24,54 @@ class MeshData
private:
const MeshType& m_mesh;
- NodeValuePerCell<Rd> m_Cjr;
- NodeValuePerCell<double> m_ljr;
- NodeValuePerCell<Rd> m_njr;
- Kokkos::View<Rd*> m_xj;
+ 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);
+ CellValue<Rd> xj(m_mesh.connectivity());
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)]);
+ 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);
-
+ CellValue<Rd> xj(m_mesh.connectivity());
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& 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)];
}
- 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);
@@ -90,33 +94,42 @@ 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];
- });
-
+ {
+ NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
+ 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)));
+ 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 int& j){
+ ljr[j] = l2Norm(m_Cjr[j]);
+ });
+ m_ljr = ljr;
+ }
+
+ {
+ NodeValuePerCell<Rd> njr(m_mesh.connectivity());
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){
+ njr[j] = (1./m_ljr[j])*m_Cjr[j];
+ });
+ 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
@@ -143,10 +156,6 @@ class MeshData
}
});
- 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);
@@ -156,51 +165,67 @@ class MeshData
ItemType::face);
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 int& jr){
+ Cjr[jr] = zero;
+ });
+
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(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);
+ const auto& cell_faces = cell_to_face_matrix.rowConst(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);
#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 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;
+ }
}
- }
- 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.length; ++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.length; ++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 int& 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 int& jr){
+ njr[jr] = (1./m_ljr[jr])*m_Cjr[jr];
+ });
+ m_njr = njr;
+ }
}
static_assert((dimension<=3), "only 1d, 2d and 3d are implemented");
}
@@ -211,22 +236,22 @@ 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;
}
@@ -244,25 +269,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(mesh.connectivity())
+ : 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 int& 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 int& jr){
+ ljr[jr] = 1;
});
+ m_ljr = ljr;
+ }
}
this->updateAllData();
}
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 87c537439bf1a27e8443865677cd6b58c132378c..91a0943dea549e8cf101d6045659201236ecd3b2 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -2,6 +2,8 @@
#define MESH_NODE_BOUNDARY_HPP
#include <Kokkos_Core.hpp>
+#include <ItemValue.hpp>
+
#include <Kokkos_Vector.hpp>
#include <TinyVector.hpp>
@@ -149,7 +151,7 @@ _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
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) {
const R2& x = xr[m_node_list[r]];
@@ -186,7 +188,7 @@ _getNormal(const MeshType& mesh)
static_assert(MeshType::dimension == 2);
typedef TinyVector<2,double> R2;
- 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());
@@ -239,7 +241,7 @@ _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) {
const R3& x = xr[m_node_list[r]];
@@ -314,7 +316,7 @@ _getOutgoingNormal(const MeshType& mesh)
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);
@@ -351,7 +353,7 @@ _getOutgoingNormal(const MeshType& mesh)
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);
@@ -388,7 +390,7 @@ _getOutgoingNormal(const MeshType& mesh)
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);
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index 52d91666b17076770d84b15f736fdf503134a2b8..be60e2634e7ae4ff380e36b52d15fc9e77492923 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -182,7 +182,7 @@ class SubItemValuePerItem<DataType,
= 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 = Kokkos::View<std::remove_const_t<DataType>*>("values", connectivity_matrix.numEntries());
}
~SubItemValuePerItem() = default;
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index a33fef0bfbe2abf2c73b34e7c3b564fb2cd81461..24865f43decc9ccce1514d619fa01a871df13176 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,16 +41,16 @@ 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();
+ using Rd = TinyVector<MeshType::dimension>;
+ const NodeValue<const Rd>& xr = mesh.xr();
+ if constexpr(MeshType::dimension == 1) {
for (unsigned int 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();
+ } else if constexpr (MeshType::dimension == 2) {
for (unsigned int r=0; r<mesh.numberOfNodes(); ++r) {
for (unsigned short i=0; i<2; ++i) {
fout << xr[r][i] << ' ';
@@ -56,7 +58,6 @@ class VTKWriter
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 (unsigned short i=0; i<3; ++i) {
fout << xr[r][i] << ' ';
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 6f83bdcaa229ca5b816711fba6ca4c3cb849deca..3bc6269eea2c7d74b4fee36d13e6da83db19fafc 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -36,14 +36,15 @@ class AcousticSolver
struct ReduceMin
{
private:
- const Kokkos::View<const double*> x_;
+ const CellValue<const double> x_;
public:
- typedef Kokkos::View<const double*>::non_const_value_type value_type;
+ typedef std::remove_const_t<CellValue<const double>::data_type> value_type;
- ReduceMin(const Kokkos::View<const double*>& x) : x_ (x) {}
+ ReduceMin(const CellValue<const double>& x) : x_ (x) {}
- typedef Kokkos::View<const double*>::size_type size_type;
+#warning was: "typedef Kokkos::View<const double*>::size_type size_type;""
+ typedef size_t size_type;
KOKKOS_INLINE_FUNCTION
void operator() (const size_type i, value_type& update) const
@@ -70,9 +71,9 @@ class AcousticSolver
};
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) {
m_rhocj[j] = rhoj[j]*cj[j];
@@ -81,10 +82,10 @@ class AcousticSolver
}
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) {
const size_t& nb_nodes =m_Ajr.numberOfSubValues(j);
@@ -96,8 +97,8 @@ 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);
@@ -122,11 +123,11 @@ class AcousticSolver
}
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);
@@ -191,12 +192,12 @@ 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;
+ const NodeValue<const Rdd> invAr = m_inv_Ar;
Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
m_ur[r]=invAr(r)*br(r);
});
@@ -206,10 +207,10 @@ 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,
@@ -224,8 +225,8 @@ class AcousticSolver
});
}
- 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));
@@ -233,37 +234,38 @@ class AcousticSolver
}
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,
+ 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<Rd>& Cjr,
- const NodeValuePerCell<double>& ljr,
- const NodeValuePerCell<Rd>& njr) {
- const Kokkos::View<const double*> rhocj = computeRhoCj(rhoj, cj);
+ 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,23 +275,23 @@ 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 CellValue<const double>& Vj,
- const Kokkos::View<const double*>& cj) const
+ 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);
@@ -314,31 +316,31 @@ class AcousticSolver
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();
+ CellValue<double>& rhoj = unknowns.rhoj();
+ CellValue<Rd>& uj = unknowns.uj();
+ CellValue<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();
+ CellValue<double>& ej = unknowns.ej();
+ CellValue<double>& pj = unknowns.pj();
+ CellValue<double>& gammaj = unknowns.gammaj();
+ CellValue<double>& cj = unknowns.cj();
- const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+ const CellValue<const Rd>& xj = m_mesh_data.xj();
const CellValue<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();
+ 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);
- const Kokkos::View<const double*> inv_mj = unknowns.invMj();
+ const CellValue<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);
@@ -357,13 +359,13 @@ class AcousticSolver
ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
});
+ NodeValue<Rd> mutable_xr = m_mesh.mutableXr();
Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
- xr[r] += dt*ur[r];
+ mutable_xr[r] += dt*ur[r];
});
-
m_mesh_data.updateAllData();
- const Kokkos::View<const double*> mj = unknowns.mj();
+ const CellValue<const double> mj = unknowns.mj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
rhoj[j] = mj[j]/Vj[j];
});
diff --git a/src/scheme/BlockPerfectGas.hpp b/src/scheme/BlockPerfectGas.hpp
index c018399e5d0644762a2474f3a99da6de20f82ffa..97a0256efb672b9c680c9edbecd99d96f20ee3fc 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.numberOfValues();
+ 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 size_t& 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.numberOfValues();
+ 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 size_t& 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 size_t& j){
+ m_cj[j] = std::sqrt(gamma[j]*(gamma[j]-1)*e[j]);
});
}
};
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index fd58774a6fd2902b96d6664bde4bbc98b0ab8fc2..ba0684ad2ba7bc35feccc845bb16481f6327b653 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -1,6 +1,9 @@
#ifndef FINITE_VOLUMES_EULER_UNKNOWNS_HPP
#define FINITE_VOLUMES_EULER_UNKNOWNS_HPP
+#include <TinyVector.hpp>
+#include <ItemValue.hpp>
+
template <typename TMeshData>
class FiniteVolumesEulerUnknowns
{
@@ -15,165 +18,165 @@ 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;
- }
+ 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;
- }
+ 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;
+ m_uj[j] = zero;
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_gammaj[j] = 1.4;
+ 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]);
+ m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
});
const CellValue<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];
+ 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];
+ 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())
{
;
}