correction T initial et ajout membre source pour kidder avec k non cst

src/main.cpp

@@ -275,7 +275,7 @@ int main(int argc, char *argv[])
       
       // ETAPE 2 DU SPLITTING - DIFFUSION
       
-      double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj);
+      double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj);
       double t_diff = t-dt_euler;
       
       if (dt_euler <= dt_diff) {
@@ -283,7 +283,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, cj);
+	  dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, nuj,cj);
 	  if (t_diff+dt_diff > t) {
 	    dt_diff = t-t_diff;
 	  }
@@ -567,6 +567,7 @@ 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);
@@ -579,14 +580,18 @@ int main(int argc, char *argv[])
 
     std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset
 	      << ":  " << rang::fgB::green << error4 << rang::fg::reset << " \n";
-    */
-    /*
+    
     double error3 = 0.;
-    error3 = finite_volumes_diffusion.error_L2_E(unknowns);
+    error3 = finite_volumes_diffusion.error_L2_E(unknowns,tmax);
 
     std::cout << "* " << rang::style::underline << "Erreur L2 E" << rang::style::reset
 	      << ":  " << rang::fgB::green << error3 << rang::fg::reset << " \n";
-    */
+    
+    double error5 = 0.;
+    error5 = finite_volumes_diffusion.error_Linf_E(unknowns, tmax);
+
+    std::cout << "* " << rang::style::underline << "Erreur L infini E" << rang::style::reset
+	      << ":  " << rang::fgB::green << error5 << 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";

src/scheme/FiniteVolumesDiffusion.hpp

@@ -441,7 +441,7 @@ public:
 
 	// 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)*((2.*xj[j][0]*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))-(6*xj[j][0]+3.)/(100*(1-t*t/(50/9))*(1-t*t/(50/9))));
 
       });
 
@@ -468,13 +468,13 @@ public:
 
     const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
 
-    //double pi = 4.*std::atan(1.);
+    double pi = 4.*std::atan(1.);
     double err_u = 0.;
     double exact_u = 0.;
     for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
       //exact_u = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant
-      // exact_u = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant
-      exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder
+      exact_u = std::sin(pi*xj[j][0])*std::exp(-t); // solution exacte cas test k non constant
+      //exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder
       err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
     }
     err_u = std::sqrt(err_u);
@@ -500,8 +500,8 @@ public:
     return err_rho;
   }
 
-  /*
-  double error_L2_E(UnknownsType& unknowns) {
+  
+  double error_L2_E(UnknownsType& unknowns, const double& t) {
 
     Kokkos::View<double*> Ej = unknowns.Ej();
 
@@ -514,13 +514,13 @@ public:
     double exact_E = 0.;
     for (size_t j=0; j<m_mesh.numberOfCells(); ++j) {
       //exact_E = (-(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.;
-      exact_E =  ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*0.2)-1.) + 2.;
+      exact_E =  ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*t)-1.) + 2.;
       err_E += (exact_E - Ej[j])*(exact_E - Ej[j])*Vj(j);
     }
     err_E = std::sqrt(err_E);
     return err_E;
   }
-  */
+  
   
   // Calcul erreur entre solution analytique et solution numerique en norme L infini (max)
   // (quand la solution exacte est connue)
@@ -555,11 +555,13 @@ public:
 
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
 
-    double exacte = -(xj[0][0]*t)/((50./9.)-t*t);
+    double pi = 4.*std::atan(1.);
+    
+    double exacte = std::sin(pi*xj[0][0])*std::exp(-t);
     double erreur = std::abs(exacte - uj[0][0]);
 
     for (size_t j=1; j<m_mesh.numberOfCells(); ++j) {
-      exacte = -(xj[j][0]*t)/((50./9.)-t*t);
+      exacte = std::sin(pi*xj[j][0])*std::exp(-t);
       if (std::abs(exacte - uj[j][0]) > erreur) {
 	erreur = std::abs(exacte - uj[j][0]);
       }
@@ -568,6 +570,29 @@ public:
     return erreur;
   }
 
+  // Calcul erreur entre solution analytique et solution numerique en norme L infini (max)
+  // (quand la solution exacte est connue)
+
+  double error_Linf_E(UnknownsType& unknowns, const double& t) {
+
+    Kokkos::View<double*> Ej = unknowns.Ej();
+
+    const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+    
+    double pi = 4.*std::atan(1.);
+    
+    double exacte = ((xj[0][0]*pi*pi*0.5)*(std::sin(pi*xj[0][0])*std::sin(pi*xj[0][0]) - std::cos(xj[0][0]*pi)*std::cos(pi*xj[0][0])) - pi*0.5*std::sin(pi*xj[0][0])*std::cos(pi*xj[0][0]))*(std::exp(-2.*t)-1.) + 2.;
+    double erreur = std::abs(exacte - Ej[0]);
+
+    for (size_t j=1; j<m_mesh.numberOfCells(); ++j) {
+      exacte = ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*t)-1.) + 2.;
+      if (std::abs(exacte - Ej[j]) > erreur) {
+	erreur = std::abs(exacte - Ej[j]);
+      }
+    }
+
+    return erreur;
+  }
   
   // Verifie la conservativite de E
 

src/scheme/FiniteVolumesEulerUnknowns.hpp

@@ -456,7 +456,7 @@ public:
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-	m_Tj[j] = 2. - 0.5*std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]); // k = x
+	m_Tj[j] = m_ej[j];
       });
 
     // Conditions aux bords de Dirichlet sur T et nu