Code Fjr et Gjr

src/scheme/FiniteVolumesDiffusion.hpp

@@ -21,7 +21,7 @@
 template<typename MeshData> // MeshData est le type generique des donnees (geometriques) attachees a un maillage
 class FiniteVolumesDiffusion 
 {
-  typedef typename MeshData::MeshType MeshType; // de type du maillage
+  typedef typename MeshData::MeshType MeshType; // type du maillage
   typedef FiniteVolumesEulerUnknowns<MeshData> UnknownsType; // type des inconnues
 
   MeshData& m_mesh_data; // reference vers les donnees attachees du maillage
@@ -73,44 +73,57 @@ private:
     }
   };
 
-  Kokkos::View<Rd**>  // Fonction qui calcule F_jr // A MODIFIER // Ajout de k
-  computeFjr(const Kokkos::View<const Rdd**>& Ajr,
-	     const Kokkos::View<const Rd*>& ur,
-	     const Kokkos::View<const Rd**>& Cjr,
+  Kokkos::View<Rd**>  // Fonction qui calcule F_jr 
+  computeFjr(const Kokkos::View<const Rd**>& Cjr,
 	     const Kokkos::View<const Rd*>& uj,
-	     const Kokkos::View<const double*>& pj) {
+	     const Kokkos::View<const double*>& Vj,
+	     const Kokkos::View<const double*>& kj) {
 	       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) = ((kjr(j,r)*Cjr(j,r) + kjr(j,r)*Cjr(j,r))/2)*((ujr(j,r)*Cjr(j,r) + ujr(j,r-1)*Cjr(j,r-1))/h(r));
+		     m_Fjr(j,r) = ((kj(j,r) + kj(j,r-1))/(2*Vj(j)))*(uj(j,r)*Cjr(j,r) + uj(j,r-1)*Cjr(j,r-1)); 
 		   }
 		 });
 
 	       return m_Fjr;
 	     }
 
-  Kokkos::View<Rd**>  // Fonction qui calcule G_jr // A MODIFIER // Ajout de k
-  computeGjr(const Kokkos::View<const Rdd**>& Ajr,
-	     const Kokkos::View<const Rd*>& ur,
-	     const Kokkos::View<const Rd**>& Cjr,
-	     const Kokkos::View<const Rd*>& uj,
-	     const Kokkos::View<const double*>& pj) {
+  Kokkos::View<Rd**>  // Fonction qui calcule G_jr 
+  computeGjr(const Kokkos::View<const Rd*>& uj,
+	     const Kokkos::View<const double*>& Fjr, {
     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_Gjr(j,r) = ((ujr(j,r)*Cjr(j,r) + ujr(j,r-1)*Cjr(j,r))/2)*Fjr(j,r);
+	  m_Gjr(j,r) = ((ujr(j,r) + ujr(j,r-1))/2)*Fjr(j,r);
 	}
       });
 
     return m_Gjr;
   }
 
+  // Calcul la liste des inverses d'une liste de matrices (pour
+  // l'instant seulement $R^{1\times 1}$)
+  void 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))};
+      });
+  }
+
+  // Calcul la liste des inverses d'une liste de reels
+  void 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);
+      });
+  }
+
 
   // Enchaine les operations pour calculer les flux (Fjr et ur) pour // A MODIFIER
   // pouvoir derouler le schema
@@ -119,59 +132,36 @@ private:
 			     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**>& Cjr) {
-    const Kokkos::View<const double*> rhocj  = computeRhoCj(rhoj, cj);
-    const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, Cjr);
-
-    const Kokkos::View<const Rdd*> Ar = computeAr(Ajr);
-    const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj);
-
-    Kokkos::View<Rd*> ur = m_ur;
     Kokkos::View<Rd**> Fjr = m_Fjr; 
-    ur  = computeUr(Ar, br);
-    Fjr = computeFjr(Ajr, ur, Cjr, uj, pj);
+    Fjr = computeFjr(Cjr, uj, xj,kj);
+    Kokkos::View<Rd**> Gjr = m_Gjr; 
+    Gjr = computeGjr(uj, Fjr);
   }
 
-  Kokkos::View<Rd*> m_br;
-  Kokkos::View<Rdd**> m_Ajr;
-  Kokkos::View<Rdd*> m_Ar;
-  Kokkos::View<Rdd*> m_inv_Ar;
   Kokkos::View<Rd**> m_Fjr;
-  Kokkos::View<Rd*> m_ur;
-  Kokkos::View<double*> m_rhocj;
-  Kokkos::View<double*> m_Vj_over_cj;
+  Kokkos::View<double*> m_CFL;
 
 public:
-  AcousticSolver(MeshData& mesh_data,
+  DiffusionSolver(MeshData& mesh_data,
 		 UnknownsType& unknowns)
     : m_mesh_data(mesh_data),
       m_mesh(mesh_data.mesh()),
       m_connectivity(m_mesh.connectivity()),
-      m_br("br", m_mesh.numberOfNodes()),
-      m_Ajr("Ajr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
-      m_Ar("Ar", m_mesh.numberOfNodes()),
-      m_inv_Ar("inv_Ar", m_mesh.numberOfNodes()),
       m_Fjr("Fjr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
-      m_ur("ur", m_mesh.numberOfNodes()),
-      m_rhocj("rho_c", m_mesh.numberOfCells()),
-      m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells())
+      m_CFL("CFL", m_mesh.numberOfCells())
   {
     ;
   }
 
-  // Calcule une evaluation du pas de temps verifiant une CFL du type
-  // c*dt/dx<1. Utilise la reduction definie dans la structure // A MODIFIER
-  // ReduceMin. Ici, dx_j=V_j
+  // Calcule une evaluation du pas de temps verifiant le CFL parabolique
+  // Utilise la reduction definie dans la structure ReduceMin. Ici, dx_j=V_j
   KOKKOS_INLINE_FUNCTION
-  double acoustic_dt(const Kokkos::View<const double*>& Vj,
-		     const Kokkos::View<const double*>& cj) const {
+  double diffusion_dt(const Kokkos::View<const double*>& Vj) const {
     Kokkos::View<double*> dt_j("dt_j", m_mesh.numberOfCells());
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	m_Vj_over_cj[j] = Vj[j]/cj[j]; // La tu mets la formule pour chaque maille
+	m_CFL[j] = rhoj[j]*(xr(r); // La tu mets la formule pour chaque maille // CONFUSION j et r !!! 
       });
 
     double dt = std::numeric_limits<double>::max();
@@ -189,9 +179,7 @@ public:
     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();
@@ -199,10 +187,9 @@ public:
     Kokkos::View<Rd*> xr = m_mesh.xr();
 
     // Calcule les flux
-    computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr);
+    computeExplicitFluxes(xr, xj, rhoj, uj, Vj, Cjr);
 
     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
@@ -227,15 +214,6 @@ public:
 	ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
       });
 
-    // deplace le maillage (ses sommets) en utilisant la vitesse
-    // donnee par le schema
-    Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
-	xr[r] += dt*ur[r];
-      });
-
-    // met a jour les quantites (geometriques) associees au maillage
-    m_mesh_data.updateAllData();
-
     // Calcul de rho avec la formule Mj = Vj rhoj
     const Kokkos::View<const double*> mj = unknowns.mj();
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){