From f8c0e599dbc7ebab2a093a442eb5d2af2d84e59e Mon Sep 17 00:00:00 2001
From: Axelle <axelle.drouard@orange.fr>
Date: Mon, 3 Feb 2025 15:18:31 +0100
Subject: [PATCH] Give symmetry BC to the reconstruction for scalar functions
---
src/scheme/EulerKineticSolverMultiD.cpp | 4 +-
src/scheme/EulerKineticSolverMultiDOrder2.cpp | 140 ++++++++++++------
2 files changed, 98 insertions(+), 46 deletions(-)
diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
index 010ddd67e..f299da04d 100644
--- a/src/scheme/EulerKineticSolverMultiD.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -46,8 +46,8 @@ getLambdaVector2D(const std::shared_ptr<const MeshVariant>& mesh,
for (size_t m = 1; m < 4 * M + 1; ++m) {
size_t i_wave = (l - 1) * 4 * M + m - 1;
double ll = l;
- lambda_vector[i_wave] = TinyVector<2>{(ll / L) * lambda * std::cos((m * pi) / (2 * M)),
- (ll / L) * lambda * std::sin((m * pi) / (2 * M))};
+ lambda_vector[i_wave] = TinyVector<2>{(ll / L) * lambda * std::cos((m * pi) / (2 * M) + 0.25 * pi),
+ (ll / L) * lambda * std::sin((m * pi) / (2 * M) + 0.25 * pi)};
}
}
// lambda_vector[4 * L * M] = TinyVector<2>{0, 0}; //!!!!!!!!!!!!!!!!!!!!
diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.cpp b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
index 8651c4399..1df024308 100644
--- a/src/scheme/EulerKineticSolverMultiDOrder2.cpp
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
@@ -454,7 +454,7 @@ class EulerKineticSolverMultiDOrder2
FaceArray<T>
compute_Flux_face(const DiscreteFunctionDPkVector<Dimension, const T>& DPk_Fn,
const CellArray<const T>& Fn,
- CellValue<bool>& is_wall_boundary_cell)
+ FaceValue<bool>& is_wall_boundary_face)
{
const size_t nb_waves = m_lambda_vector.size();
FaceArray<T> Fl(m_mesh.connectivity(), nb_waves);
@@ -484,7 +484,7 @@ class EulerKineticSolverMultiDOrder2
double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
if (li_nl > 0) {
- if (not is_wall_boundary_cell[cell_id]) {
+ if ((not is_wall_boundary_face[face_id])) { //} and (not is_tiny_vector_v<T>)) {
Fl[face_id][i_wave] += DPk_Fn(cell_id, i_wave)(m_xl[face_id]);
} else {
Fl[face_id][i_wave] += Fn[cell_id][i_wave];
@@ -661,7 +661,7 @@ class EulerKineticSolverMultiDOrder2
for (size_t k_wave = 0; k_wave < nb_waves; ++k_wave) {
rhol += Fn_rho(cell_id, k_wave)(m_xl[face_id]);
- rho_ul += Fn_rho_u(cell_id, k_wave)(m_xl[face_id]);
+ rho_ul += Fn_rho_u(cell_id, k_wave)(m_xl[face_id]); // Fj_rho_u[cell_id][k_wave]; //
rho_El += Fn_rho_E(cell_id, k_wave)(m_xl[face_id]);
}
@@ -1676,25 +1676,19 @@ class EulerKineticSolverMultiDOrder2
bc_v);
}
- CellValue<bool> is_wall_boundary_cell{p_mesh->connectivity()};
- is_wall_boundary_cell.fill(false);
+ FaceValue<bool> is_wall_boundary_face{p_mesh->connectivity()};
+ is_wall_boundary_face.fill(false);
for (auto&& bc_v : bc_list) {
std::visit(
[=](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
- auto node_list = bc.nodeList();
-
- is_wall_boundary_cell.fill(0);
+ auto face_list = bc.faceList();
- auto node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+ is_wall_boundary_face.fill(0);
- 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_wall_boundary_cell[cell_id] = 1;
- }
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ is_wall_boundary_face[face_list[i_face]] = 1;
}
}
},
@@ -1729,23 +1723,59 @@ class EulerKineticSolverMultiDOrder2
std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
// Reconstruction uniquement à l'intérieur donc pas utile
- // for (auto&& bc_descriptor : bc_descriptor_list) {
- // if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
- // symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
- // }
- // }
-
- PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
- symmetry_boundary_descriptor_list);
-
- auto reconstruction =
- PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho),
- DiscreteFunctionP0Vector(p_mesh, Fn_rho_u),
- DiscreteFunctionP0Vector(p_mesh, Fn_rho_E));
- auto DPk_Fn_rho = reconstruction[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+ for (auto&& bc_descriptor : bc_descriptor_list) {
+ if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
+ symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+ }
+ }
+
+ std::vector<std::shared_ptr<const IBoundaryDescriptor>> boundary_descriptor_list;
+
+ PolynomialReconstructionDescriptor reconstruction_descriptor_scalar(IntegrationMethodType::cell_center, 1,
+ symmetry_boundary_descriptor_list);
+
+ PolynomialReconstructionDescriptor reconstruction_descriptor_vector(IntegrationMethodType::cell_center, 1,
+ boundary_descriptor_list);
+
+ auto reconstruction_scalar =
+ PolynomialReconstruction{reconstruction_descriptor_scalar}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho),
+ // DiscreteFunctionP0Vector(p_mesh, Fn_rho_u),
+ DiscreteFunctionP0Vector(p_mesh, Fn_rho_E));
+
+ auto reconstruction_vector =
+ PolynomialReconstruction{reconstruction_descriptor_vector}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho_u));
+
+ auto DPk_Fn_rho = reconstruction_scalar[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+ auto DPk_Fn_rho_E = reconstruction_scalar[1]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+
auto DPk_Fn_rho_u =
- reconstruction[1]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
- auto DPk_Fn_rho_E = reconstruction[2]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+ reconstruction_vector[0]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
+
+ const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
+ const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
+
+ for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+ const auto cell_to_face = cell_to_face_matrix[cell_id];
+ if ((cell_id >= 500) and (cell_id < 510)) {
+ size_t i_face = 2; // for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+
+ FaceId face_id = cell_to_face[i_face];
+
+ const double sign = face_cell_is_reversed(cell_id, i_face) ? -1 : 1;
+ const TinyVector<Dimension> n = sign * nl[face_id];
+
+ // for (size_t i_wave = 0; i_wave < m_lambda_vector.size(); ++i_wave) {
+ size_t i_wave = 0;
+ std::cout << "n = " << n << ", x = " << m_xj[cell_id] << ", lambda_1 = " << m_lambda_vector[i_wave] << "\n";
+ std::cout << "F^rho_" << i_wave << "(" << cell_id << "," << face_id
+ << ") = " << DPk_Fn_rho(cell_id, i_wave)(m_xl[face_id]) << "; F^rho_u_" << i_wave << "(" << cell_id
+ << "," << face_id << ") = " << DPk_Fn_rho_u(cell_id, i_wave)(m_xl[face_id]) << "; F^rho_E_" << i_wave
+ << "(" << cell_id << "," << face_id << ") = " << DPk_Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) << "\n";
+ // }
+ //}
+ }
+ }
+ exit(0);
DiscreteFunctionDPkVector<Dimension, double> Fn_rho_lim = copy(DPk_Fn_rho);
DiscreteFunctionDPkVector<Dimension, double> Fn_rho_E_lim = copy(DPk_Fn_rho_E);
@@ -1762,10 +1792,10 @@ class EulerKineticSolverMultiDOrder2
CellArray<double> deltaLFn_rho_E{p_mesh->connectivity(), nb_waves};
if (m_scheme == 1) {
- FaceArray<double> Fl_rho = compute_Flux_face<double>(Fn_rho_lim, Fn_rho, is_wall_boundary_cell);
+ FaceArray<double> Fl_rho = compute_Flux_face<double>(Fn_rho_lim, Fn_rho, is_wall_boundary_face);
FaceArray<TinyVector<Dimension>> Fl_rho_u =
- compute_Flux_face<TinyVector<Dimension>>(Fn_rho_u_lim, Fn_rho_u, is_wall_boundary_cell);
- FaceArray<double> Fl_rho_E = compute_Flux_face<double>(Fn_rho_E_lim, Fn_rho_E, is_wall_boundary_cell);
+ compute_Flux_face<TinyVector<Dimension>>(Fn_rho_u_lim, Fn_rho_u, is_wall_boundary_face);
+ FaceArray<double> Fl_rho_E = compute_Flux_face<double>(Fn_rho_E_lim, Fn_rho_E, is_wall_boundary_face);
apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, Fn_rho, Fn_rho_u, Fn_rho_E,
bc_list);
@@ -1820,18 +1850,40 @@ class EulerKineticSolverMultiDOrder2
const CellArray<const double> M_rho_E = compute_scalar_M(rho_E, A_rho_E);
for (size_t i_subtimestep = 0; i_subtimestep < m_TimeOrder; i_subtimestep++) {
- auto reconstruction_np1 =
- PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0Vector(p_mesh, Fnp1_rho),
- DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_u),
- DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_E));
- auto DPk_Fnp1_rho = reconstruction[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+ auto reconstruction_np1_scalar =
+ PolynomialReconstruction{reconstruction_descriptor_scalar}.build(DiscreteFunctionP0Vector(p_mesh, Fnp1_rho),
+ DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_E));
+ auto reconstruction_np1_vector =
+ PolynomialReconstruction{reconstruction_descriptor_vector}.build(DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_u));
+
+ auto DPk_Fnp1_rho = reconstruction_scalar[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+ auto DPk_Fnp1_rho_E =
+ reconstruction_scalar[1]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+
auto DPk_Fnp1_rho_u =
- reconstruction[1]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
- auto DPk_Fnp1_rho_E = reconstruction[2]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+ reconstruction_vector[0]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_lim = copy(DPk_Fnp1_rho);
DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_E_lim = copy(DPk_Fnp1_rho_E);
+ // const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
+ // for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+ // const auto cell_to_face = cell_to_face_matrix[cell_id];
+ // if ((cell_id >= 500) and (cell_id < 510)) {
+ // for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+ // FaceId face_id = cell_to_face[i_face];
+
+ // std::cout << "F^rho_1(" << cell_id << "," << face_id << ") = " << DPk_Fnp1_rho(cell_id,
+ // 0)(m_xl[face_id]); std::cout << "F^rho_u_1(" << cell_id << "," << face_id
+ // << ") = " << DPk_Fnp1_rho_u(cell_id, 0)(m_xl[face_id]);
+
+ // std::cout << "F^rho_E_1(" << cell_id << "," << face_id
+ // << ") = " << DPk_Fnp1_rho_E(cell_id, 0)(m_xl[face_id]);
+ // }
+ // }
+ // }
+ // exit(0);
+
scalar_limiter(Fnp1_rho, Fnp1_rho_lim);
scalar_limiter(Fnp1_rho_E, Fnp1_rho_E_lim);
@@ -1844,10 +1896,10 @@ class EulerKineticSolverMultiDOrder2
CellArray<double> deltaLFnp1_rho_E{p_mesh->connectivity(), nb_waves};
if (m_scheme == 1) {
- FaceArray<double> Fl_rho_np1 = compute_Flux_face<double>(Fnp1_rho_lim, Fn_rho, is_wall_boundary_cell);
+ FaceArray<double> Fl_rho_np1 = compute_Flux_face<double>(Fnp1_rho_lim, Fn_rho, is_wall_boundary_face);
FaceArray<TinyVector<Dimension>> Fl_rho_u_np1 =
- compute_Flux_face<TinyVector<Dimension>>(Fnp1_rho_u_lim, Fn_rho_u, is_wall_boundary_cell);
- FaceArray<double> Fl_rho_E_np1 = compute_Flux_face<double>(Fnp1_rho_E_lim, Fn_rho_E, is_wall_boundary_cell);
+ compute_Flux_face<TinyVector<Dimension>>(Fnp1_rho_u_lim, Fn_rho_u, is_wall_boundary_face);
+ FaceArray<double> Fl_rho_E_np1 = compute_Flux_face<double>(Fnp1_rho_E_lim, Fn_rho_E, is_wall_boundary_face);
apply_face_bc(Fl_rho_np1, Fl_rho_u_np1, Fl_rho_E_np1, Fnp1_rho_lim, Fnp1_rho_u_lim, Fnp1_rho_E_lim, Fn_rho,
Fn_rho_u, Fn_rho_E, bc_list);
--
GitLab