From 887f8151fb5b742474069e25273661b7c04e67f9 Mon Sep 17 00:00:00 2001
From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr>
Date: Tue, 15 May 2018 09:16:41 +0200
Subject: [PATCH] ajout fonction calcul erreur L2 de E

---
 src/main.cpp                          | 16 +++++--
 src/scheme/FiniteVolumesDiffusion.hpp | 61 +++++++++++++++------------
 2 files changed, 47 insertions(+), 30 deletions(-)

diff --git a/src/main.cpp b/src/main.cpp
index 8dc0a24db..15c973fff 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -178,16 +178,24 @@ int main(int argc, char *argv[])
 	      << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
     double error = 0.;
-    error = finite_volumes_diffusion.error_L2(unknowns);
+    error = finite_volumes_diffusion.error_L2_u(unknowns);
 
-    std::cout << "* " << rang::style::underline << "Erreur L2" << rang::style::reset
+    std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset
 	      << ":  " << rang::fgB::green << error << rang::fg::reset << " \n";
 
+    /*
     double error2 = 0.;
     error2 = finite_volumes_diffusion.error_Linf(unknowns);
 
-    std::cout << "* " << rang::style::underline << "Erreur L infini" << rang::style::reset
+    std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset
 	      << ":  " << rang::fgB::green << error2 << rang::fg::reset << " \n";
+    */
+
+    double error3 = 0.;
+    error3 = finite_volumes_diffusion.error_L2_E(unknowns);
+
+    std::cout << "* " << rang::style::underline << "Erreur L2 E" << rang::style::reset
+	      << ":  " << rang::fgB::green << error3 << rang::fg::reset << " \n";
 
     double cons = 0.;
     cons = finite_volumes_diffusion.conservatif(unknowns);
@@ -219,7 +227,7 @@ int main(int argc, char *argv[])
      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] << ' ' << (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
+       fout << xj[j][0] << ' ' << Ej[j] << ' ' << (0.5*pi*pi*xj[j][0]*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(pi*xj[j][0])*std::cos(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
      }
      }
 
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 3f3fc0d24..9a1f422b8 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -107,6 +107,7 @@ private:
 	    //sum3 += Vj(cell_here);
 	    sum2 += (1./Vj(cell_here))*kj(cell_here);
 	    sum3 += 1./Vj(cell_here);
+	    // les deux moyennes donnent les memes resultats
 	}
 
 	m_Fl(l) = ((sum2/sum3)/Vl(l))*sum;
@@ -115,19 +116,6 @@ private:
 
     // Conditions aux bords
 
-    /* // 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);
-    */ 
-
     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))*(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)));
@@ -317,8 +305,9 @@ 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
        	Ej[j] +=  (dt*inv_mj[j]) * energy_fluxes;
 
       });
@@ -342,7 +331,7 @@ public:
   // Calcul erreur entre solution analytique et solution numerique en norme L2
   // (quand la solution exacte est connue)
 
-  double error_L2(UnknownsType& unknowns) {
+  double error_L2_u(UnknownsType& unknowns) {
 
     Kokkos::View<Rd*> uj = unknowns.uj();
 
@@ -351,15 +340,35 @@ public:
     const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
 
     double pi = 4.*std::atan(1.);
-    double erreur = 0.;
-    double exacte = 0.;
+    double err_u = 0.;
+    double exact_u = 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
-      erreur += (exacte - uj[j][0])*(exacte - uj[j][0])*Vj(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
+      err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
     }
-    erreur = std::sqrt(erreur);
-    return erreur;
+    err_u = std::sqrt(err_u);
+    return err_u;
+  }
+
+  double error_L2_E(UnknownsType& unknowns) {
+
+    Kokkos::View<double*> Ej = unknowns.Ej();
+
+    const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
+
+    const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
+
+    double pi = 4.*std::atan(1.);
+    double err_E = 0.;
+    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 = (0.5*pi*pi*xj[j][0]*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(pi*xj[j][0])*std::cos(pi*xj[j][0])) - cos(xj[j][0])*sin(xj[j][0])*pi*0.5)*(std::exp(-2.*0.2) - 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)
@@ -372,13 +381,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]);
-- 
GitLab