diff --git a/src/main.cpp b/src/main.cpp index babd25dcb30af0137ea287c5ebff1f37e20f0511..e5adc9e7bcdf93ca3918fd1d47313c7a3d9a6117 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -189,17 +189,34 @@ int main(int argc, char *argv[]) std::cout << "* " << rang::style::underline << "Erreur L infini" << rang::style::reset << ": " << rang::fgB::green << error2 << rang::fg::reset << " \n"; + double cons = 0.; + cons = finite_volumes_diffusion.conservatif(unknowns); + + std::cout << "* " << rang::style::underline << "Resultat conservativite" << rang::style::reset + << ": " << rang::fgB::green << cons << rang::fg::reset << " \n"; + + //method_cost_map["AcousticSolverWithMesh"] = timer.seconds(); method_cost_map["FiniteVolumesDiffusionWithMesh"] = timer.seconds(); - { // gnuplot output for density + { // gnuplot output for vitesse 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("comparaison_u_k_non_cst_maill_non_unif_320"); + std::ofstream fout("comparaison u"); + 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 + } + } + + { // gnuplot output for energy + const Kokkos::View<const Rd*> xj = mesh_data.xj(); + const Kokkos::View<const double*> Ej = unknowns.Ej(); + double pi = 4.*std::atan(1.); + std::ofstream fout("comparaison E"); 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] << ' ' << Ej[j] << ' ' << (std::cos(pi*xj[j][0]))*(std::cos(pi*xj[j][0]))*std::sin(pi*xj[j][0])*((pi*pi*pi)*0.5)*std::exp(-4.*pi*0.2) <<'\n'; } } diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index c01f1883c6bc98a86a6b5973d5d1036667844ab3..ecac5a30b3324952978186ba9efbdb4772118a14 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -93,10 +93,10 @@ private: Rdd sum = zero; double sum2 = 0.; for (int j=0; j<face_nb_cells(l); ++j) { - int cell_here = face_cells(l,j); - int local_face_number_in_cell = face_cell_local_face(l,j); - sum -= tensorProduct(uj(cell_here, local_face_number_in_cell), Cjr(cell_here, local_face_number_in_cell)); - sum2 += kj(cell_here); + int cell_here = face_cells(l,j); + int local_face_number_in_cell = face_cell_local_face(l,j); + sum -= tensorProduct(uj(cell_here, local_face_number_in_cell), Cjr(cell_here, local_face_number_in_cell)); + sum2 += kj(cell_here); } m_Fl(l) = ((sum2/face_nb_cells(l))/Vl(l))*sum; @@ -127,7 +127,7 @@ private: sum += uj(cell_here, local_face_number_in_cell); } - m_Gl(l) = 0.5*Fl(l)*sum; + m_Gl(l) = (1./face_nb_cells(l))*Fl(l)*sum; }); @@ -209,10 +209,6 @@ public: minVl = std::min(minVl, Vl(cell_faces(j, ll))); } - //k=2 => (kj(j+1) + 2*kj(j) + kj(j-1)) = 8 - //dt_j[j]= 0.5*rhoj(j)*Vj(j)*(2./8.)*minVl; - - // k pas forcement constant double sum = 0.; for (int m = 0; m < cell_nb_nodes(j); ++m) { sum += kj(cell_nodes(j,m)); @@ -274,8 +270,8 @@ public: momentum_fluxes += Fl(l)*Cjr(j,R); energy_fluxes += (Gl(l), Cjr(j,R)); } - uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*Sj(j) + (dt*inv_mj[j]) * momentum_fluxes; // test avec k non constant - //uj[j] += (dt*inv_mj[j]) * momentum_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 + uj[j] += (dt*inv_mj[j]) * momentum_fluxes; Ej[j] += (dt*inv_mj[j]) * energy_fluxes; }); @@ -310,8 +306,8 @@ public: double erreur = 0.; double exacte = 0.; for (size_t j=0; j<m_mesh.numberOfCells(); ++j) { - //exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant - exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant + exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant + //exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant erreur += (exacte - uj[j][0])*(exacte - uj[j][0])*Vj(j); } erreur = std::sqrt(erreur); @@ -328,13 +324,13 @@ public: const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); double pi = 4.*std::atan(1.); - //double exacte = std::sin(pi*xj[0][0])*std::exp(-2.*pi*pi*0.2); // k constant - double exacte = std::sin(pi*xj[0][0])*std::exp(-0.2); // k non constant + double exacte = std::sin(pi*xj[0][0])*std::exp(-2.*pi*pi*0.2); // k constant + //double exacte = std::sin(pi*xj[0][0])*std::exp(-0.2); // k non constant double erreur = std::abs(exacte - uj[0][0]); for (size_t j=1; j<m_mesh.numberOfCells(); ++j) { - //exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); - exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); + exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); + //exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); if (std::abs(exacte - uj[j][0]) > erreur) { erreur = std::abs(exacte - uj[j][0]); @@ -345,6 +341,24 @@ public: return erreur; } + // Verifie la conservativite de E + + double conservatif(UnknownsType& unknowns) { + + Kokkos::View<double*> Ej = unknowns.Ej(); + + const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); + + double sum = 0.; + + for (size_t j=0; j<m_mesh.numberOfCells(); ++j) { + sum += Ej[j]*Vj[j]; + } + + return sum; + + } + }; diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp index 5f8e2fceefe406e0a7effe78c12cc57d8255ae7c..cf2d371586dd46a73a0f1028fb6478d9e8b95943 100644 --- a/src/scheme/FiniteVolumesEulerUnknowns.hpp +++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp @@ -204,7 +204,18 @@ void initializeSod() }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - m_uj[j] = std::sin(pi*xj[j][0]); + if (j == m_mesh.numberOfCells()) { + m_uj[j] = zero; + } + else { + if (j == 0) { + m_uj[j] = zero; + } + else { + m_uj[j] = std::sin(pi*xj[j][0]); + } + } + }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ @@ -215,7 +226,9 @@ void initializeSod() block_eos.updateEandCFromRhoP(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]); + //m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]); + m_Ej[j] = 0.5*(m_uj[j],m_uj[j]); + //m_Ej[j] = 2.; }); const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); @@ -228,8 +241,8 @@ void initializeSod() }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - //m_kj[j] = 2.; // k constant - m_kj[j] = xj[j][0]; // k non constant, k = x + m_kj[j] = 2.; // k constant + //m_kj[j] = xj[j][0]; // k non constant, k = x }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){