Fixed acoustic scheme for d>1 (introducing ljr=|Cjr| and njr = 1/ljr Cjr)

src/mesh/MeshData.hpp

@@ -20,6 +20,8 @@ public:
 private:
   const MeshType& m_mesh;
   Kokkos::View<Rd**> m_Cjr;
+  Kokkos::View<double**>  m_ljr;
+  Kokkos::View<Rd**> m_njr;
   Kokkos::View<Rd*>  m_xj;
   Kokkos::View<double*>   m_Vj;
 
@@ -74,7 +76,7 @@ private:
   KOKKOS_INLINE_FUNCTION
   void _updateCjr() {
     if(dimension == 1) {
-      // Cjr are constant overtime
+      // Cjr/njr/ljr are constant overtime
     }
     static_assert(dimension==1, "only 1d is implemented");
   }
@@ -90,6 +92,16 @@ public:
     return m_Cjr;
   }
 
+  const Kokkos::View<const double**> ljr() const
+  {
+    return m_ljr;
+  }
+
+  const Kokkos::View<const Rd**> njr() const
+  {
+    return m_njr;
+  }
+
   const Kokkos::View<const Rd*> xj() const
   {
     return m_xj;
@@ -110,6 +122,8 @@ public:
   MeshData(const MeshType& mesh)
     : m_mesh(mesh),
       m_Cjr("Cjr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()),
+      m_ljr("ljr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()),
+      m_njr("njr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()),
       m_xj("xj", mesh.numberOfCells()),
       m_Vj("Vj", mesh.numberOfCells())
   {
@@ -119,6 +133,12 @@ public:
 	  m_Cjr(j,0)=-1;
 	  m_Cjr(j,1)= 1;
 	});
+      // in 1d njr=Cjr
+      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;
+	});
     }
     this->updateAllData();
   }

src/scheme/AcousticSolver.hpp

@@ -79,13 +79,14 @@ private:
   KOKKOS_INLINE_FUNCTION
   const Kokkos::View<const Rdd**>
   computeAjr(const Kokkos::View<const double*>& rhocj,
-	     const Kokkos::View<const Rd**>& Cjr) {
+	     const Kokkos::View<const double**>& ljr,
+	     const Kokkos::View<const Rd**>& njr) {
     const Kokkos::View<const unsigned short*> cell_nb_nodes
       = m_connectivity.cellNbNodes();
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
 	for (int r=0; r<cell_nb_nodes[j]; ++r) {
-	  m_Ajr(j,r) = tensorProduct(rhocj(j)*Cjr(j,r), Cjr(j,r));
+	  m_Ajr(j,r) = tensorProduct((rhocj(j)*ljr(j,r))*njr(j,r), njr(j,r));
 	}
       });
 
@@ -190,9 +191,11 @@ private:
 			     const Kokkos::View<const double*>& pj,
 			     const Kokkos::View<const double*>& cj,
 			     const Kokkos::View<const double*>& Vj,
-			     const Kokkos::View<const Rd**>& Cjr) {
+			     const Kokkos::View<const Rd**>& Cjr,
+			     const Kokkos::View<const double**>& ljr,
+			     const Kokkos::View<const Rd**>& njr) {
     const Kokkos::View<const double*> rhocj  = computeRhoCj(rhoj, cj);
-    const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, Cjr);
+    const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, ljr, njr);
 
     const Kokkos::View<const Rdd*> Ar = computeAr(Ajr);
     const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj);
@@ -259,9 +262,11 @@ public:
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
     const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
     const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
+    const Kokkos::View<const double**> ljr = m_mesh_data.ljr();
+    const Kokkos::View<const Rd**> njr = m_mesh_data.njr();
     Kokkos::View<Rd*> xr = m_mesh.xr();
 
-    computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr);
+    computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr, ljr, njr);
 
     const Kokkos::View<const Rd**> Fjr = m_Fjr;
     const Kokkos::View<const Rd*> ur = m_ur;