From fb508bc554fb23d1f5d8846dcc7f1940ffa927a8 Mon Sep 17 00:00:00 2001
From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr>
Date: Thu, 24 May 2018 17:17:50 +0200
Subject: [PATCH] =?UTF-8?q?mise=20=C3=A0=20jour?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/main.cpp                              | 29 ++++++-----
 src/mesh/Mesh.hpp                         |  4 +-
 src/scheme/AcousticSolver.hpp             | 10 ++++
 src/scheme/FiniteVolumesDiffusion.hpp     | 62 +++++++++++++----------
 src/scheme/FiniteVolumesEulerUnknowns.hpp |  5 +-
 5 files changed, 65 insertions(+), 45 deletions(-)

diff --git a/src/main.cpp b/src/main.cpp
index b57d75e46..7a116ff15 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -156,9 +156,6 @@ int main(int argc, char *argv[])
     double c = 0.;
     c = finite_volumes_diffusion.conservatif(unknowns);
 
-    std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset
-	      << ":  " << rang::fgB::green << c << rang::fg::reset << " \n";
-
     while((t<tmax) and (iteration<itermax)) {
      
       // ETAPE 1 DU SPLITTING - EULER
@@ -172,7 +169,7 @@ int main(int argc, char *argv[])
       t += dt_euler;
       
       // ETAPE 2 DU SPLITTING - DIFFUSION
-      /*
+      
       double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj);
  
       if (dt_euler <= dt_diff) {
@@ -190,7 +187,7 @@ int main(int argc, char *argv[])
 	}
 	std::cout << "diff : " << t_diff << std::endl;
       }
-      */
+      
 
       // DIFFUSION PURE
       
@@ -215,25 +212,28 @@ int main(int argc, char *argv[])
 
     
     double error = 0.;
-    error = finite_volumes_diffusion.error_L2_u(unknowns);
+    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 error2 = 0.;
-    error2 = finite_volumes_diffusion.error_Linf(unknowns);
+    error2 = finite_volumes_diffusion.error_Linf(unknowns, tmax);
 
-    std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset
+    std::cout << "* " << rang::style::underline << "Erreur L infini rho" << 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";
     */
+
+    std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset
+	      << ":  " << rang::fgB::green << c << rang::fg::reset << " \n";
     
     double cons = 0.;
     cons = finite_volumes_diffusion.conservatif(unknowns);
@@ -247,15 +247,16 @@ int main(int argc, char *argv[])
     
     { // gnuplot output for density
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
-     const Kokkos::View<const Rd*> uj = unknowns.uj();
+     const Kokkos::View<const double*> rhoj = unknowns.rhoj();
      double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
      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((4.*((xj[j][0]*xj[j][0])/(h*h)) + 100.-(xj[j][0]*xj[j][0])/(h*h))/100.)/h << '\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();
@@ -275,13 +276,13 @@ int main(int argc, char *argv[])
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const double*> Ej = unknowns.Ej();
      //double pi = 4.*std::atan(1.);
-     double h = std::sqrt(1. - (0.2*0.2)/(50./9.));
+     double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
      std::ofstream fout("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] << ' ' << ((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. <<'\n' ; // cas k non constant
-       fout << xj[j][0] << ' ' << Ej[j] << ' ' << (std::sqrt((4.*((xj[j][0]*xj[j][0])/(h*h)) + 100.-(xj[j][0]*xj[j][0])/(h*h))/100.)/h)*(std::sqrt((4.*((xj[j][0]*xj[j][0])/(h*h)) + 100.-(xj[j][0]*xj[j][0])/(h*h))/100.)/h) + (-(xj[j][0]*0.2)/((50./9.)-0.2*0.2))*(-(xj[j][0]*0.2)/((50./9.)-0.2*0.2))*0.5 << '\n';
+       fout << xj[j][0] << ' ' << Ej[j] << ' ' << (std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h) + (-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*(-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*0.5 << '\n';
      }
      }
 
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index d73852a87..c19db1412 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -55,9 +55,9 @@ public:
     : m_connectivity(connectivity),
       m_xr("xr", connectivity.numberOfNodes())
   {
-    const double delta_x = 1./connectivity.numberOfCells();
+    const double delta_x = (1.-0.1)/connectivity.numberOfCells();
     Kokkos::parallel_for(connectivity.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
-  	m_xr[r][0] = r*delta_x;
+  	m_xr[r][0] = 0.1+r*delta_x;
       });
   }
   
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index c22fca6c5..9abc8b463 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -349,6 +349,16 @@ 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;
+    
+    uj[0] = -(t/((50./9.)-t*t))*xj[0];
+    uj[m_mesh.numberOfCells()-1] = -(t/((50./9.)-t*t))*xj[m_mesh.numberOfCells()-1];
+    
+    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;
+    
   }
 };
 
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index ba3491a69..0b49be7ef 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -123,8 +123,10 @@ 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)));
     */
-    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];
+    //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];
+    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;
     
     return m_Fl ;
   }
@@ -164,8 +166,8 @@ 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)*xr(0);
-    m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*Fl(m_mesh.numberOfFaces()-1)*xr(m_mesh.numberOfFaces()-1);
+    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);
     
     return m_Gl ;
 
@@ -287,12 +289,6 @@ public:
     Kokkos::View<double*> kR = unknowns.kR();
 
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
-    // Premiere possibilite CL (le mieux a mon avis et le plus logique)
-    //uR[0] = (-t/((50./9.)-t*t))*xj[m_mesh.numberOfCells()-1]; 
-
-    // Deuxieme possiblitie CL
-    Kokkos::View<const Rd*> xr = m_mesh.xr();
-    //uR[0] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1];
     
     const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
     const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
@@ -326,31 +322,40 @@ public:
 	//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 = 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)));
+	//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)));
 
       });
-
+    
+    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]);
       });
 
     // met a jour les quantites (geometriques) associees au maillage
-    m_mesh_data.updateAllData();
+    //m_mesh_data.updateAllData();
 
+    /*
     // Calcul de rho avec la formule Mj = Vj rhoj
     const Kokkos::View<const double*> mj = unknowns.mj();
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
 	rhoj[j] = mj[j]/Vj[j];
-
+	
       });
+    */
   }
 
   // Calcul erreur entre solution analytique et solution numerique en norme L2
   // (quand la solution exacte est connue)
 
-  double error_L2_u(UnknownsType& unknowns) {
+  double error_L2_u(UnknownsType& unknowns, const double& t) {
 
     Kokkos::View<Rd*> uj = unknowns.uj();
 
@@ -364,13 +369,14 @@ public:
     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]*0.2)/((50./9.)-0.2*0.2);
+      exact_u = -(xj[j][0]*t)/((50./9.)-t*t);
       err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
     }
     err_u = std::sqrt(err_u);
     return err_u;
   }
 
+  /*
   double error_L2_E(UnknownsType& unknowns) {
 
     Kokkos::View<double*> Ej = unknowns.Ej();
@@ -390,29 +396,31 @@ public:
     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)
 
-  double error_Linf(UnknownsType& unknowns) {
+  double error_Linf(UnknownsType& unknowns, const double& t) {
 
-    Kokkos::View<Rd*> uj = unknowns.uj();
+    Kokkos::View<double*> rhoj = unknowns.rhoj();
 
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
 
-    double pi = 4.*std::atan(1.);
+    //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 erreur = std::abs(exacte - uj[0][0]);
+    //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); 
-
-      if (std::abs(exacte - uj[j][0]) > erreur) {
-	erreur = std::abs(exacte - uj[j][0]);
+      //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]);
       }
-
     }
 
     return erreur;
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 52e66cbcc..8804ce335 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -197,7 +197,7 @@ public:
   	  m_rhoj[j]=0.125;
   	}
 	*/
-	m_rhoj[j] = std::sqrt((4.*(xj[j][0]*xj[j][0]) + 100.-xj[j][0]*xj[j][0])/100.);
+	m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.);
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -236,7 +236,8 @@ public:
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-  	m_kj[j] =  xj[j][0]; 
+  	//m_kj[j] =  xj[j][0];
+	m_kj[j] = 0.5;
       });
 
      // Conditions aux bords de Dirichlet sur u et k
-- 
GitLab