From 593f6bf0961daee8e50fbf79b9fa3669b2e34cbc Mon Sep 17 00:00:00 2001
From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr>
Date: Fri, 7 Sep 2018 10:53:52 +0200
Subject: [PATCH] test avec kidder, convergence ok

---
 src/main.cpp                              | 22 +++++++++---------
 src/mesh/Mesh.hpp                         |  6 ++---
 src/scheme/FiniteVolumesDiffusion.hpp     |  4 ++--
 src/scheme/FiniteVolumesEulerUnknowns.hpp | 24 ++++++++++----------
 src/scheme/NoSplitting.hpp                | 27 ++++++++++++-----------
 5 files changed, 42 insertions(+), 41 deletions(-)

diff --git a/src/main.cpp b/src/main.cpp
index b7f3bb494..29e54959c 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -144,7 +144,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.8;
+    const double tmax=1.5;
     double t=0.;
 
     int itermax=std::numeric_limits<int>::max();
@@ -594,7 +594,7 @@ int main(int argc, char *argv[])
     
 
     // Erreurs sous differents normes
-    /*
+    
 
     // Density
     double error1 = 0.;
@@ -612,7 +612,7 @@ int main(int argc, char *argv[])
     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 error4 = 0.;
     error4 = finite_volumes_diffusion.error_Linf_u(unknowns, tmax);
     std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset
@@ -658,18 +658,18 @@ int main(int argc, char *argv[])
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const double*> rhoj = unknowns.rhoj();
      double h = std::sqrt(1. - (tmax*tmax)/(50./9.));
-     std::ofstream fout("rho_no_split");
+     std::ofstream fout("rho");
      fout.precision(15);
      for (size_t j=0; j<mesh.numberOfCells(); ++j) {
-       //fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; // kidder
-       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'; // kidder
+       //fout << xj[j][0] << ' ' << rhoj[j] << '\n';
      }
      }
     
     { // gnuplot output for pression
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const double*> pj = unknowns.pj();
-     std::ofstream fout("p_no_split");
+     std::ofstream fout("p");
      fout.precision(15);
      for (size_t j=0; j<mesh.numberOfCells(); ++j) {
        fout << xj[j][0] << ' ' << pj[j] << '\n';
@@ -679,7 +679,7 @@ int main(int argc, char *argv[])
      { // gnuplot output for internal energy
      const Kokkos::View<const Rd*> xj   = mesh_data.xj();
      const Kokkos::View<const double*> ej = unknowns.ej();
-     std::ofstream fout("e_no_split");
+     std::ofstream fout("e");
      fout.precision(15);
      for (size_t j=0; j<mesh.numberOfCells(); ++j) {
        fout << xj[j][0] << ' ' << ej[j] << '\n';
@@ -690,16 +690,16 @@ 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("u_no_split");
+     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(-tmax) <<'\n'; // cas k non constant
-       //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder
+       fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder
        //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << xj[j][0] << std::endl;
 
-       fout << xj[j][0] << ' ' << uj[j][0] << '\n';
+       //fout << xj[j][0] << ' ' << uj[j][0] << '\n';
      }
      }
 
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 660632ef0..f39a069a0 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -82,8 +82,8 @@ public:
       m_x0("x0", 1),
       m_xmax("xmax", 1)
   {
-    double a = -1.;
-    double b = 2.;
+    double a = 0.;
+    double b = 1.;
     m_x0[0][0] = a;
     m_xmax[0][0] = b;
     const double delta_x = (b-a)/connectivity.numberOfCells();
@@ -160,7 +160,7 @@ public:
     //std::exit(0);
    
   }
- */
+  */
   
 
   ~Mesh()
diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp
index a25682075..fe5479516 100644
--- a/src/scheme/FiniteVolumesDiffusion.hpp
+++ b/src/scheme/FiniteVolumesDiffusion.hpp
@@ -474,8 +474,8 @@ public:
 
     // Mise a jour de la vitesse et de l'energie totale specifique
     const Kokkos::View<const double*> inv_mj = unknowns.invMj();
-    Kokkos::parallel_for(m_mesh.numberOfCells()-2, KOKKOS_LAMBDA(const int& j0) {
-	const int j = j0+1;
+    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
+	//const int j = j0+1;
 	Rd momentum_fluxes = zero;
 	double energy_fluxes = 0.;
 	Rd trich = zero;
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index b60ee816d..727f22897 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -299,28 +299,28 @@ public:
     
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
 	// TAC
-	
+	/*
 	if (xj[j][0]<0.5) {
   	  m_rhoj[j]=1.;
   	} else {
   	  m_rhoj[j]=0.125;
 	}
-	
+	*/
 	// Kidder
-	//m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/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){
 	// TAC
-	
+	/*
   	if (xj[j][0]<0.5) {
   	  m_pj[j]=1.;
   	} else {
   	  m_pj[j]=0.1;
 	}
-	
+	*/
 	// Kidder
-	//m_pj[j] = 2.*std::pow(m_rhoj[j],3);
+	m_pj[j] = 2.*std::pow(m_rhoj[j],3);
       });
     
     double pi = 4.*std::atan(1.);
@@ -330,9 +330,9 @@ public:
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
   	// TAC
-	m_gammaj[j] = 1.4;
+	//m_gammaj[j] = 1.4;
 	// Kidder
-	//m_gammaj[j] = 3.;
+	m_gammaj[j] = 3.;
       });
 
     BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
@@ -352,8 +352,8 @@ public:
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
 	// Differents k (xi)
-  	//m_kj[j] =  xj[j][0];
-	m_kj[j] = 0.014;
+  	m_kj[j] =  xj[j][0];
+	//m_kj[j] = 0.5;
 
 	// TAC
 	
@@ -408,8 +408,8 @@ public:
 
     m_uL[0] = zero;
     m_uR[0] = zero;
-    m_kL[0] = 0.014;
-    m_kR[0] = 0.014;
+    m_kL[0] = 0.;
+    m_kR[0] = 1.;
 
     
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){  	
diff --git a/src/scheme/NoSplitting.hpp b/src/scheme/NoSplitting.hpp
index 92ffe8ce9..3864d773c 100644
--- a/src/scheme/NoSplitting.hpp
+++ b/src/scheme/NoSplitting.hpp
@@ -190,19 +190,19 @@ private:
 	m_ur[r]=invAr(r)*br(r);
       });
 
-    
+    /*
     m_ur[0]=zero;
     m_ur[m_mesh.numberOfNodes()-1]=zero;
-    
+    */
     //m_ur[0] = x0;
     //m_ur[m_mesh.numberOfNodes()-1] = xmax[0];
     
     // CL Kidder
-    /*
+    
     double h = std::sqrt(1. - (t*t)/(50./9.));
     m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0];
     m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0];
-    */
+    
 
     return m_ur;
   }
@@ -350,8 +350,8 @@ public:
 	// ajout second membre pour kidder (k cst)
 	//Ej[j] -= (dt*inv_mj[j])*Vj(j)*((kj(j)*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)*((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)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); 
       });
 
     // Calcul de e par la formule e = E-0.5 u^2
@@ -384,26 +384,27 @@ public:
 	}
 
 	if (j == 0) {
-	  PTj(j) = pj(j) - kj(j)*(uj[j][0]-uL[0][0])/Vl(0);
-	  //PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
+	  //PTj(j) = pj(j) - kj(j)*(uj[j][0]-uL[0][0])/Vl(0);
+	  PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
 	} else if (j == m_mesh.numberOfCells()-1) {
-	  //PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
-	  PTj(j) = pj(j) - kj(j)*(uR[0][0]-uj[j][0])/Vl(m_mesh.numberOfFaces()-1);
+	  PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t));
+	  //PTj(j) = pj(j) - kj(j)*(uR[0][0]-uj[j][0])/Vl(m_mesh.numberOfFaces()-1);
 	} else {
-	  //PTj(j) = pj(j) - kj(j)*sum/Vl(j);
 	  PTj(j) = pj(j) - kj(j)*2.*sum/sum1;
 	}
 
       });
     
     // Mise a jour de k
-    /*
+    
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
 	kj(j) = xj[j][0];
     
       });
-    */
+    
   }
 };
 
+
+
 #endif // NO_SPLITTING_HPP
-- 
GitLab