From 5e55b997f679af8d8cb3ea49513512b29121beb9 Mon Sep 17 00:00:00 2001 From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr> Date: Fri, 1 Jun 2018 15:59:45 +0200 Subject: [PATCH] ajout test sur k dans Euler_unknowns et diffusion --- src/main.cpp | 29 +++++++++--------- src/mesh/Mesh.hpp | 16 +++++++--- src/scheme/AcousticSolver.hpp | 10 +++--- src/scheme/FiniteVolumesDiffusion.hpp | 37 +++++++++++++++-------- src/scheme/FiniteVolumesEulerUnknowns.hpp | 36 +++++++++++----------- 5 files changed, 72 insertions(+), 56 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index dff800896..d81da865e 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.8; + const double tmax=0.5; double t=0.; int itermax=std::numeric_limits<int>::max(); @@ -158,7 +158,7 @@ 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 +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 @@ -226,9 +226,9 @@ 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_u", std::ios::trunc); + std::ofstream fout("film", std::ios::trunc); fout.precision(15); std::string ligne; for (size_t j = 0; j<mesh.numberOfCells(); ++j) { @@ -238,7 +238,7 @@ int main(int argc, char *argv[]) fint.close(); fout.close(); - std::ifstream rint("film_u"); + std::ifstream rint("film"); std::ofstream rout("inter", std::ios::trunc); fout.precision(15); for (size_t j = 0; j<mesh.numberOfCells(); ++j) { @@ -253,12 +253,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); @@ -276,7 +276,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,13 +301,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] << '\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'; // kidder + fout << xj[j][0] << ' ' << rhoj[j] << '\n'; } } @@ -320,8 +319,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] << ' ' << -(xj[j][0]*tmax)/((50./9.)-tmax*tmax) << '\n'; // kidder + fout << xj[j][0] << ' ' << uj[j][0] << '\n'; } } diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp index 51c7f7d00..6dbde47b6 100644 --- a/src/mesh/Mesh.hpp +++ b/src/mesh/Mesh.hpp @@ -73,7 +73,7 @@ public: // pas constant - + Mesh(const Connectivity& connectivity) : m_connectivity(connectivity), m_xr("xr", connectivity.numberOfNodes()), @@ -81,7 +81,7 @@ public: m_xmax("xmax", 1) { double a = 0.; - double b = 1.; + double b = 2.; m_x0[0][0] = a; m_xmax[0][0] = b; const double delta_x = (b-a)/connectivity.numberOfCells(); @@ -93,12 +93,18 @@ public: // pas non constant - /* + /* Mesh(const Connectivity& connectivity) : m_connectivity(connectivity), - m_xr("xr", connectivity.numberOfNodes()) + m_xr("xr", connectivity.numberOfNodes()), + m_x0("x0", 1), + m_xmax("xmax", 1) { - const double delta_x = 1./connectivity.numberOfCells(); + double a = 0.; + double b = 2.; + m_x0[0][0] = a; + m_xmax[0][0] = b; + const double delta_x = (b-a)/connectivity.numberOfCells(); Kokkos::parallel_for(connectivity.numberOfNodes(), KOKKOS_LAMBDA(const int& r){ if (r%4 == 1) { m_xr[r][0] = (r*2+1)*0.5*delta_x; diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index a082f5654..c5533a3dc 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; } diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index 18a509ad3..0857c7cdb 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))); - */ + // Kidder // 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,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 ; @@ -346,8 +348,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))); }); @@ -358,7 +360,16 @@ 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 { + kj[j]=0.05 ; + //} + } }); // met a jour les quantites (geometriques) associees au maillage diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp index c7a624a6a..e50b95fea 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); + } + //Kidder + //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.05; - } else { - m_kj[j]=0. ; - } - }*/ + //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,7 +254,7 @@ public: m_uL[0] = zero; m_uR[0] = zero; m_kL[0] = 0.; - m_kR[0] = 1.; + m_kR[0] = 0.05; } -- GitLab