nettoyage code et presentation plus claire

src/main.cpp

@@ -153,8 +153,8 @@ int main(int argc, char *argv[])
 
     BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
-    double c = 0.;
-    c = finite_volumes_diffusion.conservatif(unknowns);
+    //double c = 0.;
+    //c = finite_volumes_diffusion.conservatif(unknowns);
 
     while((t<tmax) and (iteration<itermax)) {
      
@@ -210,13 +210,11 @@ int main(int argc, char *argv[])
     std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
 	      << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
+    double error1 = 0.;
+    error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax);
 
-    double error = 0.;
-    error = finite_volumes_diffusion.error_L2_u(unknowns, tmax);
-
-    std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset
-	      << ":  " << rang::fgB::green << error << rang::fg::reset << " \n";
-
+    std::cout << "* " << rang::style::underline << "Erreur L2 rho" << rang::style::reset
+	      << ":  " << rang::fgB::green << error1 << rang::fg::reset << " \n";
     
     double error2 = 0.;
     error2 = finite_volumes_diffusion.error_Linf(unknowns, tmax);
@@ -224,6 +222,12 @@ int main(int argc, char *argv[])
     std::cout << "* " << rang::style::underline << "Erreur L infini rho" << rang::style::reset
 	      << ":  " << rang::fgB::green << error2 << rang::fg::reset << " \n";
 
+    double error = 0.;
+    error = finite_volumes_diffusion.error_L2_u(unknowns, tmax);
+
+    std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset
+	      << ":  " << rang::fgB::green << error << rang::fg::reset << " \n";
+    
     /*
     double error3 = 0.;
     error3 = finite_volumes_diffusion.error_L2_E(unknowns);
@@ -231,7 +235,7 @@ int main(int argc, char *argv[])
     std::cout << "* " << rang::style::underline << "Erreur L2 E" << rang::style::reset
 	      << ":  " << rang::fgB::green << error3 << rang::fg::reset << " \n";
     */
-
+    /*
     std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset
 	      << ":  " << rang::fgB::green << c << rang::fg::reset << " \n";
     
@@ -240,7 +244,7 @@ int main(int argc, char *argv[])
 
     std::cout << "* " << rang::style::underline << "Resultat conservativite E" << rang::style::reset
 	      << ":  " << rang::fgB::green << cons << rang::fg::reset << " \n";
-    
+    */
 
     //method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
     method_cost_map["FiniteVolumesDiffusionWithMesh"] = timer.seconds();
@@ -252,12 +256,10 @@ int main(int argc, char *argv[])
      std::ofstream fout("rho");
      fout.precision(15);
      for (size_t j=0; j<mesh.numberOfCells(); ++j) {
-       //fout << xj[j][0] << ' ' << rhoj[j] << '\n';
        fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n';
      }
      }
 
-    
      { // gnuplot output for vitesse
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const Rd*> uj = unknowns.uj();

src/scheme/AcousticSolver.hpp

@@ -349,6 +349,8 @@ public:
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
 	rhoj[j] = mj[j]/Vj[j];
       });
+
+    /*
     double h = std::sqrt(1. - (t*t)/(50./9.));
     rhoj[0] = std::sqrt((3.*((xj[0][0]*xj[0][0])/(h*h)) + 100.)/100.)/h;
     rhoj[m_mesh.numberOfCells()-1] = std::sqrt((3.*((xj[m_mesh.numberOfCells()-1][0]*xj[m_mesh.numberOfCells()-1][0])/(h*h)) + 100.)/100.)/h;
@@ -358,7 +360,7 @@ public:
     
     Ej[0] = (std::sqrt((3.*((xj[0][0]*xj[0][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[0][0]*xj[0][0])/(h*h)) + 100.)/100.)/h) + (-(xj[0][0]*t)/((50./9.)-t*t))*(-(xj[0][0]*t)/((50./9.)-t*t))*0.5;
     Ej[m_mesh.numberOfCells()-1] = (std::sqrt((3.*((xj[m_mesh.numberOfCells()-1][0]*xj[m_mesh.numberOfCells()-1][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[m_mesh.numberOfCells()-1][0]*xj[m_mesh.numberOfCells()-1][0])/(h*h)) + 100.)/100.)/h) + (-(xj[m_mesh.numberOfCells()-1][0]*t)/((50./9.)-t*t))*(-(xj[m_mesh.numberOfCells()-1][0]*t)/((50./9.)-t*t))*0.5;
-    
+    */
   }
 };
 

src/scheme/FiniteVolumesDiffusion.hpp

@@ -123,8 +123,12 @@ private:
     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)));
     //m_Fl(m_mesh.numberOfFaces()-1) = -xr[m_mesh.numberOfNodes()-1][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell)));
     */
+
+    // k = x
     //m_Fl(0,0) = -(t/((50./9.)-t*t))*xr[0][0];
     //m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*xr[m_mesh.numberOfFaces()-1][0];
+    
+    // k = 0.5
     m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5;
     m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*0.5;
     
@@ -166,6 +170,7 @@ private:
     // Conditions aux bords
     // m_Gl(0) = Fl(0)*uL(0);
     //m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
+
     m_Gl(0) = -(t/((50./9.)-t*t))*Fl(0,0)*xr(0);
     m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*Fl(m_mesh.numberOfFaces()-1,0)*xr(m_mesh.numberOfFaces()-1);
     
@@ -321,19 +326,15 @@ public:
 	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
 
+	// 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)));
+
 	// ajout second membre pour kidder (k = x)
 	//uj[j][0] += (dt*inv_mj[j])*Vj(j)*(t/((50./9.)-t*t)); 
-	Ej[j] -= (dt*inv_mj[j])*Vj(j)*((0.5*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
+	//Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
 
       });
 
-    uj[0] = -(t/((50./9.)-t*t))*xj[0];
-    uj[m_mesh.numberOfCells()-1] = -(t/((50./9.)-t*t))*xj[m_mesh.numberOfCells()-1];
-   
-    double h = std::sqrt(1. - (t*t)/(50./9.));
-    Ej[0] = (std::sqrt((3.*((xj[0][0]*xj[0][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[0][0]*xj[0][0])/(h*h)) + 100.)/100.)/h) + (-(xj[0][0]*t)/((50./9.)-t*t))*(-(xj[0][0]*t)/((50./9.)-t*t))*0.5;
-    Ej[m_mesh.numberOfCells()-1] = (std::sqrt((3.*((xj[m_mesh.numberOfCells()-1][0]*xj[m_mesh.numberOfCells()-1][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[m_mesh.numberOfCells()-1][0]*xj[m_mesh.numberOfCells()-1][0])/(h*h)) + 100.)/100.)/h) + (-(xj[m_mesh.numberOfCells()-1][0]*t)/((50./9.)-t*t))*(-(xj[m_mesh.numberOfCells()-1][0]*t)/((50./9.)-t*t))*0.5;
-    
     // Calcul de e par la formule e = E-0.5 u^2 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
 	ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
@@ -376,6 +377,25 @@ public:
     return err_u;
   }
 
+  double error_L2_rho(UnknownsType& unknowns, const double& t) {
+
+    Kokkos::View<double*> rhoj = unknowns.rhoj();
+
+    const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+
+    const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
+
+    double h = std::sqrt(1. - (t*t)/(50./9.));
+    double err_rho = 0.;
+    double exact_rho = 0.;
+    for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
+      exact_rho = std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h;
+      err_rho += (exact_rho - rhoj[j])*(exact_rho - rhoj[j])*Vj(j);
+    }
+    err_rho = std::sqrt(err_rho);
+    return err_rho;
+  }
+
   /*
   double error_L2_E(UnknownsType& unknowns) {
 
@@ -409,14 +429,10 @@ public:
 
     //double pi = 4.*std::atan(1.);
     double h = std::sqrt(1. - (t*t)/(50./9.));
-    //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::sqrt((4.*((xj[0][0]*xj[0][0])/(h*h)) + 100.-(xj[0][0]*xj[0][0])/(h*h))/100.)/h;
     double erreur = std::abs(exacte - rhoj[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::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h;
       if (std::abs(exacte - rhoj[j]) > erreur) {
 	erreur = std::abs(exacte - rhoj[j]);