diff --git a/src/main.cpp b/src/main.cpp index 7065306c5d2ede69346ee0392613397ea355ba65..a25a9b1b833c269e1e87d804765836bd9dda5a4f 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -180,7 +180,7 @@ int main(int argc, char *argv[]) // ETAPE 1 DU SPLITTING - EULER - double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj); + double dt_euler = 0.1*acoustic_solver.acoustic_dt(Vj, cj); if (t+dt_euler > tmax) { dt_euler = tmax-t; @@ -192,7 +192,7 @@ int main(int argc, char *argv[]) double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); double t_diff = t-dt_euler; - + if (dt_euler <= dt_diff) { dt_diff = dt_euler; finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns); @@ -207,6 +207,25 @@ int main(int argc, char *argv[]) } } + + // AUTRE APPROCHE DU SPLITTING (PLUS LONG) + /* + double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj); + double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); + double dt = 0.; + if (dt_euler < dt_diff) { + dt = dt_euler; + } else { + dt = dt_diff; + } + if (t+dt > tmax) { + dt = tmax-t; + } + acoustic_solver.computeNextStep(t,dt, unknowns); + finite_volumes_diffusion.computeNextStep(t, dt, unknowns); + t += dt; + */ + // DIFFUSION PURE /* diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index 0bfa990ca8b4049de0a71c00504ef9e04dce8b81..ab50c2ea4a012981b868f749d23d955c58841b55 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -129,16 +129,16 @@ private: // Kidder // k = 0.5 - //m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5; - //m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*0.5; + m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5; + m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*0.5; // k = x //m_Fl(0,0) = -(t/((50./9.)-t*t))*xr[0][0]; //m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*xr[m_mesh.numberOfFaces()-1][0]; - double h = std::sqrt(1. - (t*t)/(50./9.)); - 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]; + //double h = std::sqrt(1. - (t*t)/(50./9.)); + //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 ; } @@ -346,11 +346,11 @@ public: //uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*Sj(j) + (dt*inv_mj[j]) * momentum_fluxes; // test avec k non constant // ajout second membre pour kidder (k = 0.5) - //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((0.5*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); + Ej[j] -= (dt*inv_mj[j])*Vj(j)*((0.5*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); // ajout second membre pour kidder (k = x) - uj[j][0] += (dt*inv_mj[j])*Vj(j)*(t/((50./9.)-t*t)); - Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); + //uj[j][0] += (dt*inv_mj[j])*Vj(j)*(t/((50./9.)-t*t)); + //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); }); @@ -362,7 +362,8 @@ public: // Mise a jour de k Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { - kj(j) = xj[j][0]; + //kj(j) = xj[j][0]; + kj(j) = 0.5; /* if (xj[j][0]<0.7) { kj[j]=0.; @@ -534,7 +535,8 @@ public: Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ entropy_new(j) = std::log(pj(j)*std::pow(rhoj(j),-gammaj(j))); if (entropy_new(j) - entropy(j) < 0) { - std::cout << "ATTENTION ENTROPIE NEGATIVE !"; + std::cout << "ATTENTION ENTROPIE NEGATIVE !" << std::endl; + std::cout << "j = " << j << " difference = " << entropy_new(j) - entropy(j) << std::endl; } entropy(j) = std::log(pj(j)*std::pow(rhoj(j),-gammaj(j))); }); diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp index a5455c79938d323b0d90cf9bdc6d29b837355816..f9bccfe27aa3dbcd948cb367db98d54bd2428942 100644 --- a/src/scheme/FiniteVolumesEulerUnknowns.hpp +++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp @@ -261,8 +261,8 @@ public: }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - m_kj[j] = xj[j][0]; - //m_kj[j] = 0.5; + //m_kj[j] = xj[j][0]; + m_kj[j] = 0.5; /* if (xj[j][0]<0.7) { m_kj[j]=0.; @@ -280,8 +280,8 @@ public: m_uL[0] = zero; m_uR[0] = zero; - m_kL[0] = 0.; - m_kR[0] = 1.; + m_kL[0] = 0.5; + m_kR[0] = 0.5; Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ m_entropy(j) = std::log(m_pj[j]*std::pow(m_rhoj[j],-m_gammaj[j]));