diff --git a/src/main.cpp b/src/main.cpp index d81da865ed75c5ca855485fbd54faff72559962f..7065306c5d2ede69346ee0392613397ea355ba65 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -158,7 +158,8 @@ int main(int argc, char *argv[]) double c = 0.; c = finite_volumes_diffusion.conservatif(unknowns); - + + /* // Ecriture des valeurs initiales de rho dans un fichier const Kokkos::View<const Rd*> xj = mesh_data.xj(); std::ofstream fout("inter", std::ios::trunc); @@ -173,7 +174,8 @@ int main(int argc, char *argv[]) tempo.precision(5); tempo << std::fixed << t << '\n'; tempo.close(); - + */ + while((t<tmax) and (iteration<itermax)) { // ETAPE 1 DU SPLITTING - EULER @@ -201,10 +203,6 @@ int main(int argc, char *argv[]) dt_diff = t-t_diff; } finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns); - //dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); - //if (t_diff+dt_diff > t) { - // dt_diff = t-t_diff; - //} t_diff += dt_diff; } } @@ -226,7 +224,7 @@ int main(int argc, char *argv[]) std::cout << "temps t : " << t << std::endl; // ECRITURE DANS UN FICHIER - + /* std::ifstream fint("inter"); std::ofstream fout("film", std::ios::trunc); fout.precision(15); @@ -253,12 +251,16 @@ int main(int argc, char *argv[]) tempo.precision(5); tempo << std::fixed << t << '\n'; tempo.close(); + */ + + // ENTROPY TEST + finite_volumes_diffusion.entropie(unknowns); } std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n"; - /* + double error1 = 0.; error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax); @@ -276,7 +278,7 @@ int main(int argc, char *argv[]) std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset << ": " << rang::fgB::green << error << rang::fg::reset << " \n"; - */ + /* double error3 = 0.; error3 = finite_volumes_diffusion.error_L2_E(unknowns); @@ -301,12 +303,12 @@ int main(int argc, char *argv[]) { // gnuplot output for density 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.)); + double h = std::sqrt(1. - (tmax*tmax)/(50./9.)); std::ofstream fout("rho"); 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'; } } @@ -319,8 +321,8 @@ int main(int argc, char *argv[]) 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] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder - fout << xj[j][0] << ' ' << uj[j][0] << '\n'; + fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder + //fout << xj[j][0] << ' ' << uj[j][0] << '\n'; } } @@ -328,14 +330,14 @@ 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.); - //double h = std::sqrt(1. - (tmax*tmax)/(50./9.)); + double h = std::sqrt(1. - (tmax*tmax)/(50./9.)); std::ofstream fout("E"); 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 //fout << xj[j][0] << ' ' << Ej[j] << ' ' << ((xj[j][0]*pi*pi*0.5)*(std::sin(pi*xj[j][0])*std::sin(pi*xj[j][0]) - std::cos(xj[j][0]*pi)*std::cos(pi*xj[j][0])) - pi*0.5*std::sin(pi*xj[j][0])*std::cos(pi*xj[j][0]))*(std::exp(-2.*0.2)-1.) + 2. <<'\n' ; // cas k non constant - //fout << xj[j][0] << ' ' << Ej[j] << ' ' << (std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h) + (-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*(-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*0.5 << '\n'; // kidder - fout << xj[j][0] << ' ' << Ej[j] << '\n'; + fout << xj[j][0] << ' ' << Ej[j] << ' ' << (std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h) + (-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*(-(xj[j][0]*tmax)/((50./9.)-tmax*tmax))*0.5 << '\n'; // kidder + //fout << xj[j][0] << ' ' << Ej[j] << '\n'; } } diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp index 6dbde47b61f84ab8dba0706b10700115977f1797..1f26bcac1ae9559bfe59af8a5b34be4904490654 100644 --- a/src/mesh/Mesh.hpp +++ b/src/mesh/Mesh.hpp @@ -81,7 +81,7 @@ public: m_xmax("xmax", 1) { double a = 0.; - double b = 2.; + double b = 1.; m_x0[0][0] = a; m_xmax[0][0] = b; const double delta_x = (b-a)/connectivity.numberOfCells(); diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index c5533a3dcd298899b2501b5fb92be56454bbf05a..109bac2a8b389271bc33e5b03bbb696030b93f55 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -197,8 +197,8 @@ private: m_ur[r]=invAr(r)*br(r); }); - m_ur[0]=zero; - m_ur[m_mesh.numberOfNodes()-1]=zero; + //m_ur[0]=zero; + //m_ur[m_mesh.numberOfNodes()-1]=zero; // Kidder @@ -207,9 +207,9 @@ private: //m_ur[m_mesh.numberOfNodes()-1] = (-t/((50./9.)-t*t))*xr[m_mesh.numberOfNodes()-1]; //R(t) = x*h(t) a la place de x(t) - //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]; + 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; } @@ -330,6 +330,7 @@ public: Kokkos::View<double*> pj = unknowns.pj(); Kokkos::View<double*> gammaj = unknowns.gammaj(); Kokkos::View<double*> cj = unknowns.cj(); + Kokkos::View<double*> kj = unknowns.kj(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); @@ -380,17 +381,6 @@ public: rhoj[j] = mj[j]/Vj[j]; }); - /* - double h = std::sqrt(1. - (t*t)/(50./9.)); - rhoj[0] = std::sqrt((3.*((xj[0][0]*xj[0][0])/(h*h)) + 100.)/100.)/h; - rhoj[m_mesh.numberOfCells()-1] = std::sqrt((3.*((xj[m_mesh.numberOfCells()-1][0]*xj[m_mesh.numberOfCells()-1][0])/(h*h)) + 100.)/100.)/h; - - uj[0] = -(t/((50./9.)-t*t))*xj[0]; - uj[m_mesh.numberOfCells()-1] = -(t/((50./9.)-t*t))*xj[m_mesh.numberOfCells()-1]; - - Ej[0] = (std::sqrt((3.*((xj[0][0]*xj[0][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[0][0]*xj[0][0])/(h*h)) + 100.)/100.)/h) + (-(xj[0][0]*t)/((50./9.)-t*t))*(-(xj[0][0]*t)/((50./9.)-t*t))*0.5; - Ej[m_mesh.numberOfCells()-1] = (std::sqrt((3.*((xj[m_mesh.numberOfCells()-1][0]*xj[m_mesh.numberOfCells()-1][0])/(h*h)) + 100.)/100.)/h)*(std::sqrt((3.*((xj[m_mesh.numberOfCells()-1][0]*xj[m_mesh.numberOfCells()-1][0])/(h*h)) + 100.)/100.)/h) + (-(xj[m_mesh.numberOfCells()-1][0]*t)/((50./9.)-t*t))*(-(xj[m_mesh.numberOfCells()-1][0]*t)/((50./9.)-t*t))*0.5; - */ } }; diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index 0857c7cdb88e94e01b008a7ae365a24a5947a6f3..0bfa990ca8b4049de0a71c00504ef9e04dce8b81 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -115,7 +115,7 @@ private: }); // Conditions aux bords - + /* int cell_here = face_cells(0,0); int local_face_number_in_cell = face_cell_local_face(0,0); m_Fl(0) = -(kL(0) + kj(cell_here))*(1./(2*Vl(0)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell))); @@ -124,6 +124,8 @@ private: local_face_number_in_cell = face_cell_local_face(m_mesh.numberOfFaces()-1,0); m_Fl(m_mesh.numberOfFaces()-1) = -(kR(0) + kj(cell_here))*(1/(2.*Vl(m_mesh.numberOfFaces()-1)))*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell))); //m_Fl(m_mesh.numberOfFaces()-1) = -xr[m_mesh.numberOfNodes()-1][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uR(0), Cjr(cell_here, local_face_number_in_cell))); + */ + // Kidder // k = 0.5 @@ -134,9 +136,9 @@ private: //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 ; } @@ -176,17 +178,17 @@ private: }); // Conditions aux bords - m_Gl(0) = Fl(0)*uL(0); - m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0); + //m_Gl(0) = Fl(0)*uL(0); + //m_Gl(m_mesh.numberOfFaces()-1) = Fl(m_mesh.numberOfFaces()-1)*uR(0); // Kidder //m_Gl(0) = -(t/((50./9.)-t*t))*Fl(0,0)*xr(0); //m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*Fl(m_mesh.numberOfFaces()-1,0)*xr(m_mesh.numberOfFaces()-1); - //double h = std::sqrt(1. - (t*t)/(50./9.)); - //m_Gl(0) = -(t/((50./9.)-t*t))*h*Fl(0,0)*x0(0); - //m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*h*Fl(m_mesh.numberOfFaces()-1,0)*xmax(0); + double h = std::sqrt(1. - (t*t)/(50./9.)); + m_Gl(0) = -(t/((50./9.)-t*t))*h*Fl(0,0)*x0(0); + m_Gl(m_mesh.numberOfFaces()-1) = -(t/((50./9.)-t*t))*h*Fl(m_mesh.numberOfFaces()-1,0)*xmax(0); return m_Gl ; @@ -280,7 +282,6 @@ public: if (sum == 0.) { dt_j[j] = Vj[j]/cj[j]; } else { - // dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl; dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./sum)*minVl; } @@ -348,8 +349,8 @@ public: //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))); }); @@ -359,17 +360,20 @@ 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]; + /* if (xj[j][0]<0.7) { kj[j]=0.; } else { - //if (xj[j][0]<0.9){ - //kj[j]=0.05; - //} else { + if (xj[j][0]<0.9){ + kj[j]=0.05; + } else { kj[j]=0.05 ; - //} + } } + */ }); // met a jour les quantites (geometriques) associees au maillage @@ -402,7 +406,7 @@ public: for (size_t j=0; j<m_mesh.numberOfCells(); ++j) { //exact_u = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant // exact_u = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant - exact_u = -(xj[j][0]*t)/((50./9.)-t*t); + exact_u = -(xj[j][0]*t)/((50./9.)-t*t); // kidder err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j); } err_u = std::sqrt(err_u); @@ -515,6 +519,27 @@ public: return sum; } + + // Teste la croissance de l entropie d un gaz parfait + + void entropie(UnknownsType& unknowns) { + + Kokkos::View<double*> rhoj = unknowns.rhoj(); + Kokkos::View<double*> pj = unknowns.pj(); + Kokkos::View<double*> gammaj = unknowns.gammaj(); + + Kokkos::View<double*> entropy = unknowns.entropy(); + Kokkos::View<double*> entropy_new = unknowns.entropy_new(); + + 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 !"; + } + 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 e50b95fea886512d8dfaf9b20db09e08e8ad0266..a5455c79938d323b0d90cf9bdc6d29b837355816 100644 --- a/src/scheme/FiniteVolumesEulerUnknowns.hpp +++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp @@ -31,6 +31,8 @@ private: Kokkos::View<Rd*> m_uR; Kokkos::View<double*> m_kL; Kokkos::View<double*> m_kR; + Kokkos::View<double*> m_entropy; + Kokkos::View<double*> m_entropy_new; public: Kokkos::View<double*> rhoj() @@ -183,7 +185,26 @@ public: return m_kR; } - + Kokkos::View<double*> entropy() + { + return m_entropy; + } + + const Kokkos::View<const double*> entropy() const + { + return m_entropy; + } + + Kokkos::View<double*> entropy_new() + { + return m_entropy_new; + } + + const Kokkos::View<const double*> entropy_new() const + { + return m_entropy_new; + } + // --- Acoustic Solver --- void initializeSod() @@ -191,23 +212,27 @@ public: const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + /* 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){ + /* if (xj[j][0]<0.5) { 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); }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ @@ -215,8 +240,8 @@ public: }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - m_gammaj[j] = 1.4; - //m_gammaj[j] = 3.; + //m_gammaj[j] = 1.4; + m_gammaj[j] = 3.; }); BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj); @@ -236,17 +261,19 @@ public: }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - //m_kj[j] = xj[j][0]; + m_kj[j] = xj[j][0]; //m_kj[j] = 0.5; + /* if (xj[j][0]<0.7) { m_kj[j]=0.; } else { - //if (xj[j][0]<0.9){ - // m_kj[j]=0.05; - //} else { + if (xj[j][0]<0.9){ + m_kj[j]=0.05; + } else { m_kj[j]=0.05 ; - //} + } } + */ }); // Conditions aux bords de Dirichlet sur u et k @@ -254,8 +281,11 @@ public: m_uL[0] = zero; m_uR[0] = zero; m_kL[0] = 0.; - m_kR[0] = 0.05; + m_kR[0] = 1.; + 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])); + }); } @@ -336,10 +366,12 @@ public: m_inv_mj("inv_mj",m_mesh.numberOfCells()), m_kj("kj",m_mesh.numberOfCells()), m_Sj("Sj",m_mesh.numberOfCells()), - m_uL("uL", 1), - m_uR("uR", 1), - m_kL("kL", 1), - m_kR("kR", 1) + m_uL("uL", 1), + m_uR("uR", 1), + m_kL("kL", 1), + m_kR("kR", 1), + m_entropy("entropy", m_mesh.numberOfCells()), + m_entropy_new("entropy_new", m_mesh.numberOfCells()) { ;