diff --git a/src/main.cpp b/src/main.cpp index 2225d80d2faae3c6da7565c74f5041a0c17d16c0..728770738495e6b3756e5d949cc3d356e4577420 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -141,7 +141,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(); @@ -158,14 +158,13 @@ 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); fout.precision(15); - //fout << std::fixed << 0. << ' ' << std::fixed << t << '\n'; for (size_t j=0; j<mesh.numberOfCells(); ++j) { - fout << std::fixed << xj[j][0] << ' ' << std::fixed << rhoj[j] << '\n'; + fout << std::fixed << xj[j][0] << ' ' << std::fixed << uj[j][0] << '\n'; } fout.close(); @@ -174,7 +173,7 @@ 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 @@ -188,7 +187,7 @@ int main(int argc, char *argv[]) t += dt_euler; // ETAPE 2 DU SPLITTING - DIFFUSION - /* + double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); if (dt_euler <= dt_diff) { @@ -205,7 +204,7 @@ int main(int argc, char *argv[]) t_diff += dt_diff; } } - */ + // DIFFUSION PURE /* @@ -223,27 +222,21 @@ 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_rho", std::ios::trunc); + std::ofstream fout("film_u", std::ios::trunc); fout.precision(15); - //fout.clear(); - //fout.seekp(0, std::ios::beg); std::string ligne; - //getline(fint,ligne); - //fout << ligne << ' ' << std::fixed << t << '\n'; for (size_t j = 0; j<mesh.numberOfCells(); ++j) { getline(fint, ligne); - fout << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << rhoj[j] << '\n'; + fout << ligne << ' ' << std::fixed << xj[j][0] << ' ' << std::fixed << uj[j][0] << '\n'; } fint.close(); fout.close(); - std::ifstream rint("film_rho"); + std::ifstream rint("film_u"); std::ofstream rout("inter", std::ios::trunc); fout.precision(15); - //getline(rint, ligne); - //rout << ligne << '\n'; for (size_t j = 0; j<mesh.numberOfCells(); ++j) { getline(rint, ligne); rout << ligne << '\n'; @@ -256,12 +249,12 @@ int main(int argc, char *argv[]) tempo.precision(5); tempo << std::fixed << t << '\n'; tempo.close(); - + */ } 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); @@ -279,7 +272,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); @@ -304,13 +297,13 @@ 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] << '\n'; + //fout << xj[j][0] << ' ' << rhoj[j] << '\n'; // Kidder - // fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; + fout << xj[j][0] << ' ' << rhoj[j] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; } } @@ -323,8 +316,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] << '\n'; - //fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; + //fout << xj[j][0] << ' ' << uj[j][0] << '\n'; + fout << xj[j][0] << ' ' << uj[j][0] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; } } diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index 883c088f8528d237652dc47bcbd36f077bf29b78..a082f56543579447b17c25a79e7153e7b0698713 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 @@ -206,10 +206,10 @@ private: //m_ur[0]=(-t/((50./9.)-t*t))*xr[0]; //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]; + //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]; return m_ur; } diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index 2c1c83ffcffaddfdf219686c1f7ef5ebc9cb31d1..18a509ad3c73aa6a611287c2c430660022478464 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,7 +124,7 @@ 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))); - + */ // k = 0.5 //m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5; @@ -134,9 +134,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,15 +176,15 @@ 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); //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 ; @@ -302,7 +302,7 @@ public: Kokkos::View<double*> ej = unknowns.ej(); Kokkos::View<double*> gammaj = unknowns.gammaj(); - const Kokkos::View<const double*> kj = unknowns.kj(); + Kokkos::View<double*> kj = unknowns.kj(); Kokkos::View<Rd*> uL = unknowns.uL(); Kokkos::View<Rd*> uR = unknowns.uR(); @@ -346,8 +346,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))); }); @@ -356,6 +356,11 @@ public: ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]); }); + // Mise a jour de k + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { + kj(j) = xj[j][0]; + }); + // met a jour les quantites (geometriques) associees au maillage //m_mesh_data.updateAllData(); diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp index 3579c7d7a82f8a1c7fb3c16c39f229b2ce2a86f0..c7a624a6a451aad8608f90aaf80c03212cb7e618 100644 --- a/src/scheme/FiniteVolumesEulerUnknowns.hpp +++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp @@ -191,23 +191,23 @@ 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) { + /*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) { + /*if (xj[j][0]<0.5) { m_pj[j]=1; } else { m_pj[j]=0.1; - } + }*/ - //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 +215,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 +236,17 @@ 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) { + /*if (xj[j][0]<0.7) { m_kj[j]=0.; } else { if (xj[j][0]<0.9){ - m_kj[j]=0.; + m_kj[j]=0.05; } else { m_kj[j]=0. ; } - } + }*/ }); // Conditions aux bords de Dirichlet sur u et k @@ -254,7 +254,7 @@ public: m_uL[0] = zero; m_uR[0] = zero; m_kL[0] = 0.; - m_kR[0] = 0.; + m_kR[0] = 1.; }