corrections bords, et tests sur flux Gl et Fl

src/main.cpp

@@ -153,7 +153,7 @@ int main(int argc, char *argv[])
 
     BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
-    while((t<tmax) and (iteration<itermax)) {
+    while((t<tmax) and (iteration<1)) {
       //double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
       double dt = 0.8*finite_volumes_diffusion.diffusion_dt(rhoj, kj);
 
@@ -204,6 +204,7 @@ int main(int argc, char *argv[])
      const Kokkos::View<const Rd*> uj = unknowns.uj();
      double pi = 4.*std::atan(1.);
      std::ofstream fout("comparaison 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
@@ -215,6 +216,7 @@ int main(int argc, char *argv[])
      const Kokkos::View<const double*> Ej = unknowns.Ej();
      double pi = 4.*std::atan(1.);
      std::ofstream fout("comparaison 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
      }

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
 
-  /*      
+          
   Mesh(const Connectivity& connectivity)
   : m_connectivity(connectivity),
     m_xr("xr", connectivity.numberOfNodes())
@@ -85,7 +85,7 @@ public:
   	}
       });
   }
-  */
+  
 
   ~Mesh()
   {

src/scheme/FiniteVolumesDiffusion.hpp

@@ -93,6 +93,7 @@ private:
 
     const Kokkos::View<const double*>& Vl = m_mesh_data.Vl();
     const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
+    const Kokkos::View<const Rd*> xr = m_mesh.xr();
 
     Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
         Rdd sum = zero;
@@ -112,19 +113,32 @@ private:
       });
 
     // Conditions aux bords
-    /*
-    // Essai 1
+
+    /* // EN TRAVAUX
+    Rdd sum = zero;
+    double sum2 = 0.;
+    double sum3 = 0.;
+    for (int j=0; j<face_nb_cells(0); ++j) {
+      int cell_here = face_cells(l,j);
+      int local_face_number_in_cell = face_cell_local_face(l,j);
+      if (cell_here >= 0) {
+	sum -= tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell));
+	    sum2 += kj(cell_here)*Vj(cell_here);
+	    sum3 += Vj(cell_here);
+    */ 
+    // Essai 1 OK
     int cell_here = face_cells(0,0);
     int local_face_number_in_cell = face_cell_local_face(0,0);
-    m_Fl(0) = -(kL(0) + kj(cell_here))*0.5*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell)));
+    m_Fl(0) = -(kL(0) + kj(cell_here))*(1./(2*Vl(0)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell)));
     cell_here = face_cells(m_mesh.numberOfFaces()-1,0);
     local_face_number_in_cell = face_cell_local_face(m_mesh.numberOfFaces()-1,0);
-    m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(cell_here))*0.5*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
+    m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(cell_here))*(1/(2.*Vl(m_mesh.numberOfFaces()-1)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
 
-    // Essai 2
-    m_Fl(0) = -(kL(0) + kj(0))*0.5*(tensorProduct(uj(0), Cjr(0, 0)));
-    m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(m_mesh.numberOfCells()-1))*0.5*(tensorProduct(uj(m_mesh.numberOfCells()-1), Cjr(m_mesh.numberOfCells()-1,1))); 
-    */	
+    double pi = 4.*std::atan(1.);;
+    std::ofstream fout("comparaison flux u");
+    for (int l=0; l<m_mesh.numberOfFaces(); ++l) {
+      fout << xr[l][0]  << ' ' << m_Fl[l](0,0) << ' ' << 2.*pi*cos(pi*xr[l][0])*std::exp(-2*pi*pi*0.000015625) << std::endl;
+    }
 
     return m_Fl ;
   }
@@ -146,6 +160,8 @@ private:
 
     const Kokkos::View<const double*>& Vj = m_mesh_data.Vj();
 
+    const Kokkos::View<const Rd*> xr = m_mesh.xr();
+
     Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
         Rd sum = zero;
 	double sum2 = 0.;
@@ -163,6 +179,13 @@ private:
       m_Gl(0) = Fl(0)*uL(0);
       m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
 
+      // affichage
+      double pi = 4.*std::atan(1.);
+      std::ofstream fout("comparaison flux E");
+      for (int l=0; l<m_mesh.numberOfFaces(); ++l) {
+	fout << xr[l][0] << ' ' << m_Gl[l][0] << ' ' << sin(pi*xr[l][0])*2.*pi*cos(pi*xr[l][0])*std::exp(-4*pi*pi*0.000015625) << std::endl;
+      }
+      
     return m_Gl ;
   }
 
@@ -288,9 +311,6 @@ public:
     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);