tentative d'injection de la solution exacte de kidder pour euler

src/main.cpp

@@ -176,7 +176,7 @@ int main(int argc, char *argv[])
       }
       acoustic_solver.computeNextStep(t,dt_euler, unknowns);
       t += dt_euler;
-      
+      /*
       // ETAPE 2 DU SPLITTING - DIFFUSION
       
       double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj);
@@ -196,7 +196,7 @@ int main(int argc, char *argv[])
 	}
 	std::cout << "diff : " << t_diff << std::endl;
       }
-      
+      */
 
       // DIFFUSION PURE
       
@@ -260,10 +260,12 @@ 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./6.)*(50./6.)));
      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] << '\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';
      }
      }
     
@@ -276,20 +278,23 @@ int main(int argc, char *argv[])
      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] << '\n';
+       fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*0.2)/(1.-0.2*0.2) << '\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.);
-     std::ofstream fout("E");
+     double h = std::sqrt(1. - (0.2*0.2)/((50./6.)*(50./6.)));
+     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] << '\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';
      }
      }
 

src/scheme/FiniteVolumesEulerUnknowns.hpp

@@ -191,19 +191,23 @@ public:
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-  	if (xj[j][0]<0.5) {
+  	/*if (xj[j][0]<0.5) {
   	  m_rhoj[j]=1;
   	} else {
   	  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.);
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-  	if (xj[j][0]<0.5) {
+  	/*if (xj[j][0]<0.5) {
   	  m_pj[j]=1;
   	} else {
   	  m_pj[j]=0.1;
   	}
+	*/
+	m_pj[j] = 2.*std::pow(m_rhoj[j],3);
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
@@ -211,7 +215,8 @@ public:
       });
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
-  	m_gammaj[j] = 1.4;
+  	//m_gammaj[j] = 1.4;
+	m_gammaj[j] = 3.;
       });
 
     BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);