diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index ee7b44b7da15f559f0d57dc4d5a5545157253e97..135500f6389fd1e8fcf6b6d982a42b294cc43a65 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -38,5 +38,6 @@ add_library(
target_link_libraries(
PugsScheme
+ PugsAnalysis
${HIGHFIVE_TARGET}
)
diff --git a/src/scheme/DiscreteFunctionDPkVector.hpp b/src/scheme/DiscreteFunctionDPkVector.hpp
index 2fefc078b192304ab0a80b63782986867e13c3d8..e8ac729f4b3e002f8c1c3b2c1cdcbd5f4aef3f77 100644
--- a/src/scheme/DiscreteFunctionDPkVector.hpp
+++ b/src/scheme/DiscreteFunctionDPkVector.hpp
@@ -139,6 +139,16 @@ class DiscreteFunctionDPkVector
return m_by_cell_coefficients[cell_id];
}
+ PUGS_INLINE auto
+ componentCoefficients(const CellId cell_id, size_t i_component) const noexcept(NO_ASSERT)
+ {
+ Assert(m_mesh_v.use_count() > 0, "DiscreteFunctionDPkVector is not built");
+ Assert(i_component < m_number_of_components, "incorrect component number");
+
+ return ComponentCoefficientView{&m_by_cell_coefficients[cell_id][i_component * m_nb_coefficients_per_component],
+ m_nb_coefficients_per_component};
+ }
+
PUGS_FORCEINLINE BasisView
operator()(const CellId cell_id, size_t i_component) const noexcept(NO_ASSERT)
{
diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
index 3ea60272a74641a139f9acfb34174554e770c4bd..4022bfda3422af3d18f9281b254d108b9ede1a51 100644
--- a/src/scheme/EulerKineticSolverMultiD.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -275,7 +275,7 @@ class EulerKineticSolverMultiD
}
auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
- Array<size_t> is_boundary_cell{p_mesh->numberOfCells()};
+ CellValue<size_t> is_boundary_cell{p_mesh->connectivity()};
is_boundary_cell.fill(0);
@@ -506,7 +506,7 @@ class EulerKineticSolverMultiD
template <typename T>
NodeArray<T>
- compute_Flux1(CellArray<T>& Fn)
+ compute_Flux1_glace(const CellArray<const T>& Fn, NodeValue<bool> is_symmetry_boundary_node)
{
if constexpr (Dimension > 1) {
const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
@@ -534,12 +534,12 @@ class EulerKineticSolverMultiD
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;
+ Fr[node_id][i_wave] += li_njr * Fn[cell_id][i_wave];
sum_li_njr += li_njr;
}
}
- if (sum_li_njr != 0) {
- Fr[node_id][i_wave] /= sum_li_njr;
+ if ((sum_li_njr != 0) and (not is_symmetry_boundary_node[node_id])) {
+ Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave];
}
}
});
@@ -552,7 +552,7 @@ class EulerKineticSolverMultiD
template <typename T>
FaceArrayPerNode<T>
- compute_Flux1_eucchlyd(CellArray<T> Fn)
+ compute_Flux1_eucclhyd(const CellArray<const T>& Fn, NodeValue<bool> is_symmetry_boundary_node)
{
if constexpr (Dimension > 1) {
const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
@@ -594,12 +594,12 @@ class EulerKineticSolverMultiD
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;
+ Flr[node_id][i_face][i_wave] += li_nlr * Fn[cell_id][i_wave];
sum += li_nlr;
}
}
- if (sum != 0) {
- Flr[node_id][i_face][i_wave] /= sum;
+ if ((sum != 0) and (not is_symmetry_boundary_node[node_id])) {
+ Flr[node_id][i_face][i_wave] = 1. / sum * Flr[node_id][i_face][i_wave];
}
}
}
@@ -653,16 +653,16 @@ class EulerKineticSolverMultiD
}
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,
+ 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();
@@ -822,27 +822,297 @@ class EulerKineticSolverMultiD
}
}
- 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,
+ const CellValue<const double> rho,
+ const CellValue<const TinyVector<Dimension>> rho_u,
+ const CellValue<const double> rho_E,
+ const BoundaryConditionList& bc_list)
+ {
+ 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();
+
+ 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];
+ }
+
+ NodeValue<bool> node_is_done{p_mesh->connectivity()};
+ node_is_done.fill(false);
+
+ 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();
+ auto n = bc.outgoingNormal();
+ auto nxn = tensorProduct(n, n);
+ TinyMatrix<Dimension> I = identity;
+ auto txt = I - nxn;
+
+ // 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;
+ // }
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+ if (not node_is_done[node_id]) {
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ 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];
+
+ 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 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) {
+ sum_li_njr += 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 rhoj_prime = rho[cell_id];
+ TinyVector<Dimension> rho_uj_prime = txt * rho_u[cell_id] - nxn * rho_u[cell_id];
+ 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];
+ }
+
+ Fr_rho[node_id][i_wave] += Fn_rho_prime * li_njr_prime;
+ Fr_rho_u[node_id][i_wave] += 1. / li_njr_prime * Fn_rho_u_prime;
+ Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_njr_prime;
+
+ sum_li_njr += li_njr_prime;
+ }
+ }
+ if (sum_li_njr != 0) {
+ Fr_rho[node_id][i_wave] /= sum_li_njr;
+ Fr_rho_u[node_id][i_wave] = 1. / sum_li_njr * Fr_rho_u[node_id][i_wave];
+ Fr_rho_E[node_id][i_wave] /= sum_li_njr;
+ }
+ }
+ node_is_done[node_id] = true;
+ }
+ }
+ }
+ },
+ bc_v);
+ }
+ }
+
+ void
+ apply_eucclhyd_bc(FaceArrayPerNode<double> Flr_rho,
+ FaceArrayPerNode<TinyVector<Dimension>> Flr_rho_u,
+ FaceArrayPerNode<double> Flr_rho_E,
+ const CellValue<const double> rho,
+ const CellValue<const TinyVector<Dimension>> rho_u,
+ const CellValue<const double> rho_E,
+ const BoundaryConditionList& bc_list)
+ {
+ const size_t nb_waves = m_lambda_vector.size();
+
+ // const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+ const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
+ const auto& face_local_numbers_in_their_nodes = p_mesh->connectivity().faceLocalNumbersInTheirNodes();
+ const auto& face_to_node_matrix = p_mesh->connectivity().faceToNodeMatrix();
+ const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
+ const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
+ const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
+
+ 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);
+
+ 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];
+
+ const auto& face_to_node = face_to_node_matrix[face_id];
+ const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
+ const auto& face_to_cell = face_to_cell_matrix[face_id];
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ for (size_t i_node = 0; i_node < face_to_node.size(); ++i_node) {
+ const NodeId node_id = face_to_node[i_node];
+
+ const auto& face_local_number_in_its_node = face_local_numbers_in_their_nodes.itemArray(face_id);
+ const size_t i_face_node = face_local_number_in_its_node[i_node];
+
+ for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
+ const CellId cell_id = face_to_cell[i_cell];
+ const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
+
+ TinyVector<Dimension> n_face = nlr(face_id, i_node);
+ const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+
+ TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
+ double li_nlr_prime = sign * dot(m_lambda_vector[i_wave], nlr_prime);
+ double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
+
+ if (li_nlr > 0) {
+ sum_li_nlr[node_id][i_face_node][i_wave] += li_nlr;
+ }
+
+ if (li_nlr_prime > 0) {
+ double rhoj_prime = rho[cell_id];
+ TinyVector<Dimension> rho_uj_prime = txt * rho_u[cell_id] - nxn * rho_u[cell_id];
+ 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];
+ }
+ Flr_rho[node_id][i_face_node][i_wave] += Fn_rho_prime * li_nlr_prime;
+ Flr_rho_u[node_id][i_face_node][i_wave] += li_nlr_prime * Fn_rho_u_prime;
+ Flr_rho_E[node_id][i_face_node][i_wave] += Fn_rho_E_prime * li_nlr_prime;
+
+ sum_li_nlr[node_id][i_face_node][i_wave] += li_nlr_prime;
+ }
+ }
+ }
+ }
+ }
+ }
+ },
+ bc_v);
+ }
+
+ for (auto&& bc_v : bc_list) {
+ std::visit(
+ [=](auto&& bc) {
+ using BCType = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+ auto face_list = bc.faceList();
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+
+ const auto& face_to_node = face_to_node_matrix[face_id];
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ for (size_t i_node = 0; i_node < face_to_node.size(); ++i_node) {
+ const NodeId node_id = face_to_node[i_node];
+
+ const auto& face_local_number_in_its_node = face_local_numbers_in_their_nodes.itemArray(face_id);
+ const size_t i_face_node = face_local_number_in_its_node[i_node];
+
+ if (sum_li_nlr[node_id][i_face_node][i_wave] != 0) {
+ Flr_rho[node_id][i_face_node][i_wave] =
+ 1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho[node_id][i_face_node][i_wave];
+ Flr_rho_u[node_id][i_face_node][i_wave] =
+ 1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho_u[node_id][i_face_node][i_wave];
+ Flr_rho_E[node_id][i_face_node][i_wave] =
+ 1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho[node_id][i_face_node][i_wave];
+ }
+ }
+ }
+ }
+ }
+ },
+ bc_v);
+ }
+ }
+
+ template <typename T>
+ CellArray<T>
+ compute_deltaLFn_glace(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();
- DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
+ 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) {
- 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;
+ deltaLFn[cell_id][i_wave] += li_Cjr * Fr[node_id][i_wave];
}
}
});
@@ -852,12 +1122,19 @@ class EulerKineticSolverMultiD
}
}
- DiscreteFunctionP0Vector<double>
- compute_deltaLFn_eucclhyd(FaceArrayPerNode<double> Fr)
+ template <typename T>
+ CellArray<T>
+ compute_deltaLFn_eucclhyd(FaceArrayPerNode<T> Flr)
{
if constexpr (Dimension > 1) {
const size_t nb_waves = m_lambda_vector.size();
- DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
+ 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();
@@ -871,8 +1148,6 @@ class EulerKineticSolverMultiD
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];
@@ -896,7 +1171,7 @@ class EulerKineticSolverMultiD
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;
+ deltaLFn[cell_id][i_wave] += li_Nlr * Flr[node_id][i_face_node][i_wave];
}
}
}
@@ -1010,6 +1285,22 @@ class EulerKineticSolverMultiD
bc_v);
}
+ NodeValue<bool> is_symmetry_boundary_node{p_mesh->connectivity()};
+ is_symmetry_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_symmetry_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);
@@ -1032,59 +1323,40 @@ 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);
- 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_E));
- // auto DPk_Fn_rho = reconstruction[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
- // // auto DPk_Fn_rho_u = reconstruction[1]->get<DiscreteFunctionDPkVector<Dimension, const TinyVector < Dimension
- // // >>> ();
- // auto DPk_Fn_rho_E = reconstruction[2]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
-
- // NodeArray<double> Fr_rho = compute_Flux1(Fn_rho);
- // NodeArray<double> Fr_rho_u1 = compute_Flux1(Fn_rho_u1);
- // NodeArray<double> Fr_rho_u2 = compute_Flux1(Fn_rho_u2);
- // NodeArray<double> Fr_rho_E = compute_Flux1(Fn_rho_E);
-
- // FaceArrayPerNode<double> Flr_rho = compute_Flux1_eucchlyd(Fn_rho);
- // FaceArrayPerNode<double> Flr_rho_u1 = compute_Flux1_eucchlyd(Fn_rho_u1);
- // FaceArrayPerNode<double> Flr_rho_u2 = compute_Flux1_eucchlyd(Fn_rho_u2);
- // FaceArrayPerNode<double> Flr_rho_E = compute_Flux1_eucchlyd(Fn_rho_E);
-
- FaceArray<double> Fl_rho = compute_Flux1_face<double>(Fn_rho);
- FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux1_face<TinyVector<Dimension>>(Fn_rho_u);
- FaceArray<double> Fl_rho_E = compute_Flux1_face<double>(Fn_rho_E);
-
- apply_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);
-
- // DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn(Fr_rho);
- // DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn(Fr_rho_u1);
- // DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn(Fr_rho_u2);
- // DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn(Fr_rho_E);
-
- // DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn_eucclhyd(Flr_rho);
- // DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn_eucclhyd(Flr_rho_u1);
- // DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn_eucclhyd(Flr_rho_u2);
- // DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn_eucclhyd(Flr_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);
+ // NodeArray<double> Fr_rho = compute_Flux1_glace<double>(Fn_rho, is_symmetry_boundary_node);
+ // NodeArray<TinyVector<Dimension>> Fr_rho_u =
+ // compute_Flux1_glace<TinyVector<Dimension>>(Fn_rho_u, is_symmetry_boundary_node);
+ // NodeArray<double> Fr_rho_E = compute_Flux1_glace<double>(Fn_rho_E, is_symmetry_boundary_node);
+
+ FaceArrayPerNode<double> Flr_rho = compute_Flux1_eucclhyd<double>(Fn_rho, is_symmetry_boundary_node);
+ FaceArrayPerNode<TinyVector<Dimension>> Flr_rho_u =
+ compute_Flux1_eucclhyd<TinyVector<Dimension>>(Fn_rho_u, is_symmetry_boundary_node);
+ FaceArrayPerNode<double> Flr_rho_E = compute_Flux1_eucclhyd<double>(Fn_rho_E, is_symmetry_boundary_node);
+
+ // FaceArray<double> Fl_rho = compute_Flux1_face<double>(Fn_rho);
+ // FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux1_face<TinyVector<Dimension>>(Fn_rho_u);
+ // FaceArray<double> Fl_rho_E = compute_Flux1_face<double>(Fn_rho_E);
+
+ // apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
+ // apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
+ apply_eucclhyd_bc(Flr_rho, Flr_rho_u, Flr_rho_E, rho, rho_u, rho_E, bc_list);
+
+ // synchronize(Flr_rho);
+ // synchronize(Flr_rho_u);
+ // synchronize(Flr_rho_E);
+
+ // const CellArray<const double> deltaLFn_rho = compute_deltaLFn_glace(Fr_rho);
+ // const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = compute_deltaLFn_glace(Fr_rho_u);
+ // const CellArray<const double> deltaLFn_rho_E = compute_deltaLFn_glace(Fr_rho_E);
+
+ const CellArray<const double> deltaLFn_rho = compute_deltaLFn_eucclhyd(Flr_rho);
+ const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = compute_deltaLFn_eucclhyd(Flr_rho_u);
+ const CellArray<const double> deltaLFn_rho_E = compute_deltaLFn_eucclhyd(Flr_rho_E);
+
+ // const CellArray<const double> deltaLFn_rho = compute_deltaLFn_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);
diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.cpp b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
index 8fade4f0387b8b4d4ff77a7f10abb56973b89463..82776378ca7350104635c83f47621f00d48cebf4 100644
--- a/src/scheme/EulerKineticSolverMultiDOrder2.cpp
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
@@ -9,6 +9,7 @@
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshVariant.hpp>
+#include <mesh/StencilManager.hpp>
#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionDPkVector.hpp>
#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
@@ -280,8 +281,8 @@ class EulerKineticSolverMultiDOrder2
}
template <typename T>
- const CellValue<const T>
- compute_PFnp1(const CellArray<const T>& Fn, const CellArray<const T>& deltaLFn)
+ CellValue<const T>
+ compute_PFnp1(const CellArray<const T>& Fn, const CellArray<const T>& deltaLFn, const CellArray<const T>& deltaLFnp1)
{
CellValue<T> PFnp1{p_mesh->connectivity()};
@@ -292,15 +293,30 @@ class EulerKineticSolverMultiDOrder2
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;
+ if (m_TimeOrder == 1) {
+ 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];
- }
- });
+ 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];
+ }
+ });
+ } else if (m_TimeOrder == 2) {
+ 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] - 0.5 * dt_over_Vj * (deltaLFn[cell_id][i_wave] + deltaLFnp1[cell_id][i_wave]);
+ }
+ });
+ } else {
+ throw NotImplementedError("Eulerian scheme in multi D not implemented for order higher than 2");
+ }
return PFnp1;
}
@@ -374,7 +390,7 @@ class EulerKineticSolverMultiDOrder2
template <typename T>
FaceArrayPerNode<T>
- compute_Flux1_eucchlyd(CellArray<T> Fn)
+ compute_Flux1_eucclhyd(CellArray<T> Fn)
{
if constexpr (Dimension > 1) {
const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
@@ -794,6 +810,137 @@ class EulerKineticSolverMultiDOrder2
return is_boundary_cell;
}
+ void
+ scalar_limiter(const CellArray<const double>& f, DiscreteFunctionDPkVector<Dimension, double>& DPk_fh) const
+ {
+ StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+ StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+ auto stencil = StencilManager::instance().getCellToCellStencilArray(m_mesh.connectivity(), stencil_descriptor,
+ symmetry_boundary_descriptor_list);
+ const size_t nb_waves = m_lambda_vector.size();
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ const double fj = f[cell_id][i_wave];
+
+ double f_min = fj;
+ double f_max = fj;
+
+ const auto cell_stencil = stencil[cell_id];
+ for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+ f_min = std::min(f_min, f[cell_stencil[i_cell]][i_wave]);
+ 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 coef2 = 1.;
+ if (std::abs(f_bar_min - fj) > eps) {
+ coef2 = (f_min - fj) / ((f_bar_min - fj));
+ }
+
+ 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) {
+ coefficients[i] *= lambda;
+ }
+ }
+ });
+ }
+
+ void
+ vector_limiter(const CellArray<const TinyVector<Dimension>>& f,
+ DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>>& DPk_fh) const
+ {
+ StencilManager::BoundaryDescriptorList symmetry_boundary_descriptor_list;
+ StencilDescriptor stencil_descriptor{1, StencilDescriptor::ConnectionType::by_nodes};
+ auto stencil = StencilManager::instance().getCellToCellStencilArray(m_mesh.connectivity(), stencil_descriptor,
+ symmetry_boundary_descriptor_list);
+ const size_t nb_waves = m_lambda_vector.size();
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ const TinyVector<Dimension> fj = f[cell_id][i_wave];
+
+ TinyVector<Dimension> f_min = fj;
+ TinyVector<Dimension> f_max = fj;
+
+ const auto cell_stencil = stencil[cell_id];
+ for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ f_min[dim] = std::min(f_min[dim], f[cell_stencil[i_cell]][i_wave][dim]);
+ f_max[dim] = std::max(f_max[dim], f[cell_stencil[i_cell]][i_wave][dim]);
+ }
+ }
+
+ TinyVector<Dimension> f_bar_min = fj;
+ TinyVector<Dimension> 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 TinyVector<Dimension> f_xk = DPk_fh(cell_id, i_wave)(m_xj[cell_k_id]);
+
+ 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]));
+ }
+ }
+
+ 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];
+ for (size_t i = 1; i < coefficients.size(); ++i) {
+ coefficients[i] *= lambda;
+ }
+ }
+ });
+ }
+
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
@@ -844,9 +991,9 @@ class EulerKineticSolverMultiDOrder2
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);
+ CellArray<double> Fnp1_rho = copy(Fn_rho);
+ CellArray<TinyVector<Dimension>> Fnp1_rho_u = copy(Fn_rho_u);
+ CellArray<double> Fnp1_rho_E = copy(Fn_rho_E);
// Compute first order F
@@ -854,6 +1001,10 @@ class EulerKineticSolverMultiDOrder2
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);
+ CellValue<const double> rho_np1 = copy(rho);
+ CellValue<const TinyVector<Dimension>> rho_u_np1 = copy(rho_u);
+ CellValue<const double> rho_E_np1 = copy(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) {
@@ -874,9 +1025,19 @@ class EulerKineticSolverMultiDOrder2
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);
+ DiscreteFunctionDPkVector<Dimension, double> Fn_rho_lim = copy(DPk_Fn_rho);
+ DiscreteFunctionDPkVector<Dimension, double> Fn_rho_E_lim = copy(DPk_Fn_rho_E);
+
+ scalar_limiter(Fn_rho, Fn_rho_lim);
+ scalar_limiter(Fn_rho_E, Fn_rho_E_lim);
+
+ DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fn_rho_u_lim = copy(DPk_Fn_rho_u);
+
+ 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);
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);
@@ -897,46 +1058,121 @@ class EulerKineticSolverMultiDOrder2
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);
+ for (size_t i_subtimestep = 0; i_subtimestep < m_TimeOrder; i_subtimestep++) {
+ auto reconstruction_np1 =
+ PolynomialReconstruction{reconstruction_descriptor}.build(DiscreteFunctionP0Vector(p_mesh, Fnp1_rho),
+ DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_u),
+ DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_E));
+ auto DPk_Fnp1_rho = reconstruction[0]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
+ auto DPk_Fnp1_rho_u =
+ reconstruction[1]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
+ auto DPk_Fnp1_rho_E = reconstruction[2]->template get<DiscreteFunctionDPkVector<Dimension, const double>>();
- 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);
+ DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_lim = copy(DPk_Fnp1_rho);
+ DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_E_lim = copy(DPk_Fnp1_rho_E);
- 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);
+ scalar_limiter(Fnp1_rho, Fnp1_rho_lim);
+ scalar_limiter(Fnp1_rho_E, Fnp1_rho_E_lim);
- 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];
- }
- }
- });
+ DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fnp1_rho_u_lim = copy(DPk_Fnp1_rho_u);
+
+ 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);
+
+ 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);
+ rho_E_np1 = compute_PFnp1<double>(Fn_rho_E, deltaLFn_rho_E, deltaLFnp1_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);
+
+ if (m_TimeOrder == 1) {
+ 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];
+ }
+ }
+ });
+ } else if (m_TimeOrder == 2) {
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ if (not is_inflow_boundary_cell[cell_id]) {
+ const double c1 = 1. / (2 * m_eps + m_dt);
+ const double c2 = 2 * 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 *
+ (2 * m_eps * m_Fn_rho[cell_id][i_wave] -
+ c2 * (deltaLFn_rho[cell_id][i_wave] + deltaLFnp1_rho[cell_id][i_wave]) +
+ m_dt * M_rho_np1[cell_id][i_wave] + m_dt * (M_rho[cell_id][i_wave] - m_Fn_rho[cell_id][i_wave]));
+ Fnp1_rho_u[cell_id][i_wave] =
+ c1 * (2 * m_eps * m_Fn_rho_u[cell_id][i_wave] -
+ c2 * (deltaLFn_rho_u[cell_id][i_wave] + deltaLFnp1_rho_u[cell_id][i_wave]) +
+ m_dt * M_rho_u_np1[cell_id][i_wave] +
+ m_dt * (M_rho_u[cell_id][i_wave] - m_Fn_rho_u[cell_id][i_wave]));
+ Fnp1_rho_E[cell_id][i_wave] =
+ c1 * (2 * m_eps * m_Fn_rho_E[cell_id][i_wave] -
+ c2 * (deltaLFn_rho_E[cell_id][i_wave] + deltaLFnp1_rho_E[cell_id][i_wave]) +
+ m_dt * M_rho_E_np1[cell_id][i_wave] +
+ m_dt * (M_rho_E[cell_id][i_wave] - m_Fn_rho_E[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];
+ }
+ }
+ });
+ } else {
+ throw NotImplementedError("Lagrangian multi D scheme not implemented for time order higher than 2");
+ }
+ }
+
+ DiscreteFunctionP0Vector<const double> DFV_Fnp1_rho = DiscreteFunctionP0Vector(p_mesh, Fnp1_rho);
+ DiscreteFunctionP0Vector<const TinyVector<Dimension>> DFV_Fnp1_rho_u = DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_u);
+ DiscreteFunctionP0Vector<const double> DFV_Fnp1_rho_E = DiscreteFunctionP0Vector(p_mesh, Fnp1_rho_E);
- return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(Fnp1_rho),
- std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_u),
- std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_E));
+ return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(DFV_Fnp1_rho),
+ std::make_shared<DiscreteFunctionVariant>(DFV_Fnp1_rho_u),
+ std::make_shared<DiscreteFunctionVariant>(DFV_Fnp1_rho_E));
}
EulerKineticSolverMultiDOrder2(std::shared_ptr<const MeshType> mesh,
@@ -1183,8 +1419,8 @@ eulerKineticSolverMultiDOrder2(
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");
+ if (SpaceOrder > 2 or TimeOrder > 2) {
+ throw NotImplementedError("Euler kinetic solver in Multi D not implemented at order > 2");
}
return std::visit(