From 36b237b98b52953778a1c247f74f72ed218b6962 Mon Sep 17 00:00:00 2001 From: Axelle <axelle.drouard@orange.fr> Date: Thu, 16 Jan 2025 17:48:36 +0100 Subject: [PATCH] Test new Eucclhyd scheme and bc --- src/scheme/EulerKineticSolverMultiD.cpp | 232 ++++++++++++------------ 1 file changed, 116 insertions(+), 116 deletions(-) diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp index 9aeb37209..bc2fe33ba 100644 --- a/src/scheme/EulerKineticSolverMultiD.cpp +++ b/src/scheme/EulerKineticSolverMultiD.cpp @@ -506,7 +506,7 @@ class EulerKineticSolverMultiD template <typename T> NodeArray<T> - compute_Flux1_glace(const CellArray<const T>& Fn, NodeValue<bool> is_symmetry_boundary_node) + compute_Flux1_glace(const CellArray<const T>& Fn, NodeValue<bool> is_boundary_node) { if constexpr (Dimension > 1) { const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr(); @@ -523,7 +523,7 @@ class EulerKineticSolverMultiD parallel_for( p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - if (not is_symmetry_boundary_node[node_id]) { + if (not is_boundary_node[node_id]) { const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id); const auto& node_to_cell = node_to_cell_matrix[node_id]; for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { @@ -553,13 +553,13 @@ class EulerKineticSolverMultiD } template <typename T> - FaceArrayPerNode<T> - compute_Flux1_eucclhyd(const CellArray<const T>& Fn, NodeValue<bool> is_symmetry_boundary_node) + NodeArray<T> + compute_Flux1_eucclhyd(const CellArray<const T>& Fn, NodeValue<bool> is_boundary_node) { if constexpr (Dimension > 1) { const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr(); const size_t nb_waves = m_lambda_vector.size(); - FaceArrayPerNode<T> Flr(m_mesh.connectivity(), nb_waves); + NodeArray<T> 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(); @@ -567,47 +567,49 @@ class EulerKineticSolverMultiD const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed(); if constexpr (is_tiny_vector_v<T>) { - Flr.fill(zero); + Fr.fill(zero); } else { - Flr.fill(0.); + Fr.fill(0.); } parallel_for( p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id); - const auto& node_to_face = node_to_face_matrix[node_id]; - - 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_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]; + if (not is_boundary_node[node_id]) { + const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id); + const auto& node_to_face = node_to_face_matrix[node_id]; + for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { 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 size_t i_face_cell = face_local_number_in_its_cell[i_cell]; + 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_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]; - TinyVector<Dimension> n_face = nlr(face_id, i_node_face); - const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1; + 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]; - double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face); + TinyVector<Dimension> n_face = nlr(face_id, i_node_face); + 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 (li_nlr > 0) { - Flr[node_id][i_face][i_wave] += li_nlr * Fn[cell_id][i_wave]; - sum += li_nlr; + if (li_nlr > 0) { + Fr[node_id][i_wave] += li_nlr * Fn[cell_id][i_wave]; + sum += li_nlr; + } } } - if ((sum != 0) and (not is_symmetry_boundary_node[node_id])) { - Flr[node_id][i_face][i_wave] = 1. / sum * Flr[node_id][i_face][i_wave]; + if (sum != 0) { + Fr[node_id][i_wave] = 1. / sum * Fr[node_id][i_wave]; } } } }); - return Flr; + return Fr; } else { throw NormalError("Eucclhyd not defined in 1d"); } @@ -1030,23 +1032,28 @@ class EulerKineticSolverMultiD } void - apply_eucclhyd_bc(FaceArrayPerNode<double> Flr_rho, - FaceArrayPerNode<TinyVector<Dimension>> Flr_rho_u, - FaceArrayPerNode<double> Flr_rho_E, + apply_eucclhyd_bc(NodeArray<double> Fr_rho, + NodeArray<TinyVector<Dimension>> Fr_rho_u, + NodeArray<double> Fr_rho_E, const CellValue<const double> rho, const CellValue<const TinyVector<Dimension>> rho_u, const CellValue<const double> rho_E, + const CellArray<const double> Fn_rho, + const CellArray<const TinyVector<Dimension>> Fn_rho_u, + const CellArray<const double> Fn_rho_E, const BoundaryConditionList& bc_list) { const size_t nb_waves = m_lambda_vector.size(); - // const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr(); - const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells(); - const auto& face_local_numbers_in_their_nodes = p_mesh->connectivity().faceLocalNumbersInTheirNodes(); - const auto& face_to_node_matrix = p_mesh->connectivity().faceToNodeMatrix(); - const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix(); const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr(); const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed(); + const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells(); + 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& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix(); + 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(); TinyVector<MeshType::Dimension> inv_S = zero; for (size_t d = 0; d < MeshType::Dimension; ++d) { @@ -1055,7 +1062,6 @@ class EulerKineticSolverMultiD } inv_S[d] = 1. / inv_S[d]; } - FaceArrayPerNode<double> sum_li_nlr(m_mesh.connectivity(), nb_waves); sum_li_nlr.fill(0); @@ -1063,42 +1069,64 @@ class EulerKineticSolverMultiD 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); + if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) { + auto node_list = bc.nodeList(); + 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]; + NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()}; + nr.fill(zero); - const auto& face_to_node = face_to_node_matrix[face_id]; - 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_node = 0; i_node < node_list.size(); ++i_node) { + 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); + + 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] = 1. / sqrt(dot(nr[node_id], nr[node_id])) * nr[node_id]; + auto nxn = tensorProduct(nr[node_id], nr[node_id]); + auto txt = I - nxn; for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { - for (size_t i_node = 0; i_node < face_to_node.size(); ++i_node) { - const NodeId node_id = face_to_node[i_node]; + const auto& node_to_face = node_to_face_matrix[node_id]; + const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id); + + double sum = 0; + Fr_rho[node_id][i_wave] = 0; + Fr_rho_u[node_id][i_wave] = zero; + Fr_rho_E[node_id][i_wave] = 0; - const auto& face_local_number_in_its_node = face_local_numbers_in_their_nodes.itemArray(face_id); - const size_t i_face_node = face_local_number_in_its_node[i_node]; + 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 size_t i_node_face = node_local_number_in_its_face[i_face]; + const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(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 = nlr(face_id, i_node); const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1; + TinyVector<Dimension> n_face = sign * nlr(face_id, i_node_face); - TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face; - double li_nlr_prime = sign * dot(m_lambda_vector[i_wave], nlr_prime); - double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face); + double li_nlr = dot(m_lambda_vector[i_wave], n_face); if (li_nlr > 0) { - sum_li_nlr[node_id][i_face_node][i_wave] += li_nlr; + Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_nlr; + Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u[cell_id][i_wave]; + Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_nlr; + + sum += li_nlr; } + TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face; + double li_nlr_prime = dot(m_lambda_vector[i_wave], nlr_prime); + if (li_nlr_prime > 0) { double rhoj_prime = rho[cell_id]; TinyVector<Dimension> rho_uj_prime = txt * rho_u[cell_id] - nxn * rho_u[cell_id]; @@ -1123,48 +1151,19 @@ class EulerKineticSolverMultiD } Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1]; } - Flr_rho[node_id][i_face_node][i_wave] += Fn_rho_prime * li_nlr_prime; - Flr_rho_u[node_id][i_face_node][i_wave] += li_nlr_prime * Fn_rho_u_prime; - Flr_rho_E[node_id][i_face_node][i_wave] += Fn_rho_E_prime * li_nlr_prime; - sum_li_nlr[node_id][i_face_node][i_wave] += li_nlr_prime; + Fr_rho[node_id][i_wave] += Fn_rho_prime * li_nlr_prime; + Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u_prime; + Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_nlr_prime; + + sum += li_nlr_prime; } } } - } - } - } - }, - bc_v); - } - - for (auto&& bc_v : bc_list) { - std::visit( - [=](auto&& bc) { - using BCType = std::decay_t<decltype(bc)>; - if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) { - auto face_list = bc.faceList(); - - for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { - const FaceId face_id = face_list[i_face]; - - const auto& face_to_node = face_to_node_matrix[face_id]; - - for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { - for (size_t i_node = 0; i_node < face_to_node.size(); ++i_node) { - const NodeId node_id = face_to_node[i_node]; - - const auto& face_local_number_in_its_node = face_local_numbers_in_their_nodes.itemArray(face_id); - const size_t i_face_node = face_local_number_in_its_node[i_node]; - - if (sum_li_nlr[node_id][i_face_node][i_wave] != 0) { - Flr_rho[node_id][i_face_node][i_wave] = - 1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho[node_id][i_face_node][i_wave]; - Flr_rho_u[node_id][i_face_node][i_wave] = - 1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho_u[node_id][i_face_node][i_wave]; - Flr_rho_E[node_id][i_face_node][i_wave] = - 1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho[node_id][i_face_node][i_wave]; - } + if (sum != 0) { + Fr_rho[node_id][i_wave] /= sum; + Fr_rho_u[node_id][i_wave] = 1. / sum * Fr_rho_u[node_id][i_wave]; + Fr_rho_E[node_id][i_wave] /= sum; } } } @@ -1176,7 +1175,7 @@ class EulerKineticSolverMultiD template <typename T> CellArray<T> - compute_deltaLFn_glace(NodeArray<T> Fr) + compute_deltaLFn_node(NodeArray<T> Fr) { if constexpr (Dimension > 1) { const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr(); @@ -1206,7 +1205,7 @@ class EulerKineticSolverMultiD }); return deltaLFn; } else { - throw NormalError("Glace not defined in 1d"); + throw NormalError("Nodal solver not defined in 1d"); } } @@ -1373,8 +1372,8 @@ class EulerKineticSolverMultiD bc_v); } - NodeValue<bool> is_symmetry_boundary_node{p_mesh->connectivity()}; - is_symmetry_boundary_node.fill(false); + NodeValue<bool> is_boundary_node{p_mesh->connectivity()}; + is_boundary_node.fill(false); for (auto&& bc_v : bc_list) { std::visit( [=](auto&& bc) { @@ -1382,7 +1381,12 @@ class EulerKineticSolverMultiD if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) { auto node_list = bc.nodeList(); for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { - is_symmetry_boundary_node[node_list[i_node]] = 1; + is_boundary_node[node_list[i_node]] = 1; + } + } else if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) { + auto node_list = bc.nodeList(); + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + is_boundary_node[node_list[i_node]] = 1; } } }, @@ -1411,35 +1415,31 @@ 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_symmetry_boundary_node); - NodeArray<TinyVector<Dimension>> Fr_rho_u = - compute_Flux1_glace<TinyVector<Dimension>>(Fn_rho_u, is_symmetry_boundary_node); - NodeArray<double> Fr_rho_E = compute_Flux1_glace<double>(Fn_rho_E, is_symmetry_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); - // FaceArrayPerNode<double> Flr_rho = compute_Flux1_eucclhyd<double>(Fn_rho, is_symmetry_boundary_node); - // FaceArrayPerNode<TinyVector<Dimension>> Flr_rho_u = - // compute_Flux1_eucclhyd<TinyVector<Dimension>>(Fn_rho_u, is_symmetry_boundary_node); - // FaceArrayPerNode<double> Flr_rho_E = compute_Flux1_eucclhyd<double>(Fn_rho_E, is_symmetry_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(Flr_rho, Flr_rho_u, Flr_rho_E, rho, rho_u, 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); synchronize(Fr_rho_E); - const CellArray<const double> deltaLFn_rho = compute_deltaLFn_glace(Fr_rho); - const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = compute_deltaLFn_glace(Fr_rho_u); - const CellArray<const double> deltaLFn_rho_E = compute_deltaLFn_glace(Fr_rho_E); - - // const CellArray<const double> deltaLFn_rho = compute_deltaLFn_eucclhyd(Flr_rho); - // const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = compute_deltaLFn_eucclhyd(Flr_rho_u); - // const CellArray<const double> deltaLFn_rho_E = compute_deltaLFn_eucclhyd(Flr_rho_E); + const CellArray<const double> deltaLFn_rho = compute_deltaLFn_node(Fr_rho); + const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = compute_deltaLFn_node(Fr_rho_u); + const CellArray<const double> deltaLFn_rho_E = compute_deltaLFn_node(Fr_rho_E); // const CellArray<const double> deltaLFn_rho = compute_deltaLFn_face<double>(Fl_rho); // const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = -- GitLab