diff --git a/src/main.cpp b/src/main.cpp index c4ff74d54433a9d1820941e2036c24e604e2bf35..a87619a5d40ced673f5093a00105591168ee70fc 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -13,7 +13,7 @@ #include <Connectivity1D.hpp> #include <Mesh.hpp> -#include <BoundaryCondition.hpp> +#include <FacesBoundaryCondition.hpp> #include <AcousticSolver.hpp> #include <TinyVector.hpp> @@ -132,7 +132,9 @@ int main(int argc, char *argv[]) FacesBoundaryCondition bc0(FacesBoundaryCondition::symmetry, 0, std::vector<unsigned int>({0u})); faces_bc_list.push_back(bc0); - FacesBoundaryCondition bc1(FacesBoundaryCondition::symmetry, 0, + // FacesBoundaryCondition bc1(FacesBoundaryCondition::symmetry, 0, + // std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())})); + FacesBoundaryCondition bc1(FacesBoundaryCondition::pressure, 0.125, std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())})); faces_bc_list.push_back(bc1); } @@ -141,7 +143,7 @@ int main(int argc, char *argv[]) unknowns.initializeSod(); - AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns); + AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, faces_bc_list); typedef TinyVector<MeshType::dimension> Rd; diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index 6f6d9a86492e11b2bb80d40f9c3f4108637099ff..362352eb5bb7b5456f727134025613a2f91f2209 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -11,9 +11,10 @@ #include <Mesh.hpp> #include <MeshData.hpp> #include <FiniteVolumesEulerUnknowns.hpp> +#include <FacesBoundaryCondition.hpp> template<typename MeshData> -class AcousticSolver +class AcousticSolver { typedef typename MeshData::MeshType MeshType; typedef FiniteVolumesEulerUnknowns<MeshData> UnknownsType; @@ -21,6 +22,7 @@ class AcousticSolver MeshData& m_mesh_data; const MeshType& m_mesh; const typename MeshType::Connectivity& m_connectivity; + const std::vector<FacesBoundaryCondition>& m_faces_boundary_condition_list; constexpr static size_t dimension = MeshType::dimension; @@ -136,6 +138,68 @@ private: return m_br; } + void applyBoundaryConditions() + { + for (const auto& face_bc : m_faces_boundary_condition_list) { + switch (face_bc.type()) { + case FacesBoundaryCondition::normal_velocity: { + std::cerr << "normal_velocity NIY\n"; + std::exit(0); + break; + } + case FacesBoundaryCondition::velocity: { + std::cerr << "velocity NIY\n"; + std::exit(0); + break; + } + case FacesBoundaryCondition::pressure: { + const Kokkos::View<const unsigned int**> face_cells + = m_connectivity.faceCells(); + const Kokkos::View<const unsigned short**> face_cell_local_face + = m_connectivity.faceCellLocalFace(); + const Kokkos::View<const Rd**> Cjr + = m_mesh_data.Cjr(); + + Kokkos::parallel_for(face_bc.numberOfFaces(), KOKKOS_LAMBDA(const int& l_number) { + // quite ugly: node/faces are melt in a bad way... at least works in 1d... + const int l = face_bc.faceList()[l_number]; + assert(m_connectivity.faceNbCells() == 1); + const unsigned int j = face_cells(l,0); + const unsigned int L = face_cell_local_face(l,0); + + m_br(l) -= face_bc.value()*Cjr(j,L); + }); + break; + } + case FacesBoundaryCondition::symmetry: { + const Rdd I = identity; + const Kokkos::View<const unsigned int**> face_cells + = m_connectivity.faceCells(); + const Kokkos::View<const unsigned short**> face_cell_local_face + = m_connectivity.faceCellLocalFace(); + const Kokkos::View<const Rd**> njr + = m_mesh_data.njr(); + + Kokkos::parallel_for(face_bc.numberOfFaces(), KOKKOS_LAMBDA(const int& l_number) { + // quite ugly: node/faces are melt in a bad way... at least works in 1d... + const int l = face_bc.faceList()[l_number]; + assert(m_connectivity.faceNbCells() == 1); + const unsigned int j = face_cells(l,0); + const unsigned int L = face_cell_local_face(l,0); + + const Rd& n = njr(j,L); + const Rdd nxn = tensorProduct(n,n); + const Rdd P = I-nxn; + + m_Ar(l) = P*m_Ar(l)*P + nxn; + m_br(l) = P*m_br(l); + }); + break; + } + } + } + } + Kokkos::View<Rd*> computeUr(const Kokkos::View<const Rdd*>& Ar, const Kokkos::View<const Rd*>& br) { @@ -144,8 +208,6 @@ private: Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) { m_ur[r]=invAr(r)*br(r); }); - m_ur[0]=zero; - m_ur[m_mesh.numberOfNodes()-1]=zero; return m_ur; } @@ -200,6 +262,8 @@ private: const Kokkos::View<const Rdd*> Ar = computeAr(Ajr); const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj); + this->applyBoundaryConditions(); + Kokkos::View<Rd*> ur = m_ur; Kokkos::View<Rd**> Fjr = m_Fjr; ur = computeUr(Ar, br); @@ -217,10 +281,12 @@ private: public: AcousticSolver(MeshData& mesh_data, - UnknownsType& unknowns) + UnknownsType& unknowns, + const std::vector<FacesBoundaryCondition>& faces_bc_list) : m_mesh_data(mesh_data), m_mesh(mesh_data.mesh()), m_connectivity(m_mesh.connectivity()), + m_faces_boundary_condition_list(faces_bc_list), m_br("br", m_mesh.numberOfNodes()), m_Ajr("Ajr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()), m_Ar("Ar", m_mesh.numberOfNodes()), diff --git a/src/scheme/BoundaryCondition.hpp b/src/scheme/FacesBoundaryCondition.hpp similarity index 100% rename from src/scheme/BoundaryCondition.hpp rename to src/scheme/FacesBoundaryCondition.hpp