From 3abe3d582bfec43bfc26f71c9e3ad57835563ed9 Mon Sep 17 00:00:00 2001 From: Fanny CHOPOT <fanny.chopot.ocre@cea.fr> Date: Wed, 6 Jun 2018 10:41:45 +0200 Subject: [PATCH] mise a jour pour stephane --- src/main.cpp | 27 ++++++---- src/scheme/AcousticSolver.hpp | 20 +++++++ src/scheme/FiniteVolumesDiffusion.hpp | 64 +++++++++++++---------- src/scheme/FiniteVolumesEulerUnknowns.hpp | 10 ++-- 4 files changed, 77 insertions(+), 44 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index e86fc7901..3826df03d 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.5; + const double tmax=1.5; double t=0.; int itermax=std::numeric_limits<int>::max(); @@ -190,7 +190,7 @@ int main(int argc, char *argv[]) // ETAPE 2 DU SPLITTING - DIFFUSION - double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); + double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); double t_diff = t-dt_euler; if (dt_euler <= dt_diff) { @@ -198,7 +198,7 @@ int main(int argc, char *argv[]) finite_volumes_diffusion.computeNextStep(t_diff, dt_diff, unknowns); } else { while (t > t_diff) { - dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); + dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); if (t_diff+dt_diff > t) { dt_diff = t-t_diff; } @@ -211,7 +211,7 @@ int main(int argc, char *argv[]) // AUTRE APPROCHE DU SPLITTING (PLUS LONG) /* double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj); - double dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); + double dt_diff = 0.9*finite_volumes_diffusion.diffusion_dt(rhoj, kj, cj); double dt = 0.; if (dt_euler < dt_diff) { dt = dt_euler; @@ -273,13 +273,14 @@ int main(int argc, char *argv[]) */ // ENTROPY TEST - //finite_volumes_diffusion.entropie(unknowns); + 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); @@ -287,7 +288,7 @@ int main(int argc, char *argv[]) << ": " << rang::fgB::green << error1 << rang::fg::reset << " \n"; double error2 = 0.; - error2 = finite_volumes_diffusion.error_Linf(unknowns, tmax); + error2 = finite_volumes_diffusion.error_Linf_rho(unknowns, tmax); std::cout << "* " << rang::style::underline << "Erreur L infini rho" << rang::style::reset << ": " << rang::fgB::green << error2 << rang::fg::reset << " \n"; @@ -297,7 +298,13 @@ 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 error4 = 0.; + error4 = finite_volumes_diffusion.error_Linf_u(unknowns, tmax); + + std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset + << ": " << rang::fgB::green << error4 << rang::fg::reset << " \n"; + /* double error3 = 0.; error3 = finite_volumes_diffusion.error_L2_E(unknowns); @@ -349,14 +356,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/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index 109bac2a8..ed4c99543 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -207,9 +207,11 @@ 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]; + return m_ur; } @@ -381,6 +383,24 @@ public: rhoj[j] = mj[j]/Vj[j]; }); + // Mise a jour de k + + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { + kj(j) = xj[j][0]; + //kj(j) = 0.5; + /* + if (xj[j][0]<0.7) { + kj[j]=0.; + } else { + if (xj[j][0]<0.9){ + kj[j]=0.05; + } else { + kj[j]=0. ; + } + } + */ + }); + } }; diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index ab50c2ea4..339a15e33 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -129,16 +129,16 @@ private: // Kidder // 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; + //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; // 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]; - //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 ; } @@ -185,11 +185,12 @@ private: //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); + return m_Gl ; } @@ -280,7 +281,7 @@ public: } if (sum == 0.) { - dt_j[j] = Vj[j]/cj[j]; + dt_j[j] = std::numeric_limits<double>::max();//Vj[j]/cj[j]; } else { dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./sum)*minVl; } @@ -346,11 +347,11 @@ public: //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))); + //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,24 +359,6 @@ public: Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { 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]; - kj(j) = 0.5; - /* - 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 //m_mesh_data.updateAllData(); @@ -458,7 +441,7 @@ public: // Calcul erreur entre solution analytique et solution numerique en norme L infini (max) // (quand la solution exacte est connue) - double error_Linf(UnknownsType& unknowns, const double& t) { + double error_Linf_rho(UnknownsType& unknowns, const double& t) { Kokkos::View<double*> rhoj = unknowns.rhoj(); @@ -479,6 +462,29 @@ public: return erreur; } + // Calcul erreur entre solution analytique et solution numerique en norme L infini (max) + // (quand la solution exacte est connue) + + double error_Linf_u(UnknownsType& unknowns, const double& t) { + + Kokkos::View<Rd*> uj = unknowns.uj(); + + const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); + + double exacte = -(xj[0][0]*t)/((50./9.)-t*t); + double erreur = std::abs(exacte - uj[0][0]); + + for (size_t j=1; j<m_mesh.numberOfCells(); ++j) { + exacte = -(xj[j][0]*t)/((50./9.)-t*t); + if (std::abs(exacte - uj[j][0]) > erreur) { + erreur = std::abs(exacte - uj[j][0]); + } + } + + return erreur; + } + + // Verifie la conservativite de E double conservatif(UnknownsType& unknowns) { diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp index be16504b6..405e94c1e 100644 --- a/src/scheme/FiniteVolumesEulerUnknowns.hpp +++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp @@ -272,8 +272,8 @@ public: }); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - //m_kj[j] = xj[j][0]; - m_kj[j] = 0.5; + m_kj[j] = xj[j][0]; + //m_kj[j] = 0.5; /* if (xj[j][0]<0.7) { m_kj[j]=0.; @@ -281,7 +281,7 @@ public: if (xj[j][0]<0.9){ m_kj[j]=0.05; } else { - m_kj[j]=0.05 ; + m_kj[j]=0. ; } } */ @@ -291,8 +291,8 @@ public: m_uL[0] = zero; m_uR[0] = zero; - m_kL[0] = 0.5; - m_kR[0] = 0.5; + m_kL[0] = 0.; + 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])); -- GitLab