From f4218dce4355a2470e89851ab79906a6a8f2e60c Mon Sep 17 00:00:00 2001
From: Axelle <axelle.drouard@orange.fr>
Date: Tue, 3 Dec 2024 09:35:33 +0100
Subject: [PATCH] Set up second order in space
---
src/language/modules/SchemeModule.cpp | 39 +
src/scheme/CMakeLists.txt | 1 +
src/scheme/EulerKineticSolverMultiDOrder2.cpp | 1239 +++++++++++++++++
src/scheme/EulerKineticSolverMultiDOrder2.hpp | 25 +
4 files changed, 1304 insertions(+)
create mode 100644 src/scheme/EulerKineticSolverMultiDOrder2.cpp
create mode 100644 src/scheme/EulerKineticSolverMultiDOrder2.hpp
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index ef3e5dbc9..e1729dace 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -44,6 +44,7 @@
#include <scheme/EulerKineticSolverMoodFD.hpp>
#include <scheme/EulerKineticSolverMoodFV.hpp>
#include <scheme/EulerKineticSolverMultiD.hpp>
+#include <scheme/EulerKineticSolverMultiDOrder2.hpp>
#include <scheme/EulerKineticSolverOneFluxMood.hpp>
#include <scheme/EulerKineticSolverThirdOrderFD.hpp>
#include <scheme/EulerKineticSolverThirdOrderFV.hpp>
@@ -1041,6 +1042,44 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("euler_kinetic_MultiD_Order2",
+ std::function(
+
+ [](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 std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ 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,
+ bc_descriptor_list);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("euler_kinetic_MultiD_Order2",
+ std::function(
+
+ [](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 std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ 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,
+ bc_descriptor_list);
+ }
+
+ ));
+
this->_addBuiltinFunction("euler_kinetic_first_order_FV",
std::function(
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index e243b5249..ee7b44b7d 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -22,6 +22,7 @@ add_library(
EulerKineticSolverMoodFD.cpp
EulerKineticSolverMoodFV.cpp
EulerKineticSolverMultiD.cpp
+ EulerKineticSolverMultiDOrder2.cpp
EulerKineticSolverThirdOrderMoodFV.cpp
EulerKineticSolverThirdOrderMoodFD.cpp
EulerKineticSolverThirdOrderFV.cpp
diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.cpp b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
new file mode 100644
index 000000000..8fade4f03
--- /dev/null
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
@@ -0,0 +1,1239 @@
+#include <scheme/EulerKineticSolverMultiDOrder2.hpp>
+
+#include <language/utils/EvaluateAtPoints.hpp>
+#include <language/utils/InterpolateItemArray.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatFaceBoundary.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionDPkVector.hpp>
+#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <scheme/PolynomialReconstructionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <scheme/WallBoundaryConditionDescriptor.hpp>
+
+#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/QuadratureManager.hpp>
+#include <geometry/LineTransformation.hpp>
+#include <tuple>
+
+template <MeshConcept MeshType>
+class EulerKineticSolverMultiDOrder2
+{
+ private:
+ constexpr static size_t Dimension = MeshType::Dimension;
+ using Rd = TinyVector<Dimension>;
+
+ class SymmetryBoundaryCondition;
+ class WallBoundaryCondition;
+ class InflowListBoundaryCondition;
+ class OutflowBoundaryCondition;
+
+ using BoundaryCondition = std::
+ variant<SymmetryBoundaryCondition, WallBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition>;
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+ const double m_dt;
+ const double m_eps;
+ const std::vector<TinyVector<Dimension>>& m_lambda_vector;
+ const size_t m_SpaceOrder;
+ const size_t m_TimeOrder;
+ const CellArray<const double> m_Fn_rho;
+ const CellArray<const TinyVector<Dimension>> m_Fn_rho_u;
+ const CellArray<const double> m_Fn_rho_E;
+ const double m_gamma;
+ const std::shared_ptr<const MeshType> p_mesh;
+ const MeshType& m_mesh;
+ const CellValue<const double> m_Vj = MeshDataManager::instance().getMeshData(m_mesh).Vj();
+
+ const FaceValue<const TinyVector<Dimension>> nl = MeshDataManager::instance().getMeshData(m_mesh).nl();
+
+ CellValue<const double> m_inv_Vj;
+ const CellValue<const TinyVector<Dimension>> m_xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
+ const NodeValue<const TinyVector<Dimension>> m_xr = m_mesh.xr();
+ const FaceValue<const TinyVector<Dimension>> m_xl = MeshDataManager::instance().getMeshData(m_mesh).xl();
+
+ BoundaryConditionList
+ _getBCList(const MeshType& mesh,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+ {
+ BoundaryConditionList bc_list;
+
+ for (const auto& bc_descriptor : bc_descriptor_list) {
+ bool is_valid_boundary_condition = true;
+
+ switch (bc_descriptor->type()) {
+ case IBoundaryConditionDescriptor::Type::symmetry: {
+ bc_list.emplace_back(
+ SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::wall: {
+ bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::inflow_list: {
+ const InflowListBoundaryConditionDescriptor& inflow_list_bc_descriptor =
+ dynamic_cast<const InflowListBoundaryConditionDescriptor&>(*bc_descriptor);
+ if (inflow_list_bc_descriptor.functionSymbolIdList().size() != 2 + Dimension) {
+ std::ostringstream error_msg;
+ error_msg << "invalid number of functions for inflow boundary "
+ << inflow_list_bc_descriptor.boundaryDescriptor() << ", found "
+ << inflow_list_bc_descriptor.functionSymbolIdList().size() << ", expecting " << 2 + Dimension;
+ throw NormalError(error_msg.str());
+ }
+ auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+ Array<size_t> is_boundary_cell{p_mesh->numberOfCells()};
+
+ is_boundary_cell.fill(0);
+
+ auto node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+
+ auto node_list = node_boundary.nodeList();
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ auto cell_list = node_to_cell_matrix[node_list[i_node]];
+ for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+ const CellId cell_id = cell_list[i_cell];
+ is_boundary_cell[cell_id] = 1;
+ }
+ }
+
+ size_t nb_boundary_cells = sum(is_boundary_cell);
+
+ Array<CellId> cell_list{nb_boundary_cells};
+
+ size_t index = 0;
+ for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+ if (is_boundary_cell[cell_id]) {
+ cell_list[index++] = cell_id;
+ }
+ }
+ auto xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
+
+ Table<const double> cell_array_list =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::cell>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ xj, cell_list);
+ bc_list.emplace_back(InflowListBoundaryCondition(cell_array_list, cell_list));
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::outflow: {
+ bc_list.emplace_back(OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ break;
+ }
+ default: {
+ is_valid_boundary_condition = false;
+ }
+ }
+ if (not is_valid_boundary_condition) {
+ std::ostringstream error_msg;
+ error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ return bc_list;
+ }
+
+ public:
+ CellArray<TinyVector<MeshType::Dimension>>
+ getA(const DiscreteFunctionP0<const double>& rho,
+ const DiscreteFunctionP0<const double>& rho_u1,
+ const DiscreteFunctionP0<const double>& rho_u2,
+ const DiscreteFunctionP0<const double>& rho_E) const
+ {
+ // const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+
+ CellArray<TinyVector<Dimension>> vec_A{m_mesh.connectivity(), 4};
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) *
+ (rho_u1[cell_id] * rho_u1[cell_id] + rho_u2[cell_id] * rho_u2[cell_id]);
+ double p = (m_gamma - 1) * rho_e;
+
+ vec_A[cell_id][0][0] = rho_u1[cell_id];
+ vec_A[cell_id][1][0] = rho_u1[cell_id] * rho_u1[cell_id] / rho[cell_id] + p;
+ vec_A[cell_id][2][0] = rho_u1[cell_id] * rho_u2[cell_id] / rho[cell_id];
+ vec_A[cell_id][3][0] = (rho_E[cell_id] + p) * rho_u1[cell_id] / rho[cell_id];
+
+ vec_A[cell_id][0][1] = rho_u2[cell_id];
+ vec_A[cell_id][1][1] = rho_u1[cell_id] * rho_u2[cell_id] / rho[cell_id];
+ vec_A[cell_id][2][1] = rho_u2[cell_id] * rho_u2[cell_id] / rho[cell_id] + p;
+ vec_A[cell_id][3][1] = (rho_E[cell_id] + p) * rho_u2[cell_id] / rho[cell_id];
+ });
+ return vec_A;
+ }
+
+ const CellValue<const TinyVector<Dimension>>
+ getA_rho(const CellValue<const TinyVector<Dimension>>& rho_u) const
+ {
+ return rho_u;
+ }
+
+ const CellValue<const TinyMatrix<Dimension>>
+ getA_rho_u(const CellValue<const double>& rho,
+ const CellValue<const TinyVector<Dimension>>& rho_u,
+ const CellValue<const double>& rho_E) const
+ {
+ CellValue<TinyMatrix<Dimension>> vec_A{m_mesh.connectivity()};
+ const TinyMatrix<Dimension> I = identity;
+ CellValue<double> rho_e{m_mesh.connectivity()};
+ CellValue<double> p{m_mesh.connectivity()};
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ rho_e[cell_id] = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
+ p[cell_id] = (m_gamma - 1) * rho_e[cell_id];
+
+ vec_A[cell_id] = 1. / rho[cell_id] * tensorProduct(rho_u[cell_id], rho_u[cell_id]) + p[cell_id] * I;
+ });
+
+ return vec_A;
+ }
+
+ const CellValue<const TinyVector<Dimension>>
+ getA_rho_E(const CellValue<const double>& rho,
+ const CellValue<const TinyVector<Dimension>>& rho_u,
+ const CellValue<const double>& rho_E) const
+ {
+ CellValue<TinyVector<Dimension>> vec_A{m_mesh.connectivity()};
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ 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;
+
+ vec_A[cell_id] = (rho_E[cell_id] + p) / rho[cell_id] * rho_u[cell_id];
+ });
+
+ return vec_A;
+ }
+
+ const CellArray<const double>
+ compute_scalar_M(const CellValue<const double>& u, const CellValue<const TinyVector<Dimension>>& A_u)
+ {
+ const size_t nb_waves = m_lambda_vector.size();
+ CellArray<double> M{p_mesh->connectivity(), nb_waves};
+ TinyVector<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(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ for (size_t i = 0; i < nb_waves; ++i) {
+ M[cell_id][i] = (1. / nb_waves) * u[cell_id];
+
+ for (size_t d = 0; d < Dimension; ++d) {
+ M[cell_id][i] += inv_S[d] * m_lambda_vector[i][d] * A_u[cell_id][d];
+ }
+ }
+ });
+
+ return M;
+ }
+
+ const CellArray<const TinyVector<Dimension>>
+ compute_vector_M(const CellValue<const TinyVector<Dimension>>& u, const CellValue<const TinyMatrix<Dimension>>& A_u)
+ {
+ const size_t nb_waves = m_lambda_vector.size();
+ CellArray<TinyVector<Dimension>> M{p_mesh->connectivity(), nb_waves};
+ 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(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ M[cell_id][i_wave] = 1. / nb_waves * u[cell_id];
+ for (size_t d1 = 0; d1 < Dimension; ++d1) {
+ for (size_t d2 = 0; d2 < Dimension; ++d2) {
+ M[cell_id][i_wave][d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_u[cell_id](d2, d1);
+ }
+ }
+ }
+ });
+ return M;
+ }
+
+ template <typename T>
+ const CellValue<const T>
+ compute_PFnp1(const CellArray<const T>& Fn, const CellArray<const T>& deltaLFn)
+
+ {
+ CellValue<T> PFnp1{p_mesh->connectivity()};
+
+ if constexpr (is_tiny_vector_v<T>) {
+ PFnp1.fill(zero);
+ } else {
+ PFnp1.fill(0.);
+ }
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double inv_Vj = m_inv_Vj[cell_id];
+ const double dt_over_Vj = m_dt * inv_Vj;
+
+ for (size_t i_wave = 0; i_wave < m_lambda_vector.size(); ++i_wave) {
+ PFnp1[cell_id] += Fn[cell_id][i_wave] - dt_over_Vj * deltaLFn[cell_id][i_wave];
+ }
+ });
+
+ return PFnp1;
+ }
+
+ template <typename T>
+ const CellValue<const T>
+ compute_PFn(const CellArray<const T>& F)
+ {
+ const size_t nb_waves = m_lambda_vector.size();
+ CellValue<T> PFn{p_mesh->connectivity()};
+ if constexpr (is_tiny_vector_v<T>) {
+ PFn.fill(zero);
+ } else {
+ PFn.fill(0.);
+ }
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ for (size_t i = 0; i < nb_waves; ++i) {
+ PFn[cell_id] += F[cell_id][i];
+ }
+ });
+ return PFn;
+ }
+
+ template <typename T>
+ NodeArray<T>
+ compute_Flux1(CellArray<T>& Fn)
+ {
+ 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) {
+ 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] += Fn[cell_id][i_wave] * li_njr;
+ sum_li_njr += li_njr;
+ }
+ }
+ if (sum_li_njr != 0) {
+ Fr[node_id][i_wave] /= sum_li_njr;
+ }
+ }
+ });
+
+ return Fr;
+ } else {
+ throw NormalError("Glace not defined in 1d");
+ }
+ }
+
+ template <typename T>
+ FaceArrayPerNode<T>
+ compute_Flux1_eucchlyd(CellArray<T> Fn)
+ {
+ 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);
+ 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>) {
+ Flr.fill(zero);
+ } else {
+ Flr.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];
+
+ 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];
+
+ 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] += Fn[cell_id][i_wave] * li_nlr;
+ sum += li_nlr;
+ }
+ }
+ if (sum != 0) {
+ Flr[node_id][i_face][i_wave] /= sum;
+ }
+ }
+ }
+ });
+
+ return Flr;
+ } else {
+ throw NormalError("Eucclhyd not defined in 1d");
+ }
+ }
+
+ template <typename T>
+ FaceArray<T>
+ compute_Flux1_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);
+ const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
+ 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>) {
+ Fl.fill(zero);
+ } else {
+ Fl.fill(0.);
+ }
+
+ parallel_for(
+ p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ 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_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 = 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) {
+ Fl[face_id][i_wave] += DPk_Fn(cell_id, i_wave)(m_xl[face_id]);
+ }
+ }
+ }
+ });
+
+ return Fl;
+ }
+
+ 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)
+ {
+ const size_t nb_waves = m_lambda_vector.size();
+ const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
+ const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
+ 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];
+ }
+
+ 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 face_list = bc.faceList();
+ auto n = bc.outgoingNormal();
+ auto nxn = tensorProduct(n, n);
+ 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];
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ 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_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 = 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 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;
+
+ TinyVector<Dimension> A_rho_prime = rho_uj_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;
+
+ 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;
+
+ 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];
+ }
+
+ Fl_rho[face_id][i_wave] += Fn_rho_prime;
+ Fl_rho_u[face_id][i_wave] += Fn_rho_u_prime;
+ Fl_rho_E[face_id][i_wave] += Fn_rho_E_prime;
+ }
+ }
+ }
+ }
+ } else if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
+ 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];
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ 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_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;
+
+ auto n = sign * nl[face_id];
+ auto nxn = tensorProduct(n, n);
+ TinyMatrix<Dimension> I = identity;
+ auto txt = I - nxn;
+
+ 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 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;
+
+ TinyVector<Dimension> A_rho_prime = rho_uj_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;
+
+ 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;
+
+ 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];
+ }
+
+ Fl_rho[face_id][i_wave] += Fn_rho_prime;
+ Fl_rho_u[face_id][i_wave] += Fn_rho_u_prime;
+ Fl_rho_E[face_id][i_wave] += Fn_rho_E_prime;
+ }
+ }
+ }
+ }
+ } else if constexpr (std::is_same_v<BCType, OutflowBoundaryCondition>) {
+ 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];
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ 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_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;
+
+ auto n = sign * nl[face_id];
+
+ 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];
+ }
+ }
+ }
+ }
+ }
+ },
+ bc_v);
+ }
+ }
+
+ DiscreteFunctionP0Vector<double>
+ compute_deltaLFn(NodeArray<double> 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();
+ DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
+
+ const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+ 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) {
+ 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];
+
+ 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;
+ }
+ }
+ });
+ return deltaLFn;
+ } else {
+ throw NormalError("Glace not defined in 1d");
+ }
+ }
+
+ DiscreteFunctionP0Vector<double>
+ compute_deltaLFn_eucclhyd(FaceArrayPerNode<double> Fr)
+ {
+ if constexpr (Dimension > 1) {
+ const size_t nb_waves = m_lambda_vector.size();
+ DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
+
+ 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) {
+ deltaLFn[cell_id][i_wave] = 0;
+
+ 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] += Fr[node_id][i_face_node][i_wave] * li_Nlr;
+ }
+ }
+ }
+ });
+ return deltaLFn;
+ } else {
+ throw NormalError("Eucclhyd not defined in 1d");
+ }
+ }
+
+ template <typename T>
+ CellArray<T>
+ compute_deltaLFn_face(FaceArray<T> Fl)
+ {
+ if constexpr (Dimension > 1) {
+ const FaceValue<const TinyVector<Dimension>> Nl = MeshDataManager::instance().getMeshData(m_mesh).Nl();
+ const size_t nb_waves = m_lambda_vector.size();
+ CellArray<T> deltaLFn(p_mesh->connectivity(), nb_waves);
+
+ const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
+ const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
+
+ 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_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 double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+
+ const double li_Nl = sign * dot(m_lambda_vector[i_wave], Nl[face_id]);
+
+ deltaLFn[cell_id][i_wave] += li_Nl * Fl[face_id][i_wave];
+ }
+ }
+ });
+ return deltaLFn;
+ } else {
+ throw NormalError("Nl in meaningless in 1d");
+ }
+ }
+
+ CellValue<bool>
+ getInflowBoundaryCells(const BoundaryConditionList& bc_list)
+ {
+ CellValue<bool> is_boundary_cell{p_mesh->connectivity()};
+
+ is_boundary_cell.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, InflowListBoundaryCondition>) {
+ auto cell_list = bc.cellList();
+ for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+ const CellId cell_id = cell_list[i_cell];
+ is_boundary_cell[cell_id] = true;
+ }
+ }
+ },
+ bc_v);
+ }
+
+ return is_boundary_cell;
+ }
+
+ std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ apply(const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+ {
+ BoundaryConditionList bc_list = this->_getBCList(m_mesh, bc_descriptor_list);
+
+ const size_t nb_waves = m_lambda_vector.size();
+
+ CellValue<bool> is_inflow_boundary_cell = getInflowBoundaryCells(bc_list);
+
+ const CellValue<double> rho_ex{p_mesh->connectivity()};
+ const CellValue<TinyVector<Dimension>> rho_u_ex{p_mesh->connectivity()};
+ const CellValue<double> rho_E_ex{p_mesh->connectivity()};
+
+ rho_ex.fill(1); // !! A modifier en ne parcourant que les bords
+ rho_u_ex.fill(zero);
+ rho_E_ex.fill(0);
+
+ for (auto&& bc_v : bc_list) {
+ std::visit(
+ [=](auto&& bc) {
+ using BCType = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<BCType, InflowListBoundaryCondition>) {
+ auto cell_list = bc.cellList();
+ auto cell_array_list = bc.cellArrayList();
+ for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+ const CellId cell_id = cell_list[i_cell];
+ rho_ex[cell_id] = cell_array_list[i_cell][0];
+ rho_u_ex[cell_id][0] = cell_array_list[i_cell][1];
+ rho_u_ex[cell_id][1] = cell_array_list[i_cell][2];
+ rho_E_ex[cell_id] = cell_array_list[i_cell][3];
+ }
+ }
+ },
+ 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);
+
+ const CellArray<const double> Fn_rho_ex = compute_scalar_M(rho_ex, A_rho_ex);
+ const CellArray<const TinyVector<Dimension>> Fn_rho_u_ex = compute_vector_M(rho_u_ex, A_rho_u_ex);
+ const CellArray<const double> Fn_rho_E_ex = compute_scalar_M(rho_E_ex, A_rho_E_ex);
+
+ const CellArray<const double> Fn_rho = copy(m_Fn_rho);
+ const CellArray<const TinyVector<Dimension>> Fn_rho_u = copy(m_Fn_rho_u);
+ const CellArray<const double> Fn_rho_E = copy(m_Fn_rho_E);
+
+ DiscreteFunctionP0Vector<double> Fnp1_rho(p_mesh, nb_waves);
+ DiscreteFunctionP0Vector<TinyVector<Dimension>> Fnp1_rho_u(p_mesh, nb_waves);
+ DiscreteFunctionP0Vector<double> Fnp1_rho_E(p_mesh, nb_waves);
+
+ // Compute first order F
+
+ const CellValue<const double> rho = compute_PFn<double>(Fn_rho);
+ 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);
+
+ std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
+ // Reconstruction uniquement à l'intérieur donc pas utile
+ // for (auto&& bc_descriptor : bc_descriptor_list) {
+ // if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
+ // symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+ // }
+ // }
+
+ PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
+ symmetry_boundary_descriptor_list);
+
+ auto reconstruction =
+ PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho),
+ DiscreteFunctionP0Vector(p_mesh, Fn_rho_u),
+ DiscreteFunctionP0Vector(p_mesh, Fn_rho_E));
+ auto DPk_Fn_rho = reconstruction[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+ auto DPk_Fn_rho_u =
+ reconstruction[1]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
+ auto DPk_Fn_rho_E = reconstruction[2]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+
+ FaceArray<double> Fl_rho = compute_Flux1_face<double>(DPk_Fn_rho);
+ FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux1_face<TinyVector<Dimension>>(DPk_Fn_rho_u);
+ FaceArray<double> Fl_rho_E = compute_Flux1_face<double>(DPk_Fn_rho_E);
+
+ 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);
+
+ 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);
+
+ 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);
+ const CellValue<const TinyVector<Dimension>> A_rho_E = getA_rho_E(rho, rho_u, rho_E);
+
+ const CellArray<const double> M_rho = compute_scalar_M(rho, A_rho);
+ const CellArray<const TinyVector<Dimension>> M_rho_u = compute_vector_M(rho_u, A_rho_u);
+ const CellArray<const double> M_rho_E = compute_scalar_M(rho_E, A_rho_E);
+
+ const CellValue<const double> rho_np1 = compute_PFnp1<double>(Fn_rho, deltaLFn_rho);
+ const CellValue<const TinyVector<Dimension>> rho_u_np1 =
+ compute_PFnp1<TinyVector<Dimension>>(Fn_rho_u, deltaLFn_rho_u);
+ const CellValue<const double> rho_E_np1 = compute_PFnp1<double>(Fn_rho_E, deltaLFn_rho_E);
+
+ const CellValue<const TinyVector<Dimension>> A_rho_np1 = getA_rho(rho_u_np1);
+ const CellValue<const TinyMatrix<Dimension>> A_rho_u_np1 = getA_rho_u(rho_np1, rho_u_np1, rho_E_np1);
+ const CellValue<const TinyVector<Dimension>> A_rho_E_np1 = getA_rho_E(rho_np1, rho_u_np1, rho_E_np1);
+
+ const CellArray<const double> M_rho_np1 = compute_scalar_M(rho_np1, A_rho_np1);
+ const CellArray<const TinyVector<Dimension>> M_rho_u_np1 = compute_vector_M(rho_u_np1, A_rho_u_np1);
+ const CellArray<const double> M_rho_E_np1 = compute_scalar_M(rho_E_np1, A_rho_E_np1);
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (not is_inflow_boundary_cell[cell_id]) {
+ const double c1 = 1. / (m_eps + m_dt);
+ const double c2 = m_eps * m_dt * m_inv_Vj[cell_id];
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ Fnp1_rho[cell_id][i_wave] = c1 * (m_eps * Fn_rho[cell_id][i_wave] - c2 * deltaLFn_rho[cell_id][i_wave] +
+ m_dt * M_rho_np1[cell_id][i_wave]);
+ Fnp1_rho_u[cell_id][i_wave] =
+ c1 * (m_eps * Fn_rho_u[cell_id][i_wave] - c2 * deltaLFn_rho_u[cell_id][i_wave] +
+ m_dt * M_rho_u_np1[cell_id][i_wave]);
+ Fnp1_rho_E[cell_id][i_wave] =
+ c1 * (m_eps * Fn_rho_E[cell_id][i_wave] - c2 * deltaLFn_rho_E[cell_id][i_wave] +
+ m_dt * M_rho_E_np1[cell_id][i_wave]);
+ }
+ } else {
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ Fnp1_rho[cell_id][i_wave] = Fn_rho_ex[cell_id][i_wave];
+ Fnp1_rho_u[cell_id][i_wave] = Fn_rho_u_ex[cell_id][i_wave];
+ Fnp1_rho_E[cell_id][i_wave] = Fn_rho_E_ex[cell_id][i_wave];
+ }
+ }
+ });
+
+ return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(Fnp1_rho),
+ std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_u),
+ std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_E));
+ }
+
+ EulerKineticSolverMultiDOrder2(std::shared_ptr<const MeshType> mesh,
+ const double& dt,
+ const double& eps,
+ const double& gamma,
+ const std::vector<TinyVector<MeshType::Dimension>>& lambda_vector,
+ const size_t& SpaceOrder,
+ const size_t& TimeOrder,
+ const CellArray<const double>& Fn_rho,
+ const CellArray<const TinyVector<MeshType::Dimension>>& Fn_rho_u,
+ const CellArray<const double>& Fn_rho_E)
+ : m_dt{dt},
+ m_eps{eps},
+ m_lambda_vector{lambda_vector},
+ m_SpaceOrder{SpaceOrder},
+ m_TimeOrder{TimeOrder},
+ m_Fn_rho{Fn_rho},
+ m_Fn_rho_u{Fn_rho_u},
+ m_Fn_rho_E{Fn_rho_E},
+ m_gamma{gamma},
+ p_mesh{mesh},
+ m_mesh{*mesh}
+ {
+ if ((lambda_vector.size() != Fn_rho[CellId(0)].size()) or (lambda_vector.size() != Fn_rho_u[CellId(0)].size()) or
+ (lambda_vector.size() != Fn_rho_u[CellId(0)].size()) or
+ ((lambda_vector.size() != Fn_rho_E[CellId(0)].size()))) {
+ throw NormalError("Incompatible dimensions of lambda vector and Fn");
+ }
+
+ m_inv_Vj = [&] {
+ CellValue<double> inv_Vj{m_mesh.connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { inv_Vj[cell_id] = 1. / m_Vj[cell_id]; });
+ return inv_Vj;
+ }();
+ }
+
+ ~EulerKineticSolverMultiDOrder2() = default;
+};
+
+template <MeshConcept MeshType>
+class EulerKineticSolverMultiDOrder2<MeshType>::SymmetryBoundaryCondition
+{
+ public:
+ static constexpr const size_t Dimension = MeshType::Dimension;
+
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
+ const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;
+
+ public:
+ const Rd&
+ outgoingNormal() const
+ {
+ return m_mesh_flat_node_boundary.outgoingNormal();
+ }
+
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_flat_node_boundary.nodeList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_flat_node_boundary.nodeList();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_flat_face_boundary.faceList().size();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_flat_face_boundary.faceList();
+ }
+
+ SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
+ const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
+ : m_mesh_flat_node_boundary(mesh_flat_node_boundary), m_mesh_flat_face_boundary(mesh_flat_face_boundary)
+ {
+ ;
+ }
+
+ ~SymmetryBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class EulerKineticSolverMultiDOrder2<MeshType>::WallBoundaryCondition
+{
+ public:
+ static constexpr const size_t Dimension = MeshType::Dimension;
+
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
+
+ public:
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_node_boundary.nodeList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_face_boundary.faceList().size();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary)
+ : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary)
+ {
+ ;
+ }
+
+ ~WallBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class EulerKineticSolverMultiDOrder2<MeshType>::InflowListBoundaryCondition
+{
+ public:
+ static constexpr const size_t Dimension = MeshType::Dimension;
+
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const Table<const double> m_cell_array_list;
+ const Array<const CellId> m_cell_list;
+
+ public:
+ size_t
+ numberOfCells() const
+ {
+ return m_cell_list.size();
+ }
+
+ const Array<const CellId>&
+ cellList() const
+ {
+ return m_cell_list;
+ }
+
+ const Table<const double>&
+ cellArrayList() const
+ {
+ return m_cell_array_list;
+ }
+
+ InflowListBoundaryCondition(const Table<const double>& cell_array_list, const Array<const CellId> cell_list)
+ : m_cell_array_list(cell_array_list), m_cell_list(cell_list)
+ {
+ ;
+ }
+
+ ~InflowListBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class EulerKineticSolverMultiDOrder2<MeshType>::OutflowBoundaryCondition
+{
+ public:
+ static constexpr const size_t Dimension = MeshType::Dimension;
+
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
+
+ public:
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_node_boundary.nodeList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_face_boundary.faceList().size();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ OutflowBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary)
+ : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary)
+ {
+ ;
+ }
+
+ ~OutflowBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverMultiDOrder2(
+ const double& dt,
+ const double& gamma,
+ const std::vector<TinyVector<Dimension>>& lambda_vector,
+ const double& eps,
+ const size_t& SpaceOrder,
+ const size_t& TimeOrder,
+ 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,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+{
+ const std::shared_ptr<const MeshVariant> mesh_v = getCommonMesh({F_rho, F_rho_u, F_rho_E});
+ if (mesh_v.use_count() == 0) {
+ throw NormalError("F_rho, F_rho_u1, F_rho_u2 and F_rho_E have to be defined on the same mesh");
+ }
+ if (SpaceOrder > 1 or TimeOrder > 1) {
+ throw NotImplementedError("Euler kinetic solver in Multi D not implemented at order > 1");
+ }
+
+ return std::visit(
+ [&](auto&& p_mesh)
+ -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ using MeshType = mesh_type_t<decltype(p_mesh)>;
+ if constexpr (is_polygonal_mesh_v<MeshType> and (MeshType::Dimension == Dimension)) {
+ auto Fn_rho = F_rho->get<DiscreteFunctionP0Vector<const double>>();
+ 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,
+ Fn_rho.cellArrays(), Fn_rho_u.cellArrays(), Fn_rho_E.cellArrays());
+
+ return solver.apply(bc_descriptor_list);
+ } else {
+ throw NotImplementedError("Invalid mesh type for Multi D Euler Kinetic solver");
+ }
+ },
+ mesh_v->variant());
+}
+
+template std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverMultiDOrder2(
+ 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 std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
+
+template std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverMultiDOrder2(
+ 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 std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.hpp b/src/scheme/EulerKineticSolverMultiDOrder2.hpp
new file mode 100644
index 000000000..c533e9d12
--- /dev/null
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.hpp
@@ -0,0 +1,25 @@
+#ifndef EULER_KINETIC_SOLVER_MULTID_ORDER2_HPP
+#define EULER_KINETIC_SOLVER_MULTID_ORDER2_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+
+class IBoundaryConditionDescriptor;
+
+template <size_t Dimension>
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverMultiDOrder2(
+ const double& dt,
+ const double& gamma,
+ const std::vector<TinyVector<Dimension>>& lambda_vector,
+ const double& eps,
+ const size_t& SpaceOrder,
+ const size_t& TimeOrder,
+ 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,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
+
+#endif // EULER_KINETIC_SOLVER_MULTID_ORDER2_HPP
--
GitLab