From a3b5c00ba38166fae7e19fa2f86c7532896adf17 Mon Sep 17 00:00:00 2001
From: Axelle <axelle.drouard@orange.fr>
Date: Wed, 12 Mar 2025 17:04:28 +0100
Subject: [PATCH] Add local wavespeeds kinetic lagrangian scheme for multiD
problems
---
src/language/modules/KineticSchemeModule.cpp | 22 +
src/scheme/CMakeLists.txt | 1 +
.../EulerKineticSolverLagrangeMultiD.cpp | 252 +++-
.../EulerKineticSolverLagrangeMultiDLocal.cpp | 1061 +++++++++++++++++
.../EulerKineticSolverLagrangeMultiDLocal.hpp | 28 +
5 files changed, 1358 insertions(+), 6 deletions(-)
create mode 100644 src/scheme/EulerKineticSolverLagrangeMultiDLocal.cpp
create mode 100644 src/scheme/EulerKineticSolverLagrangeMultiDLocal.hpp
diff --git a/src/language/modules/KineticSchemeModule.cpp b/src/language/modules/KineticSchemeModule.cpp
index 9c16e6122..e5e8b929f 100644
--- a/src/language/modules/KineticSchemeModule.cpp
+++ b/src/language/modules/KineticSchemeModule.cpp
@@ -10,6 +10,7 @@
#include <scheme/EulerKineticSolverFirstOrderFV.hpp>
#include <scheme/EulerKineticSolverLagrangeFV.hpp>
#include <scheme/EulerKineticSolverLagrangeMultiD.hpp>
+#include <scheme/EulerKineticSolverLagrangeMultiDLocal.hpp>
#include <scheme/EulerKineticSolverMeanFluxMood.hpp>
#include <scheme/EulerKineticSolverMoodFD.hpp>
#include <scheme/EulerKineticSolverMoodFV.hpp>
@@ -860,6 +861,27 @@ KineticSchemeModule::KineticSchemeModule()
));
+ this->_addBuiltinFunction("euler_kinetic_lagrange_multiD_local",
+ std::function(
+
+ [](const double& dt, const size_t& L, const size_t& M, const size_t& k,
+ const double& gamma, const double& eps, const size_t& space_order,
+ const size_t& time_order, const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list) -> std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return eulerKineticSolverLagrangeMultiDLocal(dt, L, M, k, gamma, eps, space_order,
+ time_order, rho, u, E, p,
+ bc_descriptor_list);
+ }
+
+ ));
+
this->_addBuiltinFunction("euler_kinetic_acoustic_lagrange_FV",
std::function(
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index d563dcd9f..c6c6dc068 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -18,6 +18,7 @@ add_library(
WavesEquationNonUniformCellKineticSolver.cpp
EulerKineticSolver.cpp
EulerKineticSolverLagrangeMultiD.cpp
+ EulerKineticSolverLagrangeMultiDLocal.cpp
EulerKineticSolverMeanFluxMood.cpp
EulerKineticSolverOneFluxMood.cpp
EulerKineticSolverMoodFD.cpp
diff --git a/src/scheme/EulerKineticSolverLagrangeMultiD.cpp b/src/scheme/EulerKineticSolverLagrangeMultiD.cpp
index 3e8eb0b8a..15761bece 100644
--- a/src/scheme/EulerKineticSolverLagrangeMultiD.cpp
+++ b/src/scheme/EulerKineticSolverLagrangeMultiD.cpp
@@ -185,7 +185,7 @@ class EulerKineticSolverLagrangeMultiD
const size_t nb_waves = m_lambda_vector.size();
CellArray<TinyVector<Dimension>> M{p_mesh->connectivity(), nb_waves};
TinyVector<Dimension> inv_S = zero;
- TinyMatrix<Dimension> I = identity;
+ // TinyMatrix<Dimension> I = identity;
for (size_t d = 0; d < Dimension; ++d) {
for (size_t i = 0; i < nb_waves; ++i) {
inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
@@ -197,10 +197,8 @@ class EulerKineticSolverLagrangeMultiD
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] * p[cell_id] * I(d2, d1);
- }
+ for (size_t d = 0; d < Dimension; ++d) {
+ M[cell_id][i_wave][d] += 1. / Dimension * inv_S[d] * m_lambda_vector[i_wave][d] * p[cell_id];
}
}
});
@@ -251,6 +249,158 @@ class EulerKineticSolverLagrangeMultiD
return Fr;
}
+ NodeArray<double>
+ compute_Flux_Fp_lambdaj(const CellValue<const double>& p,
+ const CellValue<const TinyVector<Dimension>>& u,
+ const CellValue<const double>& lambda_j,
+ const NodeValue<const bool>& is_boundary_node)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+ const size_t nb_waves = m_lambda_vector.size();
+ NodeArray<double> 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();
+
+ Fr.fill(0.);
+
+ if (Dimension > 1) {
+ throw NotImplementedError("Non local lambda not implemented for dimension > 1");
+ }
+
+ // std::cout << lambda_j << "\n";
+
+ 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];
+
+ double lambda_r = 0.;
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+
+ lambda_r = std::max(lambda_r, lambda_j[cell_id]);
+ }
+
+ // Only in 1D
+ std::vector<TinyVector<Dimension>> lambda_r_vector(2);
+ lambda_r_vector[0][0] = lambda_r;
+ lambda_r_vector[1][0] = -lambda_r;
+ // std::cout << lambda_r_vector << "\n";
+
+ TinyVector<Dimension> inv_S = zero;
+ for (size_t d = 0; d < Dimension; ++d) {
+ for (size_t i = 0; i < nb_waves; ++i) {
+ inv_S[d] += lambda_r_vector[i][d] * lambda_r_vector[i][d];
+ }
+ inv_S[d] = 1. / inv_S[d];
+ }
+
+ 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(lambda_r_vector[i_wave], njr(cell_id, i_node));
+ if (li_njr > 1e-14) {
+ double Fj = (1. / nb_waves) * p[cell_id];
+
+ for (size_t d = 0; d < Dimension; ++d) {
+ Fj += inv_S[d] * lambda_r_vector[i_wave][d] * (lambda_r * lambda_r * u[cell_id][d]);
+ }
+
+ Fr[node_id][i_wave] += li_njr * Fj;
+ 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;
+ }
+
+ NodeArray<TinyVector<Dimension>>
+ compute_Flux_Fu_lambdaj(const CellValue<const double>& p,
+ const CellValue<const TinyVector<Dimension>>& u,
+ const CellValue<const double>& lambda_j,
+ const NodeValue<const bool>& is_boundary_node)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+ const size_t nb_waves = m_lambda_vector.size();
+ NodeArray<TinyVector<Dimension>> 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();
+
+ Fr.fill(zero);
+
+ if (Dimension > 1) {
+ throw NotImplementedError("Non local lambda not implemented for dimension > 1");
+ }
+
+ 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];
+
+ double lambda_r = 0.;
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+
+ lambda_r = std::max(lambda_r, lambda_j[cell_id]);
+ }
+
+ // Only in 1D
+ std::vector<TinyVector<Dimension>> lambda_r_vector(2);
+ lambda_r_vector[0][0] = lambda_r;
+ lambda_r_vector[1][0] = -lambda_r;
+
+ TinyVector<Dimension> inv_S = zero;
+ for (size_t d = 0; d < Dimension; ++d) {
+ for (size_t i = 0; i < nb_waves; ++i) {
+ inv_S[d] += lambda_r_vector[i][d] * lambda_r_vector[i][d];
+ }
+ inv_S[d] = 1. / inv_S[d];
+ }
+
+ 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(lambda_r_vector[i_wave], njr(cell_id, i_node));
+ if (li_njr > 1e-14) {
+ TinyVector<Dimension> Fj = (1. / nb_waves) * u[cell_id];
+
+ for (size_t d1 = 0; d1 < Dimension; ++d1) {
+ for (size_t d2 = 0; d2 < Dimension; ++d2) {
+ TinyMatrix<Dimension> I = identity;
+ Fj[d1] += inv_S[d2] * lambda_r_vector[i_wave][d2] * p[cell_id] * I(d2, d1);
+ }
+ }
+
+ Fr[node_id][i_wave] += li_njr * Fj;
+ 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;
+ }
+
void
apply_bc(NodeArray<double>& Fr_p,
NodeArray<TinyVector<Dimension>>& Fr_u,
@@ -503,6 +653,39 @@ class EulerKineticSolverLagrangeMultiD
return u;
}
+ const NodeValue<TinyVector<Dimension>>
+ compute_velocity_flux_lambdaj(NodeArray<double>& F, const CellValue<const double>& lambda_j)
+ {
+ const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
+ NodeValue<TinyVector<Dimension>> u{p_mesh->connectivity()};
+ const size_t nb_waves = m_lambda_vector.size();
+ u.fill(zero);
+
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+
+ double lambda_r = 0.;
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+
+ lambda_r = std::max(lambda_r, lambda_j[cell_id]);
+ }
+
+ const double inv_lambda_squared = (1. / (lambda_r * lambda_r));
+
+ std::vector<TinyVector<Dimension>> lambda_r_vector(2);
+ lambda_r_vector[0][0] = lambda_r;
+ lambda_r_vector[1][0] = -lambda_r;
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ u[node_id] += F[node_id][i_wave] * lambda_r_vector[i_wave];
+ }
+ u[node_id] = inv_lambda_squared * u[node_id];
+ });
+ return u;
+ }
+
const NodeValue<TinyVector<Dimension>>
compute_pressure_flux(NodeArray<TinyVector<Dimension>>& F)
{
@@ -514,7 +697,41 @@ class EulerKineticSolverLagrangeMultiD
m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
for (size_t dim = 0; dim < Dimension; ++dim) {
- p[node_id][dim] += F[node_id][i_wave][dim] * m_lambda_vector[i_wave][dim];
+ p[node_id][dim] += dot(F[node_id][i_wave], m_lambda_vector[i_wave]);
+ // p[node_id][0] += 0.5 * F[node_id][i_wave][dim] * m_lambda_vector[i_wave][dim];
+ // p[node_id][1] += 0.5 * F[node_id][i_wave][dim] * m_lambda_vector[i_wave][dim];
+ }
+ }
+ });
+ return p;
+ }
+
+ const NodeValue<TinyVector<Dimension>>
+ compute_pressure_flux_lambdaj(NodeArray<TinyVector<Dimension>>& F, const CellValue<const double>& lambda_j)
+ {
+ const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
+ NodeValue<TinyVector<Dimension>> p{p_mesh->connectivity()};
+ const size_t nb_waves = m_lambda_vector.size();
+ p.fill(zero);
+
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+
+ double lambda_r = 0.;
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+
+ lambda_r = std::max(lambda_r, lambda_j[cell_id]);
+ }
+
+ std::vector<TinyVector<Dimension>> lambda_r_vector(2);
+ lambda_r_vector[0][0] = lambda_r;
+ lambda_r_vector[1][0] = -lambda_r;
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ p[node_id][dim] += F[node_id][i_wave][dim] * lambda_r_vector[i_wave][dim];
}
}
});
@@ -638,14 +855,37 @@ class EulerKineticSolverLagrangeMultiD
CellValue<double> p_np1 = copy(p);
CellValue<double> E_np1 = copy(E);
+ const CellValue<const double> cj = [&] {
+ CellValue<double> computed_cj{m_mesh.connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { computed_cj[cell_id] = std::sqrt(m_gamma * p[cell_id] / rho[cell_id]); });
+ return computed_cj;
+ }();
+
+ const CellValue<const double> lambda_j = [&] {
+ CellValue<double> computed_lambda_j{m_mesh.connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(),
+ PUGS_LAMBDA(const CellId cell_id) { computed_lambda_j[cell_id] = std::abs(rho[cell_id] * cj[cell_id]); });
+ return computed_lambda_j;
+ }();
+
+ // std::cout << "c = " << cj << ", \n lambda = " << lambda_j << ", \n rho = " << rho << "\n";
+
const CellArray<const double> F_p = compute_M_p(p, u, lambda);
const CellArray<const TinyVector<Dimension>> F_u = compute_M_u(u, p);
NodeArray<double> Fr_p = compute_Flux<double>(F_p, is_boundary_node);
NodeArray<TinyVector<Dimension>> Fr_u = compute_Flux<TinyVector<Dimension>>(F_u, is_boundary_node);
+ // NodeArray<double> Fr_p = compute_Flux_Fp_lambdaj(p, u, lambda_j, is_boundary_node);
+ // NodeArray<TinyVector<Dimension>> Fr_u = compute_Flux_Fu_lambdaj(p, u, lambda_j, is_boundary_node);
+
apply_bc(Fr_p, Fr_u, p, u, F_p, F_u, lambda, bc_list);
+ // const NodeValue<TinyVector<Dimension>> ur = compute_velocity_flux_lambdaj(Fr_p, lambda_j);
+ // const NodeValue<TinyVector<Dimension>> pr = compute_pressure_flux_lambdaj(Fr_u, lambda_j);
const NodeValue<TinyVector<Dimension>> ur = compute_velocity_flux(Fr_p, lambda);
const NodeValue<TinyVector<Dimension>> pr = compute_pressure_flux(Fr_u);
diff --git a/src/scheme/EulerKineticSolverLagrangeMultiDLocal.cpp b/src/scheme/EulerKineticSolverLagrangeMultiDLocal.cpp
new file mode 100644
index 000000000..5672a020a
--- /dev/null
+++ b/src/scheme/EulerKineticSolverLagrangeMultiDLocal.cpp
@@ -0,0 +1,1061 @@
+#include <scheme/EulerKineticSolverLagrangeMultiDLocal.hpp>
+
+#include <language/utils/EvaluateAtPoints.hpp>
+#include <language/utils/InterpolateItemArray.hpp>
+#include <language/utils/InterpolateItemValue.hpp> //new
+#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/DiscreteFunctionDescriptorP0Vector.hpp>
+#include <scheme/DiscreteFunctionP0.hpp> //new
+#include <scheme/DiscreteFunctionP0Vector.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+#include <scheme/OutflowBoundaryConditionDescriptor.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 EulerKineticSolverLagrangeMultiDLocal
+{
+ 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 size_t m_L;
+ const size_t m_M;
+ const size_t m_k;
+ const double m_gamma;
+ const size_t m_spaceOrder;
+ const size_t m_timeOrder;
+ const CellValue<const double> m_rho;
+ const CellValue<const TinyVector<Dimension>> m_u;
+ const CellValue<const double> m_E;
+ const CellValue<const double> m_p;
+ const std::shared_ptr<const MeshType> p_mesh;
+ const MeshType& m_mesh;
+ const CellValue<const double> m_Vj = MeshDataManager::instance().getMeshData(m_mesh).Vj();
+ CellValue<const double> m_inv_Mj;
+ CellValue<const double> m_inv_Vj;
+ CellValue<const double> m_Mj;
+ const CellValue<const TinyVector<Dimension>> m_xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
+ const NodeValue<const TinyVector<Dimension>> m_xr = m_mesh.xr();
+
+ 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() != 1 + 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 " << 1 + Dimension;
+ throw NormalError(error_msg.str());
+ }
+ auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+ // CellValue<size_t> is_boundary_cell{p_mesh->connectivity()};
+
+ // 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();
+ auto xr = m_mesh.xr();
+
+ // Table<const double> cell_array_list =
+ // InterpolateItemArray<double(Rd)>::template interpolate<ItemType::cell>(inflow_list_bc_descriptor
+ // .functionSymbolIdList(),
+ // xj, cell_list);
+ Table<const double> node_array_list =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::node>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ xr, node_list);
+ // bc_list.emplace_back(InflowListBoundaryCondition(cell_array_list, cell_list));
+ bc_list.emplace_back(InflowListBoundaryCondition(node_array_list, node_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:
+ NodeArray<double>
+ compute_Flux_Fp(const CellValue<const double>& p,
+ const CellValue<const TinyVector<Dimension>>& u,
+ const NodeArray<const TinyVector<Dimension>>& lambda_vector,
+ const NodeValue<const bool>& is_boundary_node)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+ const size_t nb_waves = m_k;
+ NodeArray<double> 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();
+
+ 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];
+
+ TinyVector<Dimension> inv_S = zero;
+ for (size_t d = 0; d < Dimension; ++d) {
+ for (size_t i = 0; i < nb_waves; ++i) {
+ inv_S[d] += lambda_vector[node_id][i][d] * lambda_vector[node_id][i][d];
+ }
+ inv_S[d] = 1. / inv_S[d];
+ }
+ double lambda_r = std::sqrt(dot(lambda_vector[node_id][m_L], lambda_vector[node_id][m_L]));
+
+ 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(lambda_vector[node_id][i_wave], njr(cell_id, i_node));
+ if (li_njr > 1e-14) {
+ double Fj = (1. / nb_waves) * p[cell_id];
+
+ for (size_t d = 0; d < Dimension; ++d) {
+ Fj += inv_S[d] * lambda_vector[node_id][i_wave][d] * (lambda_r * lambda_r * u[cell_id][d]);
+ }
+
+ Fr[node_id][i_wave] += li_njr * Fj;
+ 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;
+ }
+
+ NodeArray<TinyVector<Dimension>>
+ compute_Flux_Fu(const CellValue<const double>& p,
+ const CellValue<const TinyVector<Dimension>>& u,
+ const NodeArray<const TinyVector<Dimension>>& lambda_vector,
+ const NodeValue<const bool>& is_boundary_node)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+ const size_t nb_waves = m_k;
+ NodeArray<TinyVector<Dimension>> 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();
+
+ Fr.fill(zero);
+
+ 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) {
+ TinyVector<Dimension> inv_S = zero;
+ for (size_t d = 0; d < Dimension; ++d) {
+ for (size_t i = 0; i < nb_waves; ++i) {
+ inv_S[d] += lambda_vector[node_id][i][d] * lambda_vector[node_id][i][d];
+ }
+ inv_S[d] = 1. / inv_S[d];
+ }
+
+ 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(lambda_vector[node_id][i_wave], njr(cell_id, i_node));
+ if (li_njr > 1e-14) {
+ TinyVector<Dimension> Fj = (1. / nb_waves) * u[cell_id];
+
+ for (size_t d = 0; d < Dimension; ++d) {
+ Fj[d] += 1. / Dimension * inv_S[d] * lambda_vector[node_id][i_wave][d] * p[cell_id];
+ }
+
+ Fr[node_id][i_wave] += li_njr * Fj;
+ 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;
+ }
+
+ void
+ apply_bc(NodeArray<double>& Fr_p,
+ NodeArray<TinyVector<Dimension>>& Fr_u,
+ const CellValue<const double>& p,
+ const CellValue<const TinyVector<Dimension>>& u,
+ const NodeArray<const TinyVector<Dimension>>& lambda_vector,
+ const BoundaryConditionList& bc_list)
+ {
+ const size_t nb_waves = m_k;
+ 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();
+
+ 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();
+
+ 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);
+
+ 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] += lambda_vector[node_id][i][d] * lambda_vector[node_id][i][d];
+ }
+ inv_S[d] = 1. / inv_S[d];
+ }
+ const double lambda_r = std::sqrt(dot(lambda_vector[node_id][m_L], lambda_vector[node_id][m_L]));
+
+ 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_p[node_id][i_wave] = 0;
+ Fr_u[node_id][i_wave] = zero;
+
+ 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(lambda_vector[node_id][i_wave], njr(cell_id, i_node_cell));
+
+ if (li_njr > 1e-14) {
+ const TinyVector<Dimension> A_p = lambda_r * lambda_r * u[cell_id];
+ const double A_u = p[cell_id];
+
+ double Fn_p = p[cell_id] / nb_waves;
+ TinyVector<Dimension> Fn_u = 1. / nb_waves * u[cell_id];
+
+ for (size_t d = 0; d < Dimension; ++d) {
+ Fn_p += inv_S[d] * lambda_vector[node_id][i_wave][d] * A_p[d];
+ Fn_u[d] += 1. / Dimension * inv_S[d] * lambda_vector[node_id][i_wave][d] * A_u;
+ }
+
+ Fr_p[node_id][i_wave] += Fn_p * li_njr;
+ Fr_u[node_id][i_wave] += li_njr * Fn_u;
+
+ 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(lambda_vector[node_id][i_wave], njr_prime);
+
+ if (li_njr_prime > 1e-14) {
+ double p_prime = p[cell_id];
+ TinyVector<Dimension> u_prime = txt * u[cell_id] - nxn * u[cell_id];
+
+ TinyVector<Dimension> A_p_prime = lambda_r * lambda_r * u_prime;
+ double A_u_prime = p_prime;
+
+ double Fn_p_prime = p_prime / nb_waves;
+ TinyVector<Dimension> Fn_u_prime = 1. / nb_waves * u_prime;
+
+ for (size_t d = 0; d < Dimension; ++d) {
+ Fn_p_prime += inv_S[d] * lambda_vector[node_id][i_wave][d] * A_p_prime[d];
+ Fn_u_prime[d] += 1. / Dimension * inv_S[d] * lambda_vector[node_id][i_wave][d] * A_u_prime;
+ }
+
+ Fr_p[node_id][i_wave] += Fn_p_prime * li_njr_prime;
+ Fr_u[node_id][i_wave] += li_njr_prime * Fn_u_prime;
+
+ sum += li_njr_prime;
+ }
+ }
+ if (sum != 0) {
+ Fr_p[node_id][i_wave] /= sum;
+ Fr_u[node_id][i_wave] = 1. / sum * Fr_u[node_id][i_wave];
+ }
+ }
+ }
+ } 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);
+
+ 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] += lambda_vector[node_id][i][d] * lambda_vector[node_id][i][d];
+ }
+ inv_S[d] = 1. / inv_S[d];
+ }
+ const double lambda_r = std::sqrt(dot(lambda_vector[node_id][m_L], lambda_vector[node_id][m_L]));
+
+ 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_p[node_id][i_wave] = 0;
+ Fr_u[node_id][i_wave] = zero;
+
+ 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(lambda_vector[node_id][i_wave], njr(cell_id, i_node_cell));
+
+ if (li_njr > 1e-14) {
+ const TinyVector<Dimension> A_p = lambda_r * lambda_r * u[cell_id];
+ const double A_u = p[cell_id];
+
+ double Fn_p = p[cell_id] / nb_waves;
+ TinyVector<Dimension> Fn_u = 1. / nb_waves * u[cell_id];
+
+ for (size_t d = 0; d < Dimension; ++d) {
+ Fn_p += inv_S[d] * lambda_vector[node_id][i_wave][d] * A_p[d];
+ Fn_u[d] += 1. / Dimension * inv_S[d] * lambda_vector[node_id][i_wave][d] * A_u;
+ }
+
+ Fr_p[node_id][i_wave] += Fn_p * li_njr;
+ Fr_u[node_id][i_wave] += li_njr * Fn_u;
+
+ 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(lambda_vector[node_id][i_wave], njr_prime);
+
+ if (li_njr_prime > 1e-14) {
+ double p_prime = p[cell_id];
+ TinyVector<Dimension> u_prime = txt * u[cell_id] - nxn * u[cell_id];
+
+ TinyVector<Dimension> A_p_prime = lambda_r * lambda_r * u_prime;
+ double A_u_prime = p_prime;
+
+ double Fn_p_prime = p_prime / nb_waves;
+ TinyVector<Dimension> Fn_u_prime = 1. / nb_waves * u_prime;
+
+ for (size_t d = 0; d < Dimension; ++d) {
+ Fn_p_prime += inv_S[d] * lambda_vector[node_id][i_wave][d] * A_p_prime[d];
+ Fn_u_prime[d] += 1. / Dimension * inv_S[d] * lambda_vector[node_id][i_wave][d] * A_u_prime;
+ }
+
+ Fr_p[node_id][i_wave] += Fn_p_prime * li_njr_prime;
+ Fr_u[node_id][i_wave] += li_njr_prime * Fn_u_prime;
+
+ sum += li_njr_prime;
+ }
+ }
+ if (sum != 0) {
+ Fr_p[node_id][i_wave] /= sum;
+ Fr_u[node_id][i_wave] = 1. / sum * Fr_u[node_id][i_wave];
+ }
+ }
+ }
+ }
+ },
+ bc_v);
+ }
+ }
+
+ const NodeValue<TinyVector<Dimension>>
+ compute_velocity_flux(NodeArray<double>& F, const NodeArray<const TinyVector<Dimension>>& lambda_vector)
+ {
+ NodeValue<TinyVector<Dimension>> u{p_mesh->connectivity()};
+ const size_t nb_waves = m_k;
+ u.fill(zero);
+
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ const double inv_lambda_squared = 1. / dot(lambda_vector[node_id][m_L], lambda_vector[node_id][m_L]);
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ u[node_id] += F[node_id][i_wave] * lambda_vector[node_id][i_wave];
+ }
+ u[node_id] = inv_lambda_squared * u[node_id];
+ });
+ return u;
+ }
+
+ const NodeValue<TinyVector<Dimension>>
+ compute_pressure_flux(NodeArray<TinyVector<Dimension>>& F,
+ const NodeArray<const TinyVector<Dimension>>& lambda_vector)
+ {
+ NodeValue<TinyVector<Dimension>> p{p_mesh->connectivity()};
+ const size_t nb_waves = m_k;
+ p.fill(zero);
+
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ p[node_id][dim] += dot(F[node_id][i_wave], lambda_vector[node_id][i_wave]);
+ }
+ }
+ });
+ return p;
+ }
+
+ const CellValue<const double>
+ compute_delta_velocity(const NodeValue<const TinyVector<Dimension>>& ur)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+ CellValue<double> delta{p_mesh->connectivity()};
+
+ const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+ delta.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_node = 0; i_node < cell_to_node.size(); ++i_node) {
+ const NodeId node_id = cell_to_node[i_node];
+ delta[cell_id] += dot(Cjr(cell_id, i_node), ur[node_id]);
+ }
+ });
+ return delta;
+ }
+
+ const CellValue<const TinyVector<Dimension>>
+ compute_delta_pressure(const NodeValue<const TinyVector<Dimension>>& pr)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+ CellValue<TinyVector<Dimension>> delta{p_mesh->connectivity()};
+
+ const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+ delta.fill(zero);
+
+ 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_node = 0; i_node < cell_to_node.size(); ++i_node) {
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ const NodeId node_id = cell_to_node[i_node];
+ delta[cell_id][dim] += Cjr(cell_id, i_node)[dim] * pr[node_id][dim];
+ }
+ }
+ });
+ return delta;
+ }
+
+ const CellValue<const double>
+ compute_delta_pu(const NodeValue<const TinyVector<Dimension>>& ur, const NodeValue<const TinyVector<Dimension>>& pr)
+ {
+ const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+ CellValue<double> delta{p_mesh->connectivity()};
+
+ const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ delta.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_node = 0; i_node < cell_to_node.size(); ++i_node) {
+ const NodeId node_id = cell_to_node[i_node];
+ TinyVector<Dimension> Fjr = zero;
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ Fjr[dim] = Cjr(cell_id, i_node)[dim] * pr[node_id][dim];
+ }
+ delta[cell_id] += dot(Fjr, ur[node_id]);
+ }
+ });
+ return delta;
+ }
+
+ const NodeValue<const bool>
+ getBoundaryNode(const BoundaryConditionList& bc_list)
+ {
+ 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);
+ }
+ return is_boundary_node;
+ }
+
+ const CellValue<const double>
+ build_lambda_j(const CellValue<double>& rho, const CellValue<double>& p)
+ {
+ CellValue<double> lambda_j{p_mesh->connectivity()};
+ CellValue<double> cj{m_mesh.connectivity()};
+
+ if constexpr (Dimension == 1) {
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ cj[cell_id] = std::sqrt(m_gamma * p[cell_id] / rho[cell_id]);
+ lambda_j[cell_id] = std::sqrt((3. * (m_k - 1.)) / (m_k + 1.)) * std::abs(rho[cell_id] * cj[cell_id]);
+ });
+ } else if constexpr (Dimension == 2) {
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ cj[cell_id] = std::sqrt(m_gamma * p[cell_id] / rho[cell_id]);
+ lambda_j[cell_id] =
+ (12. * m_L * m_L) / ((m_L + 1.) * (2. * m_L + 1.)) * 2. * std::abs(rho[cell_id] * cj[cell_id]);
+ });
+ } else if constexpr (Dimension == 3) {
+ if (m_k != 6) {
+ throw NotImplementedError("3D Lagrangian scheme not implemented for other than 6 waves model");
+ }
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ cj[cell_id] = std::sqrt(m_gamma * p[cell_id] / rho[cell_id]);
+ lambda_j[cell_id] = 3. * std::abs(rho[cell_id] * cj[cell_id]);
+ });
+ }
+ return lambda_j;
+ }
+
+ const NodeArray<const TinyVector<Dimension>>
+ build_nodal_wavespeeds(const CellValue<const double>& lambda_j)
+ {
+ const auto& node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+ NodeArray<TinyVector<Dimension>> lambda_vector(p_mesh->connectivity(), m_k);
+ const double pi = std::acos(-1);
+
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+
+ double lambda_r = 0.;
+ for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
+ const CellId cell_id = node_to_cell[i_cell];
+
+ lambda_r = std::max(lambda_r, lambda_j[cell_id]);
+
+ // std::cout << "node_id = " << node_id << ", lambda_r = " << lambda_r << "\n";
+ }
+
+ if constexpr (Dimension == 1) {
+ for (size_t i_wave = 0; i_wave < m_k; ++i_wave) {
+ lambda_vector[node_id][i_wave][0] = lambda_r * (1 - 2. * i_wave / (m_k - 1));
+ }
+ } else if constexpr (Dimension == 2) {
+ for (size_t l = 1; l < m_L + 1; ++l) {
+ for (size_t m = 1; m < 4 * m_M + 1; ++m) {
+ size_t i_wave = (l - 1) * 4 * m_M + m - 1;
+ double ll = l;
+ lambda_vector[node_id][i_wave] =
+ TinyVector<Dimension>{(ll / m_L) * lambda_r * std::cos((m * pi) / (2 * m_M)), // + 0.2 * pi),
+ (ll / m_L) * lambda_r * std::sin((m * pi) / (2 * m_M))}; // + 0.2 * pi)};
+ }
+ }
+ } else if constexpr (Dimension == 3) {
+ if (m_k != 6) {
+ throw NotImplementedError("3D Lagrangian scheme not implemented for other than 6 waves model");
+ }
+ lambda_vector[node_id][0] = TinyVector<Dimension>{lambda_r, 0., 0.};
+ lambda_vector[node_id][1] = TinyVector<Dimension>{-lambda_r, 0., 0.};
+ lambda_vector[node_id][2] = TinyVector<Dimension>{0., lambda_r, 0.};
+ lambda_vector[node_id][3] = TinyVector<Dimension>{0., -lambda_r, 0.};
+ lambda_vector[node_id][4] = TinyVector<Dimension>{0., 0., lambda_r};
+ lambda_vector[node_id][5] = TinyVector<Dimension>{0., 0., -lambda_r};
+ }
+ });
+ return lambda_vector;
+ }
+
+ std::tuple<std::shared_ptr<const MeshVariant>,
+ 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 auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+ const NodeValue<const bool> is_boundary_node = getBoundaryNode(bc_list);
+
+ const CellValue<double> rho = copy(m_rho);
+ const CellValue<TinyVector<Dimension>> u = copy(m_u);
+ const CellValue<double> p = copy(m_p);
+ const CellValue<double> E = copy(m_E);
+
+ CellValue<TinyVector<Dimension>> u_np1 = copy(u);
+ CellValue<double> p_np1 = copy(p);
+ CellValue<double> E_np1 = copy(E);
+
+ // std::cout << "test 1 \n";
+ const CellValue<const double> lambda_j = build_lambda_j(rho, p);
+ // std::cout << "lambda_j = " << lambda_j << "\n";
+ const NodeArray<const TinyVector<Dimension>> lambda_vector = build_nodal_wavespeeds(lambda_j);
+ // std::cout << "lambda_vector = " << lambda_vector << "\n";
+
+ NodeArray<double> Fr_p = compute_Flux_Fp(p, u, lambda_vector, is_boundary_node);
+ NodeArray<TinyVector<Dimension>> Fr_u = compute_Flux_Fu(p, u, lambda_vector, is_boundary_node);
+
+ apply_bc(Fr_p, Fr_u, p, u, lambda_vector, bc_list);
+
+ const NodeValue<TinyVector<Dimension>> ur = compute_velocity_flux(Fr_p, lambda_vector);
+ const NodeValue<TinyVector<Dimension>> pr = compute_pressure_flux(Fr_u, lambda_vector);
+
+ 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 node_list = bc.nodeList();
+ auto node_array_list = bc.nodeArrayList();
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ ur[node_id][dim] = node_array_list[i_node][dim];
+ pr[node_id][dim] = node_array_list[i_node][Dimension];
+ }
+ }
+ }
+ },
+ bc_v);
+ }
+
+ const CellValue<const double> delta_u = compute_delta_velocity(ur);
+ const CellValue<const TinyVector<Dimension>> delta_p = compute_delta_pressure(pr);
+ const CellValue<const double> delta_pu = compute_delta_pu(ur, pr);
+
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+ const double dt_over_Mj = m_dt * m_inv_Mj[cell_id];
+ u_np1[cell_id] -= dt_over_Mj * delta_p[cell_id];
+ E_np1[cell_id] -= dt_over_Mj * delta_pu[cell_id];
+ });
+
+ NodeValue<TinyVector<Dimension>> new_xr = copy(m_mesh.xr());
+ parallel_for(
+ m_mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { new_xr[node_id] += m_dt * ur[node_id]; });
+
+ std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(m_mesh.shared_connectivity(), new_xr);
+
+ CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+ CellValue<double> rho_np1{p_mesh->connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rho_np1[cell_id] = m_Mj[cell_id] / new_Vj[cell_id]; });
+
+ return std::make_tuple(std::make_shared<MeshVariant>(new_mesh),
+ std::make_shared<DiscreteFunctionVariant>(
+ DiscreteFunctionP0<const double>(new_mesh, rho_np1)),
+ std::make_shared<DiscreteFunctionVariant>(
+ DiscreteFunctionP0<const TinyVector<Dimension>>(new_mesh, u_np1)),
+ std::make_shared<DiscreteFunctionVariant>(
+ DiscreteFunctionP0<const double>(new_mesh, E_np1)));
+ }
+
+ EulerKineticSolverLagrangeMultiDLocal(std::shared_ptr<const MeshType> mesh,
+ const double& dt,
+ const size_t& L,
+ const size_t& M,
+ const size_t& k,
+ const double& gamma,
+ const size_t& space_order,
+ const size_t& time_order,
+ const CellValue<const double>& rho,
+ const CellValue<const TinyVector<MeshType::Dimension>>& u,
+ const CellValue<const double>& E,
+ const CellValue<const double>& p)
+ : m_dt{dt},
+ m_L{L},
+ m_M{M},
+ m_k{k},
+ m_gamma{gamma},
+ m_spaceOrder{space_order},
+ m_timeOrder{time_order},
+ m_rho{rho},
+ m_u{u},
+ m_E{E},
+ m_p{p},
+ p_mesh{mesh},
+ m_mesh{*mesh}
+ {
+ m_Mj = [&] {
+ CellValue<double> Mj{m_mesh.connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Mj[cell_id] = m_rho[cell_id] * m_Vj[cell_id]; });
+ return Mj;
+ }();
+
+ m_inv_Mj = [&] {
+ CellValue<double> inv_Mj{m_mesh.connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { inv_Mj[cell_id] = 1. / m_Mj[cell_id]; });
+ return inv_Mj;
+ }();
+
+ 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;
+ }();
+ }
+
+ ~EulerKineticSolverLagrangeMultiDLocal() = default;
+};
+
+template <MeshConcept MeshType>
+class EulerKineticSolverLagrangeMultiDLocal<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 EulerKineticSolverLagrangeMultiDLocal<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 EulerKineticSolverLagrangeMultiDLocal<MeshType>::InflowListBoundaryCondition
+{
+ public:
+ static constexpr const size_t Dimension = MeshType::Dimension;
+
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const Table<const double> m_node_array_list;
+ const Array<const NodeId> m_node_list;
+
+ public:
+ size_t
+ numberOfNodes() const
+ {
+ return m_node_list.size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_node_list;
+ }
+
+ const Table<const double>&
+ nodeArrayList() const
+ {
+ return m_node_array_list;
+ }
+
+ InflowListBoundaryCondition(const Table<const double>& node_array_list, const Array<const NodeId> node_list)
+ : m_node_array_list(node_array_list), m_node_list(node_list)
+ {
+ ;
+ }
+
+ ~InflowListBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class EulerKineticSolverLagrangeMultiDLocal<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;
+};
+
+std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverLagrangeMultiDLocal(
+ const double& dt,
+ const size_t& L,
+ const size_t& M,
+ const size_t& k,
+ const double& gamma,
+ const double& eps,
+ const size_t& space_order,
+ const size_t& time_order,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_n,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_n,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_n,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_n,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+{
+ const std::shared_ptr<const MeshVariant> mesh_v = getCommonMesh({rho_n, u_n, E_n, p_n});
+ if (mesh_v.use_count() == 0) {
+ throw NormalError("rho_n, u_n, E_n and p_n have to be defined on the same mesh");
+ }
+ if (space_order > 1 or time_order > 1) {
+ throw NotImplementedError("Euler kinetic solver in Multi D not implemented at order > 1");
+ }
+ double eps_bis = eps; // A retirer !!
+ eps_bis += eps_bis; // A retirer !!
+ return std::visit(
+ [&](auto&& p_mesh)
+ -> std::tuple<std::shared_ptr<const MeshVariant>, 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>) {
+ auto rho = rho_n->get<DiscreteFunctionP0<const double>>();
+ auto u = u_n->get<DiscreteFunctionP0<const TinyVector<MeshType::Dimension>>>();
+ auto E = E_n->get<DiscreteFunctionP0<const double>>();
+ auto p = p_n->get<DiscreteFunctionP0<const double>>();
+
+ EulerKineticSolverLagrangeMultiDLocal solver(p_mesh, dt, L, M, k, gamma, space_order, time_order,
+ rho.cellValues(), u.cellValues(), E.cellValues(), p.cellValues());
+
+ return solver.apply(bc_descriptor_list);
+ } else {
+ throw NotImplementedError("Invalid mesh type for multi-D Lagrangian Kinetic solver");
+ }
+ },
+ mesh_v->variant());
+}
diff --git a/src/scheme/EulerKineticSolverLagrangeMultiDLocal.hpp b/src/scheme/EulerKineticSolverLagrangeMultiDLocal.hpp
new file mode 100644
index 000000000..98989615a
--- /dev/null
+++ b/src/scheme/EulerKineticSolverLagrangeMultiDLocal.hpp
@@ -0,0 +1,28 @@
+#ifndef EULER_KINETIC_SOLVER_LAGRANGE_MULTID_LOCAL_HPP
+#define EULER_KINETIC_SOLVER_LAGRANGE_MULTID_LOCAL_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+
+class IBoundaryConditionDescriptor;
+
+std::tuple<std::shared_ptr<const MeshVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+eulerKineticSolverLagrangeMultiDLocal(
+ const double& dt,
+ const size_t& L,
+ const size_t& M,
+ const size_t& k,
+ const double& gamma,
+ const double& eps,
+ const size_t& space_order,
+ const size_t& time_order,
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
+
+#endif // EULER_KINETIC_SOLVER_LAGRANGE_MULTID_LOCAL_HPP
--
GitLab