From b4bf4a1cc12f6fac69ed7c5f89d78c73a04b37f9 Mon Sep 17 00:00:00 2001
From: HOCH PHILIPPE <philippe.hoch@gmail.com>
Date: Fri, 23 Aug 2024 05:57:15 +0200
Subject: [PATCH] True Adding second order in space for 3 current solver with
induced limitation
---
src/scheme/CellbyCellLimitation.hpp | 547 ++++
...scousFormEulerianCompositeSolver_v2_o2.cpp | 2519 +++++++++++++++++
...scousFormEulerianCompositeSolver_v2_o2.hpp | 31 +
.../RusanovEulerianCompositeSolver_o2.cpp | 2085 ++++++++++++++
.../RusanovEulerianCompositeSolver_o2.hpp | 30 +
.../RusanovEulerianCompositeSolver_v2_o2.cpp | 1995 +++++++++++++
.../RusanovEulerianCompositeSolver_v2_o2.hpp | 30 +
7 files changed, 7237 insertions(+)
create mode 100644 src/scheme/CellbyCellLimitation.hpp
create mode 100644 src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
create mode 100644 src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.hpp
create mode 100644 src/scheme/RusanovEulerianCompositeSolver_o2.cpp
create mode 100644 src/scheme/RusanovEulerianCompositeSolver_o2.hpp
create mode 100644 src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
create mode 100644 src/scheme/RusanovEulerianCompositeSolver_v2_o2.hpp
diff --git a/src/scheme/CellbyCellLimitation.hpp b/src/scheme/CellbyCellLimitation.hpp
new file mode 100644
index 000000000..5369c4fd6
--- /dev/null
+++ b/src/scheme/CellbyCellLimitation.hpp
@@ -0,0 +1,547 @@
+#ifndef CELL_BY_CELL_LIMITATION_HPP
+#define CELL_BY_CELL_LIMITATION_HPP
+
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.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/DiscreteFunctionDPk.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionDPkVector.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/PugsTraits.hpp>
+#include <variant>
+
+template <MeshConcept MeshTypeT>
+class CellByCellLimitation
+{
+ private:
+ using MeshType = MeshTypeT;
+ static constexpr size_t Dimension = MeshType::Dimension;
+
+ using Rdxd = TinyMatrix<Dimension>;
+ using Rd = TinyVector<Dimension>;
+
+ public:
+ void
+ computeLimitorVolumicScalarQuantityMinModDukowiczGradient(const MeshType& mesh,
+ const CellValue<double>& q,
+ const CellValue<Rd>& gradq,
+ CellValue<double>& Limitor_q,
+ const bool enableWeakBoundPositivityOnly = false) const
+ {
+ Assert(Dimension == 2 or Dimension == 3);
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+
+ const auto& xr = mesh.xr();
+
+ MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ // auto volumes_cell = mesh_data.Vj();
+
+ const auto& xcentroid = mesh_data.xj();
+ const auto& xface = mesh_data.xl();
+
+ if (Dimension == 2) {
+ parallel_for(
+ mesh.numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+ const double State = q[j];
+ const Rd& Cent = xcentroid[j]; // mesh.xi(j);
+ const Rd& Gradl = gradq[j];
+
+ // Min et Max des valeurs voisines
+
+ double minvalvois(State);
+ double maxvalvois(State);
+
+ // Calcul des extremas de la fonction
+
+ double minval(std::numeric_limits<double>::max());
+ double maxval(-std::numeric_limits<double>::max());
+
+ // At nodes
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ // At faces
+ const auto& cell_to_face = cell_to_face_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];
+
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const CellId& cellvois = node_to_cell[k];
+ const double statevoiscell = q[cellvois];
+ const Rd Centvois = xcentroid[cellvois]; // mesh.element()[mayvois].centroid().x();
+ const Rd Mvois(Centvois - Cent);
+
+ const double valVois = State + dot(Gradl, Mvois); // Partie MinMod
+
+ minval = std::min(minval, valVois);
+ maxval = std::max(maxval, valVois);
+
+ minvalvois = std::min(minvalvois, statevoiscell);
+ maxvalvois = std::max(maxvalvois, statevoiscell);
+ }
+
+ const Rd M(xr[node] - Cent);
+
+ const double valLineO2node = State + dot(Gradl, M); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2node);
+ maxval = std::max(maxval, valLineO2node);
+ }
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const FaceId& face = cell_to_face[l];
+
+ const Rd Mb(xface[face] - Cent);
+
+ const double valLineO2bras = State + dot(Gradl, Mb); // Idem aux demi faces
+
+ minval = std::min(minval, valLineO2bras);
+ maxval = std::max(maxval, valLineO2bras);
+ }
+
+ static const double epsZer0 = 1e-15;
+ // on applique le traitement idem de l'ordre 2..
+ double coef1 = 0, coef2 = 0;
+ if (enableWeakBoundPositivityOnly) {
+ coef1 = 1;
+ const double minGlobal = 1e-12;
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minGlobal - State) / ((minval - State));
+
+ } else {
+ if (std::fabs(maxval - State) <= epsZer0) // epsIsEqualAbs())
+ coef1 = 1.;
+ else
+ coef1 = (maxvalvois - State) / ((maxval - State));
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minvalvois - State) / ((minval - State));
+ }
+ const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+ Limitor_q[j] = alfa;
+ // grad(i) *= alfa;
+ });
+ } else {
+ // throw NormalError(" Limiteur Scal Grad non encore implemente 3D");
+ const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
+ const auto& xedge = mesh_data.xe();
+
+ parallel_for(
+ mesh.numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+ const double State = q[j];
+ const Rd& Cent = xcentroid[j]; // mesh.xi(j);
+ const Rd& Gradl = gradq[j];
+
+ // Min et Max des valeurs voisines
+
+ double minvalvois(State);
+ double maxvalvois(State);
+
+ // Calcul des extremas de la fonction
+
+ double minval(std::numeric_limits<double>::max());
+ double maxval(-std::numeric_limits<double>::max());
+
+ // At nodes
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ // At faces
+ const auto& cell_to_face = cell_to_face_matrix[j];
+ // At faces
+ const auto& cell_to_edge = cell_to_edge_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];
+
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const CellId& cellvois = node_to_cell[k];
+ const double statevoiscell = q[cellvois];
+ const Rd Centvois = xcentroid[cellvois]; // mesh.element()[mayvois].centroid().x();
+ const Rd Mvois(Centvois - Cent);
+
+ const double valVois = State + dot(Gradl, Mvois); // Partie MinMod
+
+ minval = std::min(minval, valVois);
+ maxval = std::max(maxval, valVois);
+
+ minvalvois = std::min(minvalvois, statevoiscell);
+ maxvalvois = std::max(maxvalvois, statevoiscell);
+ }
+
+ const Rd M(xr[node] - Cent);
+
+ const double valLineO2node = State + dot(Gradl, M); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2node);
+ maxval = std::max(maxval, valLineO2node);
+ }
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const FaceId& face = cell_to_face[l];
+
+ const Rd Mb(xface[face] - Cent);
+
+ const double valLineO2face = State + dot(Gradl, Mb); // Idem aux demi faces
+
+ minval = std::min(minval, valLineO2face);
+ maxval = std::max(maxval, valLineO2face);
+ }
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ const EdgeId& edge = cell_to_edge[l];
+
+ const Rd Me(xedge[edge] - Cent);
+
+ const double valLineO2edge = State + dot(Gradl, Me); // Idem aux demi faces
+
+ minval = std::min(minval, valLineO2edge);
+ maxval = std::max(maxval, valLineO2edge);
+ }
+
+ static const double epsZer0 = 1e-15;
+ // on applique le traitement idem de l'ordre 2..
+ double coef1 = 0, coef2 = 0;
+ if (enableWeakBoundPositivityOnly) {
+ coef1 = 1;
+ const double minGlobal = 1e-12;
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minGlobal - State) / ((minval - State));
+
+ } else {
+ if (std::fabs(maxval - State) <= epsZer0) // epsIsEqualAbs())
+ coef1 = 1.;
+ else
+ coef1 = (maxvalvois - State) / ((maxval - State));
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minvalvois - State) / ((minval - State));
+ }
+ const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+ Limitor_q[j] = alfa;
+ // grad(i) *= alfa;
+ });
+ }
+ }
+
+ void
+ computeLimitorInternalNrjScalarQuantityMinModDukowiczGradient(const MeshType& mesh,
+ const CellValue<double>& rho,
+ const CellValue<Rd>& gradrho,
+ const CellValue<double>& internalnrj,
+ const CellValue<Rd>& gradinternalnrj,
+ // const CellValue<Rd>& u,
+ const CellValue<Rdxd>& gradu,
+ CellValue<double>& Limitor_eps,
+ const bool enableWeakBoundPositivityOnly = false) const
+ {
+ Assert(Dimension == 2 or Dimension == 3);
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& xr = mesh.xr();
+
+ MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const auto& xcentroid = mesh_data.xj();
+ const auto& xface = mesh_data.xl();
+
+ if (Dimension == 2) {
+ parallel_for(
+ mesh.numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+ const double rhol = rho[j];
+ const Rd& Cent = xcentroid[j];
+ const Rd& grad_rhol = gradrho[j];
+
+ // const Rd& ul = u[j];
+ const Rdxd& grad_ul = gradu[j];
+
+ const double internal_nrj = internalnrj[j];
+ // const double kinetic_nrj = .5 * dot(ul, ul);
+
+ const Rd& grad_internal_nrj = gradinternalnrj[j];
+
+ // const Rd& grad_kinetic_internal_nrj = gradeps[j];
+
+ // const double qtt_limitor_density = rhol / (rhol + dot(grad,(x-xc));
+ // At nodes
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ // At faces
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ // Min et Max des valeurs voisines
+ const double State = internal_nrj;
+
+ double minvalvois(State);
+ double maxvalvois(State);
+
+ // Calcul des extremas de la fonction
+
+ double minval(std::numeric_limits<double>::max());
+ double maxval(-std::numeric_limits<double>::max());
+
+ 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];
+
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const CellId& cellvois = node_to_cell[k];
+ const Rd Centvois = xcentroid[cellvois];
+ const Rd Mvois(Centvois - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mvois));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mvois);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mvois);
+
+ const double valVois = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie MinMod
+
+ minval = std::min(minval, valVois);
+ maxval = std::max(maxval, valVois);
+
+ const double statevoiscell = internalnrj[cellvois];
+
+ minvalvois = std::min(minvalvois, statevoiscell);
+ maxvalvois = std::max(maxvalvois, statevoiscell);
+ }
+
+ const Rd Mr(xr[node] - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mr));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mr);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mr);
+
+ const double valLineO2node = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2node);
+ maxval = std::max(maxval, valLineO2node);
+ }
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const FaceId& face = cell_to_face[l];
+
+ const Rd Mb(xface[face] - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mb));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mb);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mb);
+
+ const double valLineO2face = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2face);
+ maxval = std::max(maxval, valLineO2face);
+ }
+
+ static const double epsZer0 = 1e-15;
+
+ // on applique le traitement idem de l'ordre 2..
+ double coef1 = 0, coef2 = 0;
+ if (enableWeakBoundPositivityOnly) {
+ coef1 = 1;
+ const double minGlobal = 1e-12;
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minGlobal - State) / ((minval - State));
+
+ } else {
+ if (std::fabs(maxval - State) <= epsZer0) // epsIsEqualAbs())
+ coef1 = 1.;
+ else
+ coef1 = (maxvalvois - State) / ((maxval - State));
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minvalvois - State) / ((minval - State));
+ }
+ const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+ Limitor_eps[j] = alfa;
+ });
+ } else {
+ // throw NormalError(" Limiteur Scal Grad internal nrj non encore implemente 3D");
+ const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
+
+ const auto& xedge = mesh_data.xe();
+
+ parallel_for(
+ mesh.numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+ const double rhol = rho[j];
+ const Rd& Cent = xcentroid[j];
+ const Rd& grad_rhol = gradrho[j];
+
+ // const Rd& ul = u[j];
+ const Rdxd& grad_ul = gradu[j];
+
+ const double internal_nrj = internalnrj[j];
+ // const double kinetic_nrj = .5 * dot(ul, ul);
+
+ const Rd& grad_internal_nrj = gradinternalnrj[j];
+
+ // const Rd& grad_kinetic_internal_nrj = gradeps[j];
+
+ // const double qtt_limitor_density = rhol / (rhol + dot(grad,(x-xc));
+ // At nodes
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ // At faces
+ const auto& cell_to_face = cell_to_face_matrix[j];
+ // At edges
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+ // Min et Max des valeurs voisines
+ const double State = internal_nrj;
+
+ double minvalvois(State);
+ double maxvalvois(State);
+
+ // Calcul des extremas de la fonction
+
+ double minval(std::numeric_limits<double>::max());
+ double maxval(-std::numeric_limits<double>::max());
+
+ 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];
+
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const CellId& cellvois = node_to_cell[k];
+ const Rd Centvois = xcentroid[cellvois];
+ const Rd Mvois(Centvois - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mvois));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mvois);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mvois);
+
+ const double valVois = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie MinMod
+
+ minval = std::min(minval, valVois);
+ maxval = std::max(maxval, valVois);
+
+ const double statevoiscell = internalnrj[cellvois];
+
+ minvalvois = std::min(minvalvois, statevoiscell);
+ maxvalvois = std::max(maxvalvois, statevoiscell);
+ }
+
+ const Rd Mr(xr[node] - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mr));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mr);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mr);
+
+ const double valLineO2node = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2node);
+ maxval = std::max(maxval, valLineO2node);
+ }
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const FaceId& face = cell_to_face[l];
+
+ const Rd Mb(xface[face] - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mb));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mb);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mb);
+
+ const double valLineO2face = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2face);
+ maxval = std::max(maxval, valLineO2face);
+ }
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ const EdgeId& edge = cell_to_edge[l];
+
+ const Rd Me(xedge[edge] - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Me));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Me);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Me);
+
+ const double valLineO2edge = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2edge);
+ maxval = std::max(maxval, valLineO2edge);
+ }
+
+ static const double epsZer0 = 1e-15;
+
+ // on applique le traitement idem de l'ordre 2..
+ double coef1 = 0, coef2 = 0;
+ if (enableWeakBoundPositivityOnly) {
+ coef1 = 1;
+ const double minGlobal = 1e-12;
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minGlobal - State) / ((minval - State));
+
+ } else {
+ if (std::fabs(maxval - State) <= epsZer0) // epsIsEqualAbs())
+ coef1 = 1.;
+ else
+ coef1 = (maxvalvois - State) / ((maxval - State));
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minvalvois - State) / ((minval - State));
+ }
+ const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+ Limitor_eps[j] = alfa;
+ });
+ }
+ }
+};
+#endif
diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
new file mode 100644
index 000000000..84ec1dd2b
--- /dev/null
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
@@ -0,0 +1,2519 @@
+#include <scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.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/DiscreteFunctionDPk.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionDPkVector.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/PugsTraits.hpp>
+#include <variant>
+
+template <MeshConcept MeshTypeT>
+class RoeViscousFormEulerianCompositeSolver_v2_o2
+{
+ 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>;
+
+ using Rpxd = TinyMatrix<Dimension + 2, Dimension>;
+
+ class SymmetryBoundaryCondition;
+ class InflowListBoundaryCondition;
+ class OutflowBoundaryCondition;
+ class NeumannBoundaryCondition;
+ class NeumannflatBoundaryCondition;
+
+ using BoundaryCondition = std::variant<SymmetryBoundaryCondition,
+ InflowListBoundaryCondition,
+ OutflowBoundaryCondition,
+ NeumannflatBoundaryCondition,
+ 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::neumann: {
+ if constexpr (Dimension == 2) {
+ bc_list.emplace_back(
+ NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ } else {
+ static_assert(Dimension == 3);
+ bc_list.emplace_back(
+ NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ }
+ break;
+ }
+ 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: {
+ if constexpr (Dimension == 2) {
+ bc_list.emplace_back(
+ OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ } else {
+ static_assert(Dimension == 3);
+ bc_list.emplace_back(
+ OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ }
+ break;
+ // 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:
+ CellByCellLimitation<MeshType> Limitor;
+
+ public:
+ void
+ _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>) {
+ std::cout << " Traitement Outflow \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();
+
+ const auto xj = mesh.xj();
+ const auto xr = mesh.xr();
+ const auto xf = mesh.xl();
+ const auto xe = mesh.xe();
+
+ 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] += vectorSym[dim];
+
+ 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];
+ }
+ }
+
+ // throw NormalError("Not implemented");
+ }
+ */
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _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);
+ }
+ }
+
+ void
+ _applyNeumannflatBoundaryCondition(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<NeumannflatBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement WALL \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];
+
+ vectorSym -= dot(vectorSym, normal) * normal;
+
+ 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];
+
+ vectorSym -= dot(vectorSym, normal) * normal;
+
+ 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];
+
+ vectorSym -= dot(vectorSym, normal) * normal;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _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);
+ }
+ }
+
+ public:
+ inline double
+ pression(const double rho, const double epsilon, const double gam) const
+ {
+ return (gam - 1) * rho * epsilon;
+ }
+
+ inline Rpxd
+ Flux(const double& rho, const Rd& U, const double& E, const double gam) const
+ {
+ // const R2 flux_rho = rhoU;
+ // const R22 flux_rhoU = R22(rhoU.x1() * rhoU.x1() / rho + P, rhoU.x1() * rhoU.x2() / rho, rhoU.x2() * rhoU.x1()
+ // rho,rhoU.x2() * rhoU.x2() / rho + P);
+ // const R2 flux_rhoE = ((rhoE + P) / rho) * rhoU;
+ /* 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;
+
+ */
+ const Rd& rhoU = rho * U;
+ const Rdxd rhoUTensU = tensorProduct(rhoU, U);
+ const double p = pression(rho, E - .5 * dot(U, U), gam);
+
+ Rdxd pid(identity);
+ pid *= p;
+ const Rdxd rhoUTensUPlusPid = rhoUTensU + pid;
+ const double rhoEPlusP = rho * E + p;
+
+ const Rd& rhoEPlusPtimesU = rhoEPlusP * U;
+
+ Rpxd Fluxx; // En ligne ci dessous
+
+ Fluxx[0] = rhoU;
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ // for (size_t dim2 = 0; dim2 < Dimension; ++dim2)
+ Fluxx[1 + dim] = rhoUTensUPlusPid[dim];
+ // Fluxx[1 + dim][dim2] = rhoUTensUPlusPid[dim][dim2];
+ Fluxx[1 + Dimension] = rhoEPlusPtimesU;
+ return Fluxx;
+ }
+
+ struct JacobianInformations
+ {
+ Rpxp Jacobian;
+ Rpxp LeftEigenVectors;
+ Rpxp RightEigenVectors;
+ Rp EigenValues;
+ Rpxp AbsJacobian;
+ double maxabsVpNormal;
+ // bool changingSign;
+ double MaxErrorProd;
+
+ public:
+ JacobianInformations(const Rpxp& Jacobiand,
+ const Rpxp& LeftEigenVectorsd,
+ const Rpxp& RightEigenVectorsd,
+ const Rp& EigenValuesd,
+ const Rpxp& AbsJacobiand,
+ const double maxabsVpNormald,
+ const double MaxErrorProdd)
+ : Jacobian(Jacobiand),
+ LeftEigenVectors(LeftEigenVectorsd),
+ RightEigenVectors(RightEigenVectorsd),
+ EigenValues(EigenValuesd),
+ AbsJacobian(AbsJacobiand),
+ maxabsVpNormal(maxabsVpNormald),
+ MaxErrorProd(MaxErrorProdd)
+ {}
+ ~JacobianInformations() {}
+ };
+
+ Rp
+ EvaluateEigenValuesNormal( // const double rhod, // rhoJ, uJ, EJ, gammaJ, cJ, pJ, normal
+ const Rd& Ud,
+ // const double Ed,
+ // const double gammad,
+ const double cd,
+ // const double pd,
+ const Rd& normal) const
+ {
+ Rp Eigen;
+ const double uscaln = dot(Ud, normal);
+
+ Eigen[0] = uscaln - cd;
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Eigen[dim] = uscaln;
+ Eigen[1 + Dimension] = uscaln + cd;
+
+ return Eigen;
+ }
+
+ struct RoeAverageStateStructData
+ {
+ double rho;
+ Rd U;
+ double E;
+ double H;
+ double gamma;
+ double p;
+ double c;
+ RoeAverageStateStructData(const double rhod,
+ const Rd Ud,
+ const double Ed,
+ const double Hd,
+ const double gammad,
+ const double pd,
+ const double cd)
+ : rho(rhod), U(Ud), E(Ed), H(Hd), gamma(gammad), p(pd), c(cd)
+ {}
+ };
+
+ RoeAverageStateStructData
+ RoeAverageState(const double& rhoG,
+ const Rd& UG,
+ const double& EG,
+ const double& gammaG,
+ const double& pG,
+ const double& rhoD,
+ const Rd& UD,
+ const double& ED,
+ const double gammaD,
+ const double pD) const
+
+ {
+ double gamma = .5 * (gammaG + gammaD); // ou ponderation racine roG et roD
+ double RacineRoG = sqrt(rhoG);
+ double RacineRoD = sqrt(rhoD);
+ double rho_mean = RacineRoG * RacineRoD;
+ Rd U_mean = 1. / (RacineRoG + RacineRoD) * (RacineRoG * UG + RacineRoD * UD);
+ double unorm2 = dot(U_mean, U_mean);
+ double NrjCin = .5 * unorm2;
+ const double TotEnthalpyG = (EG + pG / rhoG);
+ const double TotEnthalpyD = (ED + pD / rhoD);
+
+ double H_mean = (RacineRoG * TotEnthalpyG + RacineRoD * TotEnthalpyD) / (RacineRoG + RacineRoD);
+
+ double E_mean = H_mean / gamma + ((gamma - 1) / (gamma)) * (NrjCin);
+
+ double P_mean = rho_mean * (H_mean - E_mean);
+
+ double c2 = gamma * P_mean / rho_mean; // cspeed*cspeed;
+ // assert(fabs((gamma - 1) * rho_mean * (E_mean - .5 * (u_mean, u_mean)) - P_mean) < 1e-13); // equilibre GP
+ double c_mean = sqrt(c2); // cspeed_meandof/Area;
+
+ return RoeAverageStateStructData(rho_mean, U_mean, E_mean, H_mean, gamma, P_mean, c_mean);
+ }
+
+ bool
+ EvaluateChangingSignVpAlongNormal( // const double rhoJ,
+ const Rd& uJ,
+ // const double EJ,
+ // const double gammaJ,
+ const double cJ,
+ // const double pJ,
+
+ // const double rhoK,
+ const Rd& uK,
+ // const double EK,
+ // const double gammaK,
+ const double cK,
+ // const double pK,
+ const Rd& normal) const
+ {
+ const Rp& EigenJ = EvaluateEigenValuesNormal( // rhoJ,
+ uJ, // EJ, gammaJ,
+ cJ, // pJ,
+ normal);
+
+ const Rp& EigenK = EvaluateEigenValuesNormal( // rhoK,
+ uK, // EK, gammaK,
+ cK, // pK,
+ normal);
+
+ if (EigenJ[0] * EigenK[0] < -1e-12) // <=0) //< -1e-12)
+ return true;
+ for (size_t dim = 1; dim < Dimension + 1; ++dim) // vp multiplicite d
+ if (EigenJ[dim] * EigenK[dim] < -1e-12) // <=0) //< -1e-12)
+ return true;
+ if (EigenJ[1 + Dimension] * EigenK[1 + Dimension] < -1e-12) // <=0) //< -1e-12)
+ return true;
+
+ return false;
+ }
+
+ JacobianInformations
+ JacobianFluxAlongUnitNormal(const RoeAverageStateStructData& RoeState,
+ /*
+ const double rhoJ,
+ const Rd& uJ,
+ const double EJ,
+ const double gammaJ,
+ const double cJ,
+ const double pJ,
+
+ const double rhoK,
+ const Rd& uK,
+ const double EK,
+ const double gammaK,
+ const double cK,
+ const double pK,
+*/
+ const Rd& normal,
+ const bool check = false) const
+ {
+ Assert((l2Norm(normal) - 1) < 1e-12);
+
+ // const double& rho = RoeState.rho;
+ const Rd& u_mean = RoeState.U;
+ const double H_mean = RoeState.H;
+ const double& cspeed = RoeState.c;
+ const double& gamma = RoeState.gamma;
+ const double& uscaln = dot(u_mean, normal);
+ const double& u2 = dot(u_mean, u_mean);
+ // const R NrjCin=.5*unorm2;
+
+ const double c2 = cspeed * cspeed;
+ const double gm1 = gamma - 1;
+ const double K = c2 + gm1 * (dot(u_mean, u_mean) - H_mean);
+
+ // Le jacobien est lineaire par rapport a la normale
+ // Ref : PAr ex. le papier de D. Chauvheid sur la tension de surface
+ Rpxp Jacobian;
+
+ Rdxd CentralT = identity;
+ CentralT *= uscaln;
+ CentralT += tensorProduct(u_mean, normal) - gm1 * tensorProduct(normal, u_mean);
+
+ Jacobian(0, 0) = 0.;
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ Jacobian(0, dim + 1) = normal[dim];
+ }
+ Jacobian(0, Dimension + 1) = 0;
+
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ Jacobian(dim + 1, 0) = K * normal[dim] - uscaln * u_mean[dim];
+ for (size_t dim2 = 0; dim2 < Dimension; ++dim2)
+ Jacobian(dim + 1, dim2 + 1) = CentralT(dim, dim2);
+ Jacobian(dim + 1, Dimension + 1) = gm1 * normal[dim];
+ }
+
+ Jacobian(Dimension + 1, 0) = (K - H_mean) * uscaln;
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ Jacobian(Dimension + 1, dim + 1) = (H_mean * normal[dim] - gm1 * uscaln * u_mean[dim]);
+ Jacobian(Dimension + 1, Dimension + 1) = gamma * uscaln;
+
+ // Valeur propres..
+ Rp EigenValues;
+
+ EigenValues[0] = uscaln - cspeed;
+ for (size_t dim = 0; dim < Dimension; ++dim) // vp multiplicite d
+ EigenValues[1 + dim] = uscaln;
+ EigenValues[1 + Dimension] = uscaln + cspeed;
+
+ const double maxabsVpNormal = std::max(fabs(EigenValues[0]), fabs(EigenValues[1 + Dimension]));
+
+ // Vecteur propres a droite et gauche
+
+ // hyper plan ortho à la normale
+ std::vector<Rd> ortho(Dimension - 1);
+ if constexpr (Dimension == 2) {
+ ortho[0] = Rd{normal[1], -normal[0]}; // aussi de norme 1
+ } else {
+ const double a = normal[0];
+ const double b = normal[1];
+ const double c = normal[2];
+
+ if ((a == b) and (b == c)) {
+ static double invsqrt2 = 1. / sqrt(2.);
+ static double invsqrt6 = 1. / sqrt(6.);
+
+ ortho[0] = Rd{invsqrt2, -invsqrt2, 0};
+ ortho[1] = Rd{invsqrt6, invsqrt6, -2 * invsqrt6}; // au signe pres
+
+ } else {
+ ortho[0] = Rd{b - c, -(a - c), a - b};
+ ortho[0] *= 1. / l2Norm(ortho[0]);
+ ortho[1] = Rd{a * (b + c) - b * b - c * c, b * (a + c) - a * a - c * c, c * (a + b) - a * a - b * b};
+ ortho[1] *= 1. / l2Norm(ortho[1]); // au signe pres
+ }
+ }
+
+ Rpxp RightTligne;
+ RightTligne(0, 0) = 1;
+ for (size_t dim = 1; dim < Dimension + 1; ++dim)
+ RightTligne(0, dim) = u_mean[dim - 1] - cspeed * normal[dim - 1];
+ RightTligne(0, Dimension + 1) = H_mean - uscaln * cspeed;
+
+ RightTligne(1, 0) = 1;
+ for (size_t dim = 1; dim < Dimension + 1; ++dim)
+ RightTligne(1, dim) = u_mean[dim - 1];
+ RightTligne(1, Dimension + 1) = H_mean - c2 / gm1;
+
+ for (size_t dim = 1; dim < Dimension; ++dim) {
+ RightTligne(dim + 1, 0) = 0.;
+ for (size_t dim2 = 1; dim2 < Dimension + 1; ++dim2) {
+ RightTligne(dim + 1, dim2) = ortho[dim - 1][dim2 - 1];
+ }
+ RightTligne(dim + 1, Dimension + 1) = dot(u_mean, ortho[dim - 1]);
+ }
+
+ RightTligne(Dimension + 1, 0) = 1;
+ for (size_t dim = 1; dim < Dimension + 1; ++dim) {
+ RightTligne(Dimension + 1, dim) = u_mean[dim - 1] + cspeed * normal[dim - 1];
+ }
+ RightTligne(Dimension + 1, Dimension + 1) = H_mean + uscaln * cspeed;
+
+ Rpxp Right = transpose(RightTligne);
+
+ const double invc2 = 1. / c2;
+
+ Rpxp Left; //(zero);
+ Left(0, 0) = .5 * invc2 * (K + uscaln * cspeed);
+ for (size_t dim = 1; dim < Dimension + 1; ++dim)
+ Left(0, dim) = .5 * invc2 * (-gm1 * u_mean[dim - 1] - cspeed * normal[dim - 1]);
+ Left(0, Dimension + 1) = .5 * invc2 * gm1;
+
+ Left(1, 0) = gm1 * invc2 * (H_mean - u2);
+ for (size_t dim = 1; dim < Dimension + 1; ++dim)
+ Left(1, dim) = gm1 * invc2 * u_mean[dim - 1];
+ Left(1, 1 + Dimension) = -gm1 * invc2;
+
+ for (size_t dim = 1; dim < Dimension; ++dim) {
+ Left(1 + dim, 0) = -dot(u_mean, ortho[dim - 1]);
+ for (size_t dim2 = 0; dim2 < Dimension; ++dim2)
+ Left(1 + dim, dim2 + 1) = ortho[dim - 1][dim2];
+ Left(1 + dim, 1 + Dimension) = 0;
+ }
+
+ Left(1 + Dimension, 0) = .5 * invc2 * (K - uscaln * cspeed);
+ for (size_t dim = 1; dim < Dimension + 1; ++dim)
+ Left(1 + Dimension, dim) = .5 * invc2 * (-gm1 * u_mean[dim - 1] + cspeed * normal[dim - 1]);
+ Left(1 + Dimension, Dimension + 1) = .5 * invc2 * gm1;
+
+ // check
+ // Rpxp prod = Right * Left;
+ // std::cout << prod << "\n";
+
+ Rpxp EigenMatAbs(identity);
+ for (size_t dim = 0; dim < Dimension + 2; ++dim) {
+ EigenMatAbs(dim, dim) *= fabs(EigenValues[dim]);
+ }
+ Rpxp ProdLefAbs = Right * EigenMatAbs;
+ Rpxp AbsJacobian = ProdLefAbs * Left;
+ double maxErr = 0;
+
+ if (check) {
+ Rpxp EigenAbs(identity);
+ for (size_t dim = 0; dim < Dimension + 2; ++dim) {
+ EigenAbs(dim, dim) *= EigenValues[dim];
+ }
+ Rpxp ProdLeft = Right * EigenAbs;
+ Rpxp JacobianR = ProdLeft * Left;
+
+ // Rpxp id = identity;
+ maxErr = 0;
+ // double maxErrLeftRightId = 0;
+ for (size_t i = 0; i < Dimension + 2; ++i)
+ for (size_t j = 0; j < Dimension + 2; ++j) {
+ maxErr = std::max(maxErr, fabs(Jacobian(i, j) - JacobianR(i, j)));
+ // maxErrLeftRightId = std::max(maxErrLeftRightId, fabs(id(i, j) - prod(i, j)));
+ }
+ // std::clog << " maxErr " << maxErr << " maxErrUnit " << maxErrLeftRightId;
+ // throw NormalError("STOP");
+ }
+
+ return JacobianInformations(Jacobian, Left, Right, EigenValues, AbsJacobian, maxabsVpNormal, maxErr);
+ }
+
+ 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 size_t& degree,
+ 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);
+ std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
+
+ for (auto&& bc_descriptor : bc_descriptor_list) {
+ if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
+ symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+ }
+ }
+ static const int degre = 1;
+ PolynomialReconstructionDescriptor
+ reconstruction_descriptor(IntegrationMethodType::cell_center, // IntegrationMethodType::boundary,
+ degre, symmetry_boundary_descriptor_list);
+ //
+ //
+ //
+ CellValue<double> valrho({p_mesh->connectivity()});
+ CellValue<Rd> valu({p_mesh->connectivity()});
+ CellValue<double> valeps({p_mesh->connectivity()});
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ valrho[j] = rho[j];
+ valu[j] = u[j];
+ valeps[j] = E[j] - .5 * dot(u[j], u[j]);
+ });
+
+ DiscreteFunctionP0<double> eps(p_mesh, valeps);
+
+ // assemblage direct
+ // auto reconstructions = PolynomialReconstruction{reconstruction_descriptor}.build(rho, rho * u, rho * E);
+
+ // Pour assemblage Leibniz
+ // ca calcul au moins les gradients
+ auto reconstructions = PolynomialReconstruction{reconstruction_descriptor}.build(rho, u, eps);
+
+ DiscreteFunctionDPk rho_bar = reconstructions[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ DiscreteFunctionDPk u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ DiscreteFunctionDPk eps_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ // DiscreteFunctionDPk rho_u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ // DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ CellValue<double> Limitor_rho(p_mesh->connectivity());
+ CellValue<double> Limitor_u(p_mesh->connectivity());
+ CellValue<double> Limitor_eps(p_mesh->connectivity());
+
+ CellValue<Rd> grad_rho_cell(p_mesh->connectivity());
+ CellValue<Rdxd> grad_u_cell(p_mesh->connectivity());
+ CellValue<Rd> grad_eps_cell(p_mesh->connectivity());
+
+ // DiscreteFunctionDPk<Dimension, const double> rho_bar_lim;
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j](x) - rho[j]);
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ grad_rho_cell[j][dim] = rho_bar.coefficients(j)[1 + dim]; // si base : x, puis y , puis z
+ grad_eps_cell[j][dim] = eps_bar.coefficients(j)[1 + dim]; // idem ci dessus
+ const Rd gradu = u_bar.coefficients(j)[1 + dim];
+ for (size_t dim2 = 0; dim2 < Dimension; ++dim2) {
+ grad_u_cell[j](dim, dim2) = gradu[dim2];
+ }
+ }
+ grad_u_cell[j] = transpose(grad_u_cell[j]);
+ });
+
+ // Appels des limiteurs cell by cell
+ // Pour rho (plus ou moins classique)
+ Limitor.computeLimitorVolumicScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, Limitor_rho);
+
+ // Pour les variables d'interet (physique)
+ // Pour eps ici G.P. (avec Leibniz et toute l artillerie)
+ Limitor.computeLimitorInternalNrjScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, valeps,
+ grad_eps_cell,
+ // valu,
+ grad_u_cell, Limitor_eps);
+ // Pour u : relation deduite du principe d'induction
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Limitor_u[j] = sqrt(Limitor_eps[j]); });
+
+ //
+ // Creation des variables conservatives limitées
+ //
+ auto rho_bar_lim = copy(rho_bar);
+ auto rho_u_bar_lim = copy(u_bar); // a modifier .. du a la densite
+ auto rho_E_bar_lim = copy(eps_bar); // a modifier et a completer : densité et energie cinetique
+
+ // Assemblage des ctes et ordres élevés
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j](x) - rho[j]);
+ // rho_u_bar_lim.coefficients(j)[0] = valrho[j] * valu[j]; // * valu[j];
+ // rho_E_bar_lim.coefficients(j)[0] = valrho[j] * (valeps[j] + .5 * dot(valu[j], valu[j]));
+
+ // rho_bar_lim.coefficients(j)[0] // cte inchangée
+ rho_u_bar_lim.coefficients(j)[0] *= valrho[j]; // du coup (rho u)
+ rho_E_bar_lim.coefficients(j)[0] *= valrho[j]; // du coup (rho*epsilon)
+ rho_E_bar_lim.coefficients(j)[0] += valrho[j] * .5 * dot(valu[j], valu[j]);
+ // Rdxd grad_rhou_lim=
+ // valrho[j]*Limitor_u[j]*grad_u_cell[j]+Limitor_rho[j]*tensorProduct(valu[j],grad_rho_cell[j]);
+ //
+ Rdxd gradU_lim_transpose = Limitor_u[j] * transpose(grad_u_cell[j]);
+ const Rd gradUtu_lim = gradU_lim_transpose * valu[j];
+
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ rho_bar_lim.coefficients(j)[1 + dim] *= Limitor_rho[j];
+ rho_u_bar_lim.coefficients(j)[1 + dim] *= Limitor_u[j] * valrho[j]; // c'est la partie gradu limite
+ rho_u_bar_lim.coefficients(j)[1 + dim] += Limitor_rho[j] * grad_rho_cell[j][dim] * valu[j];
+
+ rho_E_bar_lim.coefficients(j)[1 + dim] *= Limitor_eps[j] * valrho[j];
+ rho_E_bar_lim.coefficients(j)[1 + dim] += gradUtu_lim[dim] * valrho[j];
+
+ rho_E_bar_lim.coefficients(j)[1 + dim] +=
+ (valeps[j] + .5 * dot(valu[j], valu[j])) * Limitor_rho[j] * grad_rho_cell[j][dim];
+ }
+ });
+
+ // 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];
+ });
+
+ //
+ // Remplir ici les reconstructions d'ordre élevé
+
+ //
+ 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();
+
+ const auto xr = p_mesh->xr();
+ auto xl = MeshDataManager::instance().getMeshData(*p_mesh).xl();
+ auto xe = MeshDataManager::instance().getMeshData(*p_mesh).xe();
+
+ NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()};
+ StateAtNode.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // StateAtNode[j].fill(State[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];
+
+ StateAtNode[j][l][0] = rho_bar_lim[j](xr[node]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtNode[j][l][1 + dim] = rho_u_bar_lim[j](xr[node])[dim];
+ StateAtNode[j][l][1 + Dimension] = rho_E_bar_lim[j](xr[node]);
+ }
+ });
+
+ EdgeValuePerCell<Rp> StateAtEdge{p_mesh->connectivity()};
+ StateAtEdge.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // eStateAtEdge[j].fill(State[j]);
+ 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];
+
+ StateAtEdge[j][l][0] = rho_bar_lim[j](xe[edge]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtEdge[j][l][1 + dim] = rho_u_bar_lim[j](xe[edge])[dim];
+ StateAtEdge[j][l][1 + Dimension] = rho_E_bar_lim[j](xe[edge]);
+ }
+ });
+ FaceValuePerCell<Rp> StateAtFace{p_mesh->connectivity()};
+ StateAtFace.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // StateAtFace[j].fill(State[j]);
+ 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];
+
+ StateAtFace[j][l][0] = rho_bar_lim[j](xl[face]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtFace[j][l][1 + dim] = rho_u_bar_lim[j](xl[face])[dim];
+ StateAtFace[j][l][1 + Dimension] = rho_E_bar_lim[j](xl[face]);
+ }
+ });
+
+ // 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);
+ _applyNeumannflatBoundaryCondition(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);
+ NodeValue<double> MaxErrNode{p_mesh->connectivity()};
+ MaxErrNode.fill(0);
+ FaceValue<double> MaxErrFace{p_mesh->connectivity()};
+ MaxErrFace.fill(0);
+ EdgeValue<double> MaxErrEdge{p_mesh->connectivity()};
+ MaxErrEdge.fill(0);
+
+ NodeValue<int> nbChangingSignNode{p_mesh->connectivity()};
+ nbChangingSignNode.fill(0);
+ 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);
+
+ /*
+ // Moyenne de 2 etats
+ 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;
+ */
+
+ // La moyenne de Roe entre les etats jk
+ RoeAverageStateStructData RoeS(
+ RoeAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K], gamma, p_n[K]));
+ // Viscosity j k
+ const double nrmCkr = l2Norm(Ckr_loc);
+ const Rd& normal = 1. / nrmCkr * Ckr_loc;
+ const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, normal);
+
+ const double nrmCjr = l2Norm(Cjr_loc);
+ const Rd& normalJ = 1. / nrmCjr * Cjr_loc;
+ const JacobianInformations& JacInfoJ = JacobianFluxAlongUnitNormal(RoeS, normalJ);
+
+ // Matrice de Viscosite : Valeur Abs de la jacobienne en l'etat de Roe
+ Rpxp ViscosityMatrixJK = .5 * (nrmCjr * JacInfoJ.AbsJacobian + nrmCkr * JacInfoK.AbsJacobian);
+
+ // Test Rajout Viscosité type Rusanov v2 si chgt signe vp ..
+ bool anySignChgJ = EvaluateChangingSignVpAlongNormal( // rho_n[j],
+ u_n[j], // E_n[j], gamma,
+ c_n[j], // p_n[j], rho_n[K],
+ u_n[K],
+ // E_n[K], gamma,
+ c_n[K], // p_n[K],
+ normalJ);
+ bool anySignChgK = EvaluateChangingSignVpAlongNormal( // rho_n[j],
+ u_n[j], // E_n[j], gamma,
+ c_n[j], // p_n[j], rho_n[K],
+ u_n[K],
+ // E_n[K], gamma,
+ c_n[K], // p_n[K],
+ normal);
+
+ bool anySignDif = false;
+ if (anySignChgJ or anySignChgK)
+ anySignDif = true;
+
+ if (anySignDif) {
+ Rpxp AddedViscosity(identity);
+ AddedViscosity *= std::max(JacInfoJ.maxabsVpNormal * nrmCjr, JacInfoK.maxabsVpNormal * nrmCkr);
+ ViscosityMatrixJK += AddedViscosity;
+ nbChangingSignNode[node] = 1;
+ }
+
+ MaxErrNode[node] = std::max(MaxErrNode[node], JacInfoK.MaxErrorProd);
+
+ 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) << " Max Error RoeMatrice " << max(MaxErrNode)
+ << " nb Chg Sign " << sum(nbChangingSignNode) << "\n";
+ }
+ //
+ FaceValue<int> nbChangingSignFace{p_mesh->connectivity()};
+ nbChangingSignFace.fill(0);
+ 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);
+ /*
+ // Moyenne de 2 etats
+ 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;
+ */
+ // La moyenne de Roe entre les etats jk
+ RoeAverageStateStructData RoeS(
+ RoeAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K], gamma, p_n[K]));
+ // Viscosity j k
+ const double nrmCkf = l2Norm(Ckf_loc);
+ const Rd& normal = 1. / nrmCkf * Ckf_loc;
+ const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, normal);
+ const double nrmCjf = l2Norm(Cjf_loc);
+ const Rd& normalJ = 1. / nrmCjf * Cjf_loc;
+ const JacobianInformations& JacInfoJ = JacobianFluxAlongUnitNormal(RoeS, normalJ);
+
+ // Matrice de Viscosite : Valeur Abs de la jacobienne en l'etat de Roe
+ Rpxp ViscosityMatrixJK = .5 * (nrmCjf * JacInfoJ.AbsJacobian + nrmCkf * JacInfoK.AbsJacobian);
+
+ // Test Rajout Viscosité type Rusanov v2 si chgt signe vp ..
+
+ bool anySignChgJ = EvaluateChangingSignVpAlongNormal( // rho_n[j],
+ u_n[j], // E_n[j], gamma,
+ c_n[j], // p_n[j], rho_n[K],
+ u_n[K],
+ // E_n[K], gamma,
+ c_n[K], // p_n[K],
+ normalJ);
+ bool anySignChgK = EvaluateChangingSignVpAlongNormal( // rho_n[j],
+ u_n[j], // E_n[j], gamma,
+ c_n[j], // p_n[j], rho_n[K],
+ u_n[K],
+ // E_n[K], gamma,
+ c_n[K], // p_n[K],
+ normal);
+
+ bool anySignDif = false;
+ if (anySignChgJ or anySignChgK)
+ anySignDif = true;
+
+ if (anySignDif) {
+ Rpxp AddedViscosity(identity);
+ AddedViscosity *= std::max(JacInfoJ.maxabsVpNormal * nrmCjf, JacInfoK.maxabsVpNormal * nrmCkf);
+ ViscosityMatrixJK += AddedViscosity;
+ nbChangingSignFace[face] = 1;
+ }
+
+ MaxErrFace[face] = std::max(MaxErrFace[face], JacInfoK.MaxErrorProd);
+
+ 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) << " Max Error RoeMatrice " << max(MaxErrFace)
+ << " nb Chg Sign " << sum(nbChangingSignFace) << "\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();
+
+ EdgeValue<int> nbChangingSignEdge(p_mesh->connectivity());
+ nbChangingSignEdge.fill(0);
+
+ 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);
+ /*
+ // Moyenne de 2 etats
+ 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;
+ */
+ // La moyenne de Roe entre les etats jk
+ RoeAverageStateStructData RoeS(
+ RoeAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K], gamma, p_n[K]));
+ // Viscosity j k
+ const double nrmCke = l2Norm(Cke_loc);
+ const Rd& normal = 1. / nrmCke * Cke_loc;
+ const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, normal);
+
+ const double nrmCje = l2Norm(Cje_loc);
+ const Rd& normalJ = 1. / nrmCje * Cje_loc;
+ const JacobianInformations& JacInfoJ = JacobianFluxAlongUnitNormal(RoeS, normalJ);
+
+ // Matrice de Viscosite : Valeur Abs de la jacobienne en l'etat de Roe
+ Rpxp ViscosityMatrixJK = .5 * (nrmCje * JacInfoJ.AbsJacobian + nrmCke * JacInfoK.AbsJacobian);
+
+ // Test Rajout Viscosité type Rusanov v2 si chgt signe vp ..
+
+ bool anySignChgJ = EvaluateChangingSignVpAlongNormal( // rho_n[j],
+ u_n[j], // E_n[j], gamma,
+ c_n[j], // p_n[j], rho_n[K],
+ u_n[K],
+ // E_n[K], gamma,
+ c_n[K], // p_n[K],
+ normalJ);
+ bool anySignChgK = EvaluateChangingSignVpAlongNormal( // rho_n[j],
+ u_n[j], // E_n[j], gamma,
+ c_n[j], // p_n[j], rho_n[K],
+ u_n[K],
+ // E_n[K], gamma,
+ c_n[K], // p_n[K],
+ normal);
+
+ bool anySignDif = false;
+ if (anySignChgJ or anySignChgK)
+ anySignDif = true;
+
+ if (anySignDif) {
+ Rpxp AddedViscosity(identity);
+ AddedViscosity *= std::max(JacInfoJ.maxabsVpNormal * nrmCje, JacInfoK.maxabsVpNormal * nrmCke);
+ ViscosityMatrixJK += AddedViscosity;
+ nbChangingSignEdge[edge] = 1;
+ }
+
+ MaxErrEdge[edge] = std::max(MaxErrEdge[edge], JacInfoK.MaxErrorProd);
+
+ 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) << " Max Error RoeMatrice " << max(MaxErrEdge)
+ << " nb Chg Sign " << sum(nbChangingSignEdge) << "\n";
+ }
+ } // dim 3
+
+ // Pour les assemblages
+ double theta = 2. / 3.; //.5;
+ double eta = 1. / 6.; //.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));
+ }
+
+ RoeViscousFormEulerianCompositeSolver_v2_o2() = default;
+ ~RoeViscousFormEulerianCompositeSolver_v2_o2() = default;
+};
+
+template <MeshConcept MeshType>
+class RoeViscousFormEulerianCompositeSolver_v2_o2<MeshType>::NeumannBoundaryCondition
+{
+};
+
+template <>
+class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<2>>::NeumannBoundaryCondition
+{
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
+ // const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
+ // const MeshFlatFaceBoundary<MeshType> m_mesh_flat_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 RoeViscousFormEulerianCompositeSolver_v2_o2<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 RoeViscousFormEulerianCompositeSolver_v2_o2<MeshType>::NeumannflatBoundaryCondition
+{
+};
+template <>
+class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<2>>::NeumannflatBoundaryCondition
+{
+ 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();
+ }
+
+ NeumannflatBoundaryCondition(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)
+ {
+ ;
+ }
+
+ ~NeumannflatBoundaryCondition() = default;
+};
+
+template <>
+class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<3>>::NeumannflatBoundaryCondition
+{
+ 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();
+ }
+
+ NeumannflatBoundaryCondition(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)
+ {
+ ;
+ }
+
+ ~NeumannflatBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class RoeViscousFormEulerianCompositeSolver_v2_o2<MeshType>::SymmetryBoundaryCondition
+{
+};
+
+template <>
+class RoeViscousFormEulerianCompositeSolver_v2_o2<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 RoeViscousFormEulerianCompositeSolver_v2_o2<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 RoeViscousFormEulerianCompositeSolver_v2_o2<MeshType>::InflowListBoundaryCondition
+{
+};
+
+template <>
+class RoeViscousFormEulerianCompositeSolver_v2_o2<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 RoeViscousFormEulerianCompositeSolver_v2_o2<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 RoeViscousFormEulerianCompositeSolver_v2_o2<MeshType>::OutflowBoundaryCondition
+{
+};
+
+template <>
+class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<2>>::OutflowBoundaryCondition
+{
+ 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();
+ }
+
+ OutflowBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary)
+ : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary)
+ {
+ ;
+ }
+};
+
+template <>
+class RoeViscousFormEulerianCompositeSolver_v2_o2<Mesh<3>>::OutflowBoundaryCondition
+{
+ 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();
+ }
+
+ OutflowBoundaryCondition(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)
+ {
+ ;
+ }
+};
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+roeViscousFormEulerianCompositeSolver_v2_o2(
+ 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 size_t& degree,
+ 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("RoeViscousFormEulerianCompositeSolver v2 is not available in 1D");
+ } else {
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ return RoeViscousFormEulerianCompositeSolver_v2_o2<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>>(), degree, bc_descriptor_list, dt, check);
+ } else {
+ throw NormalError("RoeViscousFormEulerianCompositeSolver v2 is only defined on polygonal meshes");
+ }
+ }
+ },
+ mesh_v->variant());
+}
diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.hpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.hpp
new file mode 100644
index 000000000..c7351a1ef
--- /dev/null
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.hpp
@@ -0,0 +1,31 @@
+#ifndef ROE_VISCOUS_FORM_EULERIAN_COMPOSITE_SOLVER_V2_O2_HPP
+#define ROE_VISCOUS_FORM_EULERIAN_COMPOSITE_SOLVER_V2_O2_HPP
+
+#include <mesh/MeshVariant.hpp>
+#include <scheme/CellbyCellLimitation.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>>
+roeViscousFormEulerianCompositeSolver_v2_o2(
+ 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 size_t& degree,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const double& dt,
+ const bool check = false);
+
+#endif // ROE_VISCOUS_FORM_EULERIAN_COMPOSITE_SOLVER_V2_O2_HPP
diff --git a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp
new file mode 100644
index 000000000..954dd8784
--- /dev/null
+++ b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp
@@ -0,0 +1,2085 @@
+#include <scheme/RusanovEulerianCompositeSolver_o2.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/DiscreteFunctionDPk.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionDPkVector.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/PugsTraits.hpp>
+#include <variant>
+
+template <MeshConcept MeshTypeT>
+class RusanovEulerianCompositeSolver_o2
+{
+ 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;
+ class NeumannflatBoundaryCondition;
+
+ using BoundaryCondition = std::variant<SymmetryBoundaryCondition,
+ InflowListBoundaryCondition,
+ OutflowBoundaryCondition,
+ NeumannflatBoundaryCondition,
+ 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::neumann: {
+ if constexpr (Dimension == 2) {
+ bc_list.emplace_back(
+ NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ } else {
+ static_assert(Dimension == 3);
+ bc_list.emplace_back(
+ NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ }
+ break;
+ }
+ 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: {
+ if constexpr (Dimension == 2) {
+ bc_list.emplace_back(
+ OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ } else {
+ static_assert(Dimension == 3);
+ bc_list.emplace_back(
+ OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ }
+ break;
+ // 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 Eulerian Composite solver";
+ throw NormalError(error_msg.str());
+ }
+ }
+
+ return bc_list;
+ }
+
+ public:
+ CellByCellLimitation<MeshType> Limitor;
+
+ public:
+ void
+ _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>) {
+ std::cout << " Traitement Outflow \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();
+
+ const auto xj = mesh.xj();
+ const auto xr = mesh.xr();
+ const auto xf = mesh.xl();
+ const auto xe = mesh.xe();
+
+ 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] += vectorSym[dim];
+
+ 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];
+ }
+ }
+
+ // throw NormalError("Not implemented");
+ }
+ */// throw NormalError("Not implemented");
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _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);
+ }
+ }
+
+ void
+ _applyNeumannflatBoundaryCondition(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<NeumannflatBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement WALL \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];
+
+ vectorSym -= dot(vectorSym, normal) * normal;
+
+ 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];
+
+ vectorSym -= dot(vectorSym, normal) * normal;
+
+ 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];
+
+ vectorSym -= dot(vectorSym, normal) * normal;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _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);
+ }
+ }
+
+ public:
+ 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 size_t& degree,
+ 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);
+ std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
+
+ for (auto&& bc_descriptor : bc_descriptor_list) {
+ if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
+ symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+ }
+ }
+ static const int degre = 1;
+ PolynomialReconstructionDescriptor
+ reconstruction_descriptor(IntegrationMethodType::cell_center, // IntegrationMethodType::boundary,
+ degre, symmetry_boundary_descriptor_list);
+ //
+ //
+ //
+ CellValue<double> valrho({p_mesh->connectivity()});
+ CellValue<Rd> valu({p_mesh->connectivity()});
+ CellValue<double> valeps({p_mesh->connectivity()});
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ valrho[j] = rho[j];
+ valu[j] = u[j];
+ valeps[j] = E[j] - .5 * dot(u[j], u[j]);
+ });
+
+ DiscreteFunctionP0<double> eps(p_mesh, valeps);
+
+ // assemblage direct
+ // auto reconstructions = PolynomialReconstruction{reconstruction_descriptor}.build(rho, rho * u, rho * E);
+
+ // Pour assemblage Leibniz
+ // ca calcul au moins les gradients
+ auto reconstructions = PolynomialReconstruction{reconstruction_descriptor}.build(rho, u, eps);
+
+ DiscreteFunctionDPk rho_bar = reconstructions[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ DiscreteFunctionDPk u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ DiscreteFunctionDPk eps_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ // DiscreteFunctionDPk rho_u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ // DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ CellValue<double> Limitor_rho(p_mesh->connectivity());
+ CellValue<double> Limitor_u(p_mesh->connectivity());
+ CellValue<double> Limitor_eps(p_mesh->connectivity());
+
+ CellValue<Rd> grad_rho_cell(p_mesh->connectivity());
+ CellValue<Rdxd> grad_u_cell(p_mesh->connectivity());
+ CellValue<Rd> grad_eps_cell(p_mesh->connectivity());
+
+ // DiscreteFunctionDPk<Dimension, const double> rho_bar_lim;
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j](x) - rho[j]);
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ grad_rho_cell[j][dim] = rho_bar.coefficients(j)[1 + dim]; // si base : x, puis y , puis z
+ grad_eps_cell[j][dim] = eps_bar.coefficients(j)[1 + dim]; // idem ci dessus
+ const Rd gradu = u_bar.coefficients(j)[1 + dim];
+ for (size_t dim2 = 0; dim2 < Dimension; ++dim2) {
+ grad_u_cell[j](dim, dim2) = gradu[dim2];
+ }
+ }
+ grad_u_cell[j] = transpose(grad_u_cell[j]);
+ });
+
+ // Appels des limiteurs cell by cell
+ // Pour rho (plus ou moins classique)
+ Limitor.computeLimitorVolumicScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, Limitor_rho);
+
+ // Pour les variables d'interet (physique)
+ // Pour eps ici G.P. (avec Leibniz et toute l artillerie)
+ Limitor.computeLimitorInternalNrjScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, valeps,
+ grad_eps_cell,
+ // valu,
+ grad_u_cell, Limitor_eps);
+ // Pour u : relation deduite du principe d'induction
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Limitor_u[j] = sqrt(Limitor_eps[j]); });
+
+ //
+ // Creation des variables conservatives limitées
+ //
+ auto rho_bar_lim = copy(rho_bar);
+ auto rho_u_bar_lim = copy(u_bar); // a modifier .. du a la densite
+ auto rho_E_bar_lim = copy(eps_bar); // a modifier et a completer : densité et energie cinetique
+
+ // Assemblage des ctes et ordres élevés
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j](x) - rho[j]);
+ // rho_u_bar_lim.coefficients(j)[0] = valrho[j] * valu[j]; // * valu[j];
+ // rho_E_bar_lim.coefficients(j)[0] = valrho[j] * (valeps[j] + .5 * dot(valu[j], valu[j]));
+
+ // rho_bar_lim.coefficients(j)[0] // cte inchangée
+ rho_u_bar_lim.coefficients(j)[0] *= valrho[j]; // du coup (rho u)
+ rho_E_bar_lim.coefficients(j)[0] *= valrho[j]; // du coup (rho*epsilon)
+ rho_E_bar_lim.coefficients(j)[0] += valrho[j] * .5 * dot(valu[j], valu[j]);
+ // Rdxd grad_rhou_lim=
+ // valrho[j]*Limitor_u[j]*grad_u_cell[j]+Limitor_rho[j]*tensorProduct(valu[j],grad_rho_cell[j]);
+ //
+ Rdxd gradU_lim_transpose = Limitor_u[j] * transpose(grad_u_cell[j]);
+ const Rd gradUtu_lim = gradU_lim_transpose * valu[j];
+
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ rho_bar_lim.coefficients(j)[1 + dim] *= Limitor_rho[j];
+ rho_u_bar_lim.coefficients(j)[1 + dim] *= Limitor_u[j] * valrho[j]; // c'est la partie gradu limite
+ rho_u_bar_lim.coefficients(j)[1 + dim] += Limitor_rho[j] * grad_rho_cell[j][dim] * valu[j];
+
+ rho_E_bar_lim.coefficients(j)[1 + dim] *= Limitor_eps[j] * valrho[j];
+ rho_E_bar_lim.coefficients(j)[1 + dim] += gradUtu_lim[dim] * valrho[j];
+
+ rho_E_bar_lim.coefficients(j)[1 + dim] +=
+ (valeps[j] + .5 * dot(valu[j], valu[j])) * Limitor_rho[j] * grad_rho_cell[j][dim];
+ }
+ });
+
+ 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];
+ });
+
+ //
+ // Remplir ici les reconstructions d'ordre élevé
+
+ //
+ 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();
+
+ const auto xr = p_mesh->xr();
+ auto xl = MeshDataManager::instance().getMeshData(*p_mesh).xl();
+ auto xe = MeshDataManager::instance().getMeshData(*p_mesh).xe();
+
+ NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()};
+ StateAtNode.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // StateAtNode[j].fill(State[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];
+
+ StateAtNode[j][l][0] = rho_bar_lim[j](xr[node]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtNode[j][l][1 + dim] = rho_u_bar_lim[j](xr[node])[dim];
+ StateAtNode[j][l][1 + Dimension] = rho_E_bar_lim[j](xr[node]);
+ }
+ });
+
+ EdgeValuePerCell<Rp> StateAtEdge{p_mesh->connectivity()};
+ StateAtEdge.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // eStateAtEdge[j].fill(State[j]);
+ 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];
+
+ StateAtEdge[j][l][0] = rho_bar_lim[j](xe[edge]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtEdge[j][l][1 + dim] = rho_u_bar_lim[j](xe[edge])[dim];
+ StateAtEdge[j][l][1 + Dimension] = rho_E_bar_lim[j](xe[edge]);
+ }
+ });
+ FaceValuePerCell<Rp> StateAtFace{p_mesh->connectivity()};
+ StateAtFace.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // StateAtFace[j].fill(State[j]);
+ 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];
+
+ StateAtFace[j][l][0] = rho_bar_lim[j](xl[face]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtFace[j][l][1 + dim] = rho_u_bar_lim[j](xl[face])[dim];
+ StateAtFace[j][l][1 + Dimension] = rho_E_bar_lim[j](xl[face]);
+ }
+ });
+
+ // 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();
+
+ NodeValue<Rpxp> ViscosityNodeMatrix{p_mesh->connectivity()};
+ ViscosityNodeMatrix.fill(identity);
+
+ parallel_for(
+ p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId l) {
+ double maxabsVpNormCjr = 0;
+ 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);
+ // double rhoNode = 0;
+ Rd uNode = zero;
+ // double ENode = 0;
+ double cNode = 0;
+
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ // rhoNode += rho_n[J];
+ uNode += u_n[J];
+ // ENode += E_n[J];
+ cNode += c_n[J];
+ }
+ // rhoNode /= node_to_cell.size();
+ uNode *= 1. / node_to_cell.size();
+ // ENode /= node_to_cell.size();
+ cNode /= node_to_cell.size();
+
+ // Boundary conditions ..
+
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cells[j];
+
+ const Rd& Cjr_loc = Cjr(J, R);
+ maxabsVpNormCjr =
+ std::max(maxabsVpNormCjr,
+ toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // rhoNode,
+ uNode, // ENode,
+ cNode, Cjr_loc));
+ }
+ ViscosityNodeMatrix[l] *= maxabsVpNormCjr;
+ });
+
+ FaceValue<Rpxp> ViscosityFaceMatrix{p_mesh->connectivity()};
+ ViscosityFaceMatrix.fill(identity);
+
+ parallel_for(
+ p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) {
+ double maxabsVpNormCjf = 0;
+ 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);
+ // double rhoFace = 0;
+ Rd uFace = zero;
+ // double EFace = 0;
+ double cFace = 0;
+
+ for (size_t j = 0; j < face_to_cell.size(); ++j) {
+ const CellId J = face_to_cell[j];
+ // rhoFace += rho_n[J];
+ uFace += u_n[J];
+ // EFace += E_n[J];
+ cFace += c_n[J];
+ }
+ // rhoFace /= face_to_cell.size();
+ uFace *= 1. / face_to_cell.size();
+ // EFace /= face_to_cell.size();
+ cFace /= face_to_cell.size();
+
+ // Boundary conditions ..
+
+ for (size_t j = 0; j < face_to_cell.size(); ++j) {
+ const CellId J = face_to_cell[j];
+ const unsigned int R = face_local_number_in_its_cells[j];
+
+ const Rd& Cjf_loc = Cjf(J, R);
+
+ 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, 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);
+ Rp statediff(zero);
+
+ 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_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 = ViscosityNodeMatrix[node] * 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] += ViscosityNodeMatrix[node] * statediff;
+
+ 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);
+ Rp statediff(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];
+
+ 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 = ViscosityFaceMatrix[face] * 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] += ViscosityFaceMatrix[face] * statediff;
+
+ 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.);
+ // 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) {
+ 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];
+ }
+ /*
+ 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) {
+ 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();
+
+ EdgeValue<Rpxp> ViscosityEdgeMatrix{p_mesh->connectivity()};
+ ViscosityEdgeMatrix.fill(identity);
+
+ parallel_for(
+ p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId l) {
+ double maxabsVpNormCje = 0;
+ 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);
+ // double rhoEdge = 0;
+ Rd uEdge = zero;
+ // double EEdge = 0;
+ double cEdge = 0;
+
+ for (size_t j = 0; j < edge_to_cell.size(); ++j) {
+ const CellId J = edge_to_cell[j];
+ // rhoEdge += rho_n[J];
+ uEdge += u_n[J];
+ // EEdge += E_n[J];
+ cEdge += c_n[J];
+ }
+ // rhoEdge /= edge_to_cell.size();
+ uEdge *= 1. / edge_to_cell.size();
+ // EEdge /= edge_to_cell.size();
+ cEdge /= edge_to_cell.size();
+
+ // Boundary conditions ..
+
+ for (size_t j = 0; j < edge_to_cell.size(); ++j) {
+ const CellId J = edge_to_cell[j];
+ const unsigned int R = edge_local_number_in_its_cells[j];
+
+ const Rd& Cje_loc = Cje(J, R);
+
+ maxabsVpNormCje =
+ std::max(maxabsVpNormCje,
+ toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // rhoEdge,
+ uEdge, // EEdge,
+ cEdge, Cje_loc));
+ }
+ ViscosityEdgeMatrix[l] *= maxabsVpNormCje;
+ });
+
+ 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);
+ Rp statediff(zero);
+
+ 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_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 = ViscosityEdgeMatrix[edge] * 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] += ViscosityEdgeMatrix[edge] * statediff;
+
+ 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 = 2. / 3.; //.5;
+ double eta = 1. / 6.; //.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_o2() = default;
+ ~RusanovEulerianCompositeSolver_o2() = default;
+};
+
+template <MeshConcept MeshType>
+class RusanovEulerianCompositeSolver_o2<MeshType>::NeumannBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver_o2<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_o2<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_o2<MeshType>::NeumannflatBoundaryCondition
+{
+};
+template <>
+class RusanovEulerianCompositeSolver_o2<Mesh<2>>::NeumannflatBoundaryCondition
+{
+ 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();
+ }
+
+ NeumannflatBoundaryCondition(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)
+ {
+ ;
+ }
+
+ ~NeumannflatBoundaryCondition() = default;
+};
+
+template <>
+class RusanovEulerianCompositeSolver_o2<Mesh<3>>::NeumannflatBoundaryCondition
+{
+ 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();
+ }
+
+ NeumannflatBoundaryCondition(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)
+ {
+ ;
+ }
+
+ ~NeumannflatBoundaryCondition() = default;
+};
+
+template <>
+class RusanovEulerianCompositeSolver_o2<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_o2<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_o2<MeshType>::InflowListBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver_o2<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_o2<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_o2<MeshType>::OutflowBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver_o2<Mesh<2>>::OutflowBoundaryCondition
+{
+ 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();
+ }
+
+ OutflowBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary)
+ : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary)
+ {
+ ;
+ }
+};
+
+template <>
+class RusanovEulerianCompositeSolver_o2<Mesh<3>>::OutflowBoundaryCondition
+{
+ 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();
+ }
+
+ OutflowBoundaryCondition(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)
+ {
+ ;
+ }
+};
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+rusanovEulerianCompositeSolver_o2(
+ 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 size_t& degree,
+ 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_o2 is not available in 1D");
+ } else {
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ return RusanovEulerianCompositeSolver_o2<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>>(), degree, bc_descriptor_list, dt, check);
+ } else {
+ throw NormalError("RusanovEulerianCompositeSolver_o2 is only defined on polygonal meshes");
+ }
+ }
+ },
+ mesh_v->variant());
+}
diff --git a/src/scheme/RusanovEulerianCompositeSolver_o2.hpp b/src/scheme/RusanovEulerianCompositeSolver_o2.hpp
new file mode 100644
index 000000000..228d7fd45
--- /dev/null
+++ b/src/scheme/RusanovEulerianCompositeSolver_o2.hpp
@@ -0,0 +1,30 @@
+#ifndef RUSANOV_EULERIAN_COMPOSITE_SOLVER_O2_HPP
+#define RUSANOV_EULERIAN_COMPOSITE_SOLVER_O2_HPP
+
+#include <memory>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/CellbyCellLimitation.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+#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_o2(
+ 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 size_t& degree,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const double& dt,
+ const bool check = false);
+
+#endif // RUSANOV_EULERIAN_COMPOSITE_SOLVER_HPP
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
new file mode 100644
index 000000000..2826f8de7
--- /dev/null
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp
@@ -0,0 +1,1995 @@
+#include <scheme/RusanovEulerianCompositeSolver_v2_o2.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/DiscreteFunctionDPk.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionDPkVector.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/PugsTraits.hpp>
+#include <variant>
+
+template <MeshConcept MeshTypeT>
+class RusanovEulerianCompositeSolver_v2_o2
+{
+ 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;
+ class NeumannflatBoundaryCondition;
+
+ using BoundaryCondition = std::variant<SymmetryBoundaryCondition,
+ InflowListBoundaryCondition,
+ OutflowBoundaryCondition,
+ NeumannflatBoundaryCondition,
+ 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::neumann: {
+ if constexpr (Dimension == 2) {
+ bc_list.emplace_back(
+ NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ } else {
+ static_assert(Dimension == 3);
+ bc_list.emplace_back(
+ NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ }
+ break;
+ }
+ 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: {
+ if constexpr (Dimension == 2) {
+ bc_list.emplace_back(
+ OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ } else {
+ static_assert(Dimension == 3);
+ bc_list.emplace_back(
+ OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+ getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ }
+ break;
+ // 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:
+ CellByCellLimitation<MeshType> Limitor;
+
+ public:
+ void
+ _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>) {
+ std::cout << " Traitement Outflow \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();
+
+ const auto xj = mesh.xj();
+ const auto xr = mesh.xr();
+ const auto xf = mesh.xl();
+ const auto xe = mesh.xe();
+
+ 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] += vectorSym[dim];
+
+ 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];
+ }
+ }
+
+ // throw NormalError("Not implemented");
+ }
+ */
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _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);
+ }
+ }
+
+ void
+ _applyNeumannflatBoundaryCondition(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<NeumannflatBoundaryCondition, T>) {
+ // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ std::cout << " Traitement WALL \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];
+
+ vectorSym -= dot(vectorSym, normal) * normal;
+
+ 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];
+
+ vectorSym -= dot(vectorSym, normal) * normal;
+
+ 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];
+
+ vectorSym -= dot(vectorSym, normal) * normal;
+
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
+ }
+ }
+ }
+ }
+ },
+ boundary_condition);
+ }
+ }
+
+ void
+ _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);
+ }
+ }
+
+ public:
+ inline double
+ pression(const double rho, const double epsilon, const double gam) const
+ {
+ return (gam - 1) * rho * epsilon;
+ }
+
+ // void
+ // computeLimitorVolumicScalarQuantityMinModDukowicz(const DiscreteFunctionDPk<Dimension, double>& q_bar,
+ // CellValue<double>& Limitor_q) const
+
+ 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 size_t& degree,
+ 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);
+
+ std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
+
+ for (auto&& bc_descriptor : bc_descriptor_list) {
+ if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
+ symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+ }
+ }
+ static const int degre = 1;
+ PolynomialReconstructionDescriptor
+ reconstruction_descriptor(IntegrationMethodType::cell_center, // IntegrationMethodType::boundary,
+ degre, symmetry_boundary_descriptor_list);
+ //
+ //
+ //
+ CellValue<double> valrho({p_mesh->connectivity()});
+ CellValue<Rd> valu({p_mesh->connectivity()});
+ CellValue<double> valeps({p_mesh->connectivity()});
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ valrho[j] = rho[j];
+ valu[j] = u[j];
+ valeps[j] = E[j] - .5 * dot(u[j], u[j]);
+ });
+
+ DiscreteFunctionP0<double> eps(p_mesh, valeps);
+
+ // assemblage direct
+ // auto reconstructions = PolynomialReconstruction{reconstruction_descriptor}.build(rho, rho * u, rho * E);
+
+ // Pour assemblage Leibniz
+ // ca calcul au moins les gradients
+ auto reconstructions = PolynomialReconstruction{reconstruction_descriptor}.build(rho, u, eps);
+
+ DiscreteFunctionDPk rho_bar = reconstructions[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ DiscreteFunctionDPk u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ DiscreteFunctionDPk eps_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ // DiscreteFunctionDPk rho_u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ // DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ CellValue<double> Limitor_rho(p_mesh->connectivity());
+ CellValue<double> Limitor_u(p_mesh->connectivity());
+ CellValue<double> Limitor_eps(p_mesh->connectivity());
+
+ CellValue<Rd> grad_rho_cell(p_mesh->connectivity());
+ CellValue<Rdxd> grad_u_cell(p_mesh->connectivity());
+ CellValue<Rd> grad_eps_cell(p_mesh->connectivity());
+
+ // DiscreteFunctionDPk<Dimension, const double> rho_bar_lim;
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j](x) - rho[j]);
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ grad_rho_cell[j][dim] = rho_bar.coefficients(j)[1 + dim]; // si base : x, puis y , puis z
+ grad_eps_cell[j][dim] = eps_bar.coefficients(j)[1 + dim]; // idem ci dessus
+ const Rd gradu = u_bar.coefficients(j)[1 + dim];
+ for (size_t dim2 = 0; dim2 < Dimension; ++dim2) {
+ grad_u_cell[j](dim, dim2) = gradu[dim2];
+ }
+ }
+ grad_u_cell[j] = transpose(grad_u_cell[j]);
+ });
+
+ // Appels des limiteurs cell by cell
+ // Pour rho (plus ou moins classique)
+ Limitor.computeLimitorVolumicScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, Limitor_rho);
+
+ // Pour les variables d'interet (physique)
+ // Pour eps ici G.P. (avec Leibniz et toute l artillerie)
+ Limitor.computeLimitorInternalNrjScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, valeps,
+ grad_eps_cell,
+ // valu,
+ grad_u_cell, Limitor_eps);
+ // Pour u : relation deduite du principe d'induction
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Limitor_u[j] = sqrt(Limitor_eps[j]); });
+
+ //
+ // Creation des variables conservatives limitées
+ //
+ auto rho_bar_lim = copy(rho_bar);
+ auto rho_u_bar_lim = copy(u_bar); // a modifier .. du a la densite
+ auto rho_E_bar_lim = copy(eps_bar); // a modifier et a completer : densité et energie cinetique
+
+ // Assemblage des ctes et ordres élevés
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j](x) - rho[j]);
+ // rho_u_bar_lim.coefficients(j)[0] = valrho[j] * valu[j]; // * valu[j];
+ // rho_E_bar_lim.coefficients(j)[0] = valrho[j] * (valeps[j] + .5 * dot(valu[j], valu[j]));
+
+ // rho_bar_lim.coefficients(j)[0] // cte inchangée
+ rho_u_bar_lim.coefficients(j)[0] *= valrho[j]; // du coup (rho u)
+ rho_E_bar_lim.coefficients(j)[0] *= valrho[j]; // du coup (rho*epsilon)
+ rho_E_bar_lim.coefficients(j)[0] += valrho[j] * .5 * dot(valu[j], valu[j]);
+ // Rdxd grad_rhou_lim=
+ // valrho[j]*Limitor_u[j]*grad_u_cell[j]+Limitor_rho[j]*tensorProduct(valu[j],grad_rho_cell[j]);
+ //
+ Rdxd gradU_lim_transpose = Limitor_u[j] * transpose(grad_u_cell[j]);
+ const Rd gradUtu_lim = gradU_lim_transpose * valu[j];
+
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ rho_bar_lim.coefficients(j)[1 + dim] *= Limitor_rho[j];
+ rho_u_bar_lim.coefficients(j)[1 + dim] *= Limitor_u[j] * valrho[j]; // c'est la partie gradu limite
+ rho_u_bar_lim.coefficients(j)[1 + dim] += Limitor_rho[j] * grad_rho_cell[j][dim] * valu[j];
+
+ rho_E_bar_lim.coefficients(j)[1 + dim] *= Limitor_eps[j] * valrho[j];
+ rho_E_bar_lim.coefficients(j)[1 + dim] += gradUtu_lim[dim] * valrho[j];
+
+ rho_E_bar_lim.coefficients(j)[1 + dim] +=
+ (valeps[j] + .5 * dot(valu[j], valu[j])) * Limitor_rho[j] * grad_rho_cell[j][dim];
+ }
+ });
+
+ // 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()};
+ //
+ // Remplir ici les reconstructions d'ordre élevé
+
+ //
+ 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();
+
+ const auto xr = p_mesh->xr();
+ auto xl = MeshDataManager::instance().getMeshData(*p_mesh).xl();
+ auto xe = MeshDataManager::instance().getMeshData(*p_mesh).xe();
+
+ NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()};
+ StateAtNode.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // StateAtNode[j].fill(State[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];
+
+ StateAtNode[j][l][0] = rho_bar_lim[j](xr[node]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtNode[j][l][1 + dim] = rho_u_bar_lim[j](xr[node])[dim];
+ StateAtNode[j][l][1 + Dimension] = rho_E_bar_lim[j](xr[node]);
+ }
+ });
+
+ EdgeValuePerCell<Rp> StateAtEdge{p_mesh->connectivity()};
+ StateAtEdge.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // eStateAtEdge[j].fill(State[j]);
+ 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];
+
+ StateAtEdge[j][l][0] = rho_bar_lim[j](xe[edge]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtEdge[j][l][1 + dim] = rho_u_bar_lim[j](xe[edge])[dim];
+ StateAtEdge[j][l][1 + Dimension] = rho_E_bar_lim[j](xe[edge]);
+ }
+ });
+ FaceValuePerCell<Rp> StateAtFace{p_mesh->connectivity()};
+ StateAtFace.fill(zero);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // StateAtFace[j].fill(State[j]);
+ 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];
+
+ StateAtFace[j][l][0] = rho_bar_lim[j](xl[face]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtFace[j][l][1 + dim] = rho_u_bar_lim[j](xl[face])[dim];
+ StateAtFace[j][l][1 + Dimension] = rho_E_bar_lim[j](xl[face]);
+ }
+ });
+
+ // 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);
+ _applyNeumannflatBoundaryCondition(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()};
+
+ 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 = 2. / 3.; //.5;
+ double eta = 1. / 6.; //.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_o2() = default;
+ ~RusanovEulerianCompositeSolver_v2_o2() = default;
+};
+
+template <MeshConcept MeshType>
+class RusanovEulerianCompositeSolver_v2_o2<MeshType>::NeumannBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver_v2_o2<Mesh<2>>::NeumannBoundaryCondition
+{
+ using Rd = TinyVector<Dimension, double>;
+
+ private:
+ const MeshNodeBoundary m_mesh_node_boundary;
+ const MeshFaceBoundary m_mesh_face_boundary;
+ // const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
+ // const MeshFlatFaceBoundary<MeshType> m_mesh_flat_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_o2<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_o2<MeshType>::NeumannflatBoundaryCondition
+{
+};
+template <>
+class RusanovEulerianCompositeSolver_v2_o2<Mesh<2>>::NeumannflatBoundaryCondition
+{
+ 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();
+ }
+
+ NeumannflatBoundaryCondition(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)
+ {
+ ;
+ }
+
+ ~NeumannflatBoundaryCondition() = default;
+};
+
+template <>
+class RusanovEulerianCompositeSolver_v2_o2<Mesh<3>>::NeumannflatBoundaryCondition
+{
+ 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();
+ }
+
+ NeumannflatBoundaryCondition(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)
+ {
+ ;
+ }
+
+ ~NeumannflatBoundaryCondition() = default;
+};
+
+template <MeshConcept MeshType>
+class RusanovEulerianCompositeSolver_v2_o2<MeshType>::SymmetryBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver_v2_o2<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_o2<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_o2<MeshType>::InflowListBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver_v2_o2<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_o2<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_o2<MeshType>::OutflowBoundaryCondition
+{
+};
+
+template <>
+class RusanovEulerianCompositeSolver_v2_o2<Mesh<2>>::OutflowBoundaryCondition
+{
+ 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();
+ }
+
+ OutflowBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary)
+ : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary)
+ {
+ ;
+ }
+};
+
+template <>
+class RusanovEulerianCompositeSolver_v2_o2<Mesh<3>>::OutflowBoundaryCondition
+{
+ 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();
+ }
+
+ OutflowBoundaryCondition(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)
+ {
+ ;
+ }
+};
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>>
+rusanovEulerianCompositeSolver_v2_o2(
+ 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 size_t& degree,
+ 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_o2<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>>(), degree, bc_descriptor_list, dt, check);
+ } else {
+ throw NormalError("RusanovEulerianCompositeSolver v2 is only defined on polygonal meshes");
+ }
+ }
+ },
+ mesh_v->variant());
+}
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.hpp b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.hpp
new file mode 100644
index 000000000..d47675e74
--- /dev/null
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.hpp
@@ -0,0 +1,30 @@
+#ifndef RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_O2_HPP
+#define RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_O2_HPP
+
+#include <memory>
+#include <mesh/MeshVariant.hpp>
+#include <scheme/CellbyCellLimitation.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+#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_o2(
+ 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 size_t& degree,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+ const double& dt,
+ const bool check = false);
+
+#endif // RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_O2_xHPP
--
GitLab