diff --git a/src/scheme/FiniteVolumesDiffusion.hpp b/src/scheme/FiniteVolumesDiffusion.hpp index 9ae88e9fae927ace7b44b77032eacedc537d7d19..c91d7b6b1929a5f1fda53274d7b5c6c51577aff7 100644 --- a/src/scheme/FiniteVolumesDiffusion.hpp +++ b/src/scheme/FiniteVolumesDiffusion.hpp @@ -78,16 +78,32 @@ private: computeFl(const Kokkos::View<const Rd**>& Cjr, const Kokkos::View<const Rd*>& uj, const Kokkos::View<const double*>& kj) { - const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes(); - const Kokkos::View<const unsigned short*> cell_nb_nodes - = m_connectivity.cellNbNodes(); + + 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(); - Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { - for (int r=0; r<cell_nb_nodes[j]; ++r) { - m_Fl(j,r) = ((kj(cell_nodes(j,r)) + kj(cell_nodes(j-1,r)))/(2*Vl(j))) - * ((uj(j,r),Cjr(j,r))+ (uj(j-1,r),Cjr(j-1,r))) ; //tensorProduct(uj(j,r),Cjr(j,r)) ? + Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) { + Rdd sum; + 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, local_face_number_in_cell), Cjr(cell_here, local_face_number_in_cell)); } + + //Fl(l) = ((kj(cell_nodes(j,r)) + kj(cell_nodes(j-1,r)))/(2*Vl(j))) + //* ((uj(j,r),Cjr(j,r))+ (uj(j-1,r),Cjr(j-1,r))) ; //tensorProduct(uj(j,r),Cjr(j,r)) ? + + // k = 1 + + m_Fl(l)= (1./(2*Vl(l)))*sum; + }); return m_Fl ; @@ -96,17 +112,29 @@ private: // Calcule un Gl Kokkos::View<Rd*> computeGl(const Kokkos::View<const Rd*>& uj, - const Kokkos::View<const Rdd*>& Fl) { - const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes(); - const Kokkos::View<const unsigned short*> cell_nb_nodes - = m_connectivity.cellNbNodes(); + const Kokkos::View<const Rdd*>& Fl) { + + 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(); + + Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) { + Rd sum = zero; + 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 += uj(cell_here, local_face_number_in_cell); + } + + m_Gl(l) = 0.5*Fl(l)*sum; - Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) { - for (int r=0; r<cell_nb_nodes[j]; ++r) { - m_Gl(j,r) = 0.5*((uj(cell_nodes(j,r)) + uj(cell_nodes(j-1,r))),Fjr(j,r)); - } }); + return m_Gl ; } @@ -136,7 +164,7 @@ private: const Kokkos::View<const double*>& kj) { Kokkos::View<Rdd*> Fl = m_Fl ; Fl = computeFl (Cjr, uj, kj); - Kokkos::View<double**> Gl = m_Gl ; + Kokkos::View<Rd*> Gl = m_Gl ; Gl = computeGl (uj, Fl ); } @@ -162,15 +190,19 @@ public: 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(); 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))); @@ -212,17 +244,23 @@ public: const Kokkos::View<const Rdd*> Fl = m_Fl ; const Kokkos::View<const Rd *> Gl = m_Gl ; - const Kokkos::View<const unsigned short*> cell_nb_nodes - = m_connectivity.cellNbNodes(); + + 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; - for (int R=0; R<cell_nb_nodes[j]; ++R) { - momentum_fluxes += Fl (j,R); - energy_fluxes += Gl (j,R); + 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;