From 1c63552656cf03543e8b1ad1f09a3c38ac1268aa Mon Sep 17 00:00:00 2001
From: Axelle <axelle.drouard@orange.fr>
Date: Thu, 14 Nov 2024 16:08:59 +0100
Subject: [PATCH] Symmetry boundary conditions
---
src/scheme/EulerKineticSolver2D.cpp | 218 ++++++++++++++++++++++++++--
1 file changed, 204 insertions(+), 14 deletions(-)
diff --git a/src/scheme/EulerKineticSolver2D.cpp b/src/scheme/EulerKineticSolver2D.cpp
index d7df345b4..470663f98 100644
--- a/src/scheme/EulerKineticSolver2D.cpp
+++ b/src/scheme/EulerKineticSolver2D.cpp
@@ -4,6 +4,7 @@
#include <language/utils/InterpolateItemArray.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/Mesh.hpp>
+#include <mesh/MeshFlatFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshVariant.hpp>
@@ -228,6 +229,8 @@ class EulerKineticSolver2D
const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
const NodeValuePerFace<const TinyVector<Dimension>> Nlr = MeshDataManager::instance().getMeshData(m_mesh).Nlr();
const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
+ const FaceValue<const TinyVector<Dimension>> nl = MeshDataManager::instance().getMeshData(m_mesh).nl();
+ const FaceValue<const TinyVector<Dimension>> Nl = MeshDataManager::instance().getMeshData(m_mesh).Nl();
CellValue<const double> m_inv_Vj;
const CellValue<const TinyVector<Dimension>> m_xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
@@ -245,7 +248,8 @@ class EulerKineticSolver2D
switch (bc_descriptor->type()) {
case IBoundaryConditionDescriptor::Type::symmetry: {
bc_list.emplace_back(
- SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
break;
}
case IBoundaryConditionDescriptor::Type::inflow_list: {
@@ -469,9 +473,6 @@ class EulerKineticSolver2D
if (li_nlr > 0) {
Fr[node_id][i_face][i_wave] += Fn[cell_id][i_wave] * li_nlr;
sum += li_nlr;
- // std::cout << "i_wave = " << i_wave << ", node_id = " << node_id << ", face_id = " << face_id << ",
- // cell_id = " << cell_id << ", lambda = " << m_lambda_vector[i_wave] << ", nlr = " << n_face << ", njr
- // = " << n_node << "\n";
}
}
if (sum != 0) {
@@ -484,6 +485,139 @@ class EulerKineticSolver2D
return Fr;
}
+ FaceArray<double>
+ compute_Flux1_face(const DiscreteFunctionP0Vector<const double> Fn)
+ {
+ const size_t nb_waves = m_lambda_vector.size();
+ FaceArray<double> Fl(m_mesh.connectivity(), nb_waves);
+ const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
+ const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
+ const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
+
+ parallel_for(
+ p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
+ const auto& face_to_cell = face_to_cell_matrix[face_id];
+
+ Fl[face_id][i_wave] = 0;
+
+ for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+ const CellId cell_id = face_to_cell[i_cell];
+ const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
+
+ TinyVector<Dimension> n_face = nl[face_id];
+ const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+
+ double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
+
+ if (li_nl > 0) {
+ Fl[face_id][i_wave] += Fn[cell_id][i_wave];
+ }
+ }
+ }
+ });
+
+ return Fl;
+ }
+
+ void
+ apply_scalar_bc(FaceArray<double> Fl_rho,
+ FaceArray<double> Fl_rho_u1,
+ FaceArray<double> Fl_rho_u2,
+ FaceArray<double> Fl_rho_E,
+ const DiscreteFunctionP0<const double> rho,
+ const DiscreteFunctionP0<const double> rho_u1,
+ const DiscreteFunctionP0<const double> rho_u2,
+ const DiscreteFunctionP0<const double> rho_E,
+ const BoundaryConditionList& bc_list)
+ {
+ const size_t nb_waves = m_lambda_vector.size();
+ const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
+ const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
+ const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
+
+ TinyVector<MeshType::Dimension> inv_S = zero;
+ for (size_t i = 0; i < nb_waves; ++i) {
+ for (size_t d = 0; d < MeshType::Dimension; ++d) {
+ inv_S[d] += std::pow(m_lambda_vector[i][d], 2);
+ }
+ }
+ for (size_t d = 0; d < MeshType::Dimension; ++d) {
+ inv_S[d] = 1. / inv_S[d];
+ }
+
+ for (auto&& bc_v : bc_list) {
+ std::visit(
+ [=, this](auto&& bc) {
+ using BCType = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+ auto face_list = bc.faceList();
+ auto n = bc.outgoingNormal();
+ auto nxn = tensorProduct(n, n);
+ TinyMatrix<Dimension> I = identity;
+ auto txt = I - nxn;
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
+ const auto& face_to_cell = face_to_cell_matrix[face_id];
+
+ for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+ const CellId cell_id = face_to_cell[i_cell];
+ const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
+
+ TinyVector<Dimension> n_face = nl[face_id];
+ const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+
+ double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
+ if (li_nl < 0) {
+ double rhoj_prime = rho[cell_id];
+ TinyVector<Dimension> rho_uj{rho_u1[cell_id], rho_u2[cell_id]};
+ TinyVector<Dimension> rho_uj_prime = txt * rho_uj - nxn * rho_uj;
+ double rho_Ej_prime = rho_E[cell_id];
+
+ double rho_e =
+ rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) *
+ (rho_u1[cell_id] * rho_u1[cell_id] + rho_u2[cell_id] * rho_u2[cell_id]);
+ double p = (m_gamma - 1) * rho_e;
+
+ TinyVector<Dimension> A_rho_prime = rho_uj_prime;
+ TinyVector<Dimension> A_rho_u1_prime{rho_uj_prime[0] * rho_uj_prime[0] / rhoj_prime + p,
+ rho_uj_prime[0] * rho_uj_prime[1] / rhoj_prime};
+ TinyVector<Dimension> A_rho_u2_prime{rho_uj_prime[0] * rho_uj_prime[1] / rhoj_prime,
+ rho_uj_prime[1] * rho_uj_prime[1] / rhoj_prime + p};
+ TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime;
+
+ double Fn_rho_prime = rhoj_prime / nb_waves +
+ inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_prime[0] +
+ inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_prime[1];
+ double Fn_rho_u1_prime = rho_uj_prime[0] / nb_waves +
+ inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_u1_prime[0] +
+ inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_u1_prime[1];
+ double Fn_rho_u2_prime = rho_uj_prime[1] / nb_waves +
+ inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_u2_prime[0] +
+ inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_u2_prime[1];
+ double Fn_rho_E_prime = rho_Ej_prime / nb_waves +
+ inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_E_prime[0] +
+ inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_E_prime[1];
+
+ Fl_rho[face_id][i_wave] += Fn_rho_prime;
+ Fl_rho_u1[face_id][i_wave] += Fn_rho_u1_prime;
+ Fl_rho_u2[face_id][i_wave] += Fn_rho_u2_prime;
+ Fl_rho_E[face_id][i_wave] += Fn_rho_E_prime;
+ }
+ }
+ }
+ }
+ }
+ },
+ bc_v);
+ }
+ }
+
DiscreteFunctionP0Vector<double>
compute_deltaLFn(NodeArray<double> Fr)
{
@@ -558,6 +692,36 @@ class EulerKineticSolver2D
return deltaLFn;
}
+ DiscreteFunctionP0Vector<double>
+ compute_deltaLFn_face(FaceArray<double> Fl)
+ {
+ const size_t nb_waves = m_lambda_vector.size();
+ DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
+
+ const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
+ const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ const auto& cell_to_face = cell_to_face_matrix[cell_id];
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ deltaLFn[cell_id][i_wave] = 0;
+
+ for (size_t i_face_cell = 0; i_face_cell < cell_to_face.size(); ++i_face_cell) {
+ const FaceId face_id = cell_to_face[i_face_cell];
+
+ const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+
+ const double li_Nl = sign * dot(m_lambda_vector[i_wave], Nl[face_id]);
+
+ deltaLFn[cell_id][i_wave] += Fl[face_id][i_wave] * li_Nl;
+ }
+ }
+ });
+ return deltaLFn;
+ }
+
CellValue<bool>
getBoundaryCells(const BoundaryConditionList& bc_list)
{
@@ -668,20 +832,32 @@ class EulerKineticSolver2D
// NodeArray<double> Fr_rho_u2 = compute_Flux1(Fn_rho_u2);
// NodeArray<double> Fr_rho_E = compute_Flux1(Fn_rho_E);
- FaceArrayPerNode<double> Fr_rho = compute_Flux1_eucchlyd(Fn_rho);
- FaceArrayPerNode<double> Fr_rho_u1 = compute_Flux1_eucchlyd(Fn_rho_u1);
- FaceArrayPerNode<double> Fr_rho_u2 = compute_Flux1_eucchlyd(Fn_rho_u2);
- FaceArrayPerNode<double> Fr_rho_E = compute_Flux1_eucchlyd(Fn_rho_E);
+ // FaceArrayPerNode<double> Fr_rho = compute_Flux1_eucchlyd(Fn_rho);
+ // FaceArrayPerNode<double> Fr_rho_u1 = compute_Flux1_eucchlyd(Fn_rho_u1);
+ // FaceArrayPerNode<double> Fr_rho_u2 = compute_Flux1_eucchlyd(Fn_rho_u2);
+ // FaceArrayPerNode<double> Fr_rho_E = compute_Flux1_eucchlyd(Fn_rho_E);
+
+ FaceArray<double> Fr_rho = compute_Flux1_face(Fn_rho);
+ FaceArray<double> Fr_rho_u1 = compute_Flux1_face(Fn_rho_u1);
+ FaceArray<double> Fr_rho_u2 = compute_Flux1_face(Fn_rho_u2);
+ FaceArray<double> Fr_rho_E = compute_Flux1_face(Fn_rho_E);
+
+ apply_scalar_bc(Fr_rho, Fr_rho_u1, Fr_rho_u2, Fr_rho_E, rho, rho_u1, rho_u2, rho_E, bc_list);
// DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn(Fr_rho);
// DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn(Fr_rho_u1);
// DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn(Fr_rho_u2);
// DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn(Fr_rho_E);
- DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn_eucclhyd(Fr_rho);
- DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn_eucclhyd(Fr_rho_u1);
- DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn_eucclhyd(Fr_rho_u2);
- DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn_eucclhyd(Fr_rho_E);
+ // DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn_eucclhyd(Fr_rho);
+ // DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn_eucclhyd(Fr_rho_u1);
+ // DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn_eucclhyd(Fr_rho_u2);
+ // DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn_eucclhyd(Fr_rho_E);
+
+ DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn_face(Fr_rho);
+ DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn_face(Fr_rho_u1);
+ DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn_face(Fr_rho_u2);
+ DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn_face(Fr_rho_E);
const CellArray<const TinyVector<Dimension>> APFn = getA(rho, rho_u1, rho_u2, rho_E);
const DiscreteFunctionP0<TinyVector<Dimension>> APFn_rho{p_mesh};
@@ -866,6 +1042,7 @@ class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition
private:
const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
+ const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;
public:
const Rd&
@@ -886,8 +1063,21 @@ class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition
return m_mesh_flat_node_boundary.nodeList();
}
- SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary)
- : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_flat_face_boundary.faceList().size();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_flat_face_boundary.faceList();
+ }
+
+ SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
+ const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
+ : m_mesh_flat_node_boundary(mesh_flat_node_boundary), m_mesh_flat_face_boundary(mesh_flat_face_boundary)
{
;
}
--
GitLab