mise a jour pour stephane

src/main.cpp

@@ -153,6 +153,12 @@ int main(int argc, char *argv[])
 
     BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
+    double c = 0.;
+    c = finite_volumes_diffusion.conservatif(unknowns);
+
+    std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset
+	      << ":  " << rang::fgB::green << c << rang::fg::reset << " \n";
+
     while((t<tmax) and (iteration<itermax)) {
      
       // ETAPE 1 DU SPLITTING - EULER
@@ -167,12 +173,11 @@ int main(int argc, char *argv[])
       
       // ETAPE 2 DU SPLITTING - DIFFUSION
       
-      double dt_diff = finite_volumes_diffusion.diffusion_dt(rhoj, kj);
+      double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj);
  
       if (dt_euler <= dt_diff) {
 	dt_diff = dt_euler;
 	finite_volumes_diffusion.computeNextStep(t, dt_diff, unknowns);
-	t += dt_diff;
       } else {
 	double t_diff = t;
 	while (t + dt_euler > t_diff) {
@@ -180,7 +185,6 @@ int main(int argc, char *argv[])
 	  dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj);
 	  t_diff += dt_diff;
 	}
-	t = t_diff;
       }
       
 
@@ -246,7 +250,7 @@ int main(int argc, char *argv[])
     { // gnuplot output for density
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const Rd*> uj = unknowns.uj();
-     std::ofstream fout("rho nounif");
+     std::ofstream fout("rho euler");
      fout.precision(15);
      for (size_t j=0; j<mesh.numberOfCells(); ++j) {
        fout << xj[j][0] << ' ' << rhoj[j] << '\n';
@@ -257,7 +261,7 @@ int main(int argc, char *argv[])
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const Rd*> uj = unknowns.uj();
      //double pi = 4.*std::atan(1.);
-     std::ofstream fout("u nounif");
+     std::ofstream fout("u euler");
      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
@@ -270,7 +274,7 @@ int main(int argc, char *argv[])
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const double*> Ej = unknowns.Ej();
      //double pi = 4.*std::atan(1.);
-     std::ofstream fout("E nounif");
+     std::ofstream fout("E euler");
      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

src/mesh/Mesh.hpp

@@ -50,7 +50,7 @@ public:
 
   // pas constant
 
-  /*
+  
   Mesh(const Connectivity& connectivity)
     : m_connectivity(connectivity),
       m_xr("xr", connectivity.numberOfNodes())
@@ -60,11 +60,11 @@ public:
   	m_xr[r][0] = r*delta_x;
       });
   }
-  */
   
-  // pas non constant
 
+  // pas non constant
   
+  /*     
   Mesh(const Connectivity& connectivity)
   : m_connectivity(connectivity),
     m_xr("xr", connectivity.numberOfNodes())
@@ -85,7 +85,7 @@ public:
   	}
       });
   }
-  
+  */
   ~Mesh()
   {
     ;

src/scheme/FiniteVolumesDiffusion.hpp

@@ -242,12 +242,15 @@ public:
 	}
 	
 	double sum = 0.;
+	
 	for (int m = 0; m < cell_nb_nodes(j); ++m) {
 	 sum += kj(cell_nodes(j,m));
 	}
 
 	dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl;
 	
+
+	
       });
     
     double dt = std::numeric_limits<double>::max();
@@ -400,6 +403,7 @@ public:
   double conservatif(UnknownsType& unknowns) {
    
     Kokkos::View<double*> Ej = unknowns.Ej();
+    Kokkos::View<double*> rhoj = unknowns.rhoj();
 
     const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
 
@@ -408,7 +412,7 @@ public:
     double sum = 0.;
 
     for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
-      sum += Ej[j]*Vj[j];
+      sum += Ej[j]*rhoj[j]*Vj[j];
     }
 
     return sum;

src/scheme/FiniteVolumesEulerUnknowns.hpp

@@ -218,8 +218,8 @@ 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]);
-	m_Ej[j] = 2.;
+  	m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
+	//m_Ej[j] = 2.;
       });
 
     const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
@@ -232,15 +232,15 @@ public:
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-  	m_kj[j] = 1.; 
+  	m_kj[j] =  0.00001*0.5; 
       });
 
      // Conditions aux bords de Dirichlet sur u et k
     
     m_uL[0] = zero;
     m_uR[0] = zero;
-    m_kL[0] = 1.;
-    m_kR[0] = 1.;
+    m_kL[0] = 0.00001*0.5;
+    m_kR[0] = 0.00001*0.5;
   }
   
   
@@ -287,8 +287,8 @@ 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
+	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){
@@ -299,8 +299,8 @@ public:
     
     m_uL[0] = zero;
     m_uR[0] = zero;
-    m_kL[0] = 0.;
-    m_kR[0] = 1.;
+    m_kL[0] = 2.;
+    m_kR[0] = 2.;
 
   }
 */