diff --git a/src/main.cpp b/src/main.cpp index 7a116ff1571f0b02d86a4688532f0d594361e6ad..895ea2577a091cb49fba9c6c591d25e4cf76e90a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -153,8 +153,8 @@ int main(int argc, char *argv[]) BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); - double c = 0.; - c = finite_volumes_diffusion.conservatif(unknowns); + //double c = 0.; + //c = finite_volumes_diffusion.conservatif(unknowns); while((t<tmax) and (iteration<itermax)) { @@ -210,19 +210,23 @@ 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 error = 0.; - error = finite_volumes_diffusion.error_L2_u(unknowns, tmax); - - std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset - << ": " << rang::fgB::green << error << rang::fg::reset << " \n"; + double error1 = 0.; + error1 = finite_volumes_diffusion.error_L2_rho(unknowns, tmax); + std::cout << "* " << rang::style::underline << "Erreur L2 rho" << rang::style::reset + << ": " << rang::fgB::green << error1 << rang::fg::reset << " \n"; double error2 = 0.; error2 = finite_volumes_diffusion.error_Linf(unknowns, tmax); std::cout << "* " << rang::style::underline << "Erreur L infini rho" << rang::style::reset << ": " << rang::fgB::green << error2 << rang::fg::reset << " \n"; + + double error = 0.; + error = finite_volumes_diffusion.error_L2_u(unknowns, tmax); + + std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset + << ": " << rang::fgB::green << error << rang::fg::reset << " \n"; /* double error3 = 0.; @@ -231,7 +235,7 @@ int main(int argc, char *argv[]) std::cout << "* " << rang::style::underline << "Erreur L2 E" << rang::style::reset << ": " << rang::fgB::green << error3 << rang::fg::reset << " \n"; */ - + /* std::cout << "* " << rang::style::underline << "Resultat conservativite rho E temps = 0" << rang::style::reset << ": " << rang::fgB::green << c << rang::fg::reset << " \n"; @@ -240,7 +244,7 @@ int main(int argc, char *argv[]) std::cout << "* " << rang::style::underline << "Resultat conservativite E" << rang::style::reset << ": " << rang::fgB::green << cons << rang::fg::reset << " \n"; - + */ //method_cost_map["AcousticSolverWithMesh"] = timer.seconds(); method_cost_map["FiniteVolumesDiffusionWithMesh"] = timer.seconds(); @@ -252,12 +256,10 @@ int main(int argc, char *argv[]) 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] << ' ' << std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h << '\n'; } } - { // gnuplot output for vitesse const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const Rd*> uj = unknowns.uj(); diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp index c19db1412272de2220ee9009dd741548f1cf9917..462ef51e6023f646c46079ff7813b955cb6e1437 100644 --- a/src/mesh/Mesh.hpp +++ b/src/mesh/Mesh.hpp @@ -85,7 +85,7 @@ public: } }); } -*/ + */ ~Mesh() { diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index 9abc8b4632584944771e2853cd615987babac737..78b6b6c62b1098d4069d2a5bbe90cc7017800fea 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -349,6 +349,8 @@ public: Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ 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; @@ -358,7 +360,7 @@ public: 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 0b49be7ef654ecc60cb95e5a714ace6977ecc466..c87a0eb3277363d7bf052abd0d3caf3619d1d4b6 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -123,8 +123,12 @@ private: 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 = x //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]; + + // k = 0.5 m_Fl(0,0) = -(t/((50./9.)-t*t))*0.5; m_Fl(m_mesh.numberOfFaces()-1,0) = -(t/((50./9.)-t*t))*0.5; @@ -166,6 +170,7 @@ 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) = -(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); @@ -321,19 +326,15 @@ public: Ej[j] += (dt*inv_mj[j]) * energy_fluxes; //uj[j] += std::exp(-t)*(dt*inv_mj[j])*Vj(j)*Sj(j) + (dt*inv_mj[j]) * momentum_fluxes; // test avec k non constant + // ajout second membre pour kidder (k = 0.5) + 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)*((0.5*t*t)/(((50./9.)-t*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[0] = -(t/((50./9.)-t*t))*xj[0]; - uj[m_mesh.numberOfCells()-1] = -(t/((50./9.)-t*t))*xj[m_mesh.numberOfCells()-1]; - - double h = std::sqrt(1. - (t*t)/(50./9.)); - 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; - + // Calcul de e par la formule e = E-0.5 u^2 Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]); @@ -376,6 +377,25 @@ public: return err_u; } + double error_L2_rho(UnknownsType& unknowns, const double& t) { + + Kokkos::View<double*> rhoj = unknowns.rhoj(); + + const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); + + const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); + + double h = std::sqrt(1. - (t*t)/(50./9.)); + double err_rho = 0.; + double exact_rho = 0.; + for (size_t j=0; j<m_mesh.numberOfCells(); ++j) { + exact_rho = std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h; + err_rho += (exact_rho - rhoj[j])*(exact_rho - rhoj[j])*Vj(j); + } + err_rho = std::sqrt(err_rho); + return err_rho; + } + /* double error_L2_E(UnknownsType& unknowns) { @@ -409,14 +429,10 @@ public: //double pi = 4.*std::atan(1.); double h = std::sqrt(1. - (t*t)/(50./9.)); - //double exacte = std::sin(pi*xj[0][0])*std::exp(-2.*pi*pi*0.2); // k constant - //double exacte = std::sin(pi*xj[0][0])*std::exp(-0.2); // k non constant double exacte = std::sqrt((4.*((xj[0][0]*xj[0][0])/(h*h)) + 100.-(xj[0][0]*xj[0][0])/(h*h))/100.)/h; double erreur = std::abs(exacte - rhoj[0]); for (size_t j=1; j<m_mesh.numberOfCells(); ++j) { - //exacte = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); - //exacte = std::sin(pi*xj[j][0])*std::exp(-0.2); exacte = std::sqrt((3.*((xj[j][0]*xj[j][0])/(h*h)) + 100.)/100.)/h; if (std::abs(exacte - rhoj[j]) > erreur) { erreur = std::abs(exacte - rhoj[j]);