diff --git a/src/main.cpp b/src/main.cpp index c86067c92979b7ec72491ff1ec08cc3ca2bb1847..b7f3bb4943d6e179dc5d8fd2dd47838baea63655 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -144,7 +144,7 @@ int main(int argc, char *argv[]) const Kokkos::View<const double*> Vj = mesh_data.Vj(); const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr(); - const double tmax=0.2; + const double tmax=0.8; double t=0.; int itermax=std::numeric_limits<int>::max(); @@ -317,7 +317,7 @@ int main(int argc, char *argv[]) finite_volumes_diffusion.computeNextStep(t, dt, unknowns); t += dt; */ - + // NAVIER-STOKES SANS SPLITTING double dt = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj,nuj, cj, nuL, nuR, kL,kR); @@ -658,18 +658,18 @@ int main(int argc, char *argv[]) const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const double*> rhoj = unknowns.rhoj(); double h = std::sqrt(1. - (tmax*tmax)/(50./9.)); - std::ofstream fout("rho"); + std::ofstream fout("rho_no_split"); fout.precision(15); for (size_t j=0; j<mesh.numberOfCells(); ++j) { - fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; // kidder - //fout << xj[j][0] << ' ' << rhoj[j] << '\n'; + //fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; // kidder + fout << xj[j][0] << ' ' << rhoj[j] << '\n'; } } { // gnuplot output for pression const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const double*> pj = unknowns.pj(); - std::ofstream fout("p"); + std::ofstream fout("p_no_split"); fout.precision(15); for (size_t j=0; j<mesh.numberOfCells(); ++j) { fout << xj[j][0] << ' ' << pj[j] << '\n'; @@ -679,7 +679,7 @@ int main(int argc, char *argv[]) { // gnuplot output for internal energy const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const double*> ej = unknowns.ej(); - std::ofstream fout("e"); + std::ofstream fout("e_no_split"); fout.precision(15); for (size_t j=0; j<mesh.numberOfCells(); ++j) { fout << xj[j][0] << ' ' << ej[j] << '\n'; @@ -690,16 +690,16 @@ 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"); + std::ofstream fout("u_no_split"); 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 //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << std::sin(pi*xj[j][0])*std::exp(-tmax) <<'\n'; // cas k non constant - fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder + //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << xj[j][0] << std::endl; - //fout << xj[j][0] << ' ' << uj[j][0] << '\n'; + fout << xj[j][0] << ' ' << uj[j][0] << '\n'; } } diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp index 007c679443e626e49b3aac3c4b5fa9d9af4a61a0..660632ef0a5c17739fd4cdc8ba5c75d12fc3a7a8 100644 --- a/src/mesh/Mesh.hpp +++ b/src/mesh/Mesh.hpp @@ -75,15 +75,15 @@ public: // pas constant - /* + Mesh(const Connectivity& connectivity) : m_connectivity(connectivity), m_xr("xr", connectivity.numberOfNodes()), m_x0("x0", 1), m_xmax("xmax", 1) { - double a = 0.; - double b = 1.; + double a = -1.; + double b = 2.; m_x0[0][0] = a; m_xmax[0][0] = b; const double delta_x = (b-a)/connectivity.numberOfCells(); @@ -91,7 +91,7 @@ public: m_xr[r][0] = a + r*delta_x; }); } - */ + // pas non constant /* @@ -123,17 +123,17 @@ public: } */ - + // pas aleatoire - + /* Mesh(const Connectivity& connectivity) : m_connectivity(connectivity), m_xr("xr", connectivity.numberOfNodes()), m_x0("x0", 1), m_xmax("xmax", 1) { - double a = 0.; - double b = 1.; + double a = -1.; + double b = 2.; m_x0[0][0] = a; m_xmax[0][0] = b; const double h = (b-a)/connectivity.numberOfCells(); @@ -145,7 +145,7 @@ public: for (int r=1; r<connectivity.numberOfNodes()-1; ++r){ const double delta_xr = dist(mt); - m_xr[r][0] = r*h + delta_xr; + m_xr[r][0] = a + r*h + delta_xr; } // creation fichier pour tracer h en fonction de x @@ -160,7 +160,7 @@ public: //std::exit(0); } - + */ ~Mesh() diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index 7351f938e0bb925acdec2838503ae00b71124885..40298dc5c01217cc5002c3f14f78e12559406ffb 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -182,18 +182,18 @@ private: }); - //m_ur[0]=zero; - //m_ur[m_mesh.numberOfNodes()-1]=zero; + m_ur[0]=zero; + m_ur[m_mesh.numberOfNodes()-1]=zero; //m_ur[0] = x0; //m_ur[m_mesh.numberOfNodes()-1] = xmax[0]; // CL Kidder - + /* double h = std::sqrt(1. - (t*t)/(50./9.)); m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0]; m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0]; - + */ return m_ur; } diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index 7dc3959daa79b2566c93ec5bca8dd0c989cd5d0c..a25682075df89fae5c16b0e50455b231f6b623cd 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -179,10 +179,10 @@ private: }); // Conditions aux bords - + m_Gl(0) = Fl(0)*uL(0); m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0); - + // Kidder /* double h = std::sqrt(1. - (t*t)/(50./9.)); diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp index 04a1d1afc13f62c77abb3348c49dadbad442431f..b60ee816d05a583dba36f7c94683471c2ed4ef74 100644 --- a/src/scheme/FiniteVolumesEulerUnknowns.hpp +++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp @@ -299,28 +299,28 @@ public: Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ // TAC - /* + if (xj[j][0]<0.5) { m_rhoj[j]=1.; } else { m_rhoj[j]=0.125; } - */ + // Kidder - m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.); + //m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.); }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ // TAC - /* + if (xj[j][0]<0.5) { - m_pj[j]=1; + m_pj[j]=1.; } else { m_pj[j]=0.1; } - */ + // Kidder - m_pj[j] = 2.*std::pow(m_rhoj[j],3); + //m_pj[j] = 2.*std::pow(m_rhoj[j],3); }); double pi = 4.*std::atan(1.); @@ -330,9 +330,9 @@ public: Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ // TAC - // m_gammaj[j] = 1.4; + m_gammaj[j] = 1.4; // Kidder - m_gammaj[j] = 3.; + //m_gammaj[j] = 3.; }); BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj); @@ -353,7 +353,7 @@ public: Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ // Differents k (xi) //m_kj[j] = xj[j][0]; - m_kj[j] = 0.005; + m_kj[j] = 0.014; // TAC @@ -408,8 +408,8 @@ public: m_uL[0] = zero; m_uR[0] = zero; - m_kL[0] = 0.005; - m_kR[0] = 0.005; + m_kL[0] = 0.014; + m_kR[0] = 0.014; Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ diff --git a/src/scheme/NoSplitting.hpp b/src/scheme/NoSplitting.hpp index b971ef440318acca9377d9d3a95a763e2a531fd6..92ffe8ce9f4180244ac2d66021266ca380c7eafd 100644 --- a/src/scheme/NoSplitting.hpp +++ b/src/scheme/NoSplitting.hpp @@ -191,18 +191,18 @@ private: }); - //m_ur[0]=zero; - //m_ur[m_mesh.numberOfNodes()-1]=zero; + m_ur[0]=zero; + m_ur[m_mesh.numberOfNodes()-1]=zero; //m_ur[0] = x0; //m_ur[m_mesh.numberOfNodes()-1] = xmax[0]; // CL Kidder - + /* double h = std::sqrt(1. - (t*t)/(50./9.)); m_ur[0]=(-t/((50./9.)-t*t))*h*x0[0]; m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*h*xmax[0]; - + */ return m_ur; } @@ -348,7 +348,7 @@ public: Ej[j] -= (dt*inv_mj[j]) * energy_fluxes; // ajout second membre pour kidder (k cst) - Ej[j] -= (dt*inv_mj[j])*Vj(j)*((kj(j)*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); + //Ej[j] -= (dt*inv_mj[j])*Vj(j)*((kj(j)*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))); @@ -384,11 +384,11 @@ public: } if (j == 0) { - //PTj(j) = pj(j) - kj(j)*(uj[j][0]-uL[0][0])/Vl(0); - PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t)); + PTj(j) = pj(j) - kj(j)*(uj[j][0]-uL[0][0])/Vl(0); + //PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t)); } else if (j == m_mesh.numberOfCells()-1) { - PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t)); - //PTj(j) = pj(j) - kj(j)*(uR[0][0]-uj[j][0])/Vl(m_mesh.numberOfFaces()-1); + //PTj(j) = pj(j) + kj(j)*(t/((50./9.)-t*t)); + PTj(j) = pj(j) - kj(j)*(uR[0][0]-uj[j][0])/Vl(m_mesh.numberOfFaces()-1); } else { //PTj(j) = pj(j) - kj(j)*sum/Vl(j); PTj(j) = pj(j) - kj(j)*2.*sum/sum1;