Feature/itemvalue

src/scheme/AcousticSolver.hpp

 #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
+  const CellValue<const double>
+  computeRhoCj(const CellValue<const double>& rhoj,
+               const CellValue<const double>& cj)
   {
-      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)
-  {
-    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;
+    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];
       });
   }