From f2d7f9f7f5431ab1c511040c2b964e23e82af5d9 Mon Sep 17 00:00:00 2001
From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr>
Date: Fri, 31 Aug 2018 16:35:22 +0200
Subject: [PATCH] autre tentative, SIGFPE

---
 src/main.cpp                              |  4 +--
 src/scheme/FiniteVolumesEulerUnknowns.hpp | 24 +++++++++++--
 src/scheme/NoSplitting.hpp                | 43 ++++++++++++++---------
 3 files changed, 49 insertions(+), 22 deletions(-)

diff --git a/src/main.cpp b/src/main.cpp
index 90feecfff..00086299a 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -144,7 +144,7 @@ int main(int argc, char *argv[])
     const Kokkos::View<const double*> Vj = mesh_data.Vj();
     const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
 
-    const double tmax=0.8;
+    const double tmax=0.2;
     double t=0.;
 
     int itermax=std::numeric_limits<int>::max();
@@ -320,7 +320,7 @@ int main(int argc, char *argv[])
 
       // NAVIER-STOKES SANS SPLITTING
 
-      double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
+      double dt = 0.05*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
       if (t+dt > tmax) {
 	dt = tmax-t;
       }
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 2b78943cd..81b26c05f 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -12,6 +12,7 @@ public:
   typedef TinyVector<dimension> Rd;
 
 private:
+
   const MeshDataType& m_mesh_data;
   const MeshType& m_mesh;
 
@@ -26,6 +27,7 @@ private:
   Kokkos::View<double*> m_mj;
   Kokkos::View<double*> m_inv_mj;
   Kokkos::View<double*> m_kj;
+  Kokkos::View<double*> m_PTj;
 
   Kokkos::View<Rd*> m_uL;
   Kokkos::View<Rd*> m_uR;
@@ -143,7 +145,17 @@ public:
   {
     return m_kj;
   }
+  
+  Kokkos::View<double*> PTj()
+  {
+    return m_PTj;
+  }
 
+  const Kokkos::View<const double*> PTj() const
+  {
+    return m_PTj;
+  }
+  
   Kokkos::View<Rd*> uL()
   {
     return m_uL;
@@ -283,6 +295,8 @@ public:
   {
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
 
+    const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
+    
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
 	// TAC
 	/*
@@ -328,7 +342,6 @@ public:
   	m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
       });
 
-    const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
   	m_mj[j] = m_rhoj[j] * Vj[j];
       });
@@ -386,7 +399,11 @@ public:
 	*/
 	
       });
-
+    
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+	m_PTj(j) = 0;
+      });
+    
      // Conditions aux bords de Dirichlet sur u et k
 
     m_uL[0] = zero;
@@ -524,7 +541,8 @@ public:
       m_nuR("nuR", 1),
       m_entropy("entropy", m_mesh.numberOfCells()),
       m_entropy_new("entropy_new", m_mesh.numberOfCells()),
-    m_S0("S0", m_mesh.numberOfCells())
+    m_S0("S0", m_mesh.numberOfCells()),
+    m_PTj("PTj",m_mesh.numberOfCells())
   {
     ;
   }
diff --git a/src/scheme/NoSplitting.hpp b/src/scheme/NoSplitting.hpp
index f046a3d86..039b535f2 100644
--- a/src/scheme/NoSplitting.hpp
+++ b/src/scheme/NoSplitting.hpp
@@ -71,7 +71,7 @@ private:
     }
   };
 
-
+  /*
   // ----
 
   KOKKOS_INLINE_FUNCTION
@@ -100,7 +100,7 @@ private:
   }
 
   // ----
-
+  */
 
   KOKKOS_INLINE_FUNCTION
   const Kokkos::View<const double*>
@@ -154,7 +154,7 @@ private:
   computeBr(const Kokkos::View<const Rdd**>& Ajr,
 	    const Kokkos::View<const Rd**>& Cjr,
 	    const Kokkos::View<const Rd*>& uj,
-	    const Kokkos::View<const double*>& pj,
+	    const Kokkos::View<const double*>& PTj,
 	    const double t) {
     const Kokkos::View<const unsigned int**>& node_cells = m_connectivity.nodeCells();
     const Kokkos::View<const unsigned short**>& node_cell_local_node = m_connectivity.nodeCellLocalNode();
@@ -167,7 +167,7 @@ private:
 	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);
+	  br += Ajr(J,R)*uj(J) + PTj(J)*Cjr(J,R);
 	}
       });
 
@@ -212,14 +212,14 @@ private:
 	     const Kokkos::View<const Rd*>& ur,
 	     const Kokkos::View<const Rd**>& Cjr,
 	     const Kokkos::View<const Rd*>& uj,
-	     const Kokkos::View<const double*>& pj) {
+	     const Kokkos::View<const double*>& PTj) {
     const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
     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_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+pj(j)*Cjr(j,r);
+	  m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+PTj(j)*Cjr(j,r);
 	}
       });
 
@@ -251,7 +251,7 @@ private:
 			     const Kokkos::View<const double*>& kj,
 			     const Kokkos::View<const double*>& rhoj,
 			     const Kokkos::View<const Rd*>& uj,
-			     const Kokkos::View<const double*>& pj,
+			     const Kokkos::View<const double*>& PTj,
 			     const Kokkos::View<const double*>& cj,
 			     const Kokkos::View<const double*>& Vj,
 			     const Kokkos::View<const Rd**>& Cjr,
@@ -259,15 +259,13 @@ private:
     const Kokkos::View<const double*> rhocj  = computeRhoCj(rhoj, cj);
     const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, Cjr);
 
-    const Kokkos::View<const double*> PTj = computePj(pj, kj, uj);
-
     const Kokkos::View<const Rdd*> Ar = computeAr(Ajr);
     const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, PTj, t);
 
     Kokkos::View<Rd*> ur = m_ur;
     Kokkos::View<Rd**> Fjr = m_Fjr;
     ur = computeUr(Ar, br, t);
-    Fjr = computeFjr(Ajr, ur, Cjr, uj, pj);
+    Fjr = computeFjr(Ajr, ur, Cjr, uj, PTj);
   }
 
   Kokkos::View<Rd*> m_br;
@@ -313,21 +311,32 @@ public:
     Kokkos::View<double*> cj = unknowns.cj();
     Kokkos::View<double*> kj = unknowns.kj();
     Kokkos::View<double*> nuj = unknowns.nuj();
+    Kokkos::View<double*> PTj = unknowns.PTj();
 
     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();
     Kokkos::View<Rd*> xr = m_mesh.xr();
 
+    const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
+
+    const Kokkos::View<const unsigned short*> cell_nb_nodes
+      = m_connectivity.cellNbNodes();
+    
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+	double sum = 0;
+	for (int m=0; m<cell_nb_nodes(j); ++m) {
+	  //sum += uj[cell_nodes(j,m)][0];
+	  sum -= (uj(cell_nodes(j,m)), Cjr(cell_nodes(j,m), m));
+	}
+	PTj(j) = pj(j) - kj(j)*sum/Vj(j);
+      });
+
     // Calcule les flux
-    computeExplicitFluxes(xr, xj, kj, rhoj, uj, pj, cj, Vj, Cjr, t);
+    computeExplicitFluxes(xr, xj, kj, rhoj, uj, PTj, cj, Vj, Cjr, t);
 
     const Kokkos::View<const Rd**> Fjr = m_Fjr;
     const Kokkos::View<const Rd*> ur = m_ur;
-    const Kokkos::View<const unsigned short*> cell_nb_nodes
-      = m_connectivity.cellNbNodes();
-    const Kokkos::View<const unsigned int**>& cell_nodes
-      = m_connectivity.cellNodes();
 
     // Mise a jour de la vitesse et de l'energie totale specifique
     const Kokkos::View<const double*> inv_mj = unknowns.invMj();
@@ -364,12 +373,12 @@ public:
       });
 
     // Mise a jour de k
-    /*
+    
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
 	kj(j) = xj[j][0];
 
       });
-    */
+    
   }
 };
 
-- 
GitLab