From a49c4dcf29e0bdb5bf2d02e0e1e97a3fff565c14 Mon Sep 17 00:00:00 2001
From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr>
Date: Fri, 22 Jun 2018 16:33:09 +0200
Subject: [PATCH] ajout transfert thermique, pb avec membre source pour
 equation de l energie pour le test diffusion pure avec k non constant

---
 src/main.cpp                              |  40 +++---
 src/mesh/Mesh.hpp                         |   4 +-
 src/scheme/AcousticSolver.hpp             |  16 ++-
 src/scheme/FiniteVolumesDiffusion.hpp     |  87 +++++++++++-
 src/scheme/FiniteVolumesEulerUnknowns.hpp | 156 ++++++++++++++++------
 5 files changed, 236 insertions(+), 67 deletions(-)

diff --git a/src/main.cpp b/src/main.cpp
index 28d36d13b..f54a78561 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -142,7 +142,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();
@@ -160,12 +160,12 @@ int main(int argc, char *argv[])
     double c = 0.;
     c = finite_volumes_diffusion.conservatif(unknowns);
 
-    
+    /*
     // Ecriture des valeurs initiales dans un fichier
     
     const Kokkos::View<const Rd*> xj = mesh_data.xj();
     const Kokkos::View<const Rd*> xr = mesh.xr();
-    /*
+    
     // rho
     std::ofstream fout1("inter1", std::ios::trunc);
     fout1.precision(15);
@@ -244,23 +244,24 @@ int main(int argc, char *argv[])
     tempo.precision(5);
     tempo << std::fixed << t << '\n';
     tempo.close();
-    */
+    
     // Fichier k indicateur
     std::ofstream diff("diffinter");
     diff.precision(5);
     for (size_t j=0; j<mesh.numberOfCells(); ++j) {
       diff << std::fixed << xj[j][0] << ' ' << std::fixed << kj[j] << '\n';
-      /*if (kj[j]>0.) {
+      if (kj[j]>0.) {
 	diff << std::fixed << xj[j][0] << ' ' << std::fixed << 4. << '\n';
       } else {
 	diff << std::fixed << xj[j][0] << ' ' << std::fixed << -0.1 << '\n';
-      }*/
+      }
     }
     diff.close();
-    
+    */
 
     while((t<tmax) and (iteration<itermax)) {
      
+      /*
       // ETAPE 1 DU SPLITTING - EULER
       
       double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
@@ -289,7 +290,7 @@ int main(int argc, char *argv[])
 	  t_diff += dt_diff;
 	}
       }
-      
+      */
       
       // AUTRE APPROCHE DU SPLITTING (PLUS LONG)
       /*
@@ -311,14 +312,13 @@ int main(int argc, char *argv[])
 
       // DIFFUSION PURE
       
-      /*
-      double dt = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj,kj);
+      double dt = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj,kj,cj);
       if (t+dt > tmax) {
 	dt = tmax-t;
       }
       finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
       t += dt;
-      */
+      
       
       block_eos.updatePandCFromRhoE();  
     
@@ -615,32 +615,37 @@ int main(int argc, char *argv[])
      { // gnuplot output for vitesse
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const Rd*> uj = unknowns.uj();
-     //double pi = 4.*std::atan(1.);
+     double pi = 4.*std::atan(1.);
      std::ofstream fout("u");
      fout.precision(15);
      for (size_t j=0; j<mesh.numberOfCells(); ++j) {
+
        //fout << xj[j][0] << ' ' << uj[j][0] <<  ' ' << std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2) <<'\n'; //cas k constant
-       //fout << xj[j][0] << ' ' << uj[j][0] <<  ' ' << std::sin(pi*xj[j][0])*std::exp(-0.2) <<'\n'; // cas k non constant
+       fout << xj[j][0] << ' ' << uj[j][0] <<  ' ' << std::sin(pi*xj[j][0])*std::exp(-tmax) <<'\n'; // cas k non constant
        //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder
-       fout << xj[j][0] << ' ' << uj[j][0] << '\n';
+       
+       //fout << xj[j][0] << ' ' << uj[j][0] << '\n';
      }
      }
 
      { // gnuplot output for energy
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const double*> Ej = unknowns.Ej();
-     //double pi = 4.*std::atan(1.);
+     double pi = 4.*std::atan(1.);
      //double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
      std::ofstream fout("E");
      fout.precision(15);
      for (size_t j=0; j<mesh.numberOfCells(); ++j) {
+
        //fout << xj[j][0] << ' ' << Ej[j] << ' ' << (-(std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0]))+(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])))*0.5*(std::exp(-4.*pi*pi*0.2)-1.) + 2. <<'\n'; // cas k constant
-       //fout << xj[j][0] << ' ' << Ej[j] << ' ' << ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*0.2)-1.) + 2. <<'\n' ; // cas k non constant
+       fout << xj[j][0] << ' ' << Ej[j] << ' ' << ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*tmax)-1.) + 2. <<'\n' ; // cas k non constant
        //fout << xj[j][0] << ' ' << Ej[j] << ' ' << (std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h) + (-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*(-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*0.5 << '\n'; // kidder
-       fout << xj[j][0] << ' ' << Ej[j] << '\n';
+
+       //fout << xj[j][0] << ' ' << Ej[j] << '\n';
      }
      }
 
+     /*
      { // gnuplot output for entropy (gaz parfait en prenant cv = 1))
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const double*> rhoj = unknowns.rhoj();
@@ -653,6 +658,7 @@ int main(int argc, char *argv[])
        fout << xj[j][0] << ' ' << std::log(pj[j]*std::pow(rhoj[j],-gammaj[j])) << ' ' << S0(j) << '\n';
      }
      }
+     */
 
   }
 
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 470236c9e..1f26bcac1 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -80,8 +80,8 @@ public:
       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();
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 1e8df86b3..f9cb525e0 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -333,6 +333,7 @@ public:
     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();
@@ -384,9 +385,11 @@ public:
       });
 
     // Mise a jour de k
-    /*
+    
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
-	//kj(j) = xj[j][0];
+	kj(j) = xj[j][0];
+	
+	/*	
 	//kj(j) = 0.5;
 	
 	if (xj[j][0]<0.7) {
@@ -398,9 +401,16 @@ public:
 	    kj[j]=0. ;
 	  }
 	}
+	*/
 	
       });
-    */
+
+    // Mise a jour de nu
+    
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+	nuj(j) = 0.5*(1.+xj[j][0]);
+      });
+    
   }
 };
 
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 5f477bd87..49f29b277 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -195,6 +195,57 @@ private:
 
   }
 
+  // Calcule un Bl
+  Kokkos::View<double*> 
+  computeBl(const Kokkos::View<const Rd**>& Cjr,
+	    const Kokkos::View<const double*>& Tj,
+	    const Kokkos::View<const double*>& nuj,
+	    const Kokkos::View<const double*>& TL,
+	    const Kokkos::View<const double*>& TR,
+	    const Kokkos::View<const double*>& nuL,
+	    const Kokkos::View<const double*>& nuR,
+	    const double& t) {
+
+    const Kokkos::View<const unsigned int**>& face_cells = m_connectivity.faceCells();
+
+    const Kokkos::View<const unsigned short*> face_nb_cells
+      = m_connectivity.faceNbCells();
+
+    const Kokkos::View<const unsigned short**> face_cell_local_face
+      = m_mesh.connectivity().faceCellLocalFace();
+
+    const Kokkos::View<const double*>& Vl = m_mesh_data.Vl();
+    const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
+
+    Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
+        Rd sum = zero;
+	double sum2 = 0.;
+	double sum3 = 0.;
+	for (int j=0; j<face_nb_cells(l); ++j) {
+	    int cell_here = face_cells(l,j);
+	    int local_face_number_in_cell = face_cell_local_face(l,j);
+	    sum -= Tj(cell_here)*Cjr(cell_here, local_face_number_in_cell);
+	    sum2 += (1./Vj(cell_here))*nuj(cell_here);
+	    sum3 += 1./Vj(cell_here);
+	}
+
+	m_Bl(l) = ((sum2/sum3)/Vl(l))*sum[0];
+	
+      });
+
+    // Conditions aux bords
+    
+    int cell_here = face_cells(0,0);
+    m_Bl(0) = (nuL(0) + nuj(cell_here))*(1./(2*Vl(0)))*(Tj(cell_here) - TL(0));
+
+    cell_here = face_cells(m_mesh.numberOfFaces()-1,0);
+    m_Bl(m_mesh.numberOfFaces()-1) = -(nuR(0) + nuj(cell_here))*(1/(2.*Vl(m_mesh.numberOfFaces()-1)))*(Tj(cell_here) - TR(0));
+
+    return m_Bl ;
+  }
+
+
+
   // Calcule la liste des inverse d'une liste de matrices 
   // (pour l'instant, juste 1x1)
   void inverse(const Kokkos::View<const Rdd*>& A,
@@ -223,15 +274,24 @@ private:
 			     const Kokkos::View<const Rd*>& uR,
 			     const Kokkos::View<const double*>& kL,
 			     const Kokkos::View<const double*>& kR,
+			     const Kokkos::View<const double*>& Tj,
+			     const Kokkos::View<const double*>& nuj,
+			     const Kokkos::View<const double*>& TL,
+			     const Kokkos::View<const double*>& TR,
+			     const Kokkos::View<const double*>& nuL,
+			     const Kokkos::View<const double*>& nuR,
 			     const double& t) { 
     Kokkos::View<Rdd*> Fl  = m_Fl ; 
     Fl  = computeFl (Cjr, uj, kj, uL, uR, kL, kR, t);
     Kokkos::View<Rd*> Gl  = m_Gl ; 
     Gl  = computeGl (uj, Fl, uL, uR, t);
+    Kokkos::View<double*> Bl  = m_Bl ; 
+    Bl  = computeBl (Cjr, Tj, nuj, TL, TR, nuL, nuR, t);
   }
 
   Kokkos::View<Rdd*> m_Fl;
   Kokkos::View<Rd*> m_Gl;
+  Kokkos::View<double*> m_Bl;
 
 public:
   FiniteVolumesDiffusion(MeshData& mesh_data,
@@ -240,7 +300,8 @@ public:
       m_mesh(mesh_data.mesh()),
       m_connectivity(m_mesh.connectivity()),
       m_Fl("Fl", m_mesh.numberOfFaces()), 
-      m_Gl("Gl", m_mesh.numberOfFaces())
+      m_Gl("Gl", m_mesh.numberOfFaces()),
+      m_Bl("Bl", m_mesh.numberOfFaces())
   {
     ;
   }
@@ -281,7 +342,7 @@ public:
 	}
 
 	if (sum == 0.) {
-	  dt_j[j] = std::numeric_limits<double>::max();//Vj[j]/cj[j]; 
+	  dt_j[j] = std::numeric_limits<double>::max(); 
 	} else {
 	  dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./sum)*minVl;
 	}
@@ -300,29 +361,36 @@ public:
   {
     Kokkos::View<double*> rhoj = unknowns.rhoj();
     Kokkos::View<Rd*> uj = unknowns.uj();
-    Kokkos::View<Rd*> Sj = unknowns.Sj();
     Kokkos::View<double*> Ej = unknowns.Ej();
+    Kokkos::View<double*> Tj = unknowns.Tj();
 
     Kokkos::View<double*> ej = unknowns.ej();
     Kokkos::View<double*> gammaj = unknowns.gammaj();
 
     Kokkos::View<double*> kj = unknowns.kj();
+    Kokkos::View<double*> nuj = unknowns.nuj();
 
     Kokkos::View<Rd*> uL = unknowns.uL();
     Kokkos::View<Rd*> uR = unknowns.uR();
     Kokkos::View<double*> kL = unknowns.kL();
     Kokkos::View<double*> kR = unknowns.kR();
 
+    Kokkos::View<double*> TL = unknowns.TL();
+    Kokkos::View<double*> TR = unknowns.TR();
+    Kokkos::View<double*> nuL = unknowns.nuL();
+    Kokkos::View<double*> nuR = unknowns.nuR();
+
     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();
 
     // Calcule les flux
-    computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, t);
+    computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, Tj, nuj, TL, TR, nuL, nuR, t);
 
     const Kokkos::View<const Rdd*> Fl  = m_Fl ;
     const Kokkos::View<const Rd *> Gl  = m_Gl ;
+    const Kokkos::View<const double*> Bl  = m_Bl ;
     
     const Kokkos::View<const unsigned short*>& cell_nb_faces
       = m_connectivity.cellNbFaces();
@@ -330,22 +398,29 @@ public:
     const Kokkos::View<const unsigned int*[2]>& cell_faces
       = m_connectivity.cellFaces();
 
+    double pi = 4.*std::atan(1.);
+
     // 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) {
 	//const int j = j0+1;
 	Rd momentum_fluxes = zero;
 	double energy_fluxes = 0.;
+	Rd trich = zero;
 	int l = 0;
 	for (int R=0; R<cell_nb_faces(j); ++R) {
 	  l = cell_faces(j,R);
 	  momentum_fluxes +=  Fl(l)*Cjr(j,R);
-	  energy_fluxes   += (Gl(l), Cjr(j,R));
+	  trich = Bl(l)*Cjr(j,R);
+	  energy_fluxes   += (Gl(l), Cjr(j,R)) + trich[0];
 	}
     
 	uj[j] += (dt*inv_mj[j]) * momentum_fluxes;
 	Ej[j] +=  (dt*inv_mj[j]) * energy_fluxes;
-	//uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*Sj(j) + (dt*inv_mj[j]) * momentum_fluxes; // test avec k non constant
+	
+	// test k non cst pour diff pure
+	uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*(std::sin(pi*xj[j][0])*(pi*pi*xj[j][0]-1.) - std::cos(xj[j][0]*pi)*pi); 
+	Ej[j] -= ((pi*pi*0.5*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])) + pi*pi*pi*std::cos(xj[j][0])*std::sin(xj[j][0]))*(std::exp(-2.*t)-1.) - pi*0.5*std::exp(-2.*t)*std::cos(xj[j][0])*std::sin(xj[j][0]) + (1.+ xj[j][0])*( (2*pi*pi*pi*std::cos(xj[j][0])*std::sin(xj[j][0])-4*std::pow(pi,4)*(std::sin(xj[j][0])*std::sin(xj[j][0])-std::cos(xj[j][0])*std::cos(xj[j][0])))*(std::exp(-2.*t)-1.) + pi*pi*0.5*std::exp(-2.*t)*(std::sin(xj[j][0])*std::sin(xj[j][0])-std::cos(xj[j][0])*std::cos(xj[j][0]))))*(dt*inv_mj[j])*Vj(j);
 
 	// ajout second membre pour kidder (k = 0.5)
 	//Ej[j] -= (dt*inv_mj[j])*Vj(j)*((0.5*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 29445d094..334450d06 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -26,15 +26,23 @@ private:
   Kokkos::View<double*> m_mj;
   Kokkos::View<double*> m_inv_mj;
   Kokkos::View<double*> m_kj;
-  Kokkos::View<Rd*> m_Sj;
+
   Kokkos::View<Rd*> m_uL;
   Kokkos::View<Rd*> m_uR;
   Kokkos::View<double*> m_kL;
   Kokkos::View<double*> m_kR;
+
   Kokkos::View<double*> m_entropy;
   Kokkos::View<double*> m_entropy_new;
   Kokkos::View<double*> m_S0;
 
+  Kokkos::View<double*> m_Tj;
+  Kokkos::View<double*> m_nuj;
+  Kokkos::View<double*> m_TL;
+  Kokkos::View<double*> m_TR;
+  Kokkos::View<double*> m_nuL;
+  Kokkos::View<double*> m_nuR;
+
 public:
   Kokkos::View<double*> rhoj()
   {
@@ -136,16 +144,6 @@ public:
     return m_kj;
   }
 
-   Kokkos::View<Rd*> Sj()
-  {
-   return m_Sj;
-  }
-
-  const Kokkos::View<const Rd*> Sj() const
-  {
-    return m_Sj;
-  }
-
   Kokkos::View<Rd*> uL()
   {
     return m_uL;
@@ -186,6 +184,67 @@ public:
     return m_kR;
   }
 
+  Kokkos::View<double*> nuj()
+  {
+    return m_nuj;
+  }
+
+  const Kokkos::View<const double*> nuj() const
+  {
+    return m_nuj;
+  }
+
+  Kokkos::View<double*> Tj()
+  {
+    return m_Tj;
+  }
+
+  const Kokkos::View<const double*> Tj() const
+  {
+    return m_Tj;
+  }
+
+  Kokkos::View<double*> TL()
+  {
+    return m_TL;
+  }
+
+  const Kokkos::View<const double*> TL() const
+  {
+    return m_TL;
+  }
+
+  Kokkos::View<double*> TR()
+  {
+    return m_TR;
+  }
+
+  const Kokkos::View<const double*> TR() const
+  {
+    return m_TR;
+  }
+
+  Kokkos::View<double*> nuL()
+  {
+    return m_nuL;
+  }
+
+  const Kokkos::View<const double*> nuL() const
+  {
+    return m_nuL;
+  }
+
+  Kokkos::View<double*> nuR()
+  {
+    return m_nuR;
+  }
+
+  const Kokkos::View<const double*> nuR() const
+  {
+    return m_nuR;
+  }
+  
+  
   Kokkos::View<double*> entropy()
   {
     return m_entropy;
@@ -216,6 +275,8 @@ public:
     return m_S0;
   }
 
+
+  /*
   // --- Acoustic Solver ---
   
   void initializeSod()
@@ -223,7 +284,6 @@ public:
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	
 	if (xj[j][0]<0.5) {
   	  m_rhoj[j]=1.;
   	} else {
@@ -235,7 +295,6 @@ public:
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	
   	if (xj[j][0]<0.5) {
   	  m_pj[j]=1;
   	} else {
@@ -245,7 +304,8 @@ public:
 	//Kidder
 	//m_pj[j] = 2.*std::pow(m_rhoj[j],3);
       });
-
+    
+    double pi = 4.*std::atan(1.);
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
   	m_uj[j] = zero;
       });
@@ -272,13 +332,14 @@ public:
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-  	//m_kj[j] =  xj[j][0];
+  	m_kj[j] =  xj[j][0];
 	//m_kj[j] = 0.5;
 
+	
 	// Sod
 
 	// k non regulier
-	/*
+	
 	if (xj[j][0]<0.7) {
   	  m_kj[j]=0.;
   	} else {
@@ -309,13 +370,14 @@ public:
 	    m_kj[j]=0. ;
 	  }
 	}
-	*/
+	
 
 	// k regulier
 	int n = 10.;
 	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];
-
+	
+	
       });
 
      // Conditions aux bords de Dirichlet sur u et k
@@ -323,17 +385,20 @@ 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){  	
 	m_entropy(j) = std::log(m_pj[j]*std::pow(m_rhoj[j],-m_gammaj[j]));
 	m_S0(j) = m_entropy(j);
       });
+    
   }
+  */
   
   
-  /* 
-
+   
+  
   // DIFFUSION PURE
 
     void initializeSod()
@@ -375,24 +440,32 @@ public:
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	m_kj[j] = 2.; // k constant
-	//m_kj[j] = xj[j][0]; // k non constant, k = x
-      });
-
-    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-    	m_Sj[j] = std::sin(pi*xj[j][0])*(pi*pi*xj[j][0]-1.) - std::cos(xj[j][0]*pi)*pi;
+	//m_kj[j] = 2.; // k constant
+	m_kj[j] = xj[j][0]; // k non constant, k = x
       });
 
     // Conditions aux bords de Dirichlet sur u et k
     
     m_uL[0] = zero;
     m_uR[0] = zero;
-    m_kL[0] = 2.;
-    m_kR[0] = 2.;
+    m_kL[0] = 0.;
+    m_kR[0] = 1.;
 
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+	m_nuj(j) = 0.5*(1.+xj[j][0]); // k = x
+      });
+
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+	m_Tj[j] = 2 - 0.5*std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]); // k = x
+      });
+
+    // Conditions aux bords de Dirichlet sur u et k
+    
+    m_TL[0] = 2.;
+    m_TR[0] = 2.;
+    m_nuL[0] = 0.5;
+    m_nuR[0] = 1.;
   }
-*/
-  
 
 
   FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)
@@ -401,6 +474,8 @@ public:
       m_rhoj("rhoj",m_mesh.numberOfCells()),
       m_uj("uj",m_mesh.numberOfCells()),
       m_Ej("Ej",m_mesh.numberOfCells()),
+      m_nuj("nuj",m_mesh.numberOfCells()),
+      m_Tj("Tj",m_mesh.numberOfCells()),
       m_ej("ej",m_mesh.numberOfCells()),
       m_pj("pj",m_mesh.numberOfCells()),
       m_gammaj("gammaj",m_mesh.numberOfCells()),
@@ -408,14 +483,17 @@ public:
       m_mj("mj",m_mesh.numberOfCells()),
       m_inv_mj("inv_mj",m_mesh.numberOfCells()),
       m_kj("kj",m_mesh.numberOfCells()),
-      m_Sj("Sj",m_mesh.numberOfCells()),
-    m_uL("uL", 1),
-    m_uR("uR", 1),
-    m_kL("kL", 1),
-    m_kR("kR", 1),
-    m_entropy("entropy", m_mesh.numberOfCells()),
-    m_entropy_new("entropy_new", m_mesh.numberOfCells()),
-    m_S0("S0", m_mesh.numberOfCells())
+      m_uL("uL", 1),
+      m_uR("uR", 1),
+      m_kL("kL", 1),
+      m_kR("kR", 1),
+      m_TL("TL", 1),
+      m_TR("TR", 1),
+      m_nuL("nuL", 1),
+      m_nuR("nuR", 1),
+      m_entropy("entropy", m_mesh.numberOfCells()),
+      m_entropy_new("entropy_new", m_mesh.numberOfCells()),
+      m_S0("S0", m_mesh.numberOfCells())
   {
     ;
   }
-- 
GitLab