From 837eca641068b4c5b9b6afe16aedec3ef698a3c0 Mon Sep 17 00:00:00 2001
From: Axelle <axelle.drouard@orange.fr>
Date: Mon, 20 Jan 2025 18:04:22 +0100
Subject: [PATCH] Fix node symmetry BC
---
src/scheme/EulerKineticSolverMultiD.cpp | 37 +++++++++++++++----------
1 file changed, 23 insertions(+), 14 deletions(-)
diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
index 791cf9a06..af1c1eb63 100644
--- a/src/scheme/EulerKineticSolverMultiD.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -842,6 +842,7 @@ class EulerKineticSolverMultiD
const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+ const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
TinyVector<MeshType::Dimension> inv_S = zero;
for (size_t d = 0; d < MeshType::Dimension; ++d) {
@@ -867,12 +868,14 @@ class EulerKineticSolverMultiD
const NodeId node_id = node_list[i_node];
const auto& node_to_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
-
+ double sum_Cjr = 0;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
const size_t i_node_cell = node_local_number_in_its_cell[i_cell];
- nr[node_id] += njr(cell_id, i_node_cell);
+ nr[node_id] += Cjr(cell_id, i_node_cell);
+
+ sum_Cjr += sqrt(dot(Cjr(cell_id, i_node_cell), Cjr(cell_id, i_node_cell)));
}
nr[node_id] = 1. / sqrt(dot(nr[node_id], nr[node_id])) * nr[node_id];
auto nxn = tensorProduct(nr[node_id], nr[node_id]);
@@ -953,11 +956,14 @@ class EulerKineticSolverMultiD
const auto& node_to_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ double sum_Cjr = 0;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
const size_t i_node_cell = node_local_number_in_its_cell[i_cell];
- nr[node_id] += njr(cell_id, i_node_cell);
+ nr[node_id] += Cjr(cell_id, i_node_cell);
+
+ sum_Cjr += sqrt(dot(Cjr(cell_id, i_node_cell), Cjr(cell_id, i_node_cell)));
}
nr[node_id] = 1. / sqrt(dot(nr[node_id], nr[node_id])) * nr[node_id];
auto nxn = tensorProduct(nr[node_id], nr[node_id]);
@@ -1054,6 +1060,7 @@ class EulerKineticSolverMultiD
const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
const auto& node_to_face_matrix = p_mesh->connectivity().nodeToFaceMatrix();
const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+ const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
TinyVector<MeshType::Dimension> inv_S = zero;
for (size_t d = 0; d < MeshType::Dimension; ++d) {
@@ -1180,11 +1187,14 @@ class EulerKineticSolverMultiD
const auto& node_to_cell = node_to_cell_matrix[node_id];
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+ double sum_Cjr = 0;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
const size_t i_node_cell = node_local_number_in_its_cell[i_cell];
- nr[node_id] += njr(cell_id, i_node_cell);
+ nr[node_id] += Cjr(cell_id, i_node_cell);
+
+ sum_Cjr += sqrt(dot(Cjr(cell_id, i_node_cell), Cjr(cell_id, i_node_cell)));
}
nr[node_id] = 1. / sqrt(dot(nr[node_id], nr[node_id])) * nr[node_id];
auto nxn = tensorProduct(nr[node_id], nr[node_id]);
@@ -1513,23 +1523,22 @@ class EulerKineticSolverMultiD
const CellValue<const TinyVector<Dimension>> rho_u = compute_PFn<TinyVector<Dimension>>(Fn_rho_u);
const CellValue<const double> rho_E = compute_PFn<double>(Fn_rho_E);
- // NodeArray<double> Fr_rho = compute_Flux1_glace<double>(Fn_rho, is_boundary_node);
- // NodeArray<TinyVector<Dimension>> Fr_rho_u = compute_Flux1_glace<TinyVector<Dimension>>(Fn_rho_u,
- // is_boundary_node); NodeArray<double> Fr_rho_E = compute_Flux1_glace<double>(Fn_rho_E,
- // is_boundary_node);
+ NodeArray<double> Fr_rho = compute_Flux1_glace<double>(Fn_rho, is_boundary_node);
+ NodeArray<TinyVector<Dimension>> Fr_rho_u = compute_Flux1_glace<TinyVector<Dimension>>(Fn_rho_u, is_boundary_node);
+ NodeArray<double> Fr_rho_E = compute_Flux1_glace<double>(Fn_rho_E, is_boundary_node);
- NodeArray<double> Fr_rho = compute_Flux1_eucclhyd<double>(Fn_rho, is_boundary_node);
- NodeArray<TinyVector<Dimension>> Fr_rho_u =
- compute_Flux1_eucclhyd<TinyVector<Dimension>>(Fn_rho_u, is_boundary_node);
- NodeArray<double> Fr_rho_E = compute_Flux1_eucclhyd<double>(Fn_rho_E, is_boundary_node);
+ // NodeArray<double> Fr_rho = compute_Flux1_eucclhyd<double>(Fn_rho, is_boundary_node);
+ // NodeArray<TinyVector<Dimension>> Fr_rho_u =
+ // compute_Flux1_eucclhyd<TinyVector<Dimension>>(Fn_rho_u, is_boundary_node);
+ // NodeArray<double> Fr_rho_E = compute_Flux1_eucclhyd<double>(Fn_rho_E, is_boundary_node);
// FaceArray<double> Fl_rho = compute_Flux1_face<double>(Fn_rho);
// FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux1_face<TinyVector<Dimension>>(Fn_rho_u);
// FaceArray<double> Fl_rho_E = compute_Flux1_face<double>(Fn_rho_E);
// apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
- // apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
- apply_eucclhyd_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
+ apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
+ // apply_eucclhyd_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
synchronize(Fr_rho);
synchronize(Fr_rho_u);
--
GitLab