Ajout fichier pour schema de NS sans splitting

src/main.cpp

@@ -19,6 +19,7 @@
 #include <Mesh.hpp>
 #include <AcousticSolver.hpp>
 #include <FiniteVolumesDiffusion.hpp>
+#include <NoSplitting.hpp>
 
 #include <TinyVector.hpp>
 #include <TinyMatrix.hpp>
@@ -136,6 +137,7 @@ int main(int argc, char *argv[])
 
     AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns);
     FiniteVolumesDiffusion<MeshDataType> finite_volumes_diffusion(mesh_data, unknowns);
+    NoSplitting<MeshDataType> no_splitting(mesh_data, unknowns);
 
     typedef TinyVector<MeshType::dimension> Rd;
 
@@ -166,6 +168,7 @@ int main(int argc, char *argv[])
     double c = 0.;
     c = finite_volumes_diffusion.conservatif(unknowns);
 
+    // Diagrammes de marche
     /*
     const Kokkos::View<const Rd*> xj = mesh_data.xj();
     int size = 3000;
@@ -175,6 +178,8 @@ int main(int argc, char *argv[])
     int i = 0;
     */
     
+    // Animations
+    /*
     // Ecriture des valeurs initiales dans un fichier
     
     const Kokkos::View<const Rd*> xj = mesh_data.xj();
@@ -271,27 +276,24 @@ int main(int argc, char *argv[])
       }
     }
     diff.close();
-    
+    */
 
 
     while((t<tmax) and (iteration<itermax)) {
      
-      
-      // ETAPE 1 DU SPLITTING - EULER
-      
-      double dt_euler = 0.9*acoustic_solver.acoustic_dt(Vj, cj);
-
+      // NAVIER-STOKES AVEC SPLITTING
+      /*
+      // Etape 1 du splitting - Euler
+      double dt_euler = 0.9*acoustic_solver.acoustic_dt(Vj, cj)
       if (t+dt_euler > tmax) {
 	dt_euler = tmax-t;
       }
       acoustic_solver.computeNextStep(t,dt_euler, unknowns);
       t += dt_euler;
       
-      // ETAPE 2 DU SPLITTING - DIFFUSION
-      
+      // Etape 2 du splitting - Diffusion
       double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
       double t_diff = t-dt_euler;
-      
       if (dt_euler <= dt_diff) {
 	dt_diff = dt_euler;
 	finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns);
@@ -305,11 +307,9 @@ int main(int argc, char *argv[])
 	  t_diff += dt_diff;
 	}
       }
-      
-
+      */
+      // Diffusion pure
       /*
-      // DIFFUSION PURE
-      
       double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj,nuL,nuR,kL,kR);
       if (t+dt > tmax) {
 	dt = tmax-t;
@@ -318,14 +318,24 @@ int main(int argc, char *argv[])
       t += dt;
       */
 
+      // NAVIER-STOKES SANS SPLITTING
+
+      double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR);
+      if (t+dt > tmax) {
+	dt = tmax-t;
+      }
+      no_splitting.computeNextStep(t,dt, unknowns);
+      t += dt;
+
+      
       block_eos.updatePandCFromRhoE();  
     
       ++iteration;
       std::cout << "temps t : " << t << std::endl;
 
       
-      // ECRITURE DANS UN FICHIER 
-      
+      // Animations
+      /*
       if ((std::fmod(t,0.001) < 0.0001) or (t == tmax)) {
     
       std::string ligne;
@@ -542,13 +552,13 @@ int main(int argc, char *argv[])
       riffout.close();
       
       }
+      */
  
- 
-      // ENTROPY TEST
+      // Entropy test
       //finite_volumes_diffusion.entropie(unknowns);
 
+      // Diagrammes de marche
       /*
-      // STOCKAGE COORDONNES ET TEMPS
       for (size_t j=0; j<mesh.numberOfCells(); ++j) {
 	x[i][j] = xj[j][0];
 	rho_marche[i][j] = rhoj[j];
@@ -564,21 +574,11 @@ int main(int argc, char *argv[])
     std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
 	      << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
-    // CREATION FICHIER x = f(t)
+
+    // Diagrammes de marche, creation fichier  x = f(t)
     /*
     std::ofstream fout("cara");
     fout.precision(15);
-    for (size_t j=0; j<mesh.numberOfCells(); ++j) {
-      for (size_t k=0; k<i; ++k) {
-	fout << tempo[k] << ' ' << x[k][j] << '\n';
-      }
-      fout << ' ' << '\n';
-      fout << ' ' << '\n';
-    }
-    */
-    /*
-    std::ofstream fout("cararho");
-    fout.precision(15);
     for (int j=0; j<mesh.numberOfCells(); ++j) {
       if (j%10 == 0) {
 	for (int k=0; k<i; ++k) {
@@ -592,56 +592,52 @@ int main(int argc, char *argv[])
     }
     */
     
+
+    // Erreurs sous differents normes
     /*
+
+    // Density
     double error1 = 0.;
     error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax);
-
     std::cout << "* " << rang::style::underline << "Erreur L2 rho" << rang::style::reset
 	      << ":  " << rang::fgB::green << error1 << rang::fg::reset << " \n";
     
     double error2 = 0.;
     error2 = finite_volumes_diffusion.error_Linf_rho(unknowns, tmax);
-
     std::cout << "* " << rang::style::underline << "Erreur L infini rho" << rang::style::reset
 	      << ":  " << rang::fgB::green << error2 << rang::fg::reset << " \n";
 
-
-    double err0 = 0.;
-    err0 = finite_volumes_diffusion.error_L1_u(unknowns, tmax);
-
-    std::cout << "* " << rang::style::underline << "Erreur L1 u" << rang::style::reset
-	      << ":  " << rang::fgB::green << err0 << rang::fg::reset << " \n";
-
+    // Vitesse
     double error = 0.;
     error = finite_volumes_diffusion.error_L2_u(unknowns, tmax);
-
     std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset
 	      << ":  " << rang::fgB::green << error << rang::fg::reset << " \n";
     
-    
     double error4 = 0.;
     error4 = finite_volumes_diffusion.error_Linf_u(unknowns, tmax);
-
     std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset
 	      << ":  " << rang::fgB::green << error4 << rang::fg::reset << " \n";
+    double err0 = 0.;
+    err0 = finite_volumes_diffusion.error_L1_u(unknowns, tmax);
+    std::cout << "* " << rang::style::underline << "Erreur L1 u" << rang::style::reset
+	      << ":  " << rang::fgB::green << err0 << rang::fg::reset << " \n";
     
-    double err1 = 0.;
-    err1 = finite_volumes_diffusion.error_L1_E(unknowns, tmax);
-
-    std::cout << "* " << rang::style::underline << "Erreur L1 E" << rang::style::reset
-	      << ":  " << rang::fgB::green << err1 << rang::fg::reset << " \n";
-
+    // Energy
     double error3 = 0.;
     error3 = finite_volumes_diffusion.error_L2_E(unknowns,tmax);
-
     std::cout << "* " << rang::style::underline << "Erreur L2 E" << rang::style::reset
 	      << ":  " << rang::fgB::green << error3 << rang::fg::reset << " \n";
     
     double error5 = 0.;
     error5 = finite_volumes_diffusion.error_Linf_E(unknowns, tmax);
-
     std::cout << "* " << rang::style::underline << "Erreur L infini E" << rang::style::reset
 	      << ":  " << rang::fgB::green << error5 << rang::fg::reset << " \n";
+
+    double err1 = 0.;
+    err1 = finite_volumes_diffusion.error_L1_E(unknowns, tmax);
+    std::cout << "* " << rang::style::underline << "Erreur L1 E" << rang::style::reset
+	      << ":  " << rang::fgB::green << err1 << rang::fg::reset << " \n";
+
     */
     	      
     std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset

src/mesh/Mesh.hpp

@@ -76,15 +76,14 @@ public:
 
   // pas constant
 
-  
   Mesh(const Connectivity& connectivity)
     : m_connectivity(connectivity),
       m_xr("xr", connectivity.numberOfNodes()),
       m_x0("x0", 1),
       m_xmax("xmax", 1)
   {
-    double a = -1.;
-    double b = 2.;
+    double a = 0.;
+    double b = 1.;
     m_x0[0][0] = a;
     m_xmax[0][0] = b;
     const double delta_x = (b-a)/connectivity.numberOfCells();
@@ -95,7 +94,6 @@ public:
   
 
   // pas non constant
-  
   /*
   Mesh(const Connectivity& connectivity)
   : m_connectivity(connectivity),
@@ -127,7 +125,6 @@ public:
 
 
   // pas aleatoire
-
   /*
   Mesh(const Connectivity& connectivity)
     : m_connectivity(connectivity),
@@ -163,7 +160,6 @@ public:
   */
 
 
-
   ~Mesh()
   {
     ;

src/scheme/AcousticSolver.hpp

@@ -182,18 +182,18 @@ private:
       });
 
     
-    m_ur[0]=zero;
-    m_ur[m_mesh.numberOfNodes()-1]=zero;
+    //m_ur[0]=zero;
+    //m_ur[m_mesh.numberOfNodes()-1]=zero;
     
     //m_ur[0] = x0;
     //m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
     
     // CL Kidder
-    /*
+    
     double h = std::sqrt(1. - (t*t)/(50./9.));
     m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0];
     m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0];
-    */
+    
 
     return m_ur;
   }
@@ -372,9 +372,9 @@ public:
 
       });
     */
-    /*
-    // Mise a jour de nu
 
+    // Mise a jour de nu
+    /*
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
 	nuj(j) = 0.5*(1.+xj[j][0]);
       });

src/scheme/FiniteVolumesEulerUnknowns.hpp

@@ -284,27 +284,29 @@ public:
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	
+	// TAC
+	/*
 	if (xj[j][0]<0.5) {
   	  m_rhoj[j]=1.;
   	} else {
   	  m_rhoj[j]=0.125;
 	}
-	
+	*/
 	// Kidder
-	//m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.);
+	m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.);
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	
+	// TAC
+	/*
   	if (xj[j][0]<0.5) {
   	  m_pj[j]=1;
   	} else {
   	  m_pj[j]=0.1;
 	}
-	
+	*/
 	// Kidder
-	//m_pj[j] = 2.*std::pow(m_rhoj[j],3);
+	m_pj[j] = 2.*std::pow(m_rhoj[j],3);
       });
     
     double pi = 4.*std::atan(1.);
@@ -313,8 +315,10 @@ public:
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-  	m_gammaj[j] = 1.4;
-	//m_gammaj[j] = 3.;
+  	// TAC
+	// m_gammaj[j] = 1.4;
+	// Kidder
+	m_gammaj[j] = 3.;
       });
 
     BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
@@ -334,15 +338,14 @@ public:
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-  	//m_kj[j] =  xj[j][0];
+	// Differents k (xi)
+  	m_kj[j] =  xj[j][0];
+	//m_kj[j] = 0.;
 
-	m_kj[j] = 0.;
-
-	/*
-	// Sod
+	// TAC
 	
 	// k non regulier
-	
+	/*
 	if (xj[j][0]<0.7) {
   	  m_kj[j]=0.;
   	} else {
@@ -374,8 +377,9 @@ public:
 	  }
 	}
 	*/
-	/*
+   
 	// k regulier
+	/*
 	int n = 1;
 	m_kj[j] = std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.7+0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.7)*(xj[j][0]<0.7+0.1/n) + std::exp(1.)*std::exp(-1./(1.-( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) )*( (xj[j][0]-(0.9-0.1/n)) / (0.1/n) ))) * (xj[j][0]>0.9-0.1/n)*(xj[j][0]<0.9) + (xj[j][0]>0.7+0.1/n)*(xj[j][0]<0.9-0.1/n);
 	m_kj[j] = 0.014*m_kj[j];
@@ -388,7 +392,7 @@ public:
     m_uL[0] = zero;
     m_uR[0] = zero;
     m_kL[0] = 0.;
-    m_kR[0] = 0.;
+    m_kR[0] = 1.;
 
     
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){  	
@@ -414,7 +418,6 @@ public:
   }
   
  
- 
   /*
   // DIFFUSION PURE
 
@@ -432,7 +435,9 @@ public:
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+        // Sans conduction thermique
 	m_uj[j] = std::sin(pi*xj[j][0]);
+	// Avec conduction thermique
 	//m_uj[j] = xj[j][0];
       });
 
@@ -444,8 +449,9 @@ public:
     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]);
+        // Sans conduction thermique
 	m_Ej[j] = 2.;
+	// Avec conduction thermique
 	//m_Ej[j] = xj[j][0]*xj[j][0]*0.5 + 2.*xj[j][0] + 1.;
       });
 

src/scheme/NoSplitting.hpp

+#ifndef NO_SPLITTING_HPP
+#define NO_SPLITTING_HPP
+
+// --- INCLUSION fichiers headers ---
+
+#include <Kokkos_Core.hpp>
+
+#include <rang.hpp>
+#include <BlockPerfectGas.hpp>
+
+#include <TinyVector.hpp>
+#include <TinyMatrix.hpp>
+#include <Mesh.hpp>
+#include <MeshData.hpp>
+#include <FiniteVolumesEulerUnknowns.hpp>
+
+// ---------------------------------
+
+// Creation classe NoSplitting
+
+template<typename MeshData> 
+class NoSplitting 
+{
+  typedef typename MeshData::MeshType MeshType; 
+  typedef FiniteVolumesEulerUnknowns<MeshData> UnknownsType; 
+
+  MeshData& m_mesh_data; 
+  MeshType& m_mesh;
+  const typename MeshType::Connectivity& m_connectivity;
+
+  constexpr static size_t dimension = MeshType::dimension; 
+
+  typedef TinyVector<dimension> Rd; 
+  typedef TinyMatrix<dimension> Rdd;
+
+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
+    {
+      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*>
+  computePj(const Kokkos::View<const double*>& pj,
+	    const Kokkos::View<const double*>& kj,
+	    const Kokkos::View<const Rd*>& uj)
+  {
+    const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
+
+    const Kokkos::View<const unsigned short*> cell_nb_nodes
+      = m_connectivity.cellNbNodes();
+
+    const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
+    
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+	double new_p = 0.;
+	double sum = 0;
+	for (int m=0; m<cell_nb_nodes(j); ++m) {
+	  sum += uj[cell_nodes(j,m)][0];
+	}
+	new_p = pj(j) - kj(j)*sum/Vj(j);
+	pj(j) = new_p;
+      });
+    return pj;
+  }
+
+  // ----
+
+
+  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) {
+	m_rhocj[j] = rhoj[j]*cj[j];
+      });
+    return m_rhocj;
+  }
+
+  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 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));
+	}
+      });
+
+    return m_Ajr;
+  }
+
+  KOKKOS_INLINE_FUNCTION 
+  const Kokkos::View<const Rdd*> 
+  computeAr(const Kokkos::View<const Rdd**>& Ajr) {
+    const Kokkos::View<const unsigned int**> node_cells = m_connectivity.nodeCells();
+    const Kokkos::View<const unsigned short**> node_cell_local_node = m_connectivity.nodeCellLocalNode();
+    const Kokkos::View<const unsigned short*> node_nb_cells = m_connectivity.nodeNbCells();
+
+    Kokkos::parallel_for(m_mesh.numberOfNodes(), 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*>
+  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 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();
+    const Kokkos::View<const unsigned short*>& node_nb_cells = m_connectivity.nodeNbCells();
+    Kokkos::View<Rd*> xr = m_mesh.xr();
+
+    Kokkos::parallel_for(m_mesh.numberOfNodes(), 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::View<Rd*>  
+  computeUr(const Kokkos::View<const Rdd*>& Ar,
+	    const Kokkos::View<const Rd*>& br,
+	    const double& t) {
+
+    inverse(Ar, m_inv_Ar);
+    const Kokkos::View<const Rdd*> invAr = m_inv_Ar;
+
+    Kokkos::View<Rd*> xr = m_mesh.xr();
+    Kokkos::View<Rd*> x0 = m_mesh.x0();
+    Kokkos::View<Rd*> xmax = m_mesh.xmax();
+
+    Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+	m_ur[r]=invAr(r)*br(r);
+      });
+
+    
+    //m_ur[0]=zero;
+    //m_ur[m_mesh.numberOfNodes()-1]=zero;
+    
+    //m_ur[0] = x0;
+    //m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
+    
+    // CL Kidder
+    
+    double h = std::sqrt(1. - (t*t)/(50./9.));
+    m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0];
+    m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0];
+    
+
+    return m_ur;
+  }
+  
+  Kokkos::View<Rd**>  
+  computeFjr(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) {
+    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);
+	}
+      });
+
+    return m_Fjr;
+  }
+
+  // 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
+  // pouvoir derouler le schema
+  KOKKOS_INLINE_FUNCTION
+  void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
+			     const Kokkos::View<const Rd*>& xj,
+			     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*>& cj,
+			     const Kokkos::View<const double*>& Vj,
+			     const Kokkos::View<const Rd**>& Cjr,
+			     const double& t) {
+    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);
+  }
+
+  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_PTj;
+  Kokkos::View<double*> m_Vj_over_cj;
+
+public:
+  NoSplitting(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_PTj("PTj", m_mesh.numberOfCells()),
+    m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells())
+  {
+    ;
+  }
+
+  // Avance la valeur des inconnues pendant un pas de temps dt
+  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();
+    Kokkos::View<double*> kj = unknowns.kj();
+    Kokkos::View<double*> nuj = unknowns.nuj();
+
+    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();
+
+    // Calcule les flux
+    computeExplicitFluxes(xr, xj, kj, rhoj, uj, pj, 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();
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+	Rd momentum_fluxes = zero;
+	double energy_fluxes = 0;
+	for (int R=0; R<cell_nb_nodes[j]; ++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;
+      });
+
+    // Calcul de e par la formule e = E-0.5 u^2
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+	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){
+	rhoj[j] = mj[j]/Vj[j];
+      });
+
+    // Mise a jour de k
+    /*
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+	kj(j) = xj[j][0];
+
+      });
+    */
+  }
+};
+
+#endif // NO_SPLITTING_HPP