#ifndef FINITE_VOLUMES_DIFFUSION_HPP #define FINITE_VOLUMES_DIFFUSION_HPP // --- INCLUSION fichiers headers --- #include <Kokkos_Core.hpp> #include <rang.hpp> #include <BlockPerfectGas.hpp> #include <TinyVector.hpp> #include <TinyMatrix.hpp> #include <Mesh.hpp> #include <MeshData.hpp> #include <FiniteVolumesEulerUnknowns.hpp> // --------------------------------- // Creation classe FiniteVolumesDiffusion template<typename MeshData> class FiniteVolumesDiffusion { typedef typename MeshData::MeshType MeshType; typedef FiniteVolumesEulerUnknowns<MeshData> UnknownsType; MeshData& m_mesh_data; const MeshType& m_mesh; const typename MeshType::Connectivity& m_connectivity; constexpr static size_t dimension = MeshType::dimension; typedef TinyVector<dimension> Rd; typedef TinyMatrix<dimension> Rdd; private: // Sert a calculer les reductions (en gros calculer le min sur des // vecteurs en parallele) struct ReduceMin { private: const Kokkos::View<const double*> x_; public: typedef Kokkos::View<const double*>::non_const_value_type value_type; ReduceMin(const Kokkos::View<const double*>& x) : x_ (x) {} typedef Kokkos::View<const double*>::size_type size_type; KOKKOS_INLINE_FUNCTION void operator() (const size_type i, value_type& update) const { if (x_(i) < update) { update = x_(i); } } KOKKOS_INLINE_FUNCTION void join (volatile value_type& dst, const volatile value_type& src) const { if (src < dst) { dst = src; } } KOKKOS_INLINE_FUNCTION void init (value_type& dst) const { // The identity under max is -Inf. dst= Kokkos::reduction_identity<value_type>::min(); } }; // Calcule un Fl Kokkos::View<Rdd*> computeFl(const Kokkos::View<const Rd**>& Cjr, const Kokkos::View<const Rd*>& uj, const Kokkos::View<const double*>& kj, const Kokkos::View<const Rd*>& uL, const Kokkos::View<const Rd*>& uR, const Kokkos::View<const double*>& kL, const Kokkos::View<const double*>& kR, const double& t) { const Kokkos::View<const unsigned int**>& face_cells = m_connectivity.faceCells(); const Kokkos::View<const unsigned short*> face_nb_cells = m_connectivity.faceNbCells(); const Kokkos::View<const unsigned short**> face_cell_local_face = m_mesh.connectivity().faceCellLocalFace(); const Kokkos::View<const double*>& Vl = m_mesh_data.Vl(); const Kokkos::View<const double*>& Vj = m_mesh_data.Vj(); const Kokkos::View<const Rd*> xr = m_mesh.xr(); Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) { Rdd sum = zero; double sum2 = 0.; double sum3 = 0.; for (int j=0; j<face_nb_cells(l); ++j) { 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 += (1./Vj(cell_here))*kj(cell_here); sum3 += 1./Vj(cell_here); } m_Fl(l) = ((sum2/sum3)/Vl(l))*sum; }); // 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))); //m_Fl(0) = -xr[0][0]*(tensorProduct(uj(cell_here), Cjr(cell_here, local_face_number_in_cell)) - tensorProduct(uL(0), Cjr(cell_here, local_face_number_in_cell))); cell_here = face_cells(m_mesh.numberOfFaces()-1,0); 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; //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]; return m_Fl ; } // Calcule un Gl Kokkos::View<Rd*> computeGl(const Kokkos::View<const Rd*>& uj, const Kokkos::View<const Rdd*>& Fl, const Kokkos::View<const Rd*>& uL, const Kokkos::View<const Rd*>& uR, const double& t) { const Kokkos::View<const unsigned int**>& face_cells = m_connectivity.faceCells(); const Kokkos::View<const unsigned short*> face_nb_cells = m_connectivity.faceNbCells(); const Kokkos::View<const unsigned short**> face_cell_local_face = m_mesh.connectivity().faceCellLocalFace(); const Kokkos::View<const double*>& Vj = m_mesh_data.Vj(); const Kokkos::View<const Rd*> xr = m_mesh.xr(); Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) { Rd sum = zero; double sum2 = 0.; for (int j=0; j<face_nb_cells(l); ++j) { int cell_here = face_cells(l,j); sum += (1./Vj(cell_here))*uj(cell_here); sum2 += 1./Vj(cell_here); } m_Gl(l) = (1./sum2)*Fl(l)*sum; }); // 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); return m_Gl ; } // Calcule la liste des inverse d'une liste de matrices // (pour l'instant, juste 1x1) void inverse(const Kokkos::View<const Rdd*>& A, Kokkos::View<Rdd*>& inv_A) const { Kokkos::parallel_for(A.size(), KOKKOS_LAMBDA(const int& r) { inv_A(r) = Rdd{1./(A(r)(0,0))}; }); } // Calcule la liste des inverses d'une liste de reels void inverse(const Kokkos::View<const double*>& x, Kokkos::View<double*>& inv_x) const { Kokkos::parallel_for(x.size(), KOKKOS_LAMBDA(const int& r) { inv_x(r) = 1./x(r); }); } // Enchaine les operations pour calculer les flux (Fjr et Gjr) pour // pouvoir derouler le schema KOKKOS_INLINE_FUNCTION void computeExplicitFluxes(const Kokkos::View<const Rd*>& uj, const Kokkos::View<const Rd**>& Cjr, const Kokkos::View<const double*>& kj, const Kokkos::View<const Rd*>& uL, const Kokkos::View<const Rd*>& uR, const Kokkos::View<const double*>& kL, const Kokkos::View<const double*>& kR, const double& t) { Kokkos::View<Rdd*> Fl = m_Fl ; Fl = computeFl (Cjr, uj, kj, uL, uR, kL, kR, t); Kokkos::View<Rd*> Gl = m_Gl ; Gl = computeGl (uj, Fl, uL, uR, t); } Kokkos::View<Rdd*> m_Fl; Kokkos::View<Rd*> m_Gl; public: FiniteVolumesDiffusion(MeshData& mesh_data, UnknownsType& unknowns) : m_mesh_data(mesh_data), m_mesh(mesh_data.mesh()), m_connectivity(m_mesh.connectivity()), m_Fl("Fl", m_mesh.numberOfFaces()), m_Gl("Gl", m_mesh.numberOfFaces()) { ; } // Calcule une evaluation du pas de temps verifiant le CFL parabolique // 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 { 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(); const Kokkos::View<const unsigned int**>& cell_faces = m_connectivity.cellFaces(); const Kokkos::View<const unsigned short*> cell_nb_faces = m_connectivity.cellNbFaces(); const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes(); const Kokkos::View<const unsigned short*> cell_nb_nodes = m_connectivity.cellNbNodes(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ double minVl = std::numeric_limits<double>::max(); for (int ll=0; ll<cell_nb_faces(j); ++ll) { minVl = std::min(minVl, Vl(cell_faces(j, ll))); } double sum = 0.; for (int m = 0; m < cell_nb_nodes(j); ++m) { 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; }); double dt = std::numeric_limits<double>::max(); Kokkos::parallel_reduce(m_mesh.numberOfCells(), ReduceMin(dt_j), dt); return dt; } // Avance la valeur des inconnues pendant un pas de temps dt void computeNextStep(const double& t, const double& dt, UnknownsType& unknowns) { Kokkos::View<double*> rhoj = unknowns.rhoj(); Kokkos::View<Rd*> uj = unknowns.uj(); Kokkos::View<Rd*> Sj = unknowns.Sj(); Kokkos::View<double*> Ej = unknowns.Ej(); Kokkos::View<double*> ej = unknowns.ej(); Kokkos::View<double*> gammaj = unknowns.gammaj(); const Kokkos::View<const double*> kj = unknowns.kj(); Kokkos::View<Rd*> uL = unknowns.uL(); Kokkos::View<Rd*> uR = unknowns.uR(); Kokkos::View<double*> kL = unknowns.kL(); Kokkos::View<double*> kR = unknowns.kR(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr(); // Calcule les flux computeExplicitFluxes(uj, Cjr, kj, uL, uR, kL, kR, t); const Kokkos::View<const Rdd*> Fl = m_Fl ; const Kokkos::View<const Rd *> Gl = m_Gl ; const Kokkos::View<const unsigned short*>& cell_nb_faces = m_connectivity.cellNbFaces(); const Kokkos::View<const unsigned int*[2]>& cell_faces = m_connectivity.cellFaces(); // Mise a jour de la vitesse et de l'energie totale specifique const Kokkos::View<const double*> inv_mj = unknowns.invMj(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { Rd momentum_fluxes = zero; double energy_fluxes = 0.; int l = 0; for (int R=0; R<cell_nb_faces(j); ++R) { l = cell_faces(j,R); momentum_fluxes += Fl(l)*Cjr(j,R); energy_fluxes += (Gl(l), Cjr(j,R)); } uj[j] += (dt*inv_mj[j]) * momentum_fluxes; 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)*((2.*xj[j][0]*t*t)/(((50./9.)-t*t)*((50./9.)-t*t))); }); // 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]); }); // met a jour les quantites (geometriques) associees au maillage //m_mesh_data.updateAllData(); /* // Calcul de rho avec la formule Mj = Vj rhoj const Kokkos::View<const double*> mj = unknowns.mj(); Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){ rhoj[j] = mj[j]/Vj[j]; }); */ } // Calcul erreur entre solution analytique et solution numerique en norme L2 // (quand la solution exacte est connue) double error_L2_u(UnknownsType& unknowns, const double& t) { Kokkos::View<Rd*> uj = unknowns.uj(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); //double pi = 4.*std::atan(1.); double err_u = 0.; double exact_u = 0.; for (size_t j=0; j<m_mesh.numberOfCells(); ++j) { //exact_u = std::sin(pi*xj[j][0])*std::exp(-2.*pi*pi*0.2); // solution exacte cas test k constant // exact_u = std::sin(pi*xj[j][0])*std::exp(-0.2); // solution exacte cas test k non constant exact_u = -(xj[j][0]*t)/((50./9.)-t*t); err_u += (exact_u - uj[j][0])*(exact_u - uj[j][0])*Vj(j); } err_u = std::sqrt(err_u); 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) { Kokkos::View<double*> Ej = unknowns.Ej(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); double pi = 4.*std::atan(1.); double err_E = 0.; double exact_E = 0.; for (size_t j=0; j<m_mesh.numberOfCells(); ++j) { //exact_E = (-(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.; exact_E = ((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.; err_E += (exact_E - Ej[j])*(exact_E - Ej[j])*Vj(j); } err_E = std::sqrt(err_E); return err_E; } */ // 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) { Kokkos::View<double*> rhoj = unknowns.rhoj(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); //double pi = 4.*std::atan(1.); double h = std::sqrt(1. - (t*t)/(50./9.)); 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::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]); } } return erreur; } // Verifie la conservativite de E double conservatif(UnknownsType& unknowns) { Kokkos::View<double*> Ej = unknowns.Ej(); Kokkos::View<double*> rhoj = unknowns.rhoj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); double sum = 0.; for (size_t j=0; j<m_mesh.numberOfCells(); ++j) { sum += Ej[j]*rhoj[j]*Vj[j]; } return sum; } // Verifie la conservativite de quantite de mvt double conservatif_mvt(UnknownsType& unknowns) { Kokkos::View<double*> rhoj = unknowns.rhoj(); Kokkos::View<Rd*> uj = unknowns.uj(); const Kokkos::View<const double*> Vj = m_mesh_data.Vj(); const Kokkos::View<const Rd*> xj = m_mesh_data.xj(); double sum = 0.; for (size_t j=0; j<m_mesh.numberOfCells(); ++j) { sum += Vj[j]*rhoj[j]*uj[j][0]; } return sum; } }; #endif // FINITE_VOLUMES_DIFFUSION_HPP