ajout solution E pour k non constant

src/main.cpp

@@ -206,8 +206,8 @@ int main(int argc, char *argv[])
      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
+       // 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
      }
      }
 
@@ -218,7 +218,8 @@ int main(int argc, char *argv[])
      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
+       //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] << ' ' << (0.5*pi*pi*xj[j][0]*(-std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0]) + std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0])) - cos(xj[j][0])*sin(xj[j][0])*pi*0.5)*(std::exp(-2.*0.2) - 1.) + 2. <<'\n' ; // cas k non 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,10 +60,10 @@ public:
   	m_xr[r][0] = r*delta_x;
       });
   }
-*/
 
-  // pas non constant
 
+  // pas non constant
+  /*
           
   Mesh(const Connectivity& connectivity)
   : m_connectivity(connectivity),
@@ -85,7 +85,7 @@ public:
   	}
       });
   }
-  
+  */
 
   ~Mesh()
   {

src/scheme/FiniteVolumesDiffusion.hpp

@@ -315,9 +315,10 @@ public:
 	  momentum_fluxes +=  Fl(l)*Cjr(j,R);
 	  energy_fluxes   += (Gl(l), Cjr(j,R));
 	}
-	//uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*Sj(j) + (dt*inv_mj[j]) * momentum_fluxes; // test avec k non constant
-	uj[j] += (dt*inv_mj[j]) * momentum_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
+	//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 
@@ -351,8 +352,8 @@ public:
     double erreur = 0.;
     double exacte = 0.;
     for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
-      exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant
-      //exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant
+      //exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant
+      exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant
       erreur += (exacte - uj[j][0])*(exacte - uj[j][0])*Vj(j);
     }
     erreur = std::sqrt(erreur);
@@ -369,13 +370,13 @@ public:
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
 
     double pi = 4.*std::atan(1.);
-    double exacte = std::sin(pi*xj[0][0])*std::exp(-2.*pi*pi*0.2); // k constant
-    //double exacte = std::sin(pi*xj[0][0])*std::exp(-0.2); // k non constant
+    //double exacte = std::sin(pi*xj[0][0])*std::exp(-2.*pi*pi*0.2); // k constant
+    double exacte = std::sin(pi*xj[0][0])*std::exp(-0.2); // k non constant
     double erreur = std::abs(exacte - uj[0][0]);
 
     for (size_t j=1; j<m_mesh.numberOfCells(); ++j) {
-      exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2);
-      //exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); 
+      //exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2);
+      exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); 
 
       if (std::abs(exacte - uj[j][0]) > erreur) {
 	erreur = std::abs(exacte - uj[j][0]);

src/scheme/FiniteVolumesEulerUnknowns.hpp

@@ -184,7 +184,6 @@ public:
   }
 
 
-
   // --- Acoustic Solver ---
 
   // void initializeSod()
@@ -277,8 +276,8 @@ void initializeSod()
       });
 
     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){