ajout rho pour splot

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.2;
+    const double tmax=0.8;
     double t=0.;
 
     int itermax=std::numeric_limits<int>::max();
@@ -169,6 +169,7 @@ int main(int argc, char *argv[])
     const Kokkos::View<const Rd*> xj = mesh_data.xj();
     int size = 3000;
     std::vector<std::vector<double>> x(size, std::vector<double>(mesh.numberOfCells()));
+    std::vector<std::vector<double>> rho_marche(size, std::vector<double>(mesh.numberOfCells()));
     std::vector<double> tempo(size);
     int i = 0;
 
@@ -274,7 +275,7 @@ int main(int argc, char *argv[])
 
     while((t<tmax) and (iteration<itermax)) {
      
-      /*
+      
       // ETAPE 1 DU SPLITTING - EULER
       
       double dt_euler = 0.9*acoustic_solver.acoustic_dt(Vj, cj);
@@ -287,7 +288,7 @@ int main(int argc, char *argv[])
       
       // ETAPE 2 DU SPLITTING - DIFFUSION
       
-      double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj);
+      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) {
@@ -295,7 +296,7 @@ int main(int argc, char *argv[])
 	finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns);
       } else {
 	while (t > t_diff) {
-	  dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, nuj,cj);
+	  dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, nuj,cj,nuL, nuR, kL,kR);
 	  if (t_diff+dt_diff > t) {
 	    dt_diff = t-t_diff;
 	  }
@@ -303,9 +304,9 @@ int main(int argc, char *argv[])
 	  t_diff += dt_diff;
 	}
       }
-      */
       
 
+      /*
       // DIFFUSION PURE
       
       double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj,kj,nuj,cj,nuL,nuR,kL,kR);
@@ -314,7 +315,7 @@ int main(int argc, char *argv[])
       }
       finite_volumes_diffusion.computeNextStep(t, dt, unknowns);
       t += dt;
-      
+      */
       
       block_eos.updatePandCFromRhoE();  
     
@@ -545,14 +546,15 @@ int main(int argc, char *argv[])
       // ENTROPY TEST
       //finite_volumes_diffusion.entropie(unknowns);
 
-      /*
+      
       // 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];
       }
       tempo[i] = t;
       i = i + 1;
-      */
+      
 
     }
 
@@ -572,21 +574,21 @@ int main(int argc, char *argv[])
       fout << ' ' << '\n';
     }
     */
-    /*
-    std::ofstream fout("cara1");
+    
+    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) {
 	  if ( (k%10 == 0) or (k == i-1) ) {
-	    fout << tempo[k] << ' ' << x[k][j] << '\n';
+	    fout << tempo[k] << ' ' << x[k][j] << ' ' << rho_marche[k][j] << '\n';
 	  }
 	}
 	fout << ' ' << '\n';
 	fout << ' ' << '\n';
       }
     }
-   */ 
+   
     
     /*
     double error1 = 0.;

src/mesh/Mesh.hpp

@@ -76,7 +76,7 @@ public:
 
   // pas constant
 
-  /*
+  
   Mesh(const Connectivity& connectivity)
     : m_connectivity(connectivity),
       m_xr("xr", connectivity.numberOfNodes()),
@@ -92,7 +92,7 @@ public:
   	m_xr[r][0] = a + r*delta_x;
       });
   }
-  */
+  
 
   // pas non constant
   
@@ -128,7 +128,7 @@ public:
 
   // pas aleatoire
 
-  
+  /*
   Mesh(const Connectivity& connectivity)
     : m_connectivity(connectivity),
       m_xr("xr", connectivity.numberOfNodes()),
@@ -160,7 +160,7 @@ public:
 
     std::exit(0);
   }
-  
+  */
 
 
 

src/scheme/FiniteVolumesEulerUnknowns.hpp

@@ -276,7 +276,7 @@ public:
   }
   
 
-  /*
+  
   // --- Acoustic Solver ---
   
   void initializeSod()
@@ -338,7 +338,7 @@ public:
 
 	//m_kj[j] = 0.;
 
-	
+	/*
 	// Sod
 
 	// k non regulier
@@ -373,12 +373,12 @@ public:
 	    m_kj[j]=0. ;
 	  }
 	}
-	
+	*/
 	
 	// 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];
+	m_kj[j] = 0.14*m_kj[j];
 	
 	
       });
@@ -397,8 +397,8 @@ public:
       });
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
 	//m_nuj(j) = 0.5*(1.+xj[j][0]); 
-	m_nuj(j) = 0.0005;
-	
+	m_nuj(j) = 0.;
+	/*
 	if (xj[j][0]<0.7) {
   	  m_nuj[j]=0.;
   	} else {
@@ -408,7 +408,7 @@ public:
 	    m_nuj[j]=0. ;
 	  }
 	}
-	
+	*/
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -419,14 +419,14 @@ public:
     
     m_TL[0] = m_ej[0];
     m_TR[0] = m_ej[m_mesh.numberOfCells()-1];
-    m_nuL[0] = 0.0005;
-    m_nuR[0] = 0.0005;
+    m_nuL[0] = 0.;
+    m_nuR[0] = 0.;
    
   }
-  */
   
   
  
+  /*
   // DIFFUSION PURE
 
     void initializeSod()
@@ -499,7 +499,7 @@ public:
     m_nuR[0] = 0.;
   
   }
- 
+  */
   
 
   FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data)