diff --git a/src/main.cpp b/src/main.cpp index b2ca417a79d2607a1a800ec6b280bb7e8c1803c5..81926430be11b975371954b7b1317f46165e4463 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -138,7 +138,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.7; + const double tmax=1.5; double t=0; int itermax=std::numeric_limits<int>::max(); @@ -170,7 +170,7 @@ int main(int argc, char *argv[]) // ETAPE 2 DU SPLITTING - DIFFUSION - double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj); + double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); if (dt_euler <= dt_diff) { dt_diff = dt_euler; @@ -179,7 +179,7 @@ int main(int argc, char *argv[]) double t_diff = t-dt_euler; while (t > t_diff) { finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns); - dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj); + dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); if (t_diff+dt_diff > t) { dt_diff = t-t_diff; } @@ -208,7 +208,7 @@ int main(int argc, char *argv[]) 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); @@ -226,7 +226,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); @@ -251,11 +251,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.)); - std::ofstream fout("rho"); + //double h = std::sqrt(1. - (tmax*tmax)/(50./9.)); + std::ofstream fout("rho1.5"); 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'; + 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'; } } @@ -268,8 +270,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'; } } @@ -277,13 +279,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'; + //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 667f1877db466a9d581834099739504c9d78429b..ef3176a0371c4a39d7f0fdecc936bace02a6222d 100644 --- a/src/mesh/Mesh.hpp +++ b/src/mesh/Mesh.hpp @@ -80,8 +80,8 @@ public: m_x0("x0", 1), m_xmax("xmax", 1) { - double a = 0.; - double b = 10.; + double a = 0.1; + 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 a4b1250e1fd3c5bc4a6707f7af7cc3055b6a8a44..883c088f8528d237652dc47bcbd36f077bf29b78 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -181,7 +181,7 @@ private: } */ - Kokkos::View<Rd*> + Kokkos::View<Rd*> // calcule u_r (vitesse au sommet du maillage et flux de vitesse) computeUr(const Kokkos::View<const Rdd*>& Ar, const Kokkos::View<const Rd*>& br, const double& t) { @@ -197,14 +197,19 @@ private: m_ur[r]=invAr(r)*br(r); }); + m_ur[0]=zero; + m_ur[m_mesh.numberOfNodes()-1]=zero; + + // Kidder + // Conditions aux limites dependant du temps //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]; + //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 44e1b5d5d96a9fe6cf7b4ced24b517dc4cad31f3..2c1c83ffcffaddfdf219686c1f7ef5ebc9cb31d1 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 ; @@ -244,12 +244,11 @@ public: // Utilise la reduction definie dans la structure ReduceMin. KOKKOS_INLINE_FUNCTION double diffusion_dt(const Kokkos::View<const double*>& rhoj, - const Kokkos::View<const double*>& kj) const { + const Kokkos::View<const double*>& kj, + const Kokkos::View<const double*>& cj) const { Kokkos::View<double*> dt_j("dt_j", m_mesh.numberOfCells()); - const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr(); - const Kokkos::View<const double*>& Vl = m_mesh_data.Vl(); const Kokkos::View<const double*>& Vj = m_mesh_data.Vj(); @@ -276,8 +275,12 @@ public: sum += kj(cell_nodes(j,m)); } - // 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; + 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; + } }); @@ -343,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))); }); diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp index 99a8f11386b5bbd813a6fe62df9fc6d2ff117e52..d7eec6457c913d6f01d2a4918672943909f2a48c 100644 --- a/src/scheme/FiniteVolumesEulerUnknowns.hpp +++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp @@ -191,24 +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; } - */ - m_rhoj[j] = std::sqrt((3.*(xj[j][0]*xj[j][0]) + 100.)/100.); + //Kidder + //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){ @@ -216,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); @@ -237,16 +236,21 @@ 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 { + m_kj[j]=0.5; + } }); // Conditions aux bords de Dirichlet sur u et k m_uL[0] = zero; m_uR[0] = zero; - m_kL[0] = xj[0][0]; - m_kR[0] = xj[m_mesh.numberOfCells()-1][0]; + m_kL[0] = 0.01; + m_kR[0] = 0.5; }