From d41c3067da42058d1b7efea654722a71e04681df Mon Sep 17 00:00:00 2001
From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr>
Date: Wed, 23 May 2018 16:59:56 +0200
Subject: [PATCH] modifs sur les CL, pas terrible

---
 src/main.cpp                              | 33 +++++++--------------
 src/scheme/AcousticSolver.hpp             |  4 +--
 src/scheme/FiniteVolumesDiffusion.hpp     | 35 ++++++++++++++---------
 src/scheme/FiniteVolumesEulerUnknowns.hpp |  6 ++--
 4 files changed, 37 insertions(+), 41 deletions(-)

diff --git a/src/main.cpp b/src/main.cpp
index 7454ef1c2..3fb41c6c3 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -138,7 +138,7 @@ int main(int argc, char *argv[])
     const Kokkos::View<const double*> Vj = mesh_data.Vj();
     const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
 
-    const double tmax=0.2;
+    const double tmax=1.;
     double t=0;
 
     int itermax=std::numeric_limits<int>::max();
@@ -159,12 +159,6 @@ int main(int argc, char *argv[])
     std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset
 	      << ":  " << rang::fgB::green << c << rang::fg::reset << " \n";
 
-    double c1 = 0.;
-    c1 = finite_volumes_diffusion.conservatif_mvt(unknowns);
-
-    std::cout << "* " << rang::style::underline << "Resultat conservativite rho u temps = 0" << rang::style::reset
-	      << ":  " << rang::fgB::green << c1 << rang::fg::reset << " \n";
-
     while((t<tmax) and (iteration<itermax)) {
      
       // ETAPE 1 DU SPLITTING - EULER
@@ -219,14 +213,14 @@ 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 error = 0.;
     error = finite_volumes_diffusion.error_L2_u(unknowns);
 
     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);
 
@@ -246,12 +240,6 @@ 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";
-
-    double cons1 = 0.;
-    cons1 = finite_volumes_diffusion.conservatif_mvt(unknowns);
-
-    std::cout << "* " << rang::style::underline << "Resultat conservativite rho u" << rang::style::reset
-	      << ":  " << rang::fgB::green << cons1 << rang::fg::reset << " \n";
     
 
     //method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
@@ -260,8 +248,8 @@ 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();
-     double h = std::sqrt(1. - (0.2*0.2)/(50./9.));
-     std::ofstream fout("rho1");
+     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';
@@ -273,28 +261,27 @@ int main(int argc, char *argv[])
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const Rd*> uj = unknowns.uj();
      //double pi = 4.*std::atan(1.);
-     std::ofstream fout("u1");
+     std::ofstream fout("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] << '\n';
-       fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*0.2)/((50./9.)-0.2*0.2) << '\n';
+       fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n';
      }
      }
 
      { // gnuplot output for energy
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
-     const Kokkos::View<const double*> ej = unknowns.ej();
+     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.));
-     std::ofstream fout("e");
+     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] << '\n';
-       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)  << '\n';
+       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';
      }
      }
 
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index b0fd38286..4962964d3 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -162,7 +162,7 @@ private:
     return m_ur;
   }
   */
-
+  
   Kokkos::View<Rd*> // calcule u_r (vitesse au sommet du maillage et flux de vitesse)
   computeUr(const Kokkos::View<const Rdd*>& Ar,
 	    const Kokkos::View<const Rd*>& br,
@@ -181,7 +181,7 @@ private:
 
     return m_ur;
   }
-
+  
   Kokkos::View<Rd**>  // Fonction qui calcule F_jr
   computeFjr(const Kokkos::View<const Rdd**>& Ajr,
 	     const Kokkos::View<const Rd*>& ur,
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index 2c9701f06..58187015f 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -116,9 +116,11 @@ private:
     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)));
+    //m_Fl(0) = -xr[0][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell)));
     cell_here = face_cells(m_mesh.numberOfFaces()-1,0);
     local_face_number_in_cell = face_cell_local_face(m_mesh.numberOfFaces()-1,0);
     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)));
 
     return m_Fl ;
   }
@@ -128,7 +130,8 @@ private:
   computeGl(const Kokkos::View<const Rd*>& uj,
 	    const Kokkos::View<const Rdd*>& Fl,
 	    const Kokkos::View<const Rd*>& uL,
-	    const Kokkos::View<const Rd*>& uR) {
+	    const Kokkos::View<const Rd*>& uR,
+	    const double& t) {
 
     const Kokkos::View<const unsigned int**>& face_cells = m_connectivity.faceCells();
 
@@ -157,6 +160,9 @@ private:
     // Conditions aux bords
     m_Gl(0) = Fl(0)*uL(0);
     m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0);
+    //const double delta_x = 1./100.;
+    //double brico = - ((xr[m_mesh.numberOfNodes()-1][0]+delta_x)*t)/((50./9.)-t*t);
+    //m_Gl[m_mesh.numberOfFaces()-1][0] = brico*Fl[m_mesh.numberOfFaces()-1](0,0);
     
     return m_Gl ;
 
@@ -189,11 +195,12 @@ private:
 			     const Kokkos::View<const Rd*>& uL,
 			     const Kokkos::View<const Rd*>& uR,
 			     const Kokkos::View<const double*>& kL,
-			     const Kokkos::View<const double*>& kR) { 
+			     const Kokkos::View<const double*>& kR,
+			     const double& t) { 
     Kokkos::View<Rdd*> Fl  = m_Fl ; 
     Fl  = computeFl (Cjr, uj, kj, uL, uR, kL, kR);
     Kokkos::View<Rd*> Gl  = m_Gl ; 
-    Gl  = computeGl (uj, Fl, uL, uR);
+    Gl  = computeGl (uj, Fl, uL, uR, t);
   }
 
   Kokkos::View<Rdd*> m_Fl;
@@ -278,20 +285,19 @@ public:
     Kokkos::View<double*> kL = unknowns.kL();
     Kokkos::View<double*> kR = unknowns.kR();
 
-    Kokkos::View<const Rd*> xr = m_mesh.xr();
+    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 Rd*> xj = m_mesh_data.xj();
-
-    uR[0] = (-t/((50./9.)-t*t))*xj[m_mesh.numberOfCells()-1];
-
     const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
     const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
 
     // Calcule les flux
-    computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR);
+    computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, t);
 
     const Kokkos::View<const Rdd*> Fl  = m_Fl ;
     const Kokkos::View<const Rd *> Gl  = m_Gl ;
@@ -318,9 +324,9 @@ 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 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)));
+	Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t)));
 
       });
 
@@ -351,12 +357,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 = 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);
       err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j);
     }
     err_u = std::sqrt(err_u);
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 00c465f1e..9b1c5f7a2 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -241,10 +241,12 @@ public:
 
      // Conditions aux bords de Dirichlet sur u et k
     
+   
+
     m_uL[0] = zero;
     m_uR[0] = zero;
-    m_kL[0] = 0.;
-    m_kR[0] = 1.;
+    m_kL[0] = xj[0][0];
+    m_kR[0] = xj[m_mesh.numberOfCells()-1][0];
   }
   
   
-- 
GitLab