From b3841d69efc888d455ac83f6256f38b888d903c9 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Mon, 23 Apr 2018 17:10:22 +0200
Subject: [PATCH] Added quite dirty BC code. Hopefully on the road to 2d/3d ...
- pressure and symmetry by now - needs a cleaner implementation before
velocity implementation
---
src/main.cpp | 8 +-
src/scheme/AcousticSolver.hpp | 74 ++++++++++++++++++-
...ndition.hpp => FacesBoundaryCondition.hpp} | 0
3 files changed, 75 insertions(+), 7 deletions(-)
rename src/scheme/{BoundaryCondition.hpp => FacesBoundaryCondition.hpp} (100%)
diff --git a/src/main.cpp b/src/main.cpp
index c4ff74d54..a87619a5d 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 6f6d9a864..362352eb5 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
--
GitLab