diff --git a/src/main.cpp b/src/main.cpp index 15f267d921445da0a9455334075ba6808d248bb1..b3e40857bb4859da61409dfb2a0fe364fde07a05 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -646,6 +646,7 @@ int main(int argc, char *argv[]) } } + /* { // gnuplot output for entropy (gaz parfait en prenant cv = 1)) const Kokkos::View<const Rd*> xj = mesh_data.xj(); diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index 0fb6058a12e75ce7055b8a3d3cc92b33651bfafe..d53830e188debc24a161f09707bd95c2b1f21536 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -140,6 +140,7 @@ private: m_Fl(0,0) = -(t/((50./9.)-t*t))*h*x0[0][0]; m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*h*xmax[0][0]; */ + return m_Fl ; } @@ -216,7 +217,7 @@ private: const Kokkos::View<const double*>& Vl = m_mesh_data.Vl(); const Kokkos::View<const double*>& Vj = m_mesh_data.Vj(); - + Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) { Rd sum = zero; double sum2 = 0.; @@ -237,10 +238,10 @@ private: int cell_here = face_cells(0,0); m_Bl(0) = (nuL(0) + nuj(cell_here))*(1./(2*Vl(0)))*(Tj(cell_here) - TL(0)); - + cell_here = face_cells(m_mesh.numberOfFaces()-1,0); m_Bl(m_mesh.numberOfFaces()-1) = -(nuR(0) + nuj(cell_here))*(1/(2.*Vl(m_mesh.numberOfFaces()-1)))*(Tj(cell_here) - TR(0)); - + return m_Bl ; } @@ -344,7 +345,7 @@ public: sum1 += nuj(cell_nodes(j,m)); } - if (sum < sum1) { + if (sum > sum1) { if (sum == 0.) { dt_j[j] = std::numeric_limits<double>::max(); } else { @@ -396,6 +397,10 @@ public: const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr(); + double pi = 4.*std::atan(1.); + TR(0) = 2-0.5*pi*pi*(std::exp(-2.*t)-1.); + // la condition au bord a droite de T depend du temps + // Calcule les flux computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, Tj, nuj, TL, TR, nuL, nuR, t); @@ -409,8 +414,6 @@ public: const Kokkos::View<const unsigned int*[2]>& cell_faces = m_connectivity.cellFaces(); - double pi = 4.*std::atan(1.); - // 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(), KOKKOS_LAMBDA(const int& j) { @@ -427,7 +430,7 @@ public: } uj[j] += (dt*inv_mj[j]) * momentum_fluxes; - Ej[j] += (dt*inv_mj[j]) * energy_fluxes; + Ej[j] += (dt*inv_mj[j]) * energy_fluxes; // test k non cst pour diff pure uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*(std::sin(pi*xj[j][0])*(pi*pi*xj[j][0]-1.) - std::cos(xj[j][0]*pi)*pi); @@ -446,18 +449,12 @@ public: 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(); - /* - // 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 de T par la formule T = E-0.5 u^2 + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { + Tj[j] = Ej[j] - 0.5 * (uj[j],uj[j]); }); - */ + } // Calcul erreur entre solution analytique et solution numerique en norme L2 diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp index b1750d2d5f39dd0a6b961a3402dc260c57275578..bcf1443d12e955dce68cd964a5e195d0cf78b5c6 100644 --- a/src/scheme/FiniteVolumesEulerUnknowns.hpp +++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp @@ -456,7 +456,7 @@ public: }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - m_Tj[j] = 2 - 0.5*std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]); // k = x + m_Tj[j] = 2. - 0.5*std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]); // k = x }); // Conditions aux bords de Dirichlet sur T et nu