diff --git a/src/main.cpp b/src/main.cpp index 9f3e6482890f4354de078896faa811bf93e9583d..397f552a84dafbf4b23d2b04ce1e8c73d2d02f6c 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -154,21 +154,35 @@ int main(int argc, char *argv[]) BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj); while((t<tmax) and (iteration<itermax)) { - //double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj); - double dt = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj); - if (t+dt>tmax) { - dt=tmax-t; - } + // ETAPE 1 DU SPLITTING - EULER - //std::cout << "t=" << t << " dt=" << dt << '\n'; - //acoustic_solver.computeNextStep(t,dt, unknowns); - finite_volumes_diffusion.computeNextStep(t, dt, unknowns); - + double dt_euler = 0.4*acoustic_solver.acoustic_dt(Vj, cj); + + if (t+dt_euler > tmax) { + dt_euler = tmax-t; + } + acoustic_solver.computeNextStep(t,dt_euler, unknowns); + + // ETAPE 2 DU SPLITTING - DIFFUSION + + double dt_diff = finite_volumes_diffusion.diffusion_dt(rhoj, kj); + std::cout << dt_euler << ' ' << dt_diff << std::endl; + if (dt_euler <= dt_diff) { + dt_diff = dt_euler; + finite_volumes_diffusion.computeNextStep(t, dt_diff, unknowns); + t += dt_euler; + } else { + while (dt_euler > dt_diff) { + finite_volumes_diffusion.computeNextStep(t, dt_diff, unknowns); + dt_diff = 0.4*finite_volumes_diffusion.diffusion_dt(rhoj, kj); + std::cout << dt_diff << '\n'; + } + t += dt_diff; + } block_eos.updatePandCFromRhoE(); - t += dt; ++iteration; std::cout << "temps t : " << t << std::endl; } @@ -177,26 +191,28 @@ 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); std::cout << "* " << rang::style::underline << "Erreur L2 u" << rang::style::reset << ": " << rang::fgB::green << error << rang::fg::reset << " \n"; - /* + double error2 = 0.; error2 = finite_volumes_diffusion.error_Linf(unknowns); std::cout << "* " << rang::style::underline << "Erreur L infini u" << rang::style::reset << ": " << rang::fgB::green << error2 << rang::fg::reset << " \n"; - */ + double error3 = 0.; error3 = finite_volumes_diffusion.error_L2_E(unknowns); std::cout << "* " << rang::style::underline << "Erreur L2 E" << rang::style::reset << ": " << rang::fgB::green << error3 << rang::fg::reset << " \n"; - + */ + double cons = 0.; cons = finite_volumes_diffusion.conservatif(unknowns); @@ -210,29 +226,29 @@ int main(int argc, char *argv[]) { // gnuplot output for vitesse const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const Rd*> uj = unknowns.uj(); - double pi = 4.*std::atan(1.); - std::ofstream fout("comparaison u"); + // double pi = 4.*std::atan(1.); + std::ofstream fout("resultat u"); fout.precision(15); 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] << ' ' << std::sin(pi*xj[j][0])*std::exp(-0.2) <<'\n'; // cas k non constant + fout << xj[j][0] << ' ' << uj[j][0] << '\n'; } } { // gnuplot output for energy const Kokkos::View<const Rd*> xj = mesh_data.xj(); const Kokkos::View<const double*> Ej = unknowns.Ej(); - double pi = 4.*std::atan(1.); - std::ofstream fout("comparaison E"); + //double pi = 4.*std::atan(1.); + std::ofstream fout("resultat 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] << ' ' << ((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] << '\n'; } } - } Kokkos::finalize(); diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp index 75239af149754a9b69a1d33f436bb84c53ef19d1..447311e00b6936a0527a042b07c2222e1e6eb43d 100644 --- a/src/mesh/Mesh.hpp +++ b/src/mesh/Mesh.hpp @@ -50,7 +50,7 @@ public: // pas constant - /* + Mesh(const Connectivity& connectivity) : m_connectivity(connectivity), m_xr("xr", connectivity.numberOfNodes()) @@ -60,11 +60,11 @@ public: m_xr[r][0] = r*delta_x; }); } - */ + // pas non constant - + /* Mesh(const Connectivity& connectivity) : m_connectivity(connectivity), m_xr("xr", connectivity.numberOfNodes()) @@ -85,7 +85,7 @@ public: } }); } - + */ ~Mesh() { diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index 275214ac9e0c580227eb2da7f8d67a3eaa6870df..833e1e50de0459ac85458fecc99d707000e6ba8c 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -103,11 +103,8 @@ private: int cell_here = face_cells(l,j); int local_face_number_in_cell = face_cell_local_face(l,j); sum -= tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)); - // sum2 += kj(cell_here)*Vj(cell_here); - //sum3 += Vj(cell_here); sum2 += (1./Vj(cell_here))*kj(cell_here); sum3 += 1./Vj(cell_here); - // les deux moyennes donnent les memes resultats } m_Fl(l) = ((sum2/sum3)/Vl(l))*sum; @@ -246,7 +243,7 @@ public: double sum = 0.; for (int m = 0; m < cell_nb_nodes(j); ++m) { - sum += kj(cell_nodes(j,m)); + sum += kj(cell_nodes(j,m)); } dt_j[j]= 0.5*rhoj(j)*Vj(j)*(1./(2.*kj(j) + sum))*minVl; @@ -306,8 +303,8 @@ public: energy_fluxes += (Gl(l), Cjr(j,R)); } - //uj[j] += (dt*inv_mj[j]) * momentum_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 + uj[j] += (dt*inv_mj[j]) * momentum_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 Ej[j] += (dt*inv_mj[j]) * energy_fluxes; }); diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp index 06610a3f2a18bb038ea4c2166b9e084ca0472c75..c19955de51305bf0ffe97e2343d4676d35865fee 100644 --- a/src/scheme/FiniteVolumesEulerUnknowns.hpp +++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp @@ -186,60 +186,69 @@ public: // --- Acoustic Solver --- - // void initializeSod() - //{ - // 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) { - // m_rhoj[j]=1; - // } else { - // m_rhoj[j]=0.125; - // } - // }); - - //Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - // if (xj[j][0]<0.5) { - // m_pj[j]=1; - // } else { - // m_pj[j]=0.1; - // } - // }); - - //Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - // m_uj[j] = zero; - // }); - - //Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - // m_gammaj[j] = 1.4; - // }); - - //BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj); - //block_eos.updateEandCFromRhoP(); - - //Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - // m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]); - // }); - - //const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); - //Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - // m_mj[j] = m_rhoj[j] * Vj[j]; - // }); - - //Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - // m_inv_mj[j] = 1./m_mj[j]; - // }); - - //Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ - // m_kj[j] = 1; - // }); - //} - - // ------------------ - -void initializeSod() + void initializeSod() { 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) { + m_rhoj[j]=1; + } else { + m_rhoj[j]=0.125; + } + }); + + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + if (xj[j][0]<0.5) { + m_pj[j]=1; + } else { + m_pj[j]=0.1; + } + }); + + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + m_uj[j] = zero; + }); + + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + m_gammaj[j] = 1.4; + }); + + BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj); + block_eos.updateEandCFromRhoP(); + + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]); + }); + + const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + m_mj[j] = m_rhoj[j] * Vj[j]; + }); + + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + m_inv_mj[j] = 1./m_mj[j]; + }); + + Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ + m_kj[j] = 1; + }); + + // Conditions aux bords de Dirichlet sur u et k + + m_uL[0] = zero; + m_uR[0] = zero; + m_kL[0] = 1.; + m_kR[0] = 1.; + } + + + + /* DIFFUSION PURE + + void initializeSod() + { + const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); double pi = 4.*std::atan(1.); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ @@ -293,6 +302,8 @@ void initializeSod() } + */ + FiniteVolumesEulerUnknowns(const MeshDataType& mesh_data) : m_mesh_data(mesh_data),