From 3d5fe299e7b08ab2bebf0fff5e51b4e7a10fce92 Mon Sep 17 00:00:00 2001
From: HOCH PHILIPPE <philippe.hoch@gmail.com>
Date: Sun, 11 Aug 2024 19:48:01 +0200
Subject: [PATCH] mise a jour Rusanov v1 et rajout Rusanov v2 correction Cjf en
2d (vu pour certains maillages)
---
src/language/modules/SchemeModule.cpp | 48 +-
src/mesh/MeshData.hpp | 61 +-
src/scheme/CMakeLists.txt | 2 +
src/scheme/RusanovEulerianCompositeSolver.cpp | 724 +++++++--
src/scheme/RusanovEulerianCompositeSolver.hpp | 6 +-
.../RusanovEulerianCompositeSolverTools.cpp | 89 ++
.../RusanovEulerianCompositeSolverTools.hpp | 24 +
.../RusanovEulerianCompositeSolver_v2.cpp | 1397 +++++++++++++++++
.../RusanovEulerianCompositeSolver_v2.hpp | 29 +
9 files changed, 2227 insertions(+), 153 deletions(-)
create mode 100644 src/scheme/RusanovEulerianCompositeSolverTools.cpp
create mode 100644 src/scheme/RusanovEulerianCompositeSolverTools.hpp
create mode 100644 src/scheme/RusanovEulerianCompositeSolver_v2.cpp
create mode 100644 src/scheme/RusanovEulerianCompositeSolver_v2.hpp
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 3960f2667..a37461b77 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -44,6 +44,7 @@
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
#include <scheme/RusanovEulerianCompositeSolver.hpp>
+#include <scheme/RusanovEulerianCompositeSolver_v2.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <scheme/VariableBCDescriptor.hpp>
#include <utils/Socket.hpp>
@@ -507,7 +508,7 @@ SchemeModule::SchemeModule()
[](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>& E, const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
@@ -515,7 +516,7 @@ SchemeModule::SchemeModule()
const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
- return rusanovEulerianCompositeSolver(rho, u, E, c, p, bc_descriptor_list, dt);
+ return rusanovEulerianCompositeSolver(rho, u, E, gamma, c, p, bc_descriptor_list, dt);
}
));
@@ -525,7 +526,43 @@ SchemeModule::SchemeModule()
[](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>& E, const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return rusanovEulerianCompositeSolver(rho, u, E, gamma, c, p, bc_descriptor_list, dt,
+ true);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("rusanov_eulerian_composite_solver_version2",
+ std::function(
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return rusanovEulerianCompositeSolver_v2(rho, u, E, gamma, c, p, bc_descriptor_list,
+ dt);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("rusanov_eulerian_composite_solver_version2_with_checks",
+ std::function(
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
@@ -533,7 +570,8 @@ SchemeModule::SchemeModule()
const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
- return rusanovEulerianCompositeSolver(rho, u, E, c, p, bc_descriptor_list, dt, true);
+ return rusanovEulerianCompositeSolver_v2(rho, u, E, gamma, c, p, bc_descriptor_list, dt,
+ true);
}
));
@@ -797,7 +835,7 @@ SchemeModule::SchemeModule()
[](const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
- return compute_dt(u, c);
+ return toolsCompositeSolver::compute_dt(u, c);
}));
this->_addBuiltinFunction("cell_volume",
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index a23d778a6..8016f370b 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -516,22 +516,7 @@ class MeshData<Mesh<Dimension>>
Cjf(j, 1) = Rd{1};
});
m_Cjf = Cjf;
- } else if constexpr (Dimension == 2) {
- const NodeValue<const Rd>& xr = m_mesh.xr();
- const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
-
- FaceValuePerCell<Rd> Cjf(m_mesh.connectivity());
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
- for (size_t R = 0; R < cell_nodes.size(); ++R) {
- int Rp = (R + 1) % cell_nodes.size();
- const Rd& xrp_xr = (xr[cell_nodes[Rp]] - xr[cell_nodes[R]]);
- Cjf(j, R) = Rd{xrp_xr[1], -xrp_xr[0]};
- }
- });
- m_Cjf = Cjf;
- } else if (Dimension == 3) {
+ } else if (Dimension >= 2) {
auto Nl = this->Nl();
const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
@@ -609,20 +594,54 @@ class MeshData<Mesh<Dimension>>
});
m_Cje = Cje;
} else if constexpr (Dimension == 2) {
+ // EdgeValuePerCell<Rd> Cje(m_mesh.connectivity());
+ // Cje = this->Cjf();
+ // auto Cjf = this->Cjf();
+ // m_xl = FaceValue<const Rd>{m_mesh.connectivity(), m_mesh.xr().arrayView()};
+ m_Cje = EdgeValuePerCell<const Rd>{m_mesh.connectivity(), this->Cjf().arrayView()};
+ /*
const NodeValue<const Rd>& xr = m_mesh.xr();
+ const auto& cell_to_edge_matrix = m_mesh.connectivity().cellToEdgeMatrix();
+ const auto& edge_to_node_matrix = m_mesh.connectivity().edgeToNodeMatrix();
const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+ // const auto& cell_edge_is_reversed = m_mesh.connectivity().cellEgdeIsReversed();
+
EdgeValuePerCell<Rd> Cje(m_mesh.connectivity());
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_nodes = cell_to_node_matrix[j];
- for (size_t R = 0; R < cell_nodes.size(); ++R) {
- int Rp = (R + 1) % cell_nodes.size();
- const Rd& xrp_xr = (xr[cell_nodes[Rp]] - xr[cell_nodes[R]]);
- Cje(j, R) = Rd{xrp_xr[1], -xrp_xr[0]};
+ const auto& cell_edge = cell_to_edge_matrix[j];
+ const auto& cell_node = cell_to_node_matrix[j];
+ assert(cell_edge.size() == cell_node.size());
+ // size_t node_locl = 0;
+ for (size_t R = 0; R < cell_edge.size(); ++R) {
+ {
+ const EdgeId& edge_id = cell_edge[R];
+ const auto& edges_nodes = edge_to_node_matrix[edge_id];
+
+ const NodeId& node0_id = edges_nodes[0];
+ const NodeId& node1_id = edges_nodes[1];
+ assert((node0_id == cell_node[R]) or (node0_id == cell_node[R]));
+ assert((node0_id == cell_node[(R + 1) % cell_node.size()]) or
+ (node0_id == cell_node[(R + 1) % cell_node.size()]));
+
+ bool reversed = true;
+ if ((node0_id == cell_node[R]) and (node1_id == cell_node[(R + 1) % cell_node.size()]))
+ reversed = false;
+
+ const Rd& xrp_xr = (xr[node1_id] - xr[node0_id]);
+
+ if (reversed)
+ Cje(j, R) = -Rd{xrp_xr[1], -xrp_xr[0]};
+ else
+ Cje(j, R) = Rd{xrp_xr[1], -xrp_xr[0]};
+
+ // node_locl++;
+ }
}
});
m_Cje = Cje;
+ */
} else if (Dimension == 3) {
auto Nlr = this->Nlr();
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 260582273..884906ec3 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -21,7 +21,9 @@ add_library(
DiscreteFunctionVectorIntegrator.cpp
DiscreteFunctionVectorInterpoler.cpp
FluxingAdvectionSolver.cpp
+ RusanovEulerianCompositeSolverTools.cpp
RusanovEulerianCompositeSolver.cpp
+ RusanovEulerianCompositeSolver_v2.cpp
)
target_link_libraries(
diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp
index e2c9f5af5..953ab6ec6 100644
--- a/src/scheme/RusanovEulerianCompositeSolver.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver.cpp
@@ -17,94 +17,6 @@
#include <variant>
-template <class Rd>
-double
-EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // const double& rho_mean,
- const Rd& U_mean,
- // const double& E_mean,
- const double& c_mean,
- const Rd& normal) // const
-{
- const double norme_normal = l2Norm(normal);
- Rd unit_normal = normal;
- unit_normal *= 1. / norme_normal;
- const double uscaln = dot(U_mean, unit_normal);
-
- return std::max(std::fabs(uscaln - c_mean) * norme_normal, std::fabs(uscaln + c_mean) * norme_normal);
-}
-
-double
-compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
-{
- const auto& c = c_v->get<DiscreteFunctionP0<const double>>();
-
- return std::visit(
- [&](auto&& p_mesh) -> double {
- const auto& mesh = *p_mesh;
-
- using MeshType = mesh_type_t<decltype(p_mesh)>;
- static constexpr size_t Dimension = MeshType::Dimension;
- using Rd = TinyVector<Dimension>;
-
- const auto& u = u_v->get<DiscreteFunctionP0<const Rd>>();
-
- if constexpr (is_polygonal_mesh_v<MeshType>) {
- const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
- // const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
-
- const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
- const NodeValuePerCell<const Rd> njr = MeshDataManager::instance().getMeshData(mesh).njr();
-
- const FaceValuePerCell<const Rd> Cjf = MeshDataManager::instance().getMeshData(mesh).Cjf();
- const FaceValuePerCell<const Rd> njf = MeshDataManager::instance().getMeshData(mesh).njf();
-
- const EdgeValuePerCell<const Rd> Cje = MeshDataManager::instance().getMeshData(mesh).Cje();
- const EdgeValuePerCell<const Rd> nje = MeshDataManager::instance().getMeshData(mesh).nje();
-
- const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
- const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
- const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
-
- CellValue<double> local_dt{mesh.connectivity()};
-
- parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
- const auto& cell_to_node = cell_to_node_matrix[j];
- const auto& cell_to_edge = cell_to_edge_matrix[j];
- const auto& cell_to_face = cell_to_face_matrix[j];
-
- double maxm(0);
- for (size_t l = 0; l < cell_to_node.size(); ++l) {
- const Rd normalCjr = Cjr(j, l);
- // maxm = std::max(maxm, EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u, c, normalCjr));
- maxm += EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u[j], c[j], normalCjr);
- }
- for (size_t l = 0; l < cell_to_face.size(); ++l) {
- const Rd normalCjf = Cjf(j, l);
- // maxm = std::max(maxm, EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u, c, normalCjr));
- maxm += EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u[j], c[j], normalCjf);
- }
-
- if constexpr (MeshType::Dimension == 3) {
- for (size_t l = 0; l < cell_to_edge.size(); ++l) {
- const Rd normalCje = Cje(j, l);
- // maxm = std::max(maxm, EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u, c, normalCjr));
- maxm += EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u[j], c[j], normalCje);
- }
- }
-
- local_dt[j] = Vj[j] / maxm; //(Sj[j] * c[j]);
- });
-
- return min(local_dt);
- } else {
- throw NormalError("unexpected mesh type");
- }
- },
- c.meshVariant()->variant());
-}
-
template <MeshConcept MeshTypeT>
class RusanovEulerianCompositeSolver
{
@@ -122,9 +34,10 @@ class RusanovEulerianCompositeSolver
class SymmetryBoundaryCondition;
class InflowListBoundaryCondition;
class OutflowBoundaryCondition;
+ class NeumannBoundaryCondition;
- using BoundaryCondition =
- std::variant<SymmetryBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition>;
+ using BoundaryCondition = std::
+ variant<SymmetryBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition, NeumannBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
@@ -226,6 +139,31 @@ class RusanovEulerianCompositeSolver
return bc_list;
}
+ public:
+ void _applyInflowBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const;
+
+ void _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const;
+
+ void _applyWallBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const;
+
+ void _applySymmetricBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const;
+
public:
double
EvaluateMaxEigenValueInGivenUnitDirection(const double& rho_mean,
@@ -242,6 +180,12 @@ class RusanovEulerianCompositeSolver
return std::max(std::fabs(uscaln - c_mean), std::fabs(uscaln + c_mean));
}
+ inline double
+ pression(const double rho, const double epsilon, const double gam) const
+ {
+ return (gam - 1) * rho * epsilon;
+ }
+
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
@@ -249,6 +193,7 @@ class RusanovEulerianCompositeSolver
const DiscreteFunctionP0<const double>& rho_n,
const DiscreteFunctionP0<const Rd>& u_n,
const DiscreteFunctionP0<const double>& E_n,
+ const double& gamma,
const DiscreteFunctionP0<const double>& c_n,
const DiscreteFunctionP0<const double>& p_n,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
@@ -280,6 +225,24 @@ class RusanovEulerianCompositeSolver
});
// CellValue<Rp> State{p_mesh->connectivity()};
+ NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()};
+ StateAtNode.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtNode[j].fill(State[j]); });
+ EdgeValuePerCell<Rp> StateAtEdge{p_mesh->connectivity()};
+ StateAtEdge.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtEdge[j].fill(State[j]); });
+ FaceValuePerCell<Rp> StateAtFace{p_mesh->connectivity()};
+ StateAtFace.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtFace[j].fill(State[j]); });
+
+ // Conditions aux limites des etats conservatifs
+
+ _applyInflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
+ //_applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
+ _applySymmetricBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
//
// Construction du flux .. ok pour l'ordre 1
@@ -301,6 +264,174 @@ class RusanovEulerianCompositeSolver
auto Flux_qtmvt = rhoUtensUPlusPid; // rhoUtensU + Pid;
auto Flux_totnrj = (rhoE + p_n) * u;
+ // Flux avec prise en compte des states at Node/Edge/Face
+
+ NodeValuePerCell<Rd> Flux_rhoAtCellNode{p_mesh->connectivity()};
+ EdgeValuePerCell<Rd> Flux_rhoAtCellEdge{p_mesh->connectivity()};
+ FaceValuePerCell<Rd> Flux_rhoAtCellFace{p_mesh->connectivity()};
+ /*
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_to_node = cell_to_node_matrix[j];
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_rhoAtCellNode[j][l][dim] = StateAtNode[j][l][0] * StateAtNode[j][l][dim];
+ }
+
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_rhoAtCellFace[j][l][dim] = StateAtFace[j][l][0] * StateAtFace[j][l][dim];
+ }
+
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_rhoAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim];
+ }
+ });
+*/
+ NodeValuePerCell<Rdxd> Flux_qtmvtAtCellNode{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
+ EdgeValuePerCell<Rdxd> Flux_qtmvtAtCellEdge{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
+ FaceValuePerCell<Rdxd> Flux_qtmvtAtCellFace{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
+ /*
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_to_node = cell_to_node_matrix[j];
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_qtmvtAtCellNode[j][l][dim] = StateAtNode[j][l][0] * StateAtNode[j][l][dim];
+ }
+
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_qtmvtAtCellFace[j][l][dim] = StateAtFace[j][l][0] * StateAtFace[j][l][dim];
+ }
+
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_qtmvtAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim];
+ }
+ });
+*/
+ // parallel_for(
+ // p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // Flux_qtmvtAtCellNode[j] = Flux_qtmvtAtCellEdge[j] = Flux_qtmvtAtCellFace[j] = Flux_qtmvt[j];
+ // });
+
+ NodeValuePerCell<Rd> Flux_totnrjAtCellNode{p_mesh->connectivity()};
+ EdgeValuePerCell<Rd> Flux_totnrjAtCellEdge{p_mesh->connectivity()};
+ FaceValuePerCell<Rd> Flux_totnrjAtCellFace{p_mesh->connectivity()};
+
+ const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ const auto& cell_to_edge_matrix = p_mesh->connectivity().cellToEdgeMatrix();
+ const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
+
+ Flux_rhoAtCellEdge.fill(zero);
+ Flux_totnrjAtCellEdge.fill(zero);
+ Flux_qtmvtAtCellEdge.fill(zero);
+
+ // Les flux aux nodes/edge/faces
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_to_node = cell_to_node_matrix[j];
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ // Etats conservatifs aux noeuds
+ const double rhonode = StateAtNode[j][l][0];
+ Rd Unode;
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Unode[dim] = StateAtNode[j][l][dim + 1] / rhonode;
+ const double rhoEnode = StateAtNode[j][l][Dimension + 1];
+ //
+ const double Enode = rhoEnode / rhonode;
+ const double epsilonnode = Enode - .5 * dot(Unode, Unode);
+ const Rd rhoUnode = rhonode * Unode;
+ const Rdxd rhoUtensUnode = tensorProduct(rhoUnode, Unode);
+
+ const double Pressionnode = pression(rhonode, epsilonnode, gamma);
+
+ const double rhoEnodePlusP = rhoEnode + Pressionnode;
+
+ Rdxd rhoUtensUPlusPidnode(identity);
+ rhoUtensUPlusPidnode *= Pressionnode;
+ rhoUtensUPlusPidnode += rhoUtensUnode;
+
+ Flux_rhoAtCellNode[j][l] = rhoUnode;
+ Flux_qtmvtAtCellNode[j][l] = rhoUtensUPlusPidnode;
+ Flux_totnrjAtCellNode[j][l] = rhoEnodePlusP * Unode;
+ }
+
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const double rhoface = StateAtFace[j][l][0];
+ Rd Uface;
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Uface[dim] = StateAtFace[j][l][dim + 1] / rhoface;
+ const double rhoEface = StateAtFace[j][l][Dimension + 1];
+ //
+ const double Eface = rhoEface / rhoface;
+ const double epsilonface = Eface - .5 * dot(Uface, Uface);
+ const Rd rhoUface = rhoface * Uface;
+ const Rdxd rhoUtensUface = tensorProduct(rhoUface, Uface);
+
+ const double Pressionface = pression(rhoface, epsilonface, gamma);
+
+ const double rhoEfacePlusP = rhoEface + Pressionface;
+
+ Rdxd rhoUtensUPlusPidface(identity);
+ rhoUtensUPlusPidface *= Pressionface;
+ rhoUtensUPlusPidface += rhoUtensUface;
+
+ Flux_rhoAtCellFace[j][l] = rhoUface;
+ Flux_qtmvtAtCellFace[j][l] = rhoUtensUPlusPidface;
+ Flux_totnrjAtCellFace[j][l] = rhoEfacePlusP * Uface;
+ }
+
+ if constexpr (Dimension == 3) {
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ const double rhoedge = StateAtEdge[j][l][0];
+ Rd Uedge;
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Uedge[dim] = StateAtEdge[j][l][dim + 1] / rhoedge;
+ const double rhoEedge = StateAtEdge[j][l][Dimension + 1];
+ //
+ const double Eedge = rhoEedge / rhoedge;
+ const double epsilonedge = Eedge - .5 * dot(Uedge, Uedge);
+ const Rd rhoUedge = rhoedge * Uedge;
+ const Rdxd rhoUtensUedge = tensorProduct(rhoUedge, Uedge);
+
+ const double Pressionedge = pression(rhoedge, epsilonedge, gamma);
+
+ const double rhoEedgePlusP = rhoEedge + Pressionedge;
+
+ Rdxd rhoUtensUPlusPidedge(identity);
+ rhoUtensUPlusPidedge *= Pressionedge;
+ rhoUtensUPlusPidedge += rhoUtensUedge;
+
+ Flux_rhoAtCellEdge[j][l] = rhoUedge;
+ Flux_qtmvtAtCellEdge[j][l] = rhoUtensUPlusPidedge;
+ Flux_totnrjAtCellEdge[j][l] = rhoEedgePlusP * Uedge;
+ }
+ }
+ });
+
+ // parallel_for(
+ // p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // Flux_totnrjAtCellNode[j] = Flux_totnrjAtCellEdge[j] = Flux_totnrjAtCellFace[j] = Flux_totnrj[j];
+ // });
+
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
auto volumes_cell = mesh_data.Vj();
@@ -351,10 +482,11 @@ class RusanovEulerianCompositeSolver
const unsigned int R = node_local_number_in_its_cells[j];
const Rd& Cjr_loc = Cjr(J, R);
- maxabsVpNormCjr = std::max(maxabsVpNormCjr,
- EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // rhoNode,
- uNode, // ENode,
- cNode, Cjr_loc));
+ maxabsVpNormCjr =
+ std::max(maxabsVpNormCjr,
+ toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // rhoNode,
+ uNode, // ENode,
+ cNode, Cjr_loc));
}
ViscosityNodeMatrix[l] *= maxabsVpNormCjr;
});
@@ -392,17 +524,18 @@ class RusanovEulerianCompositeSolver
const Rd& Cjf_loc = Cjf(J, R);
- maxabsVpNormCjf = std::max(maxabsVpNormCjf,
- EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // rhoEdge,
- uFace, // EEdge,
- cFace, Cjf_loc));
+ maxabsVpNormCjf =
+ std::max(maxabsVpNormCjf,
+ toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // rhoEdge,
+ uFace, // EEdge,
+ cFace, Cjf_loc));
}
ViscosityFaceMatrix[l] *= maxabsVpNormCjf;
});
// Computation of Viscosity Matrices at nodes and edges are done
- // We may compute Gjr, Gje
+ // We may compute Gjr, Gjf
NodeValuePerCell<Rp> Gjr{p_mesh->connectivity()};
Gjr.fill(zero);
@@ -411,16 +544,14 @@ class RusanovEulerianCompositeSolver
FaceValuePerCell<Rp> Gjf{p_mesh->connectivity()};
Gjf.fill(zero);
- const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
- const auto& cell_to_edge_matrix = p_mesh->connectivity().cellToEdgeMatrix();
- const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
-
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_to_node = cell_to_node_matrix[j];
for (size_t l = 0; l < cell_to_node.size(); ++l) {
- const NodeId& node = cell_to_node[l];
+ const NodeId& node = cell_to_node[l];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node);
+
const auto& node_to_cell = node_to_cell_matrix[node];
const Rd& Cjr_loc = Cjr(j, l);
@@ -428,17 +559,18 @@ class RusanovEulerianCompositeSolver
for (size_t k = 0; k < node_to_cell.size(); ++k) {
const CellId K = node_to_cell[k];
+ const size_t R = node_local_number_in_its_cells[k];
- const Rd& u_Cjr = Flux_qtmvt[K] * Cjr_loc;
+ const Rd& u_Cjr = Flux_qtmvtAtCellNode[K][R] * Cjr_loc; // Flux_qtmvt[K] * Cjr_loc;
// const Rp&
- statediff += State[j] - State[K];
+ statediff += StateAtNode[j][l] - StateAtNode[K][R]; // State[j] - State[K];
// const Rp& diff = ViscosityNodeMatrix[node] * statediff;
- Gjr[j][l][0] += dot(Flux_rho[K], Cjr_loc);
+ Gjr[j][l][0] += dot(Flux_rhoAtCellNode[K][R], Cjr_loc); // dot(Flux_rho[K], Cjr_loc);
for (size_t d = 0; d < Dimension; ++d)
Gjr[j][l][1 + d] += u_Cjr[d];
- Gjr[j][l][1 + Dimension] += dot(Flux_totnrj[K], Cjr_loc);
+ Gjr[j][l][1 + Dimension] += dot(Flux_totnrjAtCellNode[K][R], Cjr_loc); // dot(Flux_totnrj[K], Cjr_loc);
// Gjr[j][l] += diff;
}
@@ -479,25 +611,28 @@ class RusanovEulerianCompositeSolver
const auto& cell_to_face = cell_to_face_matrix[j];
for (size_t l = 0; l < cell_to_face.size(); ++l) {
- const FaceId& face = cell_to_face[l];
+ const FaceId& face = cell_to_face[l];
+ const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(face);
+
const auto& face_to_cell = face_to_cell_matrix[face];
const Rd& Cjf_loc = Cjf(j, l);
Rp statediff(zero);
for (size_t k = 0; k < face_to_cell.size(); ++k) {
- const CellId K = face_to_cell[k];
+ const CellId K = face_to_cell[k];
+ const unsigned int R = face_local_number_in_its_cells[k];
- const Rd& u_Cjf = Flux_qtmvt[K] * Cjf_loc;
+ const Rd& u_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc; // Flux_qtmvt[K] * Cjf_loc;
// const Rp&
- statediff += State[j] - State[K];
+ statediff += StateAtFace[j][l] - StateAtFace[K][R]; // State[j] - State[K];
// const Rp& diff = ViscosityFaceMatrix[face] * statediff;
- Gjf[j][l][0] += dot(Flux_rho[K], Cjf_loc);
+ Gjf[j][l][0] += dot(Flux_rhoAtCellFace[K][R], Cjf_loc); // dot(Flux_rho[K], Cjf_loc);
for (size_t d = 0; d < Dimension; ++d)
Gjf[j][l][1 + d] += u_Cjf[d];
- Gjf[j][l][1 + Dimension] += dot(Flux_totnrj[K], Cjf_loc);
+ Gjf[j][l][1 + Dimension] += dot(Flux_totnrjAtCellFace[K][R], Cjf_loc); // dot(Flux_totnrj[K], Cjf_loc);
// Gjf[j][l] += diff;
}
@@ -513,6 +648,9 @@ class RusanovEulerianCompositeSolver
FaceValue<double> MaxErrorConservationFace(p_mesh->connectivity());
MaxErrorConservationFace.fill(0.);
+ // const auto& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+ auto xl = mesh_data.xl();
+ auto Nl = mesh_data.Nl();
parallel_for(
p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) {
@@ -526,11 +664,21 @@ class RusanovEulerianCompositeSolver
const unsigned int R = face_local_number_in_its_cells[k];
SumGjf += Gjf[K][R];
}
+ /*
+ std::cout << " maxErr face " << xl[l] << " err " << l2Norm(SumGjf) << " size face cell"
+ << face_to_cell.size() << " Cell " << face_to_cell[0] << " / " << face_to_cell[1] << " Cjf0 "
+ << Cjf(face_to_cell[0], 0) << " / Cjf1 " << Cjf(face_to_cell[1], 1) << "\n";
+
+ if (l2Norm(SumGjf) > 1e-3)
+ throw UnexpectedError("Pb consevvation face loc");
+ */
MaxErrorConservationFace[l] = l2Norm(SumGjf);
// MaxErrorConservationFace = std::max(MaxErrorConservationFace, l2Norm(SumGjf));
}
});
std::cout << " Max Error Face " << max(MaxErrorConservationFace) << "\n";
+ if (max(MaxErrorConservationFace) > 1e-3)
+ throw UnexpectedError("Pb conservation Face");
}
if constexpr (Dimension == 3) {
@@ -573,10 +721,11 @@ class RusanovEulerianCompositeSolver
const Rd& Cje_loc = Cje(J, R);
- maxabsVpNormCje = std::max(maxabsVpNormCje,
- EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // rhoEdge,
- uEdge, // EEdge,
- cEdge, Cje_loc));
+ maxabsVpNormCje =
+ std::max(maxabsVpNormCje,
+ toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // rhoEdge,
+ uEdge, // EEdge,
+ cEdge, Cje_loc));
}
ViscosityEdgeMatrix[l] *= maxabsVpNormCje;
});
@@ -587,7 +736,9 @@ class RusanovEulerianCompositeSolver
const auto& cell_to_edge = cell_to_edge_matrix[j];
for (size_t l = 0; l < cell_to_edge.size(); ++l) {
- const EdgeId& edge = cell_to_edge[l];
+ const EdgeId& edge = cell_to_edge[l];
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge);
+
const auto& edge_to_cell = edge_to_cell_matrix[edge];
const Rd& Cje_loc = Cje(j, l);
@@ -595,17 +746,18 @@ class RusanovEulerianCompositeSolver
for (size_t k = 0; k < edge_to_cell.size(); ++k) {
const CellId K = edge_to_cell[k];
+ const size_t R = edge_local_number_in_its_cells[k];
- const Rd& u_Cje = Flux_qtmvt[K] * Cje_loc;
+ const Rd& u_Cje = Flux_qtmvtAtCellEdge[K][R] * Cje_loc; // Flux_qtmvt[K] * Cje_loc;
// const Rp&
- statediff += State[j] - State[K];
+ statediff += StateAtEdge[j][l] - StateAtEdge[K][R]; // State[j] - State[K];
// const Rp& diff = ViscosityEdgeMatrix[edge] * statediff;
- Gje[j][l][0] += dot(Flux_rho[K], Cje_loc);
+ Gje[j][l][0] += dot(Flux_rhoAtCellEdge[K][R], Cje_loc); // dot(Flux_rho[K], Cje_loc);
for (size_t d = 0; d < Dimension; ++d)
Gje[j][l][1 + d] += u_Cje[d];
- Gje[j][l][1 + Dimension] += dot(Flux_totnrj[K], Cje_loc);
+ Gje[j][l][1 + Dimension] += dot(Flux_totnrjAtCellEdge[K][R], Cje_loc); // dot(Flux_totnrj[K], Cje_loc);
// Gje[j][l] += diff;
}
@@ -700,8 +852,108 @@ class RusanovEulerianCompositeSolver
};
template <MeshConcept MeshType>
-class RusanovEulerianCompositeSolver<MeshType>::SymmetryBoundaryCondition
+class RusanovEulerianCompositeSolver<MeshType>::NeumannBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver<Mesh<2>>::NeumannBoundaryCondition
+{
+ 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();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary)
+ : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary)
+ {
+ ;
+ }
+};
+
+template <>
+class RusanovEulerianCompositeSolver<Mesh<3>>::NeumannBoundaryCondition
{
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const MeshEdgeBoundary m_mesh_edge_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
+
+ public:
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_node_boundary.nodeList().size();
+ }
+ size_t
+ numberOfEdges() const
+ {
+ return m_mesh_edge_boundary.edgeList().size();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const EdgeId>&
+ edgeList() const
+ {
+ return m_mesh_edge_boundary.edgeList();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
+ const MeshEdgeBoundary& mesh_edge_boundary,
+ const MeshFaceBoundary& mesh_face_boundary)
+ : m_mesh_node_boundary(mesh_node_boundary),
+
+ m_mesh_edge_boundary(mesh_edge_boundary),
+
+ m_mesh_face_boundary(mesh_face_boundary)
+ {
+ ;
+ }
};
template <>
@@ -991,6 +1243,7 @@ rusanovEulerianCompositeSolver(
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
+ const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
@@ -1021,7 +1274,7 @@ rusanovEulerianCompositeSolver(
return RusanovEulerianCompositeSolver<MeshType>{}.solve(p_mesh,
rho_v->get<DiscreteFunctionP0<const double>>(),
u_v->get<DiscreteFunctionP0<const Rd>>(),
- E_v->get<DiscreteFunctionP0<const double>>(),
+ E_v->get<DiscreteFunctionP0<const double>>(), gamma,
c_v->get<DiscreteFunctionP0<const double>>(),
p_v->get<DiscreteFunctionP0<const double>>(),
bc_descriptor_list, dt, check);
@@ -1032,3 +1285,224 @@ rusanovEulerianCompositeSolver(
},
mesh_v->variant());
}
+
+template <MeshConcept MeshType>
+void
+RusanovEulerianCompositeSolver<MeshType>::_applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const
+{
+ for (const auto& boundary_condition : bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<OutflowBoundaryCondition, T>) {
+ throw NormalError("Not implemented");
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <MeshConcept MeshType>
+void
+RusanovEulerianCompositeSolver<MeshType>::_applySymmetricBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const
+{
+ for (const auto& boundary_condition : bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement SYMMETRY \n";
+ const Rd& normal = bc.outgoingNormal();
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+ const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
+ // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+ const auto& face_list = bc.faceList();
+ const auto& node_list = bc.nodeList();
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+
+ const auto& node_cell_list = node_to_cell_matrix[node_id];
+ // Assert(face_cell_list.size() == 1);
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+
+ for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+ CellId node_cell_id = node_cell_list[i_cell];
+ size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim];
+ // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim];
+ }
+ }
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+
+ const auto& face_cell_list = face_to_cell_matrix[face_id];
+ Assert(face_cell_list.size() == 1);
+
+ CellId face_cell_id = face_cell_list[0];
+ size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+
+ if constexpr (Dimension == 3) {
+ const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();
+
+ const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();
+
+ const auto& edge_list = bc.edgeList();
+
+ for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
+ const EdgeId edge_id = edge_list[i_edge];
+
+ const auto& edge_cell_list = edge_to_cell_matrix[edge_id];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
+ CellId edge_cell_id = edge_cell_list[i_cell];
+ size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <MeshConcept MeshType>
+void
+RusanovEulerianCompositeSolver<MeshType>::_applyInflowBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const
+{
+ for (const auto& boundary_condition : bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<InflowListBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement INFLOW \n";
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+ const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
+ // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+ const auto& face_list = bc.faceList();
+ const auto& node_list = bc.nodeList();
+
+ const auto& face_array_list = bc.faceArrayList();
+ const 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];
+
+ const auto& node_cell_list = node_to_cell_matrix[node_id];
+ // Assert(face_cell_list.size() == 1);
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+
+ for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+ CellId node_cell_id = node_cell_list[i_cell];
+ size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateNode[node_cell_id][node_local_number_in_cell][dim] = node_array_list[i_node][dim];
+ }
+ }
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+
+ const auto& face_cell_list = face_to_cell_matrix[face_id];
+ Assert(face_cell_list.size() == 1);
+
+ CellId face_cell_id = face_cell_list[0];
+ size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][dim];
+ }
+
+ if constexpr (Dimension == 3) {
+ const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();
+
+ const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();
+ // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+ const auto& edge_list = bc.edgeList();
+
+ const auto& edge_array_list = bc.edgeArrayList();
+
+ for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
+ const EdgeId edge_id = edge_list[i_edge];
+
+ const auto& edge_cell_list = edge_to_cell_matrix[edge_id];
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ // Assert(face_cell_list.size() == 1);
+
+ for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
+ CellId edge_cell_id = edge_cell_list[i_cell];
+ size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+}
diff --git a/src/scheme/RusanovEulerianCompositeSolver.hpp b/src/scheme/RusanovEulerianCompositeSolver.hpp
index 357adf245..593b6f763 100644
--- a/src/scheme/RusanovEulerianCompositeSolver.hpp
+++ b/src/scheme/RusanovEulerianCompositeSolver.hpp
@@ -4,12 +4,13 @@
#include <mesh/MeshVariant.hpp>
#include <scheme/DiscreteFunctionVariant.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
#include <memory>
#include <vector>
-double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+// double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+// const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
@@ -18,6 +19,7 @@ rusanovEulerianCompositeSolver(
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E,
+ const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
diff --git a/src/scheme/RusanovEulerianCompositeSolverTools.cpp b/src/scheme/RusanovEulerianCompositeSolverTools.cpp
new file mode 100644
index 000000000..d4055d7e9
--- /dev/null
+++ b/src/scheme/RusanovEulerianCompositeSolverTools.cpp
@@ -0,0 +1,89 @@
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+
+template <class Rd>
+double
+toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // const double& rho_mean,
+ const Rd& U_mean,
+ // const double& E_mean,
+ const double& c_mean,
+ const Rd& normal) // const
+{
+ const double norme_normal = l2Norm(normal);
+ Rd unit_normal = normal;
+ unit_normal *= 1. / norme_normal;
+ const double uscaln = dot(U_mean, unit_normal);
+
+ return std::max(std::fabs(uscaln - c_mean) * norme_normal, std::fabs(uscaln + c_mean) * norme_normal);
+}
+
+double
+toolsCompositeSolver::compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
+{
+ const auto& c = c_v->get<DiscreteFunctionP0<const double>>();
+
+ return std::visit(
+ [&](auto&& p_mesh) -> double {
+ const auto& mesh = *p_mesh;
+
+ using MeshType = mesh_type_t<decltype(p_mesh)>;
+ static constexpr size_t Dimension = MeshType::Dimension;
+ using Rd = TinyVector<Dimension>;
+
+ const auto& u = u_v->get<DiscreteFunctionP0<const Rd>>();
+
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+ // const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
+
+ const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+ const NodeValuePerCell<const Rd> njr = MeshDataManager::instance().getMeshData(mesh).njr();
+
+ const FaceValuePerCell<const Rd> Cjf = MeshDataManager::instance().getMeshData(mesh).Cjf();
+ const FaceValuePerCell<const Rd> njf = MeshDataManager::instance().getMeshData(mesh).njf();
+
+ const EdgeValuePerCell<const Rd> Cje = MeshDataManager::instance().getMeshData(mesh).Cje();
+ const EdgeValuePerCell<const Rd> nje = MeshDataManager::instance().getMeshData(mesh).nje();
+
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+
+ CellValue<double> local_dt{mesh.connectivity()};
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ double maxm(0);
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ const Rd normalCjr = Cjr(j, l);
+ // maxm = std::max(maxm, EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u, c, normalCjr));
+ maxm += EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u[j], c[j], normalCjr);
+ }
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const Rd normalCjf = Cjf(j, l);
+ // maxm = std::max(maxm, EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u, c, normalCjr));
+ maxm += EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u[j], c[j], normalCjf);
+ }
+
+ if constexpr (MeshType::Dimension == 3) {
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ const Rd normalCje = Cje(j, l);
+ // maxm = std::max(maxm, EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u, c, normalCjr));
+ maxm += EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u[j], c[j], normalCje);
+ }
+ }
+
+ local_dt[j] = Vj[j] / maxm; //(Sj[j] * c[j]);
+ });
+
+ return min(local_dt);
+ } else {
+ throw NormalError("unexpected mesh type");
+ }
+ },
+ c.meshVariant()->variant());
+}
diff --git a/src/scheme/RusanovEulerianCompositeSolverTools.hpp b/src/scheme/RusanovEulerianCompositeSolverTools.hpp
new file mode 100644
index 000000000..cc2a4477c
--- /dev/null
+++ b/src/scheme/RusanovEulerianCompositeSolverTools.hpp
@@ -0,0 +1,24 @@
+#ifndef RUSANOV_EULERIAN_COMPOSITE_SOLVER_TOOLS_HPP
+#define RUSANOV_EULERIAN_COMPOSITE_SOLVER_TOOLS_HPP
+
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+#include <vector>
+
+namespace toolsCompositeSolver
+{
+template <class Rd>
+double EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // const double& rho_mean,
+ const Rd& U_mean,
+ // const double& E_mean,
+ const double& c_mean,
+ const Rd& normal);
+
+double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
+} // namespace toolsCompositeSolver
+#endif
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
new file mode 100644
index 000000000..3bfe9a877
--- /dev/null
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
@@ -0,0 +1,1397 @@
+#include <scheme/RusanovEulerianCompositeSolver_v2.hpp>
+
+#include <language/utils/InterpolateItemArray.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshEdgeBoundary.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatEdgeBoundary.hpp>
+#include <mesh/MeshFlatFaceBoundary.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+
+#include <variant>
+
+template <MeshConcept MeshTypeT>
+class RusanovEulerianCompositeSolver_v2
+{
+ private:
+ using MeshType = MeshTypeT;
+
+ static constexpr size_t Dimension = MeshType::Dimension;
+
+ using Rdxd = TinyMatrix<Dimension>;
+ using Rd = TinyVector<Dimension>;
+
+ using Rpxp = TinyMatrix<Dimension + 2>;
+ using Rp = TinyVector<Dimension + 2>;
+
+ class SymmetryBoundaryCondition;
+ class InflowListBoundaryCondition;
+ class OutflowBoundaryCondition;
+ class NeumannBoundaryCondition;
+
+ using BoundaryCondition = std::
+ variant<SymmetryBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition, NeumannBoundaryCondition>;
+
+ using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+ 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: {
+ if constexpr (Dimension == 2) {
+ bc_list.emplace_back(
+ SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ } else {
+ static_assert(Dimension == 3);
+ bc_list.emplace_back(
+ SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::inflow_list: {
+ const InflowListBoundaryConditionDescriptor& inflow_list_bc_descriptor =
+ dynamic_cast<const InflowListBoundaryConditionDescriptor&>(*bc_descriptor);
+ if (inflow_list_bc_descriptor.functionSymbolIdList().size() != 2 + Dimension) {
+ std::ostringstream error_msg;
+ error_msg << "invalid number of functions for inflow boundary "
+ << inflow_list_bc_descriptor.boundaryDescriptor() << ", found "
+ << inflow_list_bc_descriptor.functionSymbolIdList().size() << ", expecting " << 2 + Dimension;
+ throw NormalError(error_msg.str());
+ }
+
+ if constexpr (Dimension == 2) {
+ auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ Table<const double> node_values =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::node>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ mesh.xr(), node_boundary.nodeList());
+
+ auto xl = MeshDataManager::instance().getMeshData(mesh).xl();
+
+ auto face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ Table<const double> face_values =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::face>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ xl, face_boundary.faceList());
+
+ bc_list.emplace_back(InflowListBoundaryCondition(node_boundary, face_boundary, node_values, face_values));
+ } else {
+ static_assert(Dimension == 3);
+ auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ Table<const double> node_values =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::node>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ mesh.xr(), node_boundary.nodeList());
+
+ auto xe = MeshDataManager::instance().getMeshData(mesh).xe();
+
+ auto edge_boundary = getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ Table<const double> edge_values =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::edge>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ xe, edge_boundary.edgeList());
+
+ auto xl = MeshDataManager::instance().getMeshData(mesh).xl();
+
+ auto face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+ Table<const double> face_values =
+ InterpolateItemArray<double(Rd)>::template interpolate<ItemType::face>(inflow_list_bc_descriptor
+ .functionSymbolIdList(),
+ xl, face_boundary.faceList());
+
+ bc_list.emplace_back(InflowListBoundaryCondition(node_boundary, edge_boundary, face_boundary, node_values,
+ edge_values, face_values));
+ }
+ break;
+ }
+ case IBoundaryConditionDescriptor::Type::outflow: {
+ std::cout << "outflow not implemented yet\n";
+ 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 Rusanov v2 Eulerian Composite solver";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ return bc_list;
+ }
+
+ public:
+ void _applyInflowBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const;
+
+ void _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const;
+
+ void _applySymmetricBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const;
+
+ public:
+ double
+ EvaluateMaxEigenValueInGivenUnitDirection(const double& rho_mean,
+ const Rd& U_mean,
+ const double& E_mean,
+ const double& c_mean,
+ const Rd& normal) const
+ {
+ const double norme_normal = l2norm(normal);
+ Rd unit_normal = normal;
+ unit_normal *= 1. / norme_normal;
+ const double uscaln = dot(U_mean, unit_normal);
+
+ return std::max(std::fabs(uscaln - c_mean), std::fabs(uscaln + c_mean));
+ }
+
+ inline double
+ pression(const double rho, const double epsilon, const double gam) const
+ {
+ return (gam - 1) * rho * epsilon;
+ }
+
+ std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+ solve(const std::shared_ptr<const MeshType>& p_mesh,
+ const DiscreteFunctionP0<const double>& rho_n,
+ const DiscreteFunctionP0<const Rd>& u_n,
+ const DiscreteFunctionP0<const double>& E_n,
+ const double& gamma,
+ const DiscreteFunctionP0<const double>& c_n,
+ const DiscreteFunctionP0<const double>& p_n,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const double& dt,
+ const bool checkLocalConservation) const
+ {
+ auto rho = copy(rho_n);
+ auto u = copy(u_n);
+ auto E = copy(E_n);
+ // auto c = copy(c_n);
+ // auto p = copy(p_n);
+
+ auto bc_list = this->_getBCList(*p_mesh, bc_descriptor_list);
+
+ auto rhoE = rho * E;
+ auto rhoU = rho * u;
+
+ // Construction du vecteur des variables conservatives
+ //
+ // Ci dessous juste ordre 1
+ //
+ CellValue<Rp> State{p_mesh->connectivity()};
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ State[j][0] = rho[j];
+ for (size_t d = 0; d < Dimension; ++d)
+ State[j][1 + d] = rhoU[j][d];
+ State[j][1 + Dimension] = rhoE[j];
+ });
+
+ // CellValue<Rp> State{p_mesh->connectivity()};
+ NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()};
+ StateAtNode.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtNode[j].fill(State[j]); });
+ EdgeValuePerCell<Rp> StateAtEdge{p_mesh->connectivity()};
+ StateAtEdge.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtEdge[j].fill(State[j]); });
+ FaceValuePerCell<Rp> StateAtFace{p_mesh->connectivity()};
+ StateAtFace.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtFace[j].fill(State[j]); });
+
+ // Conditions aux limites des etats conservatifs
+
+ _applyInflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
+ //_applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
+ _applySymmetricBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
+
+ //
+ // Construction du flux .. ok pour l'ordre 1
+ //
+ CellValue<Rdxd> rhoUtensU{p_mesh->connectivity()};
+ CellValue<Rdxd> Pid(p_mesh->connectivity());
+ Pid.fill(identity);
+ CellValue<Rdxd> rhoUtensUPlusPid{p_mesh->connectivity()};
+ rhoUtensUPlusPid.fill(zero);
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ rhoUtensU[j] = tensorProduct(rhoU[j], u[j]);
+ Pid[j] *= p_n[j];
+ rhoUtensUPlusPid[j] = rhoUtensU[j] + Pid[j];
+ });
+
+ auto Flux_rho = rhoU;
+ auto Flux_qtmvt = rhoUtensUPlusPid; // rhoUtensU + Pid;
+ auto Flux_totnrj = (rhoE + p_n) * u;
+
+ // Flux avec prise en compte des states at Node/Edge/Face
+
+ NodeValuePerCell<Rd> Flux_rhoAtCellNode{p_mesh->connectivity()};
+ EdgeValuePerCell<Rd> Flux_rhoAtCellEdge{p_mesh->connectivity()};
+ FaceValuePerCell<Rd> Flux_rhoAtCellFace{p_mesh->connectivity()};
+ /*
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_to_node = cell_to_node_matrix[j];
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_rhoAtCellNode[j][l][dim] = StateAtNode[j][l][0] * StateAtNode[j][l][dim];
+ }
+
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_rhoAtCellFace[j][l][dim] = StateAtFace[j][l][0] * StateAtFace[j][l][dim];
+ }
+
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_rhoAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim];
+ }
+ });
+*/
+ NodeValuePerCell<Rdxd> Flux_qtmvtAtCellNode{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
+ EdgeValuePerCell<Rdxd> Flux_qtmvtAtCellEdge{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
+ FaceValuePerCell<Rdxd> Flux_qtmvtAtCellFace{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
+ /*
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_to_node = cell_to_node_matrix[j];
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_qtmvtAtCellNode[j][l][dim] = StateAtNode[j][l][0] * StateAtNode[j][l][dim];
+ }
+
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_qtmvtAtCellFace[j][l][dim] = StateAtFace[j][l][0] * StateAtFace[j][l][dim];
+ }
+
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Flux_qtmvtAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim];
+ }
+ });
+*/
+ // parallel_for(
+ // p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // Flux_qtmvtAtCellNode[j] = Flux_qtmvtAtCellEdge[j] = Flux_qtmvtAtCellFace[j] = Flux_qtmvt[j];
+ // });
+
+ NodeValuePerCell<Rd> Flux_totnrjAtCellNode{p_mesh->connectivity()};
+ EdgeValuePerCell<Rd> Flux_totnrjAtCellEdge{p_mesh->connectivity()};
+ FaceValuePerCell<Rd> Flux_totnrjAtCellFace{p_mesh->connectivity()};
+
+ const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ const auto& cell_to_edge_matrix = p_mesh->connectivity().cellToEdgeMatrix();
+ const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
+
+ Flux_rhoAtCellEdge.fill(zero);
+ Flux_totnrjAtCellEdge.fill(zero);
+ Flux_qtmvtAtCellEdge.fill(zero);
+
+ // Les flux aux nodes/edge/faces
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_to_node = cell_to_node_matrix[j];
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ // Etats conservatifs aux noeuds
+ const double rhonode = StateAtNode[j][l][0];
+ Rd Unode;
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Unode[dim] = StateAtNode[j][l][dim + 1] / rhonode;
+ const double rhoEnode = StateAtNode[j][l][Dimension + 1];
+ //
+ const double Enode = rhoEnode / rhonode;
+ const double epsilonnode = Enode - .5 * dot(Unode, Unode);
+ const Rd rhoUnode = rhonode * Unode;
+ const Rdxd rhoUtensUnode = tensorProduct(rhoUnode, Unode);
+
+ const double Pressionnode = pression(rhonode, epsilonnode, gamma);
+
+ const double rhoEnodePlusP = rhoEnode + Pressionnode;
+
+ Rdxd rhoUtensUPlusPidnode(identity);
+ rhoUtensUPlusPidnode *= Pressionnode;
+ rhoUtensUPlusPidnode += rhoUtensUnode;
+
+ Flux_rhoAtCellNode[j][l] = rhoUnode;
+ Flux_qtmvtAtCellNode[j][l] = rhoUtensUPlusPidnode;
+ Flux_totnrjAtCellNode[j][l] = rhoEnodePlusP * Unode;
+ }
+
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const double rhoface = StateAtFace[j][l][0];
+ Rd Uface;
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Uface[dim] = StateAtFace[j][l][dim + 1] / rhoface;
+ const double rhoEface = StateAtFace[j][l][Dimension + 1];
+ //
+ const double Eface = rhoEface / rhoface;
+ const double epsilonface = Eface - .5 * dot(Uface, Uface);
+ const Rd rhoUface = rhoface * Uface;
+ const Rdxd rhoUtensUface = tensorProduct(rhoUface, Uface);
+
+ const double Pressionface = pression(rhoface, epsilonface, gamma);
+
+ const double rhoEfacePlusP = rhoEface + Pressionface;
+
+ Rdxd rhoUtensUPlusPidface(identity);
+ rhoUtensUPlusPidface *= Pressionface;
+ rhoUtensUPlusPidface += rhoUtensUface;
+
+ Flux_rhoAtCellFace[j][l] = rhoUface;
+ Flux_qtmvtAtCellFace[j][l] = rhoUtensUPlusPidface;
+ Flux_totnrjAtCellFace[j][l] = rhoEfacePlusP * Uface;
+ }
+
+ if constexpr (Dimension == 3) {
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ const double rhoedge = StateAtEdge[j][l][0];
+ Rd Uedge;
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Uedge[dim] = StateAtEdge[j][l][dim + 1] / rhoedge;
+ const double rhoEedge = StateAtEdge[j][l][Dimension + 1];
+ //
+ const double Eedge = rhoEedge / rhoedge;
+ const double epsilonedge = Eedge - .5 * dot(Uedge, Uedge);
+ const Rd rhoUedge = rhoedge * Uedge;
+ const Rdxd rhoUtensUedge = tensorProduct(rhoUedge, Uedge);
+
+ const double Pressionedge = pression(rhoedge, epsilonedge, gamma);
+
+ const double rhoEedgePlusP = rhoEedge + Pressionedge;
+
+ Rdxd rhoUtensUPlusPidedge(identity);
+ rhoUtensUPlusPidedge *= Pressionedge;
+ rhoUtensUPlusPidedge += rhoUtensUedge;
+
+ Flux_rhoAtCellEdge[j][l] = rhoUedge;
+ Flux_qtmvtAtCellEdge[j][l] = rhoUtensUPlusPidedge;
+ Flux_totnrjAtCellEdge[j][l] = rhoEedgePlusP * Uedge;
+ }
+ }
+ });
+
+ // parallel_for(
+ // p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // Flux_totnrjAtCellNode[j] = Flux_totnrjAtCellEdge[j] = Flux_totnrjAtCellFace[j] = Flux_totnrj[j];
+ // });
+
+ MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
+
+ auto volumes_cell = mesh_data.Vj();
+
+ // Calcul des matrices de viscosité aux node/edge/face
+
+ const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+ const NodeValuePerCell<const Rd> njr = mesh_data.njr();
+
+ const FaceValuePerCell<const Rd> Cjf = mesh_data.Cjf();
+ const FaceValuePerCell<const Rd> njf = mesh_data.njf();
+
+ const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
+
+ const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
+ const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
+
+ // We compute Gjr, Gjf
+
+ NodeValuePerCell<Rp> Gjr{p_mesh->connectivity()};
+ Gjr.fill(zero);
+ EdgeValuePerCell<Rp> Gje{p_mesh->connectivity()};
+ Gje.fill(zero);
+ FaceValuePerCell<Rp> Gjf{p_mesh->connectivity()};
+ Gjf.fill(zero);
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_to_node = cell_to_node_matrix[j];
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ const NodeId& node = cell_to_node[l];
+ const auto& node_to_cell = node_to_cell_matrix[node];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node);
+
+ const Rd& Cjr_loc = Cjr(j, l);
+
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const CellId K = node_to_cell[k];
+ const size_t R = node_local_number_in_its_cells[k];
+ const Rd& Ckr_loc = Cjr(K, R);
+
+ // Une moyenne entre les etats jk
+
+ Rd uNode = .5 * (u_n[j] + u_n[K]);
+ double cNode = .5 * (c_n[j] + c_n[K]);
+
+ // Viscosity j k
+ Rpxp ViscosityMatrixJK(identity);
+ const double MaxmaxabsVpNormjk =
+ std::max(toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uNode, cNode,
+ Cjr_loc),
+ toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uNode, cNode,
+ Ckr_loc));
+
+ ViscosityMatrixJK *= MaxmaxabsVpNormjk;
+ const Rd& u_Cjr = Flux_qtmvtAtCellNode[K][R] * Cjr_loc; // Flux_qtmvt[K] * Cjr_loc;
+
+ const Rp& statediff = StateAtNode[j][l] - StateAtNode[K][R]; // State[j] - State[K];
+ const Rp& diff = ViscosityMatrixJK * statediff;
+
+ Gjr[j][l][0] += dot(Flux_rhoAtCellNode[K][R], Cjr_loc); // dot(Flux_rho[K], Cjr_loc);
+ for (size_t d = 0; d < Dimension; ++d)
+ Gjr[j][l][1 + d] += u_Cjr[d];
+ Gjr[j][l][1 + Dimension] += dot(Flux_totnrjAtCellNode[K][R], Cjr_loc); // dot(Flux_totnrj[K], Cjr_loc);
+
+ Gjr[j][l] += diff;
+ }
+
+ Gjr[j][l] *= 1. / node_to_cell.size();
+ }
+ });
+
+ if (checkLocalConservation) {
+ auto is_boundary_node = p_mesh->connectivity().isBoundaryNode();
+
+ NodeValue<double> MaxErrorConservationNode(p_mesh->connectivity());
+ MaxErrorConservationNode.fill(0.);
+ // double MaxErrorConservationNode = 0;
+ parallel_for(
+ p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId l) {
+ const auto& node_to_cell = node_to_cell_matrix[l];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(l);
+
+ if (not is_boundary_node[l]) {
+ Rp SumGjr(zero);
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const CellId K = node_to_cell[k];
+ const unsigned int R = node_local_number_in_its_cells[k];
+ SumGjr += Gjr[K][R];
+ }
+ // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr));
+ MaxErrorConservationNode[l] = l2Norm(SumGjr);
+ }
+ });
+ std::cout << " Max Error Node " << max(MaxErrorConservationNode) << "\n";
+ }
+ //
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // Edge
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const FaceId& face = cell_to_face[l];
+ const auto& face_to_cell = face_to_cell_matrix[face];
+ const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(face);
+
+ const Rd& Cjf_loc = Cjf(j, l);
+
+ for (size_t k = 0; k < face_to_cell.size(); ++k) {
+ const CellId K = face_to_cell[k];
+ const unsigned int R = face_local_number_in_its_cells[k];
+ const Rd& Ckf_loc = Cjf(K, R);
+ // Une moyenne entre les etats jk
+
+ Rd uFace = .5 * (u_n[j] + u_n[K]);
+ double cFace = .5 * (c_n[j] + c_n[K]);
+
+ // Viscosity j k
+ Rpxp ViscosityMatrixJK(identity);
+ const double MaxmaxabsVpNormjk =
+ std::max(toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uFace, cFace,
+ Cjf_loc),
+ toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uFace, cFace,
+ Ckf_loc));
+
+ ViscosityMatrixJK *= MaxmaxabsVpNormjk;
+
+ const Rd& u_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc; // Flux_qtmvt[K] * Cjf_loc;
+
+ const Rp& statediff = StateAtFace[j][l] - StateAtFace[K][R]; // State[j] - State[K];
+ const Rp& diff = ViscosityMatrixJK * statediff;
+
+ Gjf[j][l][0] += dot(Flux_rhoAtCellFace[K][R], Cjf_loc); // dot(Flux_rho[K], Cjf_loc);
+ for (size_t d = 0; d < Dimension; ++d)
+ Gjf[j][l][1 + d] += u_Cjf[d];
+ Gjf[j][l][1 + Dimension] += dot(Flux_totnrjAtCellFace[K][R], Cjf_loc); // dot(Flux_totnrj[K], Cjf_loc);
+
+ Gjf[j][l] += diff;
+ }
+
+ Gjf[j][l] *= 1. / face_to_cell.size();
+ }
+ });
+
+ if (checkLocalConservation) {
+ auto is_boundary_face = p_mesh->connectivity().isBoundaryFace();
+
+ FaceValue<double> MaxErrorConservationFace(p_mesh->connectivity());
+ MaxErrorConservationFace.fill(0.);
+
+ parallel_for(
+ p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) {
+ const auto& face_to_cell = face_to_cell_matrix[l];
+ const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(l);
+
+ if (not is_boundary_face[l]) {
+ Rp SumGjf(zero);
+ for (size_t k = 0; k < face_to_cell.size(); ++k) {
+ const CellId K = face_to_cell[k];
+ const unsigned int R = face_local_number_in_its_cells[k];
+ SumGjf += Gjf[K][R];
+ }
+ MaxErrorConservationFace[l] = l2Norm(SumGjf);
+ // MaxErrorConservationFace = std::max(MaxErrorConservationFace, l2Norm(SumGjf));
+ }
+ });
+ std::cout << " Max Error Face " << max(MaxErrorConservationFace) << "\n";
+ }
+
+ if constexpr (Dimension == 3) {
+ const auto& edge_to_cell_matrix = p_mesh->connectivity().edgeToCellMatrix();
+ const auto& edge_local_numbers_in_their_cells = p_mesh->connectivity().edgeLocalNumbersInTheirCells();
+
+ const EdgeValuePerCell<const Rd> Cje = mesh_data.Cje();
+ const EdgeValuePerCell<const Rd> nje = mesh_data.nje();
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // Edge
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ const EdgeId& edge = cell_to_edge[l];
+ const auto& edge_to_cell = edge_to_cell_matrix[edge];
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge);
+
+ const Rd& Cje_loc = Cje(j, l);
+
+ for (size_t k = 0; k < edge_to_cell.size(); ++k) {
+ const CellId K = edge_to_cell[k];
+ const size_t R = edge_local_number_in_its_cells[k];
+
+ const Rd& Cke_loc = Cje(K, R);
+
+ // Une moyenne entre les etats jk
+
+ Rd uEdge = .5 * (u_n[j] + u_n[K]);
+ double cEdge = .5 * (c_n[j] + c_n[K]);
+
+ // Viscosity j k
+ Rpxp ViscosityMatrixJK(identity);
+ const double MaxmaxabsVpNormjk =
+ std::max(toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uEdge, cEdge,
+ Cje_loc),
+ toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uEdge, cEdge,
+ Cke_loc));
+
+ ViscosityMatrixJK *= MaxmaxabsVpNormjk;
+
+ const Rd& u_Cje = Flux_qtmvtAtCellEdge[K][R] * Cje_loc; // Flux_qtmvt[K] * Cje_loc;
+
+ const Rp& statediff = StateAtEdge[j][l] - StateAtEdge[K][R]; // State[j] - State[K];
+ const Rp& diff = ViscosityMatrixJK * statediff;
+
+ Gje[j][l][0] += dot(Flux_rhoAtCellEdge[K][R], Cje_loc); // dot(Flux_rho[K], Cje_loc);
+ for (size_t d = 0; d < Dimension; ++d)
+ Gje[j][l][1 + d] += u_Cje[d];
+ Gje[j][l][1 + Dimension] += dot(Flux_totnrjAtCellEdge[K][R], Cje_loc); // dot(Flux_totnrj[K], Cje_loc);
+
+ Gje[j][l] += diff;
+ }
+
+ Gje[j][l] *= 1. / edge_to_cell.size();
+ }
+ });
+
+ if (checkLocalConservation) {
+ auto is_boundary_edge = p_mesh->connectivity().isBoundaryEdge();
+
+ EdgeValue<double> MaxErrorConservationEdge(p_mesh->connectivity());
+ MaxErrorConservationEdge.fill(0.);
+ // double MaxErrorConservationEdge = 0;
+ parallel_for(
+ p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId l) {
+ const auto& edge_to_cell = edge_to_cell_matrix[l];
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(l);
+
+ if (not is_boundary_edge[l]) {
+ Rp SumGje(zero);
+ for (size_t k = 0; k < edge_to_cell.size(); ++k) {
+ const CellId K = edge_to_cell[k];
+ const unsigned int R = edge_local_number_in_its_cells[k];
+ SumGje += Gje[K][R];
+ }
+ // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2norm(SumGje));
+ MaxErrorConservationEdge[l] = l2Norm(SumGje);
+ }
+ });
+ std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << "\n";
+ }
+ } // dim 3
+
+ // Pour les assemblages
+ double theta = .5;
+ double eta = .2;
+ if constexpr (Dimension == 2) {
+ eta = 0;
+ } else {
+ theta = 1. / 3.;
+ eta = 1. / 3.;
+
+ // theta = .5;
+ // eta = 0;
+ }
+ //
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_to_node = cell_to_node_matrix[j];
+
+ Rp SumFluxesNode(zero);
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ SumFluxesNode += Gjr[j][l];
+ }
+ // SumFluxesNode *= (1 - theta);
+
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+ Rp SumFluxesEdge(zero);
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ SumFluxesEdge += Gje[j][l];
+ }
+
+ const auto& cell_to_face = cell_to_face_matrix[j];
+ Rp SumFluxesFace(zero);
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ SumFluxesFace += Gjf[j][l];
+ }
+ // SumFluxesEdge *= (theta);
+
+ const Rp SumFluxes = (1 - theta - eta) * SumFluxesNode + eta * SumFluxesEdge + theta * SumFluxesFace;
+
+ State[j] -= dt / volumes_cell[j] * SumFluxes;
+
+ rho[j] = State[j][0];
+ for (size_t d = 0; d < Dimension; ++d)
+ u[j][d] = State[j][1 + d] / rho[j];
+ E[j] = State[j][1 + Dimension] / rho[j];
+ });
+
+ return std::make_tuple(std::make_shared<const DiscreteFunctionVariant>(rho),
+ std::make_shared<const DiscreteFunctionVariant>(u),
+ std::make_shared<const DiscreteFunctionVariant>(E));
+ }
+
+ RusanovEulerianCompositeSolver_v2() = default;
+ ~RusanovEulerianCompositeSolver_v2() = default;
+};
+
+template <MeshConcept MeshType>
+class RusanovEulerianCompositeSolver_v2<MeshType>::NeumannBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver_v2<Mesh<2>>::NeumannBoundaryCondition
+{
+ 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();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary)
+ : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary)
+ {
+ ;
+ }
+};
+
+template <>
+class RusanovEulerianCompositeSolver_v2<Mesh<3>>::NeumannBoundaryCondition
+{
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const MeshEdgeBoundary m_mesh_edge_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
+
+ public:
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_node_boundary.nodeList().size();
+ }
+ size_t
+ numberOfEdges() const
+ {
+ return m_mesh_edge_boundary.edgeList().size();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const EdgeId>&
+ edgeList() const
+ {
+ return m_mesh_edge_boundary.edgeList();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
+ const MeshEdgeBoundary& mesh_edge_boundary,
+ const MeshFaceBoundary& mesh_face_boundary)
+ : m_mesh_node_boundary(mesh_node_boundary),
+
+ m_mesh_edge_boundary(mesh_edge_boundary),
+
+ m_mesh_face_boundary(mesh_face_boundary)
+ {
+ ;
+ }
+};
+
+template <MeshConcept MeshType>
+class RusanovEulerianCompositeSolver_v2<MeshType>::SymmetryBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver_v2<Mesh<2>>::SymmetryBoundaryCondition
+{
+ public:
+ 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();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_flat_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_flat_node_boundary.nodeList();
+ }
+
+ 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 <>
+class RusanovEulerianCompositeSolver_v2<Mesh<3>>::SymmetryBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
+ const MeshFlatEdgeBoundary<MeshType> m_mesh_flat_edge_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();
+ }
+
+ size_t
+ numberOfEdges() const
+ {
+ return m_mesh_flat_edge_boundary.edgeList().size();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_flat_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_flat_node_boundary.nodeList();
+ }
+
+ const Array<const EdgeId>&
+ edgeList() const
+ {
+ return m_mesh_flat_edge_boundary.edgeList();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_flat_face_boundary.faceList();
+ }
+
+ SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
+ const MeshFlatEdgeBoundary<MeshType>& mesh_flat_edge_boundary,
+ const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
+ : m_mesh_flat_node_boundary(mesh_flat_node_boundary),
+ m_mesh_flat_edge_boundary(mesh_flat_edge_boundary),
+ m_mesh_flat_face_boundary(mesh_flat_face_boundary)
+ {
+ ;
+ }
+
+ ~SymmetryBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class RusanovEulerianCompositeSolver_v2<MeshType>::InflowListBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver_v2<Mesh<2>>::InflowListBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
+ const Table<const double> m_node_array_list;
+ const Table<const double> m_face_array_list;
+
+ public:
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_node_boundary.nodeList().size();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ const Table<const double>&
+ nodeArrayList() const
+ {
+ return m_node_array_list;
+ }
+
+ const Table<const double>&
+ faceArrayList() const
+ {
+ return m_face_array_list;
+ }
+
+ InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
+ const MeshFaceBoundary& mesh_face_boundary,
+ const Table<const double>& node_array_list,
+ const Table<const double>& face_array_list)
+ : m_mesh_node_boundary(mesh_node_boundary),
+ m_mesh_face_boundary(mesh_face_boundary),
+ m_node_array_list(node_array_list),
+ m_face_array_list(face_array_list)
+ {
+ ;
+ }
+
+ ~InflowListBoundaryCondition() = default;
+};
+
+template <>
+class RusanovEulerianCompositeSolver_v2<Mesh<3>>::InflowListBoundaryCondition
+{
+ public:
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const MeshEdgeBoundary m_mesh_edge_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
+ const Table<const double> m_node_array_list;
+ const Table<const double> m_edge_array_list;
+ const Table<const double> m_face_array_list;
+
+ public:
+ size_t
+ numberOfNodes() const
+ {
+ return m_mesh_node_boundary.nodeList().size();
+ }
+
+ size_t
+ numberOfEdges() const
+ {
+ return m_mesh_edge_boundary.edgeList().size();
+ }
+
+ size_t
+ numberOfFaces() const
+ {
+ return m_mesh_face_boundary.faceList().size();
+ }
+
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ const Array<const EdgeId>&
+ edgeList() const
+ {
+ return m_mesh_edge_boundary.edgeList();
+ }
+
+ const Array<const FaceId>&
+ faceList() const
+ {
+ return m_mesh_face_boundary.faceList();
+ }
+
+ const Table<const double>&
+ nodeArrayList() const
+ {
+ return m_node_array_list;
+ }
+
+ const Table<const double>&
+ edgeArrayList() const
+ {
+ return m_edge_array_list;
+ }
+
+ const Table<const double>&
+ faceArrayList() const
+ {
+ return m_face_array_list;
+ }
+
+ InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
+ const MeshEdgeBoundary& mesh_edge_boundary,
+ const MeshFaceBoundary& mesh_face_boundary,
+ const Table<const double>& node_array_list,
+ const Table<const double>& edge_array_list,
+ const Table<const double>& face_array_list)
+ : m_mesh_node_boundary(mesh_node_boundary),
+ m_mesh_edge_boundary(mesh_edge_boundary),
+ m_mesh_face_boundary(mesh_face_boundary),
+ m_node_array_list(node_array_list),
+ m_edge_array_list(edge_array_list),
+ m_face_array_list(face_array_list)
+ {
+ ;
+ }
+
+ ~InflowListBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class RusanovEulerianCompositeSolver_v2<MeshType>::OutflowBoundaryCondition
+{
+};
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+rusanovEulerianCompositeSolver_v2(
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
+ const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const double& dt,
+ const bool check)
+{
+ std::shared_ptr mesh_v = getCommonMesh({rho_v, u_v, E_v, c_v, p_v});
+ if (not mesh_v) {
+ throw NormalError("discrete functions are not defined on the same mesh");
+ }
+
+ if (not checkDiscretizationType({rho_v, u_v, E_v}, DiscreteFunctionType::P0)) {
+ throw NormalError("acoustic solver expects P0 functions");
+ }
+
+ return std::visit(
+ PUGS_LAMBDA(auto&& p_mesh)
+ ->std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ using MeshType = mesh_type_t<decltype(p_mesh)>;
+ static constexpr size_t Dimension = MeshType::Dimension;
+ using Rd = TinyVector<Dimension>;
+
+ if constexpr (Dimension == 1) {
+ throw NormalError("RusanovEulerianCompositeSolver v2 is not available in 1D");
+ } else {
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ return RusanovEulerianCompositeSolver_v2<MeshType>{}
+ .solve(p_mesh, rho_v->get<DiscreteFunctionP0<const double>>(), u_v->get<DiscreteFunctionP0<const Rd>>(),
+ E_v->get<DiscreteFunctionP0<const double>>(), gamma, c_v->get<DiscreteFunctionP0<const double>>(),
+ p_v->get<DiscreteFunctionP0<const double>>(), bc_descriptor_list, dt, check);
+ } else {
+ throw NormalError("RusanovEulerianCompositeSolver v2 is only defined on polygonal meshes");
+ }
+ }
+ },
+ mesh_v->variant());
+}
+
+template <MeshConcept MeshType>
+void
+RusanovEulerianCompositeSolver_v2<MeshType>::_applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const
+{
+ for (const auto& boundary_condition : bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<OutflowBoundaryCondition, T>) {
+ throw NormalError("Not implemented");
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <MeshConcept MeshType>
+void
+RusanovEulerianCompositeSolver_v2<MeshType>::_applySymmetricBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const
+{
+ for (const auto& boundary_condition : bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement SYMMETRY \n";
+ const Rd& normal = bc.outgoingNormal();
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+ const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
+ // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+ const auto& face_list = bc.faceList();
+ const auto& node_list = bc.nodeList();
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+
+ const auto& node_cell_list = node_to_cell_matrix[node_id];
+ // Assert(face_cell_list.size() == 1);
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+
+ for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+ CellId node_cell_id = node_cell_list[i_cell];
+ size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim];
+ // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim];
+ }
+ }
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+
+ const auto& face_cell_list = face_to_cell_matrix[face_id];
+ Assert(face_cell_list.size() == 1);
+
+ CellId face_cell_id = face_cell_list[0];
+ size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+
+ if constexpr (Dimension == 3) {
+ const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();
+
+ const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();
+
+ const auto& edge_list = bc.edgeList();
+
+ for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
+ const EdgeId edge_id = edge_list[i_edge];
+
+ const auto& edge_cell_list = edge_to_cell_matrix[edge_id];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
+ CellId edge_cell_id = edge_cell_list[i_cell];
+ size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];
+
+ Rd vectorSym(zero);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim];
+
+ Rdxd MatriceProj(identity);
+ MatriceProj -= tensorProduct(normal, normal);
+ vectorSym = MatriceProj * vectorSym;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+}
+
+template <MeshConcept MeshType>
+void
+RusanovEulerianCompositeSolver_v2<MeshType>::_applyInflowBoundaryCondition(const BoundaryConditionList& bc_list,
+ const MeshType& mesh,
+ NodeValuePerCell<Rp>& stateNode,
+ EdgeValuePerCell<Rp>& stateEdge,
+ FaceValuePerCell<Rp>& stateFace) const
+{
+ for (const auto& boundary_condition : bc_list) {
+ std::visit(
+ [&](auto&& bc) {
+ using T = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<InflowListBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement INFLOW \n";
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+ const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
+ // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+ const auto& face_list = bc.faceList();
+ const auto& node_list = bc.nodeList();
+
+ const auto& face_array_list = bc.faceArrayList();
+ const 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];
+
+ const auto& node_cell_list = node_to_cell_matrix[node_id];
+ // Assert(face_cell_list.size() == 1);
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+
+ for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+ CellId node_cell_id = node_cell_list[i_cell];
+ size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateNode[node_cell_id][node_local_number_in_cell][dim] = node_array_list[i_node][dim];
+ }
+ }
+
+ for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+ const FaceId face_id = face_list[i_face];
+
+ const auto& face_cell_list = face_to_cell_matrix[face_id];
+ Assert(face_cell_list.size() == 1);
+
+ CellId face_cell_id = face_cell_list[0];
+ size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][dim];
+ }
+
+ if constexpr (Dimension == 3) {
+ const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();
+
+ const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();
+ // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
+
+ const auto& edge_list = bc.edgeList();
+
+ const auto& edge_array_list = bc.edgeArrayList();
+
+ for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
+ const EdgeId edge_id = edge_list[i_edge];
+
+ const auto& edge_cell_list = edge_to_cell_matrix[edge_id];
+ // Assert(face_cell_list.size() == 1);
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
+
+ for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
+ CellId edge_cell_id = edge_cell_list[i_cell];
+ size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];
+
+ for (size_t dim = 0; dim < Dimension + 2; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+}
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.hpp b/src/scheme/RusanovEulerianCompositeSolver_v2.hpp
new file mode 100644
index 000000000..749c09a18
--- /dev/null
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2.hpp
@@ -0,0 +1,29 @@
+#ifndef RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+#define RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+
+#include <memory>
+#include <vector>
+
+// double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+// const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+rusanovEulerianCompositeSolver_v2(
+ const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E,
+ const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const double& dt,
+ const bool check = false);
+
+#endif // RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
--
GitLab