diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt index 009bb3bf1699600c6e520ec20ebf8092fddab94e..6e79e4324e5dd0a905b935ad4a616a67c8e55447 100644 --- a/src/scheme/CMakeLists.txt +++ b/src/scheme/CMakeLists.txt @@ -28,6 +28,7 @@ add_library( RusanovEulerianCompositeSolver.cpp RusanovEulerianCompositeSolver_v2.cpp RoeViscousFormEulerianCompositeSolver_v2.cpp + RoeFluxFormEulerianCompositeSolver_v2.cpp RusanovEulerianCompositeSolver_o2.cpp RusanovEulerianCompositeSolver_v2_o2.cpp RusanovEulerianCompositeSolver_v2_order_n.cpp diff --git a/src/scheme/CellbyCellLimitation.hpp b/src/scheme/CellbyCellLimitation.hpp index a443e71ef6781058e48c093516c2f2b5aaf89a68..ecd607a0f9ab81ab6b970755fcdfa8abe936f551 100644 --- a/src/scheme/CellbyCellLimitation.hpp +++ b/src/scheme/CellbyCellLimitation.hpp @@ -1,9 +1,6 @@ #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/ItemValueUtils.hpp> #include <mesh/Mesh.hpp> @@ -24,6 +21,7 @@ #include <scheme/DiscreteFunctionP0.hpp> #include <scheme/DiscreteFunctionUtils.hpp> #include <scheme/InflowListBoundaryConditionDescriptor.hpp> +#include <scheme/LimitationTools.hpp> #include <scheme/PolynomialReconstruction.hpp> #include <scheme/PolynomialReconstructionDescriptor.hpp> #include <utils/PugsTraits.hpp> @@ -44,7 +42,8 @@ class CellByCellLimitation density_limiter(const MeshType& mesh, const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list, const DiscreteFunctionP0<const double>& rho, - DiscreteFunctionDPk<Dimension, double>& rho_L) const + DiscreteFunctionDPk<Dimension, double>& rho_L, + const bool enableWeakBoundPositivityOnly = false) const { const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix(); const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix(); @@ -95,9 +94,13 @@ class CellByCellLimitation for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { const FaceId face_id = face_list[i_face]; - const Rd& x0 = xr[face_to_node_matrix[face_id][0]]; - const Rd& x2 = xr[face_to_node_matrix[face_id][1]]; - const Rd& x1 = .5 * (x0 + x2); + const auto& face_to_node = face_to_node_matrix[face_id]; + const NodeId& nodei = face_to_node[0]; + const NodeId& nodef = face_to_node[face_to_node.size() - 1]; + + const Rd x0 = xr[nodei]; + const Rd x1 = xl[face_id]; + const Rd x2 = xr[nodef]; const double rho_x0 = rho_L[cell_id](x0); const double rho_x1 = rho_L[cell_id](x1); @@ -106,16 +109,120 @@ class CellByCellLimitation rho_bar_min = std::min(rho_bar_min, std::min(rho_x0, std::min(rho_x1, rho_x2))); rho_bar_max = std::max(rho_bar_max, std::max(rho_x0, std::max(rho_x1, rho_x2))); } - // const LineParabolicTransformation<Dimension> t(x0, x1, x2); - // for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) { - // const double rho_xk = rho_L[cell_id](t(qf.point(iq))); - - // rho_bar_min = std::min(rho_bar_min, rho_xk); - // rho_bar_max = std::max(rho_bar_max, rho_xk); - //} } else { + throw NotImplementedError("not implement in 3D"); + } + /* + const double eps = 1E-14; + double coef1 = 1; + if (std::abs(rho_bar_max - rhoj) > eps) { + coef1 = (rho_max - rhoj) / ((rho_bar_max - rhoj)); + } + + double coef2 = 1.; + if (std::abs(rho_bar_min - rhoj) > eps) { + coef2 = (rho_min - rhoj) / ((rho_bar_min - rhoj)); + } +*/ + const double lambda = // std::max(0., std::min(1., std::min(coef1, coef2))); + toolsLimitation::computeCoefLimitationBarthJespersen(rhoj, rho_bar_min, rho_bar_max, rho_min, rho_max, + enableWeakBoundPositivityOnly); + + auto coefficients = rho_L.coefficients(cell_id); + + coefficients[0] = (1 - lambda) * rho[cell_id] + lambda * coefficients[0]; + + for (size_t i = 1; i < coefficients.size(); ++i) { + coefficients[i] *= lambda; + } + }); + } + + void + density_limiter(const MeshType& mesh, + const FaceArray<TinyVector<2>>& QuadratureFace, + const EdgeArray<TinyVector<2>>& QuadratureEdge, + const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list, + const DiscreteFunctionP0<const double>& rho, + DiscreteFunctionDPk<Dimension, double>& rho_L, + const bool enableWeakBoundPositivityOnly = false) const + { + const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix(); + const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix(); + + const auto& xr = mesh.xr(); + // const auto& xl = mesh.xl(); + + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + + // auto stencil = StencilManager::instance()._getStencilArray(mesh.connectivity(), 1); + auto stencil = StencilManager::instance() + .getCellToCellStencilArray(mesh.connectivity(), + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, + symmetry_boundary_descriptor_list); + auto xj = mesh_data.xj(); + const auto& xl = mesh_data.xl(); + + // const QuadratureFormula<1> qf = + // QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree)); + + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + const double rhoj = rho[cell_id]; + + double rho_min = rhoj; + double rho_max = rhoj; + + const auto cell_stencil = stencil[cell_id]; + for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) { + rho_min = std::min(rho_min, rho[cell_stencil[i_cell]]); + rho_max = std::max(rho_max, rho[cell_stencil[i_cell]]); } + double rho_bar_min = rhoj; + double rho_bar_max = rhoj; + + for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) { + const CellId cell_k_id = cell_stencil[i_cell]; + const double rho_xk = rho_L[cell_id](xj[cell_k_id]); + + rho_bar_min = std::min(rho_bar_min, rho_xk); + rho_bar_max = std::max(rho_bar_max, rho_xk); + } + + if constexpr (Dimension == 2) { + auto face_list = cell_to_face_matrix[cell_id]; + + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; + + const auto& face_to_node = face_to_node_matrix[face_id]; + const NodeId& nodei = face_to_node[0]; + const NodeId& nodef = face_to_node[face_to_node.size() - 1]; + + const Rd& x0 = xr[nodei]; + const Rd& x1 = xr[nodef]; + + const double rho_x0 = rho_L[cell_id](x0); + const double rho_x1 = rho_L[cell_id](x1); + + rho_bar_min = std::min(rho_bar_min, std::min(rho_x0, rho_x1)); + rho_bar_max = std::min(rho_bar_max, std::max(rho_x0, rho_x1)); + + for (size_t i = 0; i < QuadratureFace[face_id].size(); ++i) { + // const Rd& x_quad = xd + QuadratureFace[face][i][1] * (xf - xd); + const Rd& x_quad = xr[nodei] + QuadratureFace[face_id][i][1] * (xr[nodef] - xr[nodei]); + + const double rho_xk = rho_L[cell_id](x_quad); + + rho_bar_min = std::min(rho_bar_min, rho_xk); + rho_bar_max = std::max(rho_bar_max, rho_xk); + } + } + } else { + throw NotImplementedError("not implement in 3D"); + } + /* const double eps = 1E-14; double coef1 = 1; if (std::abs(rho_bar_max - rhoj) > eps) { @@ -126,8 +233,11 @@ class CellByCellLimitation if (std::abs(rho_bar_min - rhoj) > eps) { coef2 = (rho_min - rhoj) / ((rho_bar_min - rhoj)); } +*/ + const double lambda = // std::max(0., std::min(1., std::min(coef1, coef2))); - const double lambda = std::max(0., std::min(1., std::min(coef1, coef2))); + toolsLimitation::computeCoefLimitationBarthJespersen(rhoj, rho_bar_min, rho_bar_max, rho_min, rho_max, + enableWeakBoundPositivityOnly); auto coefficients = rho_L.coefficients(cell_id); @@ -139,26 +249,27 @@ class CellByCellLimitation }); } + // + // Idem for epsilon + // + void specific_internal_nrj_limiter( const MeshType& mesh, const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list, - - // const DiscreteFunctionP0<const double>& rho, - // const DiscreteFunctionDPk<Dimension, double>& rho_L, const DiscreteFunctionP0<const double>& epsilon, - // const DiscreteFunctionDPk<Dimension, double>& - auto epsilon_R(const CellId cell_id, const Rd& x), - CellValue<double>& lambda_epsilon) // const - // DiscreteFunctionDPk<Dimension, double>& epsilon_R) const + auto epsilon_R, //(const CellId cell_id, const Rd& x), + CellValue<double>& lambda_epsilon, + const bool enableWeakBoundPositivityOnly = false) const { const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix(); const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix(); const auto& xr = mesh.xr(); - const auto& xl = mesh.xl(); + // const auto& xl = mesh.xl(); MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + const auto& xl = mesh_data.xl(); // auto stencil = StencilManager::instance()._getStencilArray(mesh.connectivity(), 1); auto stencil = StencilManager::instance() @@ -216,7 +327,7 @@ class CellByCellLimitation // epsilon_R_max = std::max(epsilon_R_max, epsilon_xk); //} } - + /* const double eps = 1E-14; double coef1 = 1; if (std::abs(epsilon_R_max - epsilonj) > eps) { @@ -227,11 +338,114 @@ class CellByCellLimitation if (std::abs(epsilon_R_min - epsilonj) > eps) { coef2 = (epsilon_min - epsilonj) / ((epsilon_R_min - epsilonj)); } +*/ + lambda_epsilon[cell_id] = // std::max(0., std::min(1., std::min(coef1, coef2))); + toolsLimitation::computeCoefLimitationBarthJespersen(epsilonj, epsilon_R_min, epsilon_R_min, epsilon_min, + epsilon_max, enableWeakBoundPositivityOnly); + }); + } - lambda_epsilon[cell_id] = std::max(0., std::min(1., std::min(coef1, coef2))); + // + // + // + void + specific_internal_nrj_limiter( + const MeshType& mesh, + const FaceArray<TinyVector<2>>& QuadratureFace, + const EdgeArray<TinyVector<2>>& QuadratureEdge, + const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list, + const DiscreteFunctionP0<const double>& epsilon, + auto epsilon_R, //(const CellId cell_id, const Rd& x), + CellValue<double>& lambda_epsilon, + const bool enableWeakBoundPositivityOnly = false) const + { + const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix(); + const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix(); + + const auto& xr = mesh.xr(); + // const auto& xl = mesh.xl(); + + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + + // auto stencil = StencilManager::instance()._getStencilArray(mesh.connectivity(), 1); + auto stencil = StencilManager::instance() + .getCellToCellStencilArray(mesh.connectivity(), + StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes}, + symmetry_boundary_descriptor_list); + auto xj = mesh_data.xj(); + + // const QuadratureFormula<1> qf = + // QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree)); + + parallel_for( + mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { + const double epsilonj = epsilon[cell_id]; + + double epsilon_min = epsilonj; + double epsilon_max = epsilonj; + + const auto cell_stencil = stencil[cell_id]; + for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) { + epsilon_min = std::min(epsilon_min, epsilon[cell_stencil[i_cell]]); + epsilon_max = std::max(epsilon_max, epsilon[cell_stencil[i_cell]]); + } + + double epsilon_R_min = epsilonj; + double epsilon_R_max = epsilonj; + + for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) { + const CellId cell_k_id = cell_stencil[i_cell]; + const double epsilon_xk = epsilon_R(cell_id, xj[cell_k_id]); + + epsilon_R_min = std::min(epsilon_R_min, epsilon_xk); + epsilon_R_max = std::max(epsilon_R_max, epsilon_xk); + } + + auto face_list = cell_to_face_matrix[cell_id]; + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; + + const Rd& x0 = xr[face_to_node_matrix[face_id][0]]; + // const Rd& x1 = xl[face_id][0]; + const Rd& x2 = xr[face_to_node_matrix[face_id][1]]; + const double epsilon_x0 = epsilon_R(cell_id, x0); + // const double epsilon_x1 = epsilon_R(cell_id, x1); + const double epsilon_x2 = epsilon_R(cell_id, x2); + epsilon_R_min = std::min(epsilon_R_min, std::min(epsilon_x0, epsilon_x2)); + epsilon_R_max = std::max(epsilon_R_max, std::max(epsilon_x0, epsilon_x2)); + + for (size_t i = 0; i < QuadratureFace[face_id].size(); ++i) { + // const Rd& x_quad = xd + QuadratureFace[face][i][1] * (xf - xd); + const Rd& x_quad = x0 + QuadratureFace[face_id][i][1] * (x2 - x0); + + const double epsilon_xk = epsilon_R(cell_id, x_quad); + + epsilon_R_min = std::min(epsilon_R_min, epsilon_xk); + epsilon_R_max = std::max(epsilon_R_max, epsilon_xk); + } + } + /* + const double eps = 1E-14; + double coef1 = 1; + if (std::abs(epsilon_R_max - epsilonj) > eps) { + coef1 = (epsilon_max - epsilonj) / ((epsilon_R_max - epsilonj)); + } + + double coef2 = 1.; + if (std::abs(epsilon_R_min - epsilonj) > eps) { + coef2 = (epsilon_min - epsilonj) / ((epsilon_R_min - epsilonj)); + } +*/ + lambda_epsilon[cell_id] = // std::max(0., std::min(1., std::min(coef1, coef2))); + toolsLimitation::computeCoefLimitationBarthJespersen(epsilonj, epsilon_R_min, epsilon_R_max, epsilon_min, + epsilon_max, enableWeakBoundPositivityOnly); }); } + // + // Other version (CellValue Coef for limitation ) + // + void computeLimitorVolumicScalarQuantityMinModDukowiczGradient(const MeshType& mesh, const CellValue<double>& q, diff --git a/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.cpp new file mode 100644 index 0000000000000000000000000000000000000000..999d343f2b8b1d18592dce15d39bdf5d66cd9c96 --- /dev/null +++ b/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.cpp @@ -0,0 +1,2846 @@ +#include <scheme/RoeFluxFormEulerianCompositeSolver_v2.hpp> + +#include <language/utils/InterpolateItemArray.hpp> +#include <mesh/Mesh.hpp> +#include <mesh/MeshData.hpp> +#include <mesh/MeshDataManager.hpp> +#include <mesh/MeshEdgeBoundary.hpp> +#include <mesh/MeshFaceBoundary.hpp> +#include <mesh/MeshFlatEdgeBoundary.hpp> +#include <mesh/MeshFlatFaceBoundary.hpp> +#include <mesh/MeshFlatNodeBoundary.hpp> +#include <mesh/MeshNodeBoundary.hpp> +#include <mesh/MeshTraits.hpp> +#include <mesh/MeshVariant.hpp> +#include <mesh/SubItemValuePerItemUtils.hpp> +#include <scheme/DiscreteFunctionUtils.hpp> +#include <scheme/InflowListBoundaryConditionDescriptor.hpp> + +#include <variant> + +template <MeshConcept MeshTypeT> +class RoeFluxFormEulerianCompositeSolver_v2 +{ + private: + using MeshType = MeshTypeT; + + static constexpr size_t Dimension = MeshType::Dimension; + + using Rdxd = TinyMatrix<Dimension>; + using Rd = TinyVector<Dimension>; + + using Rpxp = TinyMatrix<Dimension + 2>; + using Rp = TinyVector<Dimension + 2>; + + using Rpxd = TinyMatrix<Dimension + 2, Dimension>; + + class SymmetryBoundaryCondition; + class InflowListBoundaryCondition; + class OutflowBoundaryCondition; + class WallBoundaryCondition; + class NeumannflatBoundaryCondition; + + using BoundaryCondition = std::variant<SymmetryBoundaryCondition, + InflowListBoundaryCondition, + OutflowBoundaryCondition, + NeumannflatBoundaryCondition, + WallBoundaryCondition>; + + using BoundaryConditionList = std::vector<BoundaryCondition>; + + Rpxp + SignofASquareDiagonalisableMatrix(const Rpxp& Left, const Rp& Diag_vp, const Rpxp& Right) const + { + Rpxp MatriceDiag(identity); + for (size_t i = 0; i < Left.numberOfRows(); ++i) + MatriceDiag(i, i) *= signe(Diag_vp[i]); + + return Left * (MatriceDiag * Right); + + // return (Left * Diag_vp * Right); + Rpxp M(zero); + + for (size_t i = 0; i < Left.numberOfRows(); ++i) + for (size_t j = 0; j < Left.numberOfColumns(); ++j) { + M(i, j) = signe(Diag_vp[i]) * Right(i, j); + } + return Left * M; + }; + + 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::wall: { + if constexpr (Dimension == 2) { + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); + } else { + static_assert(Dimension == 3); + bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()), + getMeshFaceBoundary(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 Roe v2 Eulerian Composite solver"; + throw NormalError(error_msg.str()); + } + } + + return bc_list; + } + + 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] = stateFace[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] = stateFace[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 + _applyWallBoundaryCondition(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<WallBoundaryCondition, T>) { + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + std::cout << " Traitement WALL local (non flat) \n"; + // const Rd& normal = bc.outgoingNormal(); + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& node_to_face_matrix = mesh.connectivity().nodeToFaceMatrix(); + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& node_local_numbers_in_their_faces = mesh.connectivity().nodeLocalNumbersInTheirFaces(); + 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 Cjr = mesh_data.Cjr(); + const auto Cjf = mesh_data.Cjf(); + + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + const NodeId node_id = node_list[i_node]; + + const auto& node_face_list = node_to_face_matrix[node_id]; + // Assert(face_cell_list.size() == 1); + // const auto& node_local_number_in_its_faces = node_local_numbers_in_their_faces.itemArray(node_id); + + // on va chercher les normale d'une face issue du noeud de CL et contenue dans le faceList + Rd normal(zero); + int nbnormal = 0; + + for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) { + FaceId node_face_id = node_face_list[i_face]; + + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (node_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_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(facebc_id, 0); + + // Normal locale approchée + Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell); + normalloc *= 1. / l2Norm(normalloc); + normal += normalloc; + ++nbnormal; + + break; + } + } + } + + if (nbnormal == 0) + continue; + normal *= 1. / nbnormal; + + normal *= 1. / l2Norm(normal); + 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); + + // Normal locale approchée + Rd normal(Cjf(face_cell_id, face_local_number_in_cell)); + normal *= 1. / l2Norm(normal); + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateFace[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_to_face_matrix = mesh.connectivity().edgeToFaceMatrix(); + + const auto& edge_local_numbers_in_their_faces = mesh.connectivity().edgeLocalNumbersInTheirFaces(); + + 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_face_list = edge_to_face_matrix[edge_id]; + + // on va chercher les normale d'une face issue du edge de CL et contenue dans le faceList + Rd normal(zero); + int nbnormal = 0; + for (size_t i_face = 0; i_face < edge_face_list.size(); ++i_face) { + FaceId edge_face_id = edge_face_list[i_face]; + + for (size_t i_facebc = 0; i_facebc < face_list.size(); ++i_facebc) { + const FaceId facebc_id = face_list[i_facebc]; + if (edge_face_id == facebc_id) { + const auto& face_cell_list = face_to_cell_matrix[facebc_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(facebc_id, 0); + + // Normal locale approchée + Rd normalloc = Cjf(face_cell_id, face_local_number_in_cell); + normalloc *= 1. / l2Norm(normalloc); + normal += normalloc; + ++nbnormal; + break; + } + } + } + + if (nbnormal == 0) + continue; + normal *= 1. / nbnormal; + + normal *= 1. / l2Norm(normal); + + 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); + + 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; + } + + enum TypeAverageState + { + RoeAverage, + MeanAverage + }; + struct AverageStateStructData + { + double rho; + Rd U; + double E; + double H; + double gamma; + double p; + double c; + AverageStateStructData(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) + {} + }; + + AverageStateStructData + 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 AverageStateStructData(rho_mean, U_mean, E_mean, H_mean, gamma, P_mean, c_mean); + } + + AverageStateStructData + MeanAverageState(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 rho_mean = .5 * (rhoG + rhoD); + Rd U_mean = .5 / (rhoG + rhoD) * (rhoG * UG + rhoD * 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 H_mean = .5 / (rhoG + rhoD) * (rhoG * TotEnthalpyG + rhoD * TotEnthalpyD); + + 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 AverageStateStructData(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 AverageStateStructData& 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); + + auto bc_list = this->_getBCList(*p_mesh, bc_descriptor_list); + + auto rhoE = rho * E; + auto rhoU = rho * u; + + // Construction du vecteur des variables conservatives + // + // Ci dessous juste ordre 1 + // + CellValue<Rp> State{p_mesh->connectivity()}; + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + State[j][0] = rho[j]; + for (size_t d = 0; d < Dimension; ++d) + State[j][1 + d] = rhoU[j][d]; + State[j][1 + Dimension] = rhoE[j]; + }); + + // CellValue<Rp> State{p_mesh->connectivity()}; + NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()}; + StateAtNode.fill(zero); + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtNode[j].fill(State[j]); }); + EdgeValuePerCell<Rp> StateAtEdge{p_mesh->connectivity()}; + StateAtEdge.fill(zero); + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtEdge[j].fill(State[j]); }); + FaceValuePerCell<Rp> StateAtFace{p_mesh->connectivity()}; + StateAtFace.fill(zero); + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtFace[j].fill(State[j]); }); + + // Conditions aux limites des etats conservatifs + + _applyInflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + //_applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + _applySymmetricBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + _applyNeumannflatBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + _applyWallBoundaryCondition(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 "signe" 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 Gr, Gf + + NodeValue<Rpxd> Gr{p_mesh->connectivity()}; + Gr.fill(zero); + EdgeValue<Rpxd> Ge{p_mesh->connectivity()}; + Ge.fill(zero); + FaceValue<Rpxd> Gf{p_mesh->connectivity()}; + Gf.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); + + // NodeValuePerCell<Rpxp> Lambda_rj{p_mesh->connectivity()}; + CellValuePerNode<Rpxp> Lambda_rj{p_mesh->connectivity()}; + Lambda_rj.fill(zero); + auto is_boundary_node = p_mesh->connectivity().isBoundaryNode(); + + parallel_for( + p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { + const auto& node_to_cell = node_to_cell_matrix[r]; + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r); + + if (1 == 0) { + // if (is_boundary_node[r]) { + for (size_t l = 0; l < node_to_cell.size(); ++l) { + const CellId& j = node_to_cell[l]; + const size_t R = node_local_number_in_its_cells[l]; + + // const Rd& Cjr_loc = Cjr(j, R); + // const double norm_Cjr_loc = l2Norm(Cjr_loc); + // const Rd& Njr_loc = 1. / norm_Cjr_loc * Cjr_loc; + + const double rho_bd = StateAtNode[j][R][0]; + Rd U_bd; + for (size_t dim = 0; dim < Dimension; ++dim) + U_bd[dim] = StateAtNode[j][R][dim + 1] / rho_bd; + const double E_bd = StateAtNode[j][R][Dimension + 1] / rho_bd; + const double p_bd = pression(rho_bd, E_bd - .5 * dot(U_bd, U_bd), gamma); + AverageStateStructData RoeS( + RoeAverageState(rho_bd, U_bd, E_bd, gamma, p_bd, rho_bd, U_bd, E_bd, gamma, p_bd)); + // MeanAverageState(rho_bd, U_bd, E_bd, gamma, p_bd, rho_bd, U_bd, E_bd, gamma, p_bd)); + + const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, njr(j, R)); + Lambda_rj[r][l] = SignofASquareDiagonalisableMatrix(JacInfoK.RightEigenVectors, JacInfoK.EigenValues, + JacInfoK.LeftEigenVectors); + } + } else { + Rpxp somme_Lambda = zero; + + for (size_t l = 0; l < node_to_cell.size(); ++l) { + const CellId& j = node_to_cell[l]; + const size_t R = node_local_number_in_its_cells[l]; + + Rpxp Lambda_risum = zero; + // for (size_t k = 0; k < l; ++k) { + for (size_t k = 0; k < node_to_cell.size(); ++k) { + if (k == l) + continue; + const CellId K = node_to_cell[k]; + + // Etat de Roe entre j et K + // Puis Eval de la Jacobienne dans la dir. njr + // Ou bien autre Moyenne + AverageStateStructData 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])); + // MeanAverageState(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])); + + const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, njr(j, R)); + Lambda_risum += SignofASquareDiagonalisableMatrix(JacInfoK.RightEigenVectors, JacInfoK.EigenValues, + JacInfoK.LeftEigenVectors); + } + + // if (node_to_cell.size() > 1) + Lambda_risum *= 1. / std::max(1., (node_to_cell.size() - 1.)); + Lambda_rj[r][l] = Lambda_risum; + + somme_Lambda += Lambda_risum; + } + + // somme_Lambda *= (-0.5 / (node_to_cell.size() - 1.)); + // for (size_t l = 0; l < node_to_cell.size(); ++l) // boucle sur les mailles voisines + // { + // Lambda_rj[r][l] *= (node_to_cell.size() * 0.5 / (node_to_cell.size() - 1.)); + // Lambda_rj[r][l] += somme_Lambda; + // } // fin boucle sur i + + somme_Lambda *= -1. / (node_to_cell.size()); + for (size_t l = 0; l < node_to_cell.size(); ++l) // boucle sur les mailles voisines + { + Lambda_rj[r][l] += somme_Lambda; + } + } + }); + + synchronize(Lambda_rj); + + if (checkLocalConservation) { + NodeValue<double> MaxFrobNormLambda{p_mesh->connectivity()}; + MaxFrobNormLambda.fill(0.); + parallel_for( + p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { + if (not is_boundary_node[r]) { + const auto& node_to_cell = node_to_cell_matrix[r]; + Rpxp somme_Lambda = zero; + for (size_t l = 0; l < node_to_cell.size(); l++) // boucle sur les mailles voisines + somme_Lambda += Lambda_rj[r][l]; + + const double FrobNorm = frobeniusNorm(somme_Lambda); + MaxFrobNormLambda[r] = FrobNorm; + if (FrobNorm > 1e-9) { + std::cout << " Somme matrice noeud " << somme_Lambda << "\n"; + throw NormalError("Problem of TLBL at node"); + } + } + }); + std::cout << " Max Frobenius Norm (node) " << max(MaxFrobNormLambda) << "\n"; + } + + // Ici on met jf .. + FaceValuePerCell<Rpxp> Lambda_jf{p_mesh->connectivity()}; + Lambda_jf.fill(zero); + auto is_boundary_face = p_mesh->connectivity().isBoundaryFace(); + + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_to_face = cell_to_face_matrix[j]; + // const auto& cell_local_number_in_its_faces = cell_local_number_in_its_faces.itemArray(j); + + for (size_t l = 0; l < cell_to_face.size(); ++l) { + const FaceId& face = cell_to_face[l]; + if (is_boundary_face[face]) { + const double rho_bd = StateAtFace[j][l][0]; + Rd U_bd; + for (size_t dim = 0; dim < Dimension; ++dim) + U_bd[dim] = StateAtFace[j][l][dim + 1] / rho_bd; + const double E_bd = StateAtFace[j][l][Dimension + 1] / rho_bd; + const double p_bd = pression(rho_bd, E_bd - .5 * dot(U_bd, U_bd), gamma); + AverageStateStructData RoeS( + RoeAverageState(rho_bd, U_bd, E_bd, gamma, p_bd, rho_bd, U_bd, E_bd, gamma, p_bd)); + // MeanAverageState(rho_bd, U_bd, E_bd, gamma, p_bd, rho_bd, U_bd, E_bd, gamma, p_bd)); + + const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, njf(j, l)); + + Lambda_jf[j][l] = SignofASquareDiagonalisableMatrix(JacInfoK.RightEigenVectors, JacInfoK.EigenValues, + JacInfoK.LeftEigenVectors); + // continue; + } else { + 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); + // const double norm_Cjf_loc = l2Norm(Cjf_loc); + // const Rd& Njf_loc = njf(j, l); + + CellId K = face_to_cell[0]; + // unsigned int R = face_local_number_in_its_cells[0]; + + if (face_to_cell.size() == 1) { + K = j; + // R = l; + } else { + const CellId K1 = face_to_cell[0]; + const CellId K2 = face_to_cell[1]; + + if (j == K1) { + K = K2; + // R = face_local_number_in_its_cells[1]; + } + } + + AverageStateStructData 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])); + // MeanAverageState(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])); + + const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, njf(j, l)); + + Lambda_jf[j][l] = SignofASquareDiagonalisableMatrix(JacInfoK.RightEigenVectors, JacInfoK.EigenValues, + JacInfoK.LeftEigenVectors); + } + } + }); + synchronize(Lambda_jf); + + // + // Like node .. node / cells + // + CellValuePerEdge<Rpxp> Lambda_ej{p_mesh->connectivity()}; + Lambda_ej.fill(zero); + auto is_boundary_edge = p_mesh->connectivity().isBoundaryEdge(); + + const EdgeValuePerCell<const Rd> Cje = mesh_data.Cje(); + const EdgeValuePerCell<const Rd> nje = mesh_data.nje(); + + 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(); + + parallel_for( + p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId e) { + const auto& edge_to_cell = edge_to_cell_matrix[e]; + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(e); + + Rpxp somme_Lambda = zero; + + for (size_t l = 0; l < edge_to_cell.size(); ++l) { + const CellId& j = edge_to_cell[l]; + const size_t R = edge_local_number_in_its_cells[l]; + // const Rd& Cje_loc = Cje(j, R); + // const double& norm_Cje_loc = l2Norm(Cje_loc); + // const Rd& Nje_loc = 1. / (norm_Cje_loc)*Cje_loc; + + Rpxp Lambda_risum = zero; + + for (size_t k = 0; k < edge_to_cell.size(); ++k) { + if (k == l) + continue; + const CellId K = edge_to_cell[k]; + // const size_t s = edge_local_number_in_its_cells[k]; + // const Rd& Cks_loc = Cje(K, s); + // const double norm_Cks_loc = l2Norm(Cks_loc); + // const Rd& Nks_loc = 1. / norm_Cks_loc * Cks_loc; + + // Etat de Roe entre j et K + // Puis Eval de la Jacobienne dans la dir. Nks_loc + // Ou bien autre Moyenne + AverageStateStructData 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])); + // MeanAverageState(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])); + + const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, nje(j, R)); + Lambda_risum += SignofASquareDiagonalisableMatrix(JacInfoK.LeftEigenVectors, JacInfoK.EigenValues, + JacInfoK.RightEigenVectors); + } + // if (edge_to_cell.size() > 1) + Lambda_risum *= 1. / std::max(1., (edge_to_cell.size() - 1.)); + Lambda_ej[e][l] = Lambda_risum; + + somme_Lambda += Lambda_risum; + } + // somme_Lambda *= (-0.5 / (node_to_cell.size() - 1.)); + // for (size_t l = 0; l < node_to_cell.size(); ++l) // boucle sur les mailles voisines + // { + // Lambda_rj[r][l] *= (node_to_cell.size() * 0.5 / (node_to_cell.size() - 1.)); + // Lambda_rj[r][l] += somme_Lambda; + // } // fin boucle sur i + + somme_Lambda *= -1. / (edge_to_cell.size()); + for (size_t l = 0; l < edge_to_cell.size(); l++) // boucle sur les mailles voisines + { + Lambda_ej[e][l] += somme_Lambda; + } + }); + synchronize(Lambda_ej); + } + if (checkLocalConservation) { + const auto& edge_to_cell_matrix = p_mesh->connectivity().edgeToCellMatrix(); + + EdgeValue<double> MaxFrobNormLambda{p_mesh->connectivity()}; + MaxFrobNormLambda.fill(0.); + // auto is_boundary_edge = p_mesh->connectivity().isBoundaryEdge(); + parallel_for( + p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId e) { + if (not is_boundary_edge[e]) { + const auto& edge_to_cell = edge_to_cell_matrix[e]; + Rpxp somme_Lambda = zero; + for (size_t l = 0; l < edge_to_cell.size(); l++) // boucle sur les mailles voisines + somme_Lambda += Lambda_ej[e][l]; + const double FrobNorm = frobeniusNorm(somme_Lambda); + MaxFrobNormLambda[e] = FrobNorm; + + if (FrobNorm > 1e-9) { + std::cout << " Somme matrice edge " << somme_Lambda << "\n"; + throw NormalError("Problem of TLBL at edge"); + } + } + }); + std::cout << " Max Frobenius Norm (edge) " << max(MaxFrobNormLambda) << "\n"; + } + + NodeValue<int> nbChangingSignNode{p_mesh->connectivity()}; + nbChangingSignNode.fill(0); + + parallel_for( + p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node) { + 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); + // if (is_boundary_node[node]) { + if (1 == 0) { + Rpxd FluxMoy; + FluxMoy = zero; + Rp Umoy = zero; + for (size_t l = 0; l < node_to_cell.size(); ++l) { + const CellId& j = node_to_cell[l]; + const size_t R = node_local_number_in_its_cells[l]; + + const Rpxp Lambda_rl = Lambda_rj[node][l]; + Rpxp MatriceIDPlusLambda_rj = identity; + MatriceIDPlusLambda_rj += Lambda_rl; + + Umoy += StateAtNode[j][R]; + + Rpxd Flux; + for (size_t d = 0; d < Dimension; ++d) { + Flux(0, d) = Flux_rhoAtCellNode[j][R][d]; // dot(Flux_rhoAtCellNode[j][R], Cjr_loc); + // dot(Flux_rho[K], Cjr_loc); + for (size_t dd = 0; dd < Dimension; ++dd) + Flux(dd + 1, d) = Flux_qtmvtAtCellNode[j][R](dd, d); // u_Cjr[d]; + + Flux(Dimension + 1, d) = Flux_totnrjAtCellNode[j][R][d]; // dot(Flux_totnrjAtCellNode[j][R], Cjr_loc); + // dot(Flux_totnrj[K], Cjr_loc); + } + FluxMoy += MatriceIDPlusLambda_rj * Flux; + } + // Gr[node] = (1. / node_to_cell.size()) * FluxMoy; // MatriceIDPlusLambda_rj * Flux; + Umoy *= 1. / node_to_cell.size(); + const double rhonode = Umoy[0]; + Rd rhoUnode; + for (size_t d = 0; d < Dimension; ++d) + rhoUnode[d] = Umoy[d + 1]; + const double rhoEnode = Umoy[Dimension + 1]; + const double Enode = rhoEnode / rhonode; + const Rd Unode = 1. / rhonode * rhoUnode; + const double epsilonnode = Enode - .5 * dot(Unode, 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; + + Rpxd Flux_Umoy; + for (size_t d = 0; d < Dimension; ++d) { + Flux_Umoy(0, d) = rhoUnode[d]; + for (size_t dd = 0; dd < Dimension; ++dd) + Flux_Umoy(dd + 1, d) = rhoUtensUPlusPidnode(dd, d); + Flux_Umoy(Dimension + 1, d) = rhoEnodePlusP * Unode[d]; + } + Gr[node] = Flux_Umoy; //(1. / node_to_cell.size()) * FluxMoy; // MatriceIDPlusLambda_rj * Flux; + } else { + for (size_t l = 0; l < node_to_cell.size(); ++l) { + const CellId& j = node_to_cell[l]; + const size_t R = node_local_number_in_its_cells[l]; + // const Rd& Cjr_loc = Cjr(j, R); + const Rpxp Lambda_rl = Lambda_rj[node][l]; + + Rpxp MatriceIDPlusLambda_rj = identity; + MatriceIDPlusLambda_rj += Lambda_rl; + + Rpxd Flux; + for (size_t d = 0; d < Dimension; ++d) { + Flux(0, d) = Flux_rhoAtCellNode[j][R][d]; + for (size_t dd = 0; dd < Dimension; ++dd) + Flux(dd + 1, d) = Flux_qtmvtAtCellNode[j][R](dd, d); + Flux(Dimension + 1, d) = Flux_totnrjAtCellNode[j][R][d]; + } + Gr[node] += MatriceIDPlusLambda_rj * Flux; + } + Gr[node] *= 1. / (node_to_cell.size()); + } + }); + synchronize(Gr); + + if (checkLocalConservation) { + 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); + double maxGjrAbs = 0; + 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 += Gr[l] * Cjr(K, R); + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gr[l] * Cjr(K, R))); + } + // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr)); + MaxErrorConservationNode[l] = l2Norm(SumGjr) / maxGjrAbs; + } + }); + 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); + const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed(); + + parallel_for( + p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId f) { + const auto& face_to_cell = face_to_cell_matrix[f]; + const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(f); + + // if (1 == 0) { //} + if (is_boundary_face[f]) { + const CellId& j = face_to_cell[0]; + const size_t R = face_local_number_in_its_cells[0]; + + Rpxd FluxJ; + for (size_t d = 0; d < Dimension; ++d) { + FluxJ(0, d) = Flux_rhoAtCellFace[j][R][d]; + for (size_t dd = 0; dd < Dimension; ++dd) + FluxJ(dd + 1, d) = Flux_qtmvtAtCellFace[j][R](dd, d); + FluxJ(Dimension + 1, d) = Flux_totnrjAtCellFace[j][R][d]; + } + + Gf[f] = FluxJ; // MatriceIdplusLambda * FluxJ; + } else { + for (size_t l = 0; l < face_to_cell.size(); ++l) { + const CellId& j = face_to_cell[l]; + const size_t R = face_local_number_in_its_cells[l]; + + const Rpxp& Lambda_fj = Lambda_jf[j][R]; + + Rpxp MatriceIDPlusLambda_fj = identity; + MatriceIDPlusLambda_fj += Lambda_fj; + + // const Rp& statediff = StateAtNode[j][l] - StateAtNode[K][R]; // State[j] - State[K]; + // const Rp& diff = ViscosityMatrixJK * statediff; + + Rpxd Flux; + for (size_t d = 0; d < Dimension; ++d) { + Flux(0, d) = Flux_rhoAtCellFace[j][R][d]; + for (size_t dd = 0; dd < Dimension; ++dd) + Flux(dd + 1, d) = Flux_qtmvtAtCellFace[j][R](dd, d); + Flux(Dimension + 1, d) = Flux_totnrjAtCellFace[j][R][d]; + } + + Gf[f] += MatriceIDPlusLambda_fj * Flux; + } + Gf[f] *= 1. / (face_to_cell.size()); + } + }); + synchronize(Gf); + + if (checkLocalConservation) { + 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); + double maxGjrAbs = 0; + 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 += Gf[l] * Cjf(K, R); + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gf[l] * Cjf(K, R))); + } + MaxErrorConservationFace[l] = l2Norm(SumGjf) / maxGjrAbs; + // 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(); + + EdgeValue<int> nbChangingSignEdge(p_mesh->connectivity()); + nbChangingSignEdge.fill(0); + + parallel_for( + p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId e) { + // Edge + + const auto& edge_to_cell = edge_to_cell_matrix[e]; + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(e); + + for (size_t l = 0; l < edge_to_cell.size(); ++l) { + const CellId& j = edge_to_cell[l]; + const size_t R = edge_local_number_in_its_cells[l]; + // const Rd& Cjr_loc = Cjr(j, R); + const Rpxp Lambda_el = Lambda_ej[e][l]; + + Rpxp MatriceIDPlusLambda_ej = identity; + MatriceIDPlusLambda_ej += Lambda_el; + + Rpxd Flux; + for (size_t d = 0; d < Dimension; ++d) { + Flux(0, d) = Flux_rhoAtCellEdge[j][R][d]; + for (size_t dd = 0; dd < Dimension; ++dd) + Flux(dd + 1, d) = Flux_qtmvtAtCellEdge[j][R](dd, d); + Flux(Dimension + 1, d) = Flux_totnrjAtCellEdge[j][R][d]; + } + Ge[e] += MatriceIDPlusLambda_ej * Flux; + } + Ge[e] *= 1. / (edge_to_cell.size()); + /* + 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 + AverageStateStructData 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])); + // MeanAverageState(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; // false; // 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(); + } + */ + }); + synchronize(Ge); + + if (checkLocalConservation) { + 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); + double maxGjrAbs = 0; + 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 += Ge[l] * Cje(K, R); + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Ge[l] * Cje(K, R))); + } + // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2Norm(SumGje)); + MaxErrorConservationEdge[l] = l2Norm(SumGje) / maxGjrAbs; + } + }); + 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.; // std::atan(1); // pi / 4.; // 2. / 3.; // 2. / 3.; // 2. / 3.; //.5;sth + double eta = 1. / 6.; //.2; + if constexpr (Dimension == 2) { + eta = 0; + // theta = 1; // pour schema aux aretes + // theta=0; //pour schema aux noeuds + } + + // 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) { + const NodeId& node = cell_to_node[l]; + SumFluxesNode += Gr[node] * Cjr(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) { + const EdgeId& edge = cell_to_edge[l]; + + SumFluxesEdge += Ge[edge] * Cje(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) { + const FaceId& face = cell_to_face[l]; + + SumFluxesFace += Gf[face] * Cjf(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)); + } + + RoeFluxFormEulerianCompositeSolver_v2() = default; + ~RoeFluxFormEulerianCompositeSolver_v2() = default; +}; + +template <MeshConcept MeshType> +class RoeFluxFormEulerianCompositeSolver_v2<MeshType>::WallBoundaryCondition +{ +}; + +template <> +class RoeFluxFormEulerianCompositeSolver_v2<Mesh<2>>::WallBoundaryCondition +{ + 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(); + } + + WallBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) + : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary) + { + ; + } +}; + +template <> +class RoeFluxFormEulerianCompositeSolver_v2<Mesh<3>>::WallBoundaryCondition +{ + 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(); + } + + WallBoundaryCondition(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 RoeFluxFormEulerianCompositeSolver_v2<MeshType>::NeumannflatBoundaryCondition +{ +}; +template <> +class RoeFluxFormEulerianCompositeSolver_v2<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 RoeFluxFormEulerianCompositeSolver_v2<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 RoeFluxFormEulerianCompositeSolver_v2<MeshType>::SymmetryBoundaryCondition +{ +}; + +template <> +class RoeFluxFormEulerianCompositeSolver_v2<Mesh<2>>::SymmetryBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary; + const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary; + + public: + const Rd& + outgoingNormal() const + { + return m_mesh_flat_node_boundary.outgoingNormal(); + } + + size_t + numberOfNodes() const + { + return m_mesh_flat_node_boundary.nodeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_flat_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_flat_node_boundary.nodeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_flat_face_boundary.faceList(); + } + + SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary, + const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary) + : m_mesh_flat_node_boundary(mesh_flat_node_boundary), m_mesh_flat_face_boundary(mesh_flat_face_boundary) + { + ; + } + + ~SymmetryBoundaryCondition() = default; +}; + +template <> +class RoeFluxFormEulerianCompositeSolver_v2<Mesh<3>>::SymmetryBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary; + const MeshFlatEdgeBoundary<MeshType> m_mesh_flat_edge_boundary; + const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary; + + public: + const Rd& + outgoingNormal() const + { + return m_mesh_flat_node_boundary.outgoingNormal(); + } + + size_t + numberOfNodes() const + { + return m_mesh_flat_node_boundary.nodeList().size(); + } + + size_t + numberOfEdges() const + { + return m_mesh_flat_edge_boundary.edgeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_flat_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_flat_node_boundary.nodeList(); + } + + const Array<const EdgeId>& + edgeList() const + { + return m_mesh_flat_edge_boundary.edgeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_flat_face_boundary.faceList(); + } + + SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary, + const MeshFlatEdgeBoundary<MeshType>& mesh_flat_edge_boundary, + const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary) + : m_mesh_flat_node_boundary(mesh_flat_node_boundary), + m_mesh_flat_edge_boundary(mesh_flat_edge_boundary), + m_mesh_flat_face_boundary(mesh_flat_face_boundary) + { + ; + } + + ~SymmetryBoundaryCondition() = default; +}; + +template <MeshConcept MeshType> +class RoeFluxFormEulerianCompositeSolver_v2<MeshType>::InflowListBoundaryCondition +{ +}; + +template <> +class RoeFluxFormEulerianCompositeSolver_v2<Mesh<2>>::InflowListBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshNodeBoundary m_mesh_node_boundary; + const MeshFaceBoundary m_mesh_face_boundary; + const Table<const double> m_node_array_list; + const Table<const double> m_face_array_list; + + public: + size_t + numberOfNodes() const + { + return m_mesh_node_boundary.nodeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_node_boundary.nodeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_face_boundary.faceList(); + } + + const Table<const double>& + nodeArrayList() const + { + return m_node_array_list; + } + + const Table<const double>& + faceArrayList() const + { + return m_face_array_list; + } + + InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, + const MeshFaceBoundary& mesh_face_boundary, + const Table<const double>& node_array_list, + const Table<const double>& face_array_list) + : m_mesh_node_boundary(mesh_node_boundary), + m_mesh_face_boundary(mesh_face_boundary), + m_node_array_list(node_array_list), + m_face_array_list(face_array_list) + { + ; + } + + ~InflowListBoundaryCondition() = default; +}; + +template <> +class RoeFluxFormEulerianCompositeSolver_v2<Mesh<3>>::InflowListBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshNodeBoundary m_mesh_node_boundary; + const MeshEdgeBoundary m_mesh_edge_boundary; + const MeshFaceBoundary m_mesh_face_boundary; + const Table<const double> m_node_array_list; + const Table<const double> m_edge_array_list; + const Table<const double> m_face_array_list; + + public: + size_t + numberOfNodes() const + { + return m_mesh_node_boundary.nodeList().size(); + } + + size_t + numberOfEdges() const + { + return m_mesh_edge_boundary.edgeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_node_boundary.nodeList(); + } + + const Array<const EdgeId>& + edgeList() const + { + return m_mesh_edge_boundary.edgeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_face_boundary.faceList(); + } + + const Table<const double>& + nodeArrayList() const + { + return m_node_array_list; + } + + const Table<const double>& + edgeArrayList() const + { + return m_edge_array_list; + } + + const Table<const double>& + faceArrayList() const + { + return m_face_array_list; + } + + InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, + const MeshEdgeBoundary& mesh_edge_boundary, + const MeshFaceBoundary& mesh_face_boundary, + const Table<const double>& node_array_list, + const Table<const double>& edge_array_list, + const Table<const double>& face_array_list) + : m_mesh_node_boundary(mesh_node_boundary), + m_mesh_edge_boundary(mesh_edge_boundary), + m_mesh_face_boundary(mesh_face_boundary), + m_node_array_list(node_array_list), + m_edge_array_list(edge_array_list), + m_face_array_list(face_array_list) + { + ; + } + + ~InflowListBoundaryCondition() = default; +}; + +template <MeshConcept MeshType> +class RoeFluxFormEulerianCompositeSolver_v2<MeshType>::OutflowBoundaryCondition +{ +}; + +template <> +class RoeFluxFormEulerianCompositeSolver_v2<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 RoeFluxFormEulerianCompositeSolver_v2<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>> +roeFluxFormEulerianCompositeSolver_v2( + const std::shared_ptr<const DiscreteFunctionVariant>& rho_v, + const std::shared_ptr<const DiscreteFunctionVariant>& u_v, + const std::shared_ptr<const DiscreteFunctionVariant>& E_v, + const double& gamma, + const std::shared_ptr<const DiscreteFunctionVariant>& c_v, + const std::shared_ptr<const DiscreteFunctionVariant>& p_v, + // const 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("RoeFluxFormFormEulerianCompositeSolver v2 is not available in 1D"); + } else { + if constexpr (is_polygonal_mesh_v<MeshType>) { + return RoeFluxFormEulerianCompositeSolver_v2<MeshType>{} + .solve(p_mesh, rho_v->get<DiscreteFunctionP0<const double>>(), u_v->get<DiscreteFunctionP0<const Rd>>(), + E_v->get<DiscreteFunctionP0<const double>>(), gamma, c_v->get<DiscreteFunctionP0<const double>>(), + p_v->get<DiscreteFunctionP0<const double>>(), // degree, + bc_descriptor_list, dt, check); + } else { + throw NormalError("RoeFluxFormEulerianCompositeSolver v2 is only defined on polygonal meshes"); + } + } + }, + mesh_v->variant()); +} diff --git a/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.hpp b/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.hpp new file mode 100644 index 0000000000000000000000000000000000000000..e7f6519c64c63b6e0374a9f4c0efa1d345a3c43b --- /dev/null +++ b/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.hpp @@ -0,0 +1,89 @@ +#ifndef ROE_FLUX_FORM_EULERIAN_COMPOSITE_SOLVER_V2_HPP +#define ROE_FLUX_FORM_EULERIAN_COMPOSITE_SOLVER_V2_HPP + +#include <mesh/MeshTraits.hpp> + +#include "JacobianAndStructuralInfoForSystemsofEquations.hpp" +#include <scheme/RusanovEulerianCompositeSolverTools.hpp> + +#include <memory> +#include <tuple> +#include <vector> + +inline double +signe(const double a) +{ + if (a < 0) + return -1.; + if (a > 0) + return 1; + + return 0; +} +class Rp; +class Rpxp; + +std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>> +roeFluxFormEulerianCompositeSolver_v2( + const std::shared_ptr<const DiscreteFunctionVariant>& rho, + const std::shared_ptr<const DiscreteFunctionVariant>& u, + const std::shared_ptr<const DiscreteFunctionVariant>& E, + const double& gamma, + const std::shared_ptr<const DiscreteFunctionVariant>& c, + const std::shared_ptr<const DiscreteFunctionVariant>& p, + // const size_t& degree, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, + const double& dt, + const bool check = false); +/* +class DiscreteFunctionVariant; +class IBoundaryConditionDescriptor; +class MeshVariant; +template <size_t Dimension> +class Mesh; + +template <size_t Dimension, int number_of_unknowns> +class RoeFluxFormEulerianCompositeSolver_v2 +{ + private: + static constexpr int p = number_of_unknowns + (Dimension - 1); + using Rp = TinyVector<p>; // number_of_unknowns + (Dimension - 1)>; + + using Rpxp = TinyMatrix<p, p, double>; + // private: + public: + JacobianStructuralInfoSystemsOfEquations<Dimension, number_of_unknowns> JacobianInfos; + + // KN<std::vector<std::vector<Rn> > > _sign_jacobian_matrix_dof_in_cell; + // KN<std::vector<std::vector<Rn> > > _vp_jacobian_matrix_dof_in_cell; + // //KN<std::vector<Rnm> > _jac_jacobian_matrix_dof_in_cell; + // KN<std::vector<std::vector<Rnm> > > _right_eigenv_matrix_dof_in_cell; + // KN<std::vector<std::vector<Rnm> > > _left_eigenv_matrix_dof_in_cell; + // //KN<std::vector<Rnm> > _rotation_matrix_dof_in_cell; + // KN<std::vector<R2nm> > _MatriceU_dof_in_cell; + + NodeValuePerCell<const Rp> _vp_jacobian_matrix_node_in_cell; + EdgeValuePerCell<const Rp> _vp_jacobian_matrix_edge_in_cell; + FaceValuePerCell<const Rp> _vp_jacobian_matrix_face_in_cell; + + NodeValuePerCell<const Rpxp> _right_eigenv_matrix_node_in_cell; + EdgeValuePerCell<const Rpxp> _right_eigenv_matrix_edge_in_cell; + FaceValuePerCell<const Rpxp> _right_eigenv_matrix_face_in_cell; + + NodeValuePerCell<const Rpxp> _left_eigenv_matrix_node_in_cell; + EdgeValuePerCell<const Rpxp> _left_eigenv_matrix_edge_in_cell; + FaceValuePerCell<const Rpxp> _left_eigenv_matrix_face_in_cell; + + NodeValuePerCell<const Rpxp> _MatriceU_node_in_cell; + EdgeValuePerCell<const Rpxp> _MatriceU_edge_in_cell; + FaceValuePerCell<const Rpxp> _MatriceU_face_in_cell; + + public: + RoeFluxFormEulerianCompositeSolver_v2() {} + + ~RoeFluxFormEulerianCompositeSolver_v2() {} +}; +*/ +#endif