diff --git a/src/main.cpp b/src/main.cpp index e5f4b0b7fb67ff43f694cfc6a3a68b9158cf7590..d53f39766d0169f58ce22e42db1b9b4bc8a9ee1d 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -153,8 +153,14 @@ int main(int argc, char *argv[]) BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); - while((t<tmax) and (iteration<itermax)) { + double c = 0.; + c = finite_volumes_diffusion.conservatif(unknowns); + + std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset + << ": " << rang::fgB::green << c << rang::fg::reset << " \n"; + while((t<tmax) and (iteration<itermax)) { + // ETAPE 1 DU SPLITTING - EULER double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj); @@ -167,12 +173,11 @@ int main(int argc, char *argv[]) // ETAPE 2 DU SPLITTING - DIFFUSION - double dt_diff = finite_volumes_diffusion.diffusion_dt(rhoj, kj); + double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj); if (dt_euler <= dt_diff) { dt_diff = dt_euler; finite_volumes_diffusion.computeNextStep(t, dt_diff, unknowns); - t += dt_diff; } else { double t_diff = t; while (t + dt_euler > t_diff) { @@ -180,12 +185,11 @@ int main(int argc, char *argv[]) dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj); t_diff += dt_diff; } - t = t_diff; } // DIFFUSION PURE - + /* double dt = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj,kj); if (t+dt > tmax) { @@ -246,7 +250,7 @@ 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(); - std::ofstream fout("rho nounif"); + std::ofstream fout("rho euler"); fout.precision(15); for (size_t j=0; j<mesh.numberOfCells(); ++j) { fout << xj[j][0] << ' ' << rhoj[j] << '\n'; @@ -257,7 +261,7 @@ 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 nounif"); + std::ofstream fout("u euler"); 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 @@ -270,7 +274,7 @@ int main(int argc, char *argv[]) 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("E nounif"); + std::ofstream fout("E euler"); 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 diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp index 41e7b5b7f47f8b7ace653901a909b2cb6650b502..f1c38049b8f92f759b5203dcc6b9e6d0a242b196 100644 --- a/src/mesh/Mesh.hpp +++ b/src/mesh/Mesh.hpp @@ -50,7 +50,7 @@ public: // pas constant - /* + Mesh(const Connectivity& connectivity) : m_connectivity(connectivity), m_xr("xr", connectivity.numberOfNodes()) @@ -60,11 +60,11 @@ public: m_xr[r][0] = r*delta_x; }); } - */ + // pas non constant - + /* Mesh(const Connectivity& connectivity) : m_connectivity(connectivity), m_xr("xr", connectivity.numberOfNodes()) @@ -85,7 +85,7 @@ public: } }); } - + */ ~Mesh() { ; diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index db4daf818a2d68979ae50f65a58548b99ff726f4..fee4ed8d60cca3ab24521f06614143116ff1c088 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -240,14 +240,17 @@ public: for (int ll=0; ll<cell_nb_faces(j); ++ll) { minVl = std::min(minVl, Vl(cell_faces(j, ll))); } - + double sum = 0.; + for (int m = 0; m < cell_nb_nodes(j); ++m) { sum += kj(cell_nodes(j,m)); } dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl; + + }); double dt = std::numeric_limits<double>::max(); @@ -400,6 +403,7 @@ public: double conservatif(UnknownsType& unknowns) { Kokkos::View<double*> Ej = unknowns.Ej(); + Kokkos::View<double*> rhoj = unknowns.rhoj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); @@ -408,7 +412,7 @@ public: double sum = 0.; for (size_t j=0; j<m_mesh.numberOfCells(); ++j) { - sum += Ej[j]*Vj[j]; + sum += Ej[j]*rhoj[j]*Vj[j]; } return sum; diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp index fe4b0a30d833b48714b1fc1b5d778ec01835f4fd..01fc749928df8af873d8074e32fab16b4fe06042 100644 --- a/src/scheme/FiniteVolumesEulerUnknowns.hpp +++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp @@ -185,7 +185,7 @@ public: // --- Acoustic Solver --- - + void initializeSod() { const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); @@ -218,8 +218,8 @@ public: 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] = 2.; + m_Ej[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(); @@ -232,19 +232,19 @@ public: }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - m_kj[j] = 1.; + m_kj[j] = 0.00001*0.5; }); // Conditions aux bords de Dirichlet sur u et k m_uL[0] = zero; m_uR[0] = zero; - m_kL[0] = 1.; - m_kR[0] = 1.; + m_kL[0] = 0.00001*0.5; + m_kR[0] = 0.00001*0.5; } - - /* + + /* // DIFFUSION PURE @@ -254,7 +254,7 @@ public: double pi = 4.*std::atan(1.); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - m_rhoj[j]=1.; + m_rhoj[j]=1.; }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ @@ -287,8 +287,8 @@ public: }); 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){ @@ -299,8 +299,8 @@ public: m_uL[0] = zero; m_uR[0] = zero; - m_kL[0] = 0.; - m_kR[0] = 1.; + m_kL[0] = 2.; + m_kR[0] = 2.; } */