diff --git a/src/language/modules/KineticSchemeModule.cpp b/src/language/modules/KineticSchemeModule.cpp index a111e052203d5a1bae5b01978546d556afc42e14..bfe0d148ee68cba722cb8e0c74b66e0069efb682 100644 --- a/src/language/modules/KineticSchemeModule.cpp +++ b/src/language/modules/KineticSchemeModule.cpp @@ -564,7 +564,7 @@ KineticSchemeModule::KineticSchemeModule() std::function( [](const double& dt, const double& gamma, const std::vector<TinyVector<2>>& lambda_vector, - const double& eps, const size_t& SpaceOrder, + const double& eps, const size_t& SpaceOrder, const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_E, @@ -572,8 +572,8 @@ KineticSchemeModule::KineticSchemeModule() bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>> { - return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho, - F_rho_u, F_E, bc_descriptor_list); + return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, scheme, + F_rho, F_rho_u, F_E, bc_descriptor_list); } )); @@ -582,7 +582,7 @@ KineticSchemeModule::KineticSchemeModule() std::function( [](const double& dt, const double& gamma, const std::vector<TinyVector<3>>& lambda_vector, - const double& eps, const size_t& SpaceOrder, + const double& eps, const size_t& SpaceOrder, const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_E, @@ -590,8 +590,8 @@ KineticSchemeModule::KineticSchemeModule() bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>> { - return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho, - F_rho_u, F_E, bc_descriptor_list); + return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, scheme, + F_rho, F_rho_u, F_E, bc_descriptor_list); } )); @@ -650,7 +650,7 @@ KineticSchemeModule::KineticSchemeModule() [](const double& dt, const double& gamma, const std::vector<TinyVector<2>>& lambda_vector, const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder, - const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, + const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_E, const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& @@ -658,7 +658,7 @@ KineticSchemeModule::KineticSchemeModule() std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>> { return eulerKineticSolverMultiDOrder2(dt, gamma, lambda_vector, eps, SpaceOrder, - TimeOrder, F_rho, F_rho_u, F_E, + TimeOrder, scheme, F_rho, F_rho_u, F_E, bc_descriptor_list); } @@ -669,7 +669,7 @@ KineticSchemeModule::KineticSchemeModule() [](const double& dt, const double& gamma, const std::vector<TinyVector<3>>& lambda_vector, const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder, - const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, + const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_E, const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& @@ -677,7 +677,7 @@ KineticSchemeModule::KineticSchemeModule() std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>> { return eulerKineticSolverMultiDOrder2(dt, gamma, lambda_vector, eps, SpaceOrder, - TimeOrder, F_rho, F_rho_u, F_E, + TimeOrder, scheme, F_rho, F_rho_u, F_E, bc_descriptor_list); } diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp index af1c1eb6387f1f396d3f8bc0c2978d04f5b8bbaf..010ddd67e32e392a8b4f233db47bc92e1cfe56e0 100644 --- a/src/scheme/EulerKineticSolverMultiD.cpp +++ b/src/scheme/EulerKineticSolverMultiD.cpp @@ -228,6 +228,7 @@ class EulerKineticSolverMultiD const double m_eps; const std::vector<TinyVector<Dimension>>& m_lambda_vector; const size_t m_SpaceOrder; + const size_t m_scheme; const CellArray<const double> m_Fn_rho; const CellArray<const TinyVector<Dimension>> m_Fn_rho_u; const CellArray<const double> m_Fn_rho_E; @@ -868,14 +869,11 @@ 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] += 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]); @@ -956,14 +954,11 @@ 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] += 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]); @@ -1031,6 +1026,66 @@ class EulerKineticSolverMultiD } } } + } else if constexpr (std::is_same_v<BCType, OutflowBoundaryCondition>) { + auto node_list = bc.nodeList(); + TinyMatrix<Dimension> I = identity; + + NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()}; + nr.fill(zero); + + 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] += 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]); + auto txt = I - nxn; + + for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { + 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; + + 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]; + + double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell)); + + if (li_njr > 0) { + Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_njr; + Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u[cell_id][i_wave]; + Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_njr; + + sum += li_njr; + } + + TinyVector<Dimension> njr_prime = txt * njr(cell_id, i_node_cell) - nxn * njr(cell_id, i_node_cell); + double li_njr_prime = dot(m_lambda_vector[i_wave], njr_prime); + + if (li_njr_prime > 0) { + Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_njr_prime; + Fr_rho_u[node_id][i_wave] += li_njr_prime * Fn_rho_u[cell_id][i_wave]; + Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_njr_prime; + + sum += li_njr_prime; + } + } + 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; + } + } + } } }, bc_v); @@ -1093,7 +1148,7 @@ class EulerKineticSolverMultiD 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); } nr[node_id] = 1. / sqrt(dot(nr[node_id], nr[node_id])) * nr[node_id]; auto nxn = tensorProduct(nr[node_id], nr[node_id]); @@ -1187,14 +1242,11 @@ 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] += 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]); @@ -1275,6 +1327,80 @@ class EulerKineticSolverMultiD } } } + } else if constexpr (std::is_same_v<BCType, OutflowBoundaryCondition>) { + auto node_list = bc.nodeList(); + + TinyMatrix<Dimension> I = identity; + + NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()}; + nr.fill(zero); + + 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] += 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]); + auto txt = I - nxn; + + for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { + 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; + + 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]; + + 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); + + double li_nlr = dot(m_lambda_vector[i_wave], n_face); + + if (li_nlr > 0) { + 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) { + Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_nlr_prime; + Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u[cell_id][i_wave]; + Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_nlr_prime; + + sum += li_nlr_prime; + } + } + } + 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; + } + } + } } }, bc_v); @@ -1317,66 +1443,6 @@ class EulerKineticSolverMultiD } } - template <typename T> - CellArray<T> - compute_deltaLFn_eucclhyd(FaceArrayPerNode<T> Flr) - { - if constexpr (Dimension > 1) { - const size_t nb_waves = m_lambda_vector.size(); - CellArray<T> deltaLFn(p_mesh->connectivity(), nb_waves); - - if constexpr (is_tiny_vector_v<T>) { - deltaLFn.fill(zero); - } else { - deltaLFn.fill(0.); - } - - const NodeValuePerFace<const TinyVector<Dimension>> Nlr = MeshDataManager::instance().getMeshData(m_mesh).Nlr(); - - 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]; - - for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { - 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 auto& face_to_node = face_to_node_matrix[face_id]; - - for (size_t i_node_face = 0; i_node_face < face_to_node.size(); ++i_node_face) { - const NodeId node_id = face_to_node[i_node_face]; - const auto& node_to_face = node_to_face_matrix[node_id]; - - auto local_face_number_in_node = [&](FaceId face_number) { - for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) { - if (face_number == node_to_face[i_face]) { - return i_face; - } - } - return std::numeric_limits<size_t>::max(); - }; - - const size_t i_face_node = local_face_number_in_node(face_id); - - const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1; - - const double li_Nlr = sign * dot(m_lambda_vector[i_wave], Nlr(face_id, i_node_face)); - - deltaLFn[cell_id][i_wave] += li_Nlr * Flr[node_id][i_face_node][i_wave]; - } - } - } - }); - return deltaLFn; - } else { - throw NormalError("Eucclhyd not defined in 1d"); - } - } - template <typename T> CellArray<T> compute_deltaLFn_face(FaceArray<T> Fl) @@ -1523,35 +1589,56 @@ 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); + CellArray<double> deltaLFn_rho{p_mesh->connectivity(), nb_waves}; + CellArray<TinyVector<Dimension>> deltaLFn_rho_u{p_mesh->connectivity(), nb_waves}; + CellArray<double> deltaLFn_rho_E{p_mesh->connectivity(), nb_waves}; + + if (m_scheme == 1) { + 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); - // 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); + 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); - // 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); + synchronize(Fl_rho); + synchronize(Fl_rho_u); + synchronize(Fl_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); + deltaLFn_rho = compute_deltaLFn_face<double>(Fl_rho); + deltaLFn_rho_u = compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u); + deltaLFn_rho_E = compute_deltaLFn_face<double>(Fl_rho_E); - synchronize(Fr_rho); - synchronize(Fr_rho_u); - synchronize(Fr_rho_E); + } else if (m_scheme == 2) { + 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); - 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); + 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); + synchronize(Fr_rho); + synchronize(Fr_rho_u); + synchronize(Fr_rho_E); - // const CellArray<const double> deltaLFn_rho = compute_deltaLFn_face<double>(Fl_rho); - // const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = - // compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u); - // const CellArray<const double> deltaLFn_rho_E = compute_deltaLFn_face<double>(Fl_rho_E); + deltaLFn_rho = compute_deltaLFn_node(Fr_rho); + deltaLFn_rho_u = compute_deltaLFn_node(Fr_rho_u); + deltaLFn_rho_E = compute_deltaLFn_node(Fr_rho_E); + + } else if (m_scheme == 3) { + 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); + + 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); + + deltaLFn_rho = compute_deltaLFn_node(Fr_rho); + deltaLFn_rho_u = compute_deltaLFn_node(Fr_rho_u); + deltaLFn_rho_E = compute_deltaLFn_node(Fr_rho_E); + } const CellValue<const TinyVector<Dimension>> A_rho = getA_rho(rho_u); const CellValue<const TinyMatrix<Dimension>> A_rho_u = getA_rho_u(rho, rho_u, rho_E); @@ -1609,6 +1696,7 @@ class EulerKineticSolverMultiD const double& gamma, const std::vector<TinyVector<MeshType::Dimension>>& lambda_vector, const size_t& SpaceOrder, + const size_t& scheme, const CellArray<const double>& Fn_rho, const CellArray<const TinyVector<MeshType::Dimension>>& Fn_rho_u, const CellArray<const double>& Fn_rho_E) @@ -1616,6 +1704,7 @@ class EulerKineticSolverMultiD m_eps{eps}, m_lambda_vector{lambda_vector}, m_SpaceOrder{SpaceOrder}, + m_scheme{scheme}, m_Fn_rho{Fn_rho}, m_Fn_rho_u{Fn_rho_u}, m_Fn_rho_E{Fn_rho_E}, @@ -1834,6 +1923,7 @@ eulerKineticSolverMultiD(const double& dt, const std::vector<TinyVector<Dimension>>& lambda_vector, const double& eps, const size_t& SpaceOrder, + const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E, @@ -1857,7 +1947,7 @@ eulerKineticSolverMultiD(const double& dt, auto Fn_rho_u = F_rho_u->get<DiscreteFunctionP0Vector<const TinyVector<MeshType::Dimension>>>(); auto Fn_rho_E = F_rho_E->get<DiscreteFunctionP0Vector<const double>>(); - EulerKineticSolverMultiD solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, Fn_rho.cellArrays(), + EulerKineticSolverMultiD solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, scheme, Fn_rho.cellArrays(), Fn_rho_u.cellArrays(), Fn_rho_E.cellArrays()); return solver.apply(bc_descriptor_list); @@ -1876,6 +1966,7 @@ eulerKineticSolverMultiD(const double& dt, const std::vector<TinyVector<2>>& lambda_vector, const double& eps, const size_t& SpaceOrder, + const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E, @@ -1889,6 +1980,7 @@ eulerKineticSolverMultiD(const double& dt, const std::vector<TinyVector<3>>& lambda_vector, const double& eps, const size_t& SpaceOrder, + const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E, diff --git a/src/scheme/EulerKineticSolverMultiD.hpp b/src/scheme/EulerKineticSolverMultiD.hpp index 987f10f8a97189ca0db782cf0545e92b0d5bb933..4053facd06336991f11a3d8addfcbdf86e49e9e6 100644 --- a/src/scheme/EulerKineticSolverMultiD.hpp +++ b/src/scheme/EulerKineticSolverMultiD.hpp @@ -32,6 +32,7 @@ eulerKineticSolverMultiD(const double& dt, const std::vector<TinyVector<Dimension>>& lambda_vector, const double& eps, const size_t& SpaceOrder, + const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E, diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.cpp b/src/scheme/EulerKineticSolverMultiDOrder2.cpp index 82776378ca7350104635c83f47621f00d48cebf4..02c995fdd75728747366ba4c1afb2563af28340d 100644 --- a/src/scheme/EulerKineticSolverMultiDOrder2.cpp +++ b/src/scheme/EulerKineticSolverMultiDOrder2.cpp @@ -50,6 +50,7 @@ class EulerKineticSolverMultiDOrder2 const std::vector<TinyVector<Dimension>>& m_lambda_vector; const size_t m_SpaceOrder; const size_t m_TimeOrder; + const size_t m_scheme; const CellArray<const double> m_Fn_rho; const CellArray<const TinyVector<Dimension>> m_Fn_rho_u; const CellArray<const double> m_Fn_rho_E; @@ -451,7 +452,7 @@ class EulerKineticSolverMultiDOrder2 template <typename T> FaceArray<T> - compute_Flux1_face(const DiscreteFunctionDPkVector<Dimension, const T> DPk_Fn) + compute_Flux_face(const DiscreteFunctionDPkVector<Dimension, const T> DPk_Fn) { const size_t nb_waves = m_lambda_vector.size(); FaceArray<T> Fl(m_mesh.connectivity(), nb_waves); @@ -490,17 +491,125 @@ class EulerKineticSolverMultiDOrder2 return Fl; } + template <typename T> + NodeArray<T> + compute_Flux_glace(const DiscreteFunctionDPkVector<Dimension, const T> DPk_Fn, NodeValue<bool> is_boundary_node) + { + if constexpr (Dimension > 1) { + const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr(); + const size_t nb_waves = m_lambda_vector.size(); + NodeArray<T> Fr(m_mesh.connectivity(), nb_waves); + const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells(); + const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix(); + + if constexpr (is_tiny_vector_v<T>) { + Fr.fill(zero); + } else { + Fr.fill(0.); + } + + parallel_for( + p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId 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) { + double sum_li_njr = 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 unsigned int i_node = node_local_number_in_its_cell[i_cell]; + + double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node)); + if (li_njr > 0) { + Fr[node_id][i_wave] += li_njr * DPk_Fn(cell_id, i_wave)(m_xr[node_id]); + sum_li_njr += li_njr; + } + } + if (sum_li_njr != 0) { + Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave]; + } + } + } + }); + + return Fr; + } else { + throw NormalError("Glace not defined in 1d"); + } + } + + template <typename T> + NodeArray<T> + compute_Flux_eucclhyd(const DiscreteFunctionDPkVector<Dimension, const T> DPk_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(); + 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(); + const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix(); + const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed(); + + if constexpr (is_tiny_vector_v<T>) { + Fr.fill(zero); + } else { + Fr.fill(0.); + } + + parallel_for( + p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { + 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_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]; + + 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_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) { + Fr[node_id][i_wave] += li_nlr * DPk_Fn(cell_id, i_wave)(m_xl[face_id]); + sum += li_nlr; + } + } + } + if (sum != 0) { + Fr[node_id][i_wave] = 1. / sum * Fr[node_id][i_wave]; + } + } + } + }); + + return Fr; + } else { + throw NormalError("Eucclhyd not defined in 1d"); + } + } + void - apply_bc(FaceArray<double> Fl_rho, - FaceArray<TinyVector<Dimension>> Fl_rho_u, - FaceArray<double> Fl_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) + apply_face_bc(FaceArray<double> Fl_rho, + FaceArray<TinyVector<Dimension>> Fl_rho_u, + FaceArray<double> Fl_rho_E, + DiscreteFunctionDPkVector<Dimension, double> Fn_rho, + DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fn_rho_u, + DiscreteFunctionDPkVector<Dimension, double> Fn_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(); @@ -537,27 +646,36 @@ class EulerKineticSolverMultiDOrder2 const CellId cell_id = face_to_cell[i_cell]; const size_t i_face_cell = face_local_number_in_its_cell[i_cell]; + double rhol = 0; + TinyVector<Dimension> rho_ul = zero; + double rho_El = 0; + + 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_El += Fn_rho_E(cell_id, k_wave)(m_xl[face_id]); + } + 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_u[cell_id]; - TinyVector<Dimension> rho_uj_prime = txt * rho_uj - nxn * rho_uj; - double rho_Ej_prime = rho_E[cell_id]; + double rhol_prime = rhol; + TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul; + double rho_El_prime = rho_El; - double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]); - double p = (m_gamma - 1) * rho_e; + double rho_el = rho_El - 0.5 * (1. / rhol) * dot(rho_ul, rho_ul); + double pl = (m_gamma - 1) * rho_el; - TinyVector<Dimension> A_rho_prime = rho_uj_prime; + TinyVector<Dimension> A_rho_prime = rho_ul_prime; TinyMatrix<Dimension> A_rho_u_prime = - 1. / rhoj_prime * tensorProduct(rho_uj_prime, rho_uj_prime) + p * I; - TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime; + 1. / rhol_prime * tensorProduct(rho_ul_prime, rho_ul_prime) + pl * I; + TinyVector<Dimension> A_rho_E_prime = (rho_El_prime + pl) / rhol_prime * rho_ul_prime; - double Fn_rho_prime = rhoj_prime / nb_waves; - TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_uj_prime; - double Fn_rho_E_prime = rho_Ej_prime / nb_waves; + double Fn_rho_prime = rhol_prime / nb_waves; + TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ul_prime; + double Fn_rho_E_prime = rho_El_prime / nb_waves; for (size_t d1 = 0; d1 < Dimension; ++d1) { Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1]; @@ -587,8 +705,18 @@ class EulerKineticSolverMultiDOrder2 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]; - const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1; + double rhol = 0; + TinyVector<Dimension> rho_ul = zero; + double rho_El = 0; + + 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_El += Fn_rho_E(cell_id, k_wave)(m_xl[face_id]); + } + + const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1; auto n = sign * nl[face_id]; auto nxn = tensorProduct(n, n); TinyMatrix<Dimension> I = identity; @@ -596,22 +724,21 @@ class EulerKineticSolverMultiDOrder2 double li_nl = dot(m_lambda_vector[i_wave], n); if (li_nl < 0) { - double rhoj_prime = rho[cell_id]; - TinyVector<Dimension> rho_uj = rho_u[cell_id]; - TinyVector<Dimension> rho_uj_prime = txt * rho_uj - nxn * rho_uj; - double rho_Ej_prime = rho_E[cell_id]; + double rhol_prime = rhol; + TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul; + double rho_El_prime = rho_El; - double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]); - double p = (m_gamma - 1) * rho_e; + double rho_el = rho_El - 0.5 * (1. / rhol) * dot(rho_ul, rho_ul); + double pl = (m_gamma - 1) * rho_el; - TinyVector<Dimension> A_rho_prime = rho_uj_prime; + TinyVector<Dimension> A_rho_prime = rho_ul_prime; TinyMatrix<Dimension> A_rho_u_prime = - 1. / rhoj_prime * tensorProduct(rho_uj_prime, rho_uj_prime) + p * I; - TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime; + 1. / rhol_prime * tensorProduct(rho_ul_prime, rho_ul_prime) + pl * I; + TinyVector<Dimension> A_rho_E_prime = (rho_El_prime + pl) / rhol_prime * rho_ul_prime; - double Fn_rho_prime = rhoj_prime / nb_waves; - TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_uj_prime; - double Fn_rho_E_prime = rho_Ej_prime / nb_waves; + double Fn_rho_prime = rhol_prime / nb_waves; + TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ul_prime; + double Fn_rho_E_prime = rho_El_prime / nb_waves; for (size_t d1 = 0; d1 < Dimension; ++d1) { Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1]; @@ -647,9 +774,9 @@ class EulerKineticSolverMultiDOrder2 double li_nl = dot(m_lambda_vector[i_wave], n); if (li_nl < 0) { - Fl_rho[face_id][i_wave] += Fn_rho[cell_id][i_wave]; - Fl_rho_u[face_id][i_wave] += Fn_rho_u[cell_id][i_wave]; - Fl_rho_E[face_id][i_wave] += Fn_rho_E[cell_id][i_wave]; + Fl_rho[face_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]); + Fl_rho_u[face_id][i_wave] += Fn_rho_u(cell_id, i_wave)(m_xl[face_id]); + Fl_rho_E[face_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]); } } } @@ -660,88 +787,654 @@ class EulerKineticSolverMultiDOrder2 } } - DiscreteFunctionP0Vector<double> - compute_deltaLFn(NodeArray<double> Fr) + void + apply_glace_bc(NodeArray<double> Fr_rho, + NodeArray<TinyVector<Dimension>> Fr_rho_u, + NodeArray<double> Fr_rho_E, + DiscreteFunctionDPkVector<Dimension, double> Fn_rho, + DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fn_rho_u, + DiscreteFunctionDPkVector<Dimension, double> Fn_rho_E, + const BoundaryConditionList& bc_list) { - if constexpr (Dimension > 1) { - const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr(); - const size_t nb_waves = m_lambda_vector.size(); - DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves); + const size_t nb_waves = m_lambda_vector.size(); + 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(); - const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix(); + TinyVector<MeshType::Dimension> inv_S = zero; + for (size_t d = 0; d < MeshType::Dimension; ++d) { + for (size_t i = 0; i < nb_waves; ++i) { + inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d]; + } + inv_S[d] = 1. / inv_S[d]; + } - parallel_for( - p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { - const auto& cell_to_node = cell_to_node_matrix[cell_id]; + 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 node_list = bc.nodeList(); - for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { - deltaLFn[cell_id][i_wave] = 0; - for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { - const NodeId node_id = cell_to_node[i_node]; + TinyMatrix<Dimension> I = identity; - double li_Cjr = dot(m_lambda_vector[i_wave], Cjr(cell_id, i_node)); - deltaLFn[cell_id][i_wave] += Fr[node_id][i_wave] * li_Cjr; + NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()}; + nr.fill(zero); + + 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] += 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]); + auto txt = I - nxn; + + for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { + 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; + + 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]; + + double rhor = 0; + TinyVector<Dimension> rho_ur = zero; + double rho_Er = 0; + + for (size_t k_wave = 0; k_wave < nb_waves; ++k_wave) { + rhor += Fn_rho(cell_id, k_wave)(m_xr[node_id]); + rho_ur += Fn_rho_u(cell_id, k_wave)(m_xr[node_id]); + rho_Er += Fn_rho_E(cell_id, k_wave)(m_xr[node_id]); + } + + double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell)); + + if (li_njr > 0) { + Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xr[node_id]) * li_njr; + Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u(cell_id, i_wave)(m_xr[node_id]); + Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xr[node_id]) * li_njr; + + sum += li_njr; + } + + TinyVector<Dimension> njr_prime = txt * njr(cell_id, i_node_cell) - nxn * njr(cell_id, i_node_cell); + double li_njr_prime = dot(m_lambda_vector[i_wave], njr_prime); + + if (li_njr_prime > 0) { + double rhor_prime = rhor; + TinyVector<Dimension> rho_ur_prime = txt * rho_ur - nxn * rho_ur; + double rho_Er_prime = rho_Er; + + double rho_er = rho_Er - 0.5 * (1. / rhor) * dot(rho_ur, rho_ur); + double pr = (m_gamma - 1) * rho_er; + + TinyVector<Dimension> A_rho_prime = rho_ur_prime; + TinyMatrix<Dimension> A_rho_u_prime = + 1. / rhor_prime * tensorProduct(rho_ur_prime, rho_ur_prime) + pr * I; + TinyVector<Dimension> A_rho_E_prime = (rho_Er_prime + pr) / rhor_prime * rho_ur_prime; + + double Fn_rho_prime = rhor_prime / nb_waves; + TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ur_prime; + double Fn_rho_E_prime = rho_Er_prime / nb_waves; + + for (size_t d1 = 0; d1 < Dimension; ++d1) { + Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1]; + for (size_t d2 = 0; d2 < Dimension; ++d2) { + Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1); + } + Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1]; + } + + Fr_rho[node_id][i_wave] += Fn_rho_prime * li_njr_prime; + Fr_rho_u[node_id][i_wave] += li_njr_prime * Fn_rho_u_prime; + Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_njr_prime; + + sum += li_njr_prime; + } + } + 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; + } + } + } + } else if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) { + auto node_list = bc.nodeList(); + + TinyMatrix<Dimension> I = identity; + + NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()}; + nr.fill(zero); + + 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] += 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]); + auto txt = I - nxn; + + for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { + 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; + + 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]; + + double rhor = 0; + TinyVector<Dimension> rho_ur = zero; + double rho_Er = 0; + + for (size_t k_wave = 0; k_wave < nb_waves; ++k_wave) { + rhor += Fn_rho(cell_id, k_wave)(m_xr[node_id]); + rho_ur += Fn_rho_u(cell_id, k_wave)(m_xr[node_id]); + rho_Er += Fn_rho_E(cell_id, k_wave)(m_xr[node_id]); + } + + double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell)); + + if (li_njr > 0) { + Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xr[node_id]) * li_njr; + Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u(cell_id, i_wave)(m_xr[node_id]); + Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xr[node_id]) * li_njr; + + sum += li_njr; + } + + TinyVector<Dimension> njr_prime = txt * njr(cell_id, i_node_cell) - nxn * njr(cell_id, i_node_cell); + double li_njr_prime = dot(m_lambda_vector[i_wave], njr_prime); + + if (li_njr_prime > 0) { + double rhor_prime = rhor; + TinyVector<Dimension> rho_ur_prime = txt * rho_ur - nxn * rho_ur; + double rho_Er_prime = rho_Er; + + double rho_er = rho_Er - 0.5 * (1. / rhor) * dot(rho_ur, rho_ur); + double pr = (m_gamma - 1) * rho_er; + + TinyVector<Dimension> A_rho_prime = rho_ur_prime; + TinyMatrix<Dimension> A_rho_u_prime = + 1. / rhor_prime * tensorProduct(rho_ur_prime, rho_ur_prime) + pr * I; + TinyVector<Dimension> A_rho_E_prime = (rho_Er_prime + pr) / rhor_prime * rho_ur_prime; + + double Fn_rho_prime = rhor_prime / nb_waves; + TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ur_prime; + double Fn_rho_E_prime = rho_Er_prime / nb_waves; + + for (size_t d1 = 0; d1 < Dimension; ++d1) { + Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1]; + for (size_t d2 = 0; d2 < Dimension; ++d2) { + Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1); + } + Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1]; + } + + Fr_rho[node_id][i_wave] += Fn_rho_prime * li_njr_prime; + Fr_rho_u[node_id][i_wave] += li_njr_prime * Fn_rho_u_prime; + Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_njr_prime; + + sum += li_njr_prime; + } + } + 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; + } + } + } + } else if constexpr (std::is_same_v<BCType, OutflowBoundaryCondition>) { + auto node_list = bc.nodeList(); + TinyMatrix<Dimension> I = identity; + + NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()}; + nr.fill(zero); + + 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] += 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]); + auto txt = I - nxn; + + for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { + 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; + + 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]; + + double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell)); + + if (li_njr > 0) { + Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xr[node_id]) * li_njr; + Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u(cell_id, i_wave)(m_xr[node_id]); + Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xr[node_id]) * li_njr; + + sum += li_njr; + } + + TinyVector<Dimension> njr_prime = txt * njr(cell_id, i_node_cell) - nxn * njr(cell_id, i_node_cell); + double li_njr_prime = dot(m_lambda_vector[i_wave], njr_prime); + + if (li_njr_prime > 0) { + Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xr[node_id]) * li_njr_prime; + Fr_rho_u[node_id][i_wave] += li_njr_prime * Fn_rho_u(cell_id, i_wave)(m_xr[node_id]); + Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xr[node_id]) * li_njr_prime; + + sum += li_njr_prime; + } + } + 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; + } + } } } - }); - return deltaLFn; - } else { - throw NormalError("Glace not defined in 1d"); + }, + bc_v); } } - DiscreteFunctionP0Vector<double> - compute_deltaLFn_eucclhyd(FaceArrayPerNode<double> Fr) + void + apply_eucclhyd_bc(NodeArray<double> Fr_rho, + NodeArray<TinyVector<Dimension>> Fr_rho_u, + NodeArray<double> Fr_rho_E, + DiscreteFunctionDPkVector<Dimension, double> Fn_rho, + DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fn_rho_u, + DiscreteFunctionDPkVector<Dimension, double> Fn_rho_E, + const BoundaryConditionList& bc_list) { - if constexpr (Dimension > 1) { - const size_t nb_waves = m_lambda_vector.size(); - DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves); + const size_t nb_waves = m_lambda_vector.size(); - 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 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(); + const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr(); - 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(); + TinyVector<MeshType::Dimension> inv_S = zero; + for (size_t d = 0; d < MeshType::Dimension; ++d) { + for (size_t i = 0; i < nb_waves; ++i) { + inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d]; + } + inv_S[d] = 1. / inv_S[d]; + } + FaceArrayPerNode<double> sum_li_nlr(m_mesh.connectivity(), nb_waves); + sum_li_nlr.fill(0); - parallel_for( - p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { - const auto& cell_to_face = cell_to_face_matrix[cell_id]; + 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, WallBoundaryCondition>) { + auto node_list = bc.nodeList(); - for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { - deltaLFn[cell_id][i_wave] = 0; + TinyMatrix<Dimension> I = identity; - 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 auto& face_to_node = face_to_node_matrix[face_id]; + NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()}; + nr.fill(zero); + + 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] += 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]); + auto txt = I - nxn; + + for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { + 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; + + 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]; + + double rhol = 0; + TinyVector<Dimension> rho_ul = zero; + double rho_El = 0; + + 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_El += Fn_rho_E(cell_id, k_wave)(m_xl[face_id]); + } + + 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); + + double li_nlr = dot(m_lambda_vector[i_wave], n_face); + + if (li_nlr > 0) { + Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr; + Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]); + Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * li_nlr; + + sum += li_nlr; + } - for (size_t i_node_face = 0; i_node_face < face_to_node.size(); ++i_node_face) { - const NodeId node_id = face_to_node[i_node_face]; - const auto& node_to_face = node_to_face_matrix[node_id]; + TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face; + double li_nlr_prime = dot(m_lambda_vector[i_wave], nlr_prime); - auto local_face_number_in_node = [&](FaceId face_number) { - for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) { - if (face_number == node_to_face[i_face]) { - return i_face; + if (li_nlr_prime > 0) { + double rhol_prime = rhol; + TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul; + double rho_El_prime = rho_El; + + double rho_el = rho_El - 0.5 * (1. / rhol) * dot(rho_ul, rho_ul); + double pl = (m_gamma - 1) * rho_el; + + TinyVector<Dimension> A_rho_prime = rho_ul_prime; + TinyMatrix<Dimension> A_rho_u_prime = + 1. / rhol_prime * tensorProduct(rho_ul_prime, rho_ul_prime) + pl * I; + TinyVector<Dimension> A_rho_E_prime = (rho_El_prime + pl) / rhol_prime * rho_ul_prime; + + double Fn_rho_prime = rhol_prime / nb_waves; + TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ul_prime; + double Fn_rho_E_prime = rho_El_prime / nb_waves; + + for (size_t d1 = 0; d1 < Dimension; ++d1) { + Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1]; + for (size_t d2 = 0; d2 < Dimension; ++d2) { + Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1); + } + Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1]; + } + + 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; + } + } + } + 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; + } + } + } + } else if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) { + auto node_list = bc.nodeList(); + + TinyMatrix<Dimension> I = identity; + + NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()}; + nr.fill(zero); + + 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] += 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]); + auto txt = I - nxn; + + for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { + 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; + + 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]; + + double rhol = 0; + TinyVector<Dimension> rho_ul = zero; + double rho_El = 0; + + 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_El += Fn_rho_E(cell_id, k_wave)(m_xl[face_id]); + } + + 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); + + double li_nlr = dot(m_lambda_vector[i_wave], n_face); + + if (li_nlr > 0) { + Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr; + Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]); + Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * 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 rhol_prime = rhol; + TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul; + double rho_El_prime = rho_El; + + double rho_el = rho_El - 0.5 * (1. / rhol) * dot(rho_ul, rho_ul); + double pl = (m_gamma - 1) * rho_el; + + TinyVector<Dimension> A_rho_prime = rho_ul_prime; + TinyMatrix<Dimension> A_rho_u_prime = + 1. / rhol_prime * tensorProduct(rho_ul_prime, rho_ul_prime) + pl * I; + TinyVector<Dimension> A_rho_E_prime = (rho_El_prime + pl) / rhol_prime * rho_ul_prime; + + double Fn_rho_prime = rhol_prime / nb_waves; + TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ul_prime; + double Fn_rho_E_prime = rho_El_prime / nb_waves; + + for (size_t d1 = 0; d1 < Dimension; ++d1) { + Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1]; + for (size_t d2 = 0; d2 < Dimension; ++d2) { + Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1); + } + Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1]; + } + + 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; } } - return std::numeric_limits<size_t>::max(); - }; + } + 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; + } + } + } + } else if constexpr (std::is_same_v<BCType, OutflowBoundaryCondition>) { + auto node_list = bc.nodeList(); + + TinyMatrix<Dimension> I = identity; + + NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()}; + nr.fill(zero); + + 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] += 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]); + auto txt = I - nxn; + + for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) { + 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; + + 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]; - const size_t i_face_node = local_face_number_in_node(face_id); + 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); - const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1; + double li_nlr = dot(m_lambda_vector[i_wave], n_face); - const double li_Nlr = sign * dot(m_lambda_vector[i_wave], Nlr(face_id, i_node_face)); + if (li_nlr > 0) { + Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr; + Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]); + Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * li_nlr; - deltaLFn[cell_id][i_wave] += Fr[node_id][i_face_node][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) { + Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr_prime; + Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]); + Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * li_nlr_prime; + + sum += li_nlr_prime; + } + } + } + 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; + } } } } + }, + bc_v); + } + } + + template <typename T> + CellArray<T> + compute_deltaLFn_node(NodeArray<T> Fr) + { + if constexpr (Dimension > 1) { + const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr(); + const size_t nb_waves = m_lambda_vector.size(); + CellArray<T> deltaLFn(p_mesh->connectivity(), nb_waves); + + const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix(); + + if constexpr (is_tiny_vector_v<T>) { + deltaLFn.fill(zero); + } else { + deltaLFn.fill(0.); + } + + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId 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) { + for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) { + const NodeId node_id = cell_to_node[i_node]; + + double li_Cjr = dot(m_lambda_vector[i_wave], Cjr(cell_id, i_node)); + deltaLFn[cell_id][i_wave] += li_Cjr * Fr[node_id][i_wave]; + } + } }); return deltaLFn; } else { - throw NormalError("Eucclhyd not defined in 1d"); + throw NormalError("Nodal solver not defined in 1d"); } } @@ -833,30 +1526,23 @@ class EulerKineticSolverMultiDOrder2 f_max = std::max(f_max, f[cell_stencil[i_cell]][i_wave]); } - double f_bar_min = fj; - double f_bar_max = fj; - - for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) { - const CellId cell_k_id = cell_stencil[i_cell]; - const double f_xk = DPk_fh(cell_id, i_wave)(m_xj[cell_k_id]); - - f_bar_min = std::min(f_bar_min, f_xk); - f_bar_max = std::max(f_bar_max, f_xk); - } - const double eps = 1E-14; - double coef1 = 1; - if (std::abs(f_bar_max - fj) > eps) { - coef1 = (f_max - fj) / ((f_bar_max - fj)); - } + double coef = 1; + double lambda = 1.; - double coef2 = 1.; - if (std::abs(f_bar_min - fj) > eps) { - coef2 = (f_min - fj) / ((f_bar_min - fj)); + for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) { + const CellId cell_i_id = cell_stencil[i_cell]; + const double f_bar_i = DPk_fh(cell_id, i_wave)(m_xj[cell_i_id]); + const double Delta = (f_bar_i - fj); + + if (Delta > eps) { + coef = std::min(1., (f_max - fj) / Delta); + } else if (Delta < -eps) { + coef = std::min(1., (f_min - fj) / Delta); + } + lambda = std::min(lambda, coef); } - const double lambda = std::max(0., std::min(1., std::min(coef1, coef2))); - auto coefficients = DPk_fh.componentCoefficients(cell_id, i_wave); coefficients[0] = (1 - lambda) * f[cell_id][i_wave] + lambda * coefficients[0]; for (size_t i = 1; i < coefficients.size(); ++i) { @@ -892,45 +1578,26 @@ class EulerKineticSolverMultiDOrder2 } } - TinyVector<Dimension> f_bar_min = fj; - TinyVector<Dimension> f_bar_max = fj; + const double eps = 1E-14; + TinyVector<Dimension> coef; + double lambda = 1.; for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) { - const CellId cell_k_id = cell_stencil[i_cell]; - const TinyVector<Dimension> f_xk = DPk_fh(cell_id, i_wave)(m_xj[cell_k_id]); + const CellId cell_i_id = cell_stencil[i_cell]; + const TinyVector<Dimension> f_bar_i = DPk_fh(cell_id, i_wave)(m_xj[cell_i_id]); + TinyVector<Dimension> Delta; for (size_t dim = 0; dim < Dimension; ++dim) { - f_bar_min[dim] = std::min(f_bar_min[dim], f_xk[dim]); - f_bar_max[dim] = std::max(f_bar_max[dim], f_xk[dim]); - } - } - - const double eps = 1E-14; - TinyVector<Dimension> coef1; - for (size_t dim = 0; dim < Dimension; ++dim) { - coef1[dim] = 1; - if (std::abs(f_bar_max[dim] - fj[dim]) > eps) { - coef1[dim] = (f_max[dim] - fj[dim]) / ((f_bar_max[dim] - fj[dim])); - } - } - - TinyVector<Dimension> coef2; - for (size_t dim = 0; dim < Dimension; ++dim) { - coef2[dim] = 1; - if (std::abs(f_bar_min[dim] - fj[dim]) > eps) { - coef2[dim] = (f_min[dim] - fj[dim]) / ((f_bar_min[dim] - fj[dim])); + Delta[dim] = (f_bar_i[dim] - fj[dim]); + coef[dim] = 1.; + if (Delta[dim] > eps) { + coef[dim] = std::min(1., (f_max[dim] - fj[dim]) / Delta[dim]); + } else if (Delta[dim] < -eps) { + coef[dim] = std::min(1., (f_min[dim] - fj[dim]) / Delta[dim]); + } + lambda = std::min(lambda, coef[dim]); } } - - double min_coef1 = coef1[0]; - double min_coef2 = coef2[0]; - for (size_t dim = 1; dim < Dimension; ++dim) { - min_coef1 = std::min(min_coef1, coef1[dim]); - min_coef2 = std::min(min_coef2, coef2[dim]); - } - - const double lambda = std::max(0., std::min(1., std::min(min_coef1, min_coef2))); - auto coefficients = DPk_fh.componentCoefficients(cell_id, i_wave); coefficients[0] = (1 - lambda) * f[cell_id][i_wave] + lambda * coefficients[0]; @@ -979,6 +1646,27 @@ class EulerKineticSolverMultiDOrder2 bc_v); } + NodeValue<bool> is_boundary_node{p_mesh->connectivity()}; + is_boundary_node.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, SymmetryBoundaryCondition>) { + 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; + } + } 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; + } + } + }, + bc_v); + } + const CellValue<const TinyVector<Dimension>> A_rho_ex = getA_rho(rho_u_ex); const CellValue<const TinyMatrix<Dimension>> A_rho_u_ex = getA_rho_u(rho_ex, rho_u_ex, rho_E_ex); const CellValue<const TinyVector<Dimension>> A_rho_E_ex = getA_rho_E(rho_ex, rho_u_ex, rho_E_ex); @@ -1035,20 +1723,57 @@ class EulerKineticSolverMultiDOrder2 vector_limiter(Fn_rho_u, Fn_rho_u_lim); - FaceArray<double> Fl_rho = compute_Flux1_face<double>(Fn_rho_lim); - FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux1_face<TinyVector<Dimension>>(Fn_rho_u_lim); - FaceArray<double> Fl_rho_E = compute_Flux1_face<double>(Fn_rho_E_lim); + CellArray<double> deltaLFn_rho{p_mesh->connectivity(), nb_waves}; + CellArray<TinyVector<Dimension>> deltaLFn_rho_u{p_mesh->connectivity(), nb_waves}; + 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); + FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux_face<TinyVector<Dimension>>(Fn_rho_u_lim); + FaceArray<double> Fl_rho_E = compute_Flux_face<double>(Fn_rho_E_lim); + + apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, bc_list); - apply_bc(Fl_rho, Fl_rho_u, Fl_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list); + synchronize(Fl_rho); + synchronize(Fl_rho_u); + synchronize(Fl_rho_E); - synchronize(Fl_rho); - synchronize(Fl_rho_u); - synchronize(Fl_rho_E); + deltaLFn_rho = compute_deltaLFn_face<double>(Fl_rho); + deltaLFn_rho_u = compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u); + deltaLFn_rho_E = compute_deltaLFn_face<double>(Fl_rho_E); - const CellArray<const double> deltaLFn_rho = compute_deltaLFn_face<double>(Fl_rho); - const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = - compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u); - const CellArray<const double> deltaLFn_rho_E = compute_deltaLFn_face<double>(Fl_rho_E); + } else if (m_scheme == 2) { + NodeArray<double> Fr_rho = compute_Flux_glace<double>(Fn_rho_lim, is_boundary_node); + NodeArray<TinyVector<Dimension>> Fr_rho_u = + compute_Flux_glace<TinyVector<Dimension>>(Fn_rho_u_lim, is_boundary_node); + NodeArray<double> Fr_rho_E = compute_Flux_glace<double>(Fn_rho_E_lim, is_boundary_node); + + apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, bc_list); + + synchronize(Fr_rho); + synchronize(Fr_rho_u); + synchronize(Fr_rho_E); + + deltaLFn_rho = compute_deltaLFn_node(Fr_rho); + deltaLFn_rho_u = compute_deltaLFn_node(Fr_rho_u); + deltaLFn_rho_E = compute_deltaLFn_node(Fr_rho_E); + + } else if (m_scheme == 3) { + NodeArray<double> Fr_rho = compute_Flux_eucclhyd<double>(Fn_rho_lim, is_boundary_node); + NodeArray<TinyVector<Dimension>> Fr_rho_u = + compute_Flux_eucclhyd<TinyVector<Dimension>>(Fn_rho_u_lim, is_boundary_node); + NodeArray<double> Fr_rho_E = compute_Flux_eucclhyd<double>(Fn_rho_E_lim, is_boundary_node); + + apply_eucclhyd_bc(Fr_rho, Fr_rho_u, Fr_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, bc_list); + + synchronize(Fr_rho); + synchronize(Fr_rho_u); + synchronize(Fr_rho_E); + + deltaLFn_rho = compute_deltaLFn_node(Fr_rho); + deltaLFn_rho_u = compute_deltaLFn_node(Fr_rho_u); + deltaLFn_rho_E = compute_deltaLFn_node(Fr_rho_E); + } const CellValue<const TinyVector<Dimension>> A_rho = getA_rho(rho_u); const CellValue<const TinyMatrix<Dimension>> A_rho_u = getA_rho_u(rho, rho_u, rho_E); @@ -1078,21 +1803,56 @@ class EulerKineticSolverMultiDOrder2 vector_limiter(Fnp1_rho_u, Fnp1_rho_u_lim); - FaceArray<double> Fl_rho_np1 = compute_Flux1_face<double>(Fnp1_rho_lim); - FaceArray<TinyVector<Dimension>> Fl_rho_u_np1 = compute_Flux1_face<TinyVector<Dimension>>(Fnp1_rho_u_lim); - FaceArray<double> Fl_rho_E_np1 = compute_Flux1_face<double>(Fnp1_rho_E_lim); - - apply_bc(Fl_rho_np1, Fl_rho_u_np1, Fl_rho_E_np1, rho_np1, rho_u_np1, rho_E_np1, Fnp1_rho, Fnp1_rho_u, Fnp1_rho_E, - bc_list); - - synchronize(Fl_rho_np1); - synchronize(Fl_rho_u_np1); - synchronize(Fl_rho_E_np1); - - const CellArray<const double> deltaLFnp1_rho = compute_deltaLFn_face<double>(Fl_rho_np1); - const CellArray<const TinyVector<Dimension>> deltaLFnp1_rho_u = - compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u_np1); - const CellArray<const double> deltaLFnp1_rho_E = compute_deltaLFn_face<double>(Fl_rho_E_np1); + CellArray<double> deltaLFnp1_rho{p_mesh->connectivity(), nb_waves}; + CellArray<TinyVector<Dimension>> deltaLFnp1_rho_u{p_mesh->connectivity(), nb_waves}; + 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); + FaceArray<TinyVector<Dimension>> Fl_rho_u_np1 = compute_Flux_face<TinyVector<Dimension>>(Fnp1_rho_u_lim); + FaceArray<double> Fl_rho_E_np1 = compute_Flux_face<double>(Fnp1_rho_E_lim); + + 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, bc_list); + + synchronize(Fl_rho_np1); + synchronize(Fl_rho_u_np1); + synchronize(Fl_rho_E_np1); + + deltaLFnp1_rho = compute_deltaLFn_face<double>(Fl_rho_np1); + deltaLFnp1_rho_u = compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u_np1); + deltaLFnp1_rho_E = compute_deltaLFn_face<double>(Fl_rho_E_np1); + } else if (m_scheme == 2) { + NodeArray<double> Fr_rho_np1 = compute_Flux_glace<double>(Fnp1_rho_lim, is_boundary_node); + NodeArray<TinyVector<Dimension>> Fr_rho_u_np1 = + compute_Flux_glace<TinyVector<Dimension>>(Fnp1_rho_u_lim, is_boundary_node); + NodeArray<double> Fr_rho_E_np1 = compute_Flux_glace<double>(Fnp1_rho_E_lim, is_boundary_node); + + apply_glace_bc(Fr_rho_np1, Fr_rho_u_np1, Fr_rho_E_np1, Fnp1_rho_lim, Fnp1_rho_u_lim, Fnp1_rho_E_lim, bc_list); + + synchronize(Fr_rho_np1); + synchronize(Fr_rho_u_np1); + synchronize(Fr_rho_E_np1); + + deltaLFnp1_rho = compute_deltaLFn_node<double>(Fr_rho_np1); + deltaLFnp1_rho_u = compute_deltaLFn_node<TinyVector<Dimension>>(Fr_rho_u_np1); + deltaLFnp1_rho_E = compute_deltaLFn_node<double>(Fr_rho_E_np1); + } else if (m_scheme == 3) { + NodeArray<double> Fr_rho_np1 = compute_Flux_eucclhyd<double>(Fnp1_rho_lim, is_boundary_node); + NodeArray<TinyVector<Dimension>> Fr_rho_u_np1 = + compute_Flux_eucclhyd<TinyVector<Dimension>>(Fnp1_rho_u_lim, is_boundary_node); + NodeArray<double> Fr_rho_E_np1 = compute_Flux_eucclhyd<double>(Fnp1_rho_E_lim, is_boundary_node); + + apply_eucclhyd_bc(Fr_rho_np1, Fr_rho_u_np1, Fr_rho_E_np1, Fnp1_rho_lim, Fnp1_rho_u_lim, Fnp1_rho_E_lim, + bc_list); + + synchronize(Fr_rho_np1); + synchronize(Fr_rho_u_np1); + synchronize(Fr_rho_E_np1); + + deltaLFnp1_rho = compute_deltaLFn_node<double>(Fr_rho_np1); + deltaLFnp1_rho_u = compute_deltaLFn_node<TinyVector<Dimension>>(Fr_rho_u_np1); + deltaLFnp1_rho_E = compute_deltaLFn_node<double>(Fr_rho_E_np1); + } rho_np1 = compute_PFnp1<double>(Fn_rho, deltaLFn_rho, deltaLFnp1_rho); rho_u_np1 = compute_PFnp1<TinyVector<Dimension>>(Fn_rho_u, deltaLFn_rho_u, deltaLFnp1_rho_u); @@ -1182,6 +1942,7 @@ class EulerKineticSolverMultiDOrder2 const std::vector<TinyVector<MeshType::Dimension>>& lambda_vector, const size_t& SpaceOrder, const size_t& TimeOrder, + const size_t& scheme, const CellArray<const double>& Fn_rho, const CellArray<const TinyVector<MeshType::Dimension>>& Fn_rho_u, const CellArray<const double>& Fn_rho_E) @@ -1190,6 +1951,7 @@ class EulerKineticSolverMultiDOrder2 m_lambda_vector{lambda_vector}, m_SpaceOrder{SpaceOrder}, m_TimeOrder{TimeOrder}, + m_scheme{scheme}, m_Fn_rho{Fn_rho}, m_Fn_rho_u{Fn_rho_u}, m_Fn_rho_E{Fn_rho_E}, @@ -1410,6 +2172,7 @@ eulerKineticSolverMultiDOrder2( const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder, + const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E, @@ -1433,7 +2196,7 @@ eulerKineticSolverMultiDOrder2( auto Fn_rho_u = F_rho_u->get<DiscreteFunctionP0Vector<const TinyVector<MeshType::Dimension>>>(); auto Fn_rho_E = F_rho_E->get<DiscreteFunctionP0Vector<const double>>(); - EulerKineticSolverMultiDOrder2 solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, TimeOrder, + EulerKineticSolverMultiDOrder2 solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, TimeOrder, scheme, Fn_rho.cellArrays(), Fn_rho_u.cellArrays(), Fn_rho_E.cellArrays()); return solver.apply(bc_descriptor_list); @@ -1454,6 +2217,7 @@ eulerKineticSolverMultiDOrder2( const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder, + const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E, @@ -1469,6 +2233,7 @@ eulerKineticSolverMultiDOrder2( const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder, + const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E, diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.hpp b/src/scheme/EulerKineticSolverMultiDOrder2.hpp index c533e9d12ffc4fe1003d814fc47b4a4d97fce004..5eebcaf5a641de327c5a21e7354e49b6b1c6bbce 100644 --- a/src/scheme/EulerKineticSolverMultiDOrder2.hpp +++ b/src/scheme/EulerKineticSolverMultiDOrder2.hpp @@ -17,6 +17,7 @@ eulerKineticSolverMultiDOrder2( const double& eps, const size_t& SpaceOrder, const size_t& TimeOrder, + const size_t& scheme, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u, const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,