diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 4e610f1459b01b2fb643dc45e891446d8c1e3241..3376c89de8da3ae851850aab70ad4db0e7ba00f0 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -942,13 +942,15 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
- -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>> {
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
return eulerKineticSolver2D(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho, F_rho_u1,
- F_rho_u2, F_E, rho_ex, rho_u_ex, rho_E_ex);
+ F_rho_u2, F_E, rho_ex, rho_u_ex, rho_E_ex,
+ bc_descriptor_list);
}
));
diff --git a/src/scheme/EulerKineticSolver2D.cpp b/src/scheme/EulerKineticSolver2D.cpp
index eabcc1f3a764fe28b0235ff55f5fdd9d40349af0..0f2513f67ce26de57d583909ba5c1696675c3336 100644
--- a/src/scheme/EulerKineticSolver2D.cpp
+++ b/src/scheme/EulerKineticSolver2D.cpp
@@ -3,11 +3,16 @@
#include <language/utils/EvaluateAtPoints.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/Mesh.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshVariant.hpp>
#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
#include <scheme/DiscreteFunctionP0Vector.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
@@ -197,6 +202,12 @@ class EulerKineticSolver2D
constexpr static size_t Dimension = MeshType::Dimension;
using Rd = TinyVector<Dimension>;
+ class FixedBoundaryCondition;
+ class SymmetryBoundaryCondition;
+
+ using BoundaryCondition = std::variant<FixedBoundaryCondition, SymmetryBoundaryCondition>;
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
const double m_dt;
const double m_eps;
const std::vector<TinyVector<Dimension>>& m_lambda_vector;
@@ -221,6 +232,39 @@ class EulerKineticSolver2D
const CellValue<const TinyVector<Dimension>> m_xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
const NodeValue<const TinyVector<Dimension>> m_xr = m_mesh.xr();
+ BoundaryConditionList
+ _getBCList(const MeshType& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ {
+ BoundaryConditionList bc_list;
+
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ bc_list.emplace_back(
+ SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::fixed: {
+ bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ break;
+ }
+ default: {
+ is_valid_boundary_condition = false;
+ }
+ }
+ if (not is_valid_boundary_condition) {
+ std::ostringstream error_msg;
+ error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ return bc_list;
+ }
+
public:
CellArray<TinyVector<MeshType::Dimension>>
getA(const DiscreteFunctionP0<const double>& rho,
@@ -473,24 +517,33 @@ class EulerKineticSolver2D
}
CellValue<bool>
- getBoundaryCells()
+ getBoundaryCells(const BoundaryConditionList& bc_list)
{
CellValue<bool> is_boundary_cell{p_mesh->connectivity()};
- const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
- auto is_boundary_node = p_mesh->connectivity().isBoundaryNode();
+ auto is_boundary_node = p_mesh->connectivity().isBoundaryNode();
is_boundary_cell.fill(false);
- parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
- const auto& cell_to_node = cell_to_node_matrix[cell_id];
- for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
- const NodeId node_id = cell_to_node[i_node];
- if (is_boundary_node[node_id]) {
- is_boundary_cell[cell_id] = true;
+ auto node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+
+ for (auto&& bc_v : bc_list) {
+ std::visit(
+ [=](auto&& bc) {
+ using BCType = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
+ auto node_list = bc.nodeList();
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ auto cell_list = node_to_cell_matrix[node_list[i_node]];
+ for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+ const CellId cell_id = cell_list[i_cell];
+ is_boundary_cell[cell_id] = true;
+ }
+ }
}
- }
- });
+ },
+ bc_v);
+ }
+
return is_boundary_cell;
}
@@ -498,11 +551,13 @@ class EulerKineticSolver2D
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
- apply()
+ apply(const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
+ BoundaryConditionList bc_list = this->_getBCList(m_mesh, bc_descriptor_list);
+
const size_t nb_waves = m_lambda_vector.size();
- CellValue<bool> is_boundary_cell = getBoundaryCells();
+ CellValue<bool> is_boundary_cell = getBoundaryCells(bc_list);
DiscreteFunctionP0<double> rho_u1_ex{p_mesh};
DiscreteFunctionP0<double> rho_u2_ex{p_mesh};
@@ -629,14 +684,6 @@ class EulerKineticSolver2D
Fnp1_rho_E[cell_id][i_wave] =
c1 * (m_eps * Fn_rho_E[cell_id][i_wave] - c2 * deltaLFn_rho_E[cell_id][i_wave] +
m_dt * MPFnp1_rho_E[cell_id][i_wave]);
- // std::cout << "inside Fnp1_rho[" << cell_id << "][" << i_wave << "] = " << Fnp1_rho[cell_id][i_wave] <<
- // "\n"; std::cout << "inside Fnp1_rho_u1[" << cell_id << "][" << i_wave << "] = " <<
- // Fnp1_rho_u1[cell_id][i_wave]
- // << "\n";
- // std::cout << "inside Fnp1_rho_u2[" << cell_id << "][" << i_wave << "] = " << Fnp1_rho_u2[cell_id][i_wave]
- // << "\n";
- // std::cout << "inside Fnp1_rho_E[" << cell_id << "][" << i_wave << "] = " << Fnp1_rho_E[cell_id][i_wave]
- // << "\n";
}
} else {
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
@@ -644,16 +691,6 @@ class EulerKineticSolver2D
Fnp1_rho_u1[cell_id][i_wave] = Fn_rho_u1_ex[cell_id][i_wave];
Fnp1_rho_u2[cell_id][i_wave] = Fn_rho_u2_ex[cell_id][i_wave];
Fnp1_rho_E[cell_id][i_wave] = Fn_rho_E_ex[cell_id][i_wave];
- // std::cout << "boundary Fnp1_rho[" << cell_id << "][" << i_wave << "] = " << Fnp1_rho[cell_id][i_wave]
- // << "\n";
- // std::cout << "boundary Fnp1_rho_u1[" << cell_id << "][" << i_wave << "] = " <<
- // Fnp1_rho_u1[cell_id][i_wave]
- // << "\n";
- // std::cout << "boundary Fnp1_rho_u2[" << cell_id << "][" << i_wave << "] = " <<
- // Fnp1_rho_u2[cell_id][i_wave]
- // << "\n";
- // std::cout << "boundary Fnp1_rho_E[" << cell_id << "][" << i_wave << "] = " << Fnp1_rho_E[cell_id][i_wave]
- // << "\n";
}
}
});
@@ -723,7 +760,8 @@ eulerKineticSolver2D(const double& dt,
const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex)
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
const std::shared_ptr<const MeshVariant> mesh_v =
getCommonMesh({F_rho, F_rho_u1, F_rho_u2, F_rho_E, rho_ex, rho_u_ex, rho_E_ex});
@@ -751,10 +789,67 @@ eulerKineticSolver2D(const double& dt,
EulerKineticSolver2D solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, Fn_rho, Fn_rho_u1, Fn_rho_u2,
Fn_rho_E, rho_n_ex, rho_u_n_ex, rho_E_n_ex);
- return solver.apply();
+ return solver.apply(bc_descriptor_list);
} else {
throw NotImplementedError("Invalid mesh type for 2D Euler Kinetic solver");
}
},
mesh_v->variant());
}
+
+template <MeshConcept MeshType>
+class EulerKineticSolver2D<MeshType>::FixedBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ FixedBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
+
+ ~FixedBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition
+{
+ public:
+ static constexpr const size_t Dimension = MeshType::Dimension;
+
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
+
+ public:
+ const Rd&
+ outgoingNormal() const
+ {
+ return m_mesh_flat_node_boundary.outgoingNormal();
+ }
+
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_flat_node_boundary.nodeList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_flat_node_boundary.nodeList();
+ }
+
+ SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary)
+ : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+ {
+ ;
+ }
+
+ ~SymmetryBoundaryCondition() = default;
+};
diff --git a/src/scheme/EulerKineticSolver2D.hpp b/src/scheme/EulerKineticSolver2D.hpp
index 437611c304b6d95e56a514b6bd1cf31cb2002511..76d024c82658527ee2fe895ada184929fabbb6ef 100644
--- a/src/scheme/EulerKineticSolver2D.hpp
+++ b/src/scheme/EulerKineticSolver2D.hpp
@@ -4,6 +4,8 @@
#include <language/utils/FunctionSymbolId.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
+class IBoundaryConditionDescriptor;
+
std::vector<TinyVector<2>> getLambdaVector(const std::shared_ptr<const MeshVariant>& mesh,
const double& lambda,
const size_t& L,
@@ -38,6 +40,7 @@ eulerKineticSolver2D(const double& dt,
const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex);
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
#endif // EULER_KINETIC_SOLVER_2D_HPP