correction solution exacte gnuplot et modification CL dans computeUr (Acoustic)

src/main.cpp

@@ -265,7 +265,7 @@ int main(int argc, char *argv[])
      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((4.*((xj[j][0]*xj[j][0])/(h*h)) + 100.-(xj[j][0]*xj[j][0])/(h*h))/100.)/h << '\n';
      }
      }
     
@@ -279,7 +279,7 @@ int main(int argc, char *argv[])
        //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)/(1.-0.2*0.2) << '\n';
+       fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*0.2)/((50./9.)-0.2*0.2) << '\n';
      }
      }
 
@@ -294,7 +294,7 @@ int main(int argc, char *argv[])
        //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)  << '\n';
      }
      }
 

src/scheme/AcousticSolver.hpp

@@ -147,6 +147,7 @@ private:
     return m_br;
   }
 
+  /*
   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) {
@@ -160,6 +161,26 @@ 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,
+	    const double& t) {
+
+    inverse(Ar, m_inv_Ar);
+    const Kokkos::View<const Rdd*> invAr = m_inv_Ar;
+
+    Kokkos::View<Rd*> xr = m_mesh.xr();
+
+    Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
+	m_ur[r]=invAr(r)*br(r);
+      });
+    m_ur[0]=zero;
+    m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1];
+
+    return m_ur;
+  }
 
   Kokkos::View<Rd**>  // Fonction qui calcule F_jr
   computeFjr(const Kokkos::View<const Rdd**>& Ajr,
@@ -207,7 +228,8 @@ private:
 			     const Kokkos::View<const double*>& pj,
 			     const Kokkos::View<const double*>& cj,
 			     const Kokkos::View<const double*>& Vj,
-			     const Kokkos::View<const Rd**>& Cjr) {
+			     const Kokkos::View<const Rd**>& Cjr,
+			     const double& t) {
     const Kokkos::View<const double*> rhocj  = computeRhoCj(rhoj, cj);
     const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, Cjr);
 
@@ -216,7 +238,8 @@ private:
 
     Kokkos::View<Rd*> ur = m_ur;
     Kokkos::View<Rd**> Fjr = m_Fjr;
-    ur  = computeUr(Ar, br);
+    //ur  = computeUr(Ar, br);
+    ur = computeUr(Ar, br, t);
     Fjr = computeFjr(Ajr, ur, Cjr, uj, pj);
   }
 
@@ -282,7 +305,7 @@ public:
     Kokkos::View<Rd*> xr = m_mesh.xr();
 
     // Calcule les flux
-    computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr);
+    computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr, t);
 
     const Kokkos::View<const Rd**> Fjr = m_Fjr;
     const Kokkos::View<const Rd*> ur = m_ur;

src/scheme/FiniteVolumesDiffusion.hpp

@@ -278,8 +278,6 @@ public:
     Kokkos::View<double*> kL = unknowns.kL();
     Kokkos::View<double*> kR = unknowns.kR();
 
-    uR[0][0] = -t/(1.-t*t);
-
     const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
     const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
     const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();