From f71b0006f2a219ea7bcc98b2f66227b42ea87938 Mon Sep 17 00:00:00 2001
From: Axelle <axelle.drouard@orange.fr>
Date: Tue, 12 Nov 2024 17:18:22 +0100
Subject: [PATCH] Fix Eucchlyd scheme for kinetic solver
---
src/scheme/EulerKineticSolver2D.cpp | 69 +++++++++--------------------
1 file changed, 21 insertions(+), 48 deletions(-)
diff --git a/src/scheme/EulerKineticSolver2D.cpp b/src/scheme/EulerKineticSolver2D.cpp
index 923f6d68f..eabcc1f3a 100644
--- a/src/scheme/EulerKineticSolver2D.cpp
+++ b/src/scheme/EulerKineticSolver2D.cpp
@@ -352,9 +352,10 @@ class EulerKineticSolver2D
const size_t nb_waves = m_lambda_vector.size();
FaceArrayPerNode<double> Fr(m_mesh.connectivity(), nb_waves);
const auto& node_local_numbers_in_their_faces = p_mesh->connectivity().nodeLocalNumbersInTheirFaces();
+ const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
const auto& node_to_face_matrix = p_mesh->connectivity().nodeToFaceMatrix();
const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
- const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
parallel_for(
p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
@@ -363,39 +364,28 @@ class EulerKineticSolver2D
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
- const FaceId face_id = node_to_face[i_face];
- const auto& face_to_cell = face_to_cell_matrix[face_id];
- const unsigned int i_node_face = node_local_number_in_its_face[i_face];
- Fr[node_id][i_face][i_wave] = 0;
- double sum = 0;
+ const FaceId face_id = node_to_face[i_face];
+ 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];
+ const size_t i_node_face = node_local_number_in_its_face[i_face];
+ Fr[node_id][i_face][i_wave] = 0;
+ double sum = 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 auto& cell_to_node = cell_to_node_matrix[cell_id];
-
- auto local_node_number_in_cell = [&](NodeId node_number) {
- for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
- if (node_number == cell_to_node[i_node]) {
- return i_node;
- }
- }
- return std::numeric_limits<size_t>::max();
- };
-
- const size_t i_node_cell = local_node_number_in_cell(node_id);
+ const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
TinyVector<Dimension> n_face = nlr(face_id, i_node_face);
- TinyVector<Dimension> n_node = njr(cell_id, i_node_cell);
+ const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+
+ double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
- if (dot(n_face, n_node) < 0) {
- for (size_t d = 0; d < Dimension; ++d) {
- n_face[d] = -n_face[d];
- }
- }
- double li_nlr = dot(m_lambda_vector[i_wave], n_face);
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) {
@@ -439,15 +429,14 @@ class EulerKineticSolver2D
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& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
- const auto& face_to_node_matrix = p_mesh->connectivity().faceToNodeMatrix();
- const auto& node_to_face_matrix = p_mesh->connectivity().nodeToFaceMatrix();
+ const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
+ const auto& face_to_node_matrix = p_mesh->connectivity().faceToNodeMatrix();
+ const auto& node_to_face_matrix = p_mesh->connectivity().nodeToFaceMatrix();
+ 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];
- const auto& cell_to_node = cell_to_node_matrix[cell_id];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
deltaLFn[cell_id][i_wave] = 0;
@@ -469,28 +458,12 @@ class EulerKineticSolver2D
return std::numeric_limits<size_t>::max();
};
- auto local_node_number_in_cell = [&](NodeId node_number) {
- for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
- if (node_number == cell_to_node[i_node]) {
- return i_node;
- }
- }
- return std::numeric_limits<size_t>::max();
- };
-
const size_t i_face_node = local_face_number_in_node(face_id);
- const size_t i_node_cell = local_node_number_in_cell(node_id);
- TinyVector<Dimension> n_face = nlr(face_id, i_node_face);
- TinyVector<Dimension> n_node = njr(cell_id, i_node_cell);
+ const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
- if (dot(n_face, n_node) < 0) {
- for (size_t d = 0; d < Dimension; ++d) {
- n_face[d] = -n_face[d];
- }
- }
+ const double li_Nlr = sign * dot(m_lambda_vector[i_wave], Nlr(face_id, i_node_face));
- const double li_Nlr = dot(m_lambda_vector[i_wave], n_face);
deltaLFn[cell_id][i_wave] += Fr[node_id][i_face_node][i_wave] * li_Nlr;
}
}
--
GitLab