diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt index 13f1117862a81296211c5cead802b49268132a70..eb8969789477ee3f46bffcdf248fc7fa7f44eaaa 100644 --- a/src/scheme/CMakeLists.txt +++ b/src/scheme/CMakeLists.txt @@ -29,6 +29,7 @@ add_library( RoeViscousFormEulerianCompositeSolver_v2.cpp RusanovEulerianCompositeSolver_o2.cpp RusanovEulerianCompositeSolver_v2_o2.cpp + RusanovEulerianCompositeSolver_v2_order_n.cpp RoeViscousFormEulerianCompositeSolver_v2_o2.cpp ) diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2_order_n.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2_order_n.cpp new file mode 100644 index 0000000000000000000000000000000000000000..430c64fcb4576587755e4bdd801a07f4a314ef30 --- /dev/null +++ b/src/scheme/RusanovEulerianCompositeSolver_v2_order_n.cpp @@ -0,0 +1,2484 @@ +#include <scheme/RusanovEulerianCompositeSolver_v2_order_n.hpp> + +#include <analysis/GaussQuadratureDescriptor.hpp> +#include <analysis/QuadratureManager.hpp> +#include <geometry/SquareTransformation.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/StencilManager.hpp> +#include <mesh/SubItemValuePerItemUtils.hpp> +#include <scheme/DiscreteFunctionDPk.hpp> +#include <scheme/DiscreteFunctionDPkVariant.hpp> +#include <scheme/DiscreteFunctionDPkVector.hpp> +#include <scheme/DiscreteFunctionUtils.hpp> +#include <scheme/InflowListBoundaryConditionDescriptor.hpp> +#include <scheme/PolynomialReconstruction.hpp> +#include <utils/PugsTraits.hpp> +#include <variant> + +template <MeshConcept MeshTypeT> +class RusanovEulerianCompositeSolver_v2_order_n +{ + private: + using MeshType = MeshTypeT; + + static constexpr size_t Dimension = MeshType::Dimension; + + using Rdxd = TinyMatrix<Dimension>; + using Rd = TinyVector<Dimension>; + + using Rpxp = TinyMatrix<Dimension + 2>; + using Rp = TinyVector<Dimension + 2>; + + class SymmetryBoundaryCondition; + class InflowListBoundaryCondition; + class OutflowBoundaryCondition; + class WallBoundaryCondition; + class NeumannflatBoundaryCondition; + + using BoundaryCondition = std::variant<SymmetryBoundaryCondition, + InflowListBoundaryCondition, + OutflowBoundaryCondition, + NeumannflatBoundaryCondition, + WallBoundaryCondition>; + + using BoundaryConditionList = std::vector<BoundaryCondition>; + + const size_t m_quadrature_degree; + + 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 Rusanov v2 order n Eulerian Composite solver"; + throw NormalError(error_msg.str()); + } + } + + return bc_list; + } + + public: + CellByCellLimitation<MeshType> Limitor; + + public: + 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& 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 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_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); + + // Normal locale approchée + Rd normal(zero); + 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]; + normal += Cjr(node_cell_id, node_local_number_in_cell); + } + normal *= 1. / node_cell_list.size(); + normal *= 1. / l2Norm(normal); + + 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] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); + + const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); + + const auto& edge_list = bc.edgeList(); + + const auto Cje = mesh_data.Cje(); + + 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); + + // Normal locale approchée + Rd normal(zero); + 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]; + normal += Cje(edge_cell_id, edge_local_number_in_cell); + } + normal *= 1. / edge_cell_list.size(); + normal *= 1. / l2Norm(normal); + + 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 + _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list, + const MeshType& mesh, + NodeValuePerCell<Rp>& stateNode, + EdgeValuePerCell<Rp>& stateEdge, + FaceValuePerCell<Rp>& stateFace) const + { + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<OutflowBoundaryCondition, T>) { + std::cout << " Traitement Outflow \n"; + // const Rd& normal = bc.outgoingNormal(); + /* + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); + // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& face_list = bc.faceList(); + const auto& node_list = bc.nodeList(); + + const auto xj = mesh.xj(); + const auto xr = mesh.xr(); + const auto xf = mesh.xl(); + const auto xe = mesh.xe(); + + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + const NodeId node_id = node_list[i_node]; + + const auto& node_cell_list = node_to_cell_matrix[node_id]; + // Assert(face_cell_list.size() == 1); + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id); + + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + + for (size_t dim = 0; dim < Dimension + 2; ++dim) + stateNode[node_cell_id][node_local_number_in_cell][dim] += vectorSym[dim]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + + Rdxd MatriceProj(identity); + MatriceProj -= tensorProduct(normal, normal); + vectorSym = MatriceProj * vectorSym; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim]; + } + } + + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; + + const auto& face_cell_list = face_to_cell_matrix[face_id]; + Assert(face_cell_list.size() == 1); + + CellId face_cell_id = face_cell_list[0]; + size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0); + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim]; + + Rdxd MatriceProj(identity); + MatriceProj -= tensorProduct(normal, normal); + vectorSym = MatriceProj * vectorSym; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); + + const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); + + const auto& edge_list = bc.edgeList(); + + for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { + const EdgeId edge_id = edge_list[i_edge]; + + const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; + // Assert(face_cell_list.size() == 1); + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim]; + + Rdxd MatriceProj(identity); + MatriceProj -= tensorProduct(normal, normal); + vectorSym = MatriceProj * vectorSym; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + } + + // throw NormalError("Not implemented"); + } + */ + } + }, + boundary_condition); + } + } + + void + _applySymmetricBoundaryCondition(const BoundaryConditionList& bc_list, + const MeshType& mesh, + NodeValuePerCell<Rp>& stateNode, + EdgeValuePerCell<Rp>& stateEdge, + FaceValuePerCell<Rp>& stateFace) const + { + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) { + // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + std::cout << " Traitement SYMMETRY \n"; + const Rd& normal = bc.outgoingNormal(); + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); + // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& face_list = bc.faceList(); + const auto& node_list = bc.nodeList(); + + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + const NodeId node_id = node_list[i_node]; + + const auto& node_cell_list = node_to_cell_matrix[node_id]; + // Assert(face_cell_list.size() == 1); + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id); + + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + + Rdxd MatriceProj(identity); + MatriceProj -= tensorProduct(normal, normal); + vectorSym = MatriceProj * vectorSym; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim]; + } + } + + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; + + const auto& face_cell_list = face_to_cell_matrix[face_id]; + Assert(face_cell_list.size() == 1); + + CellId face_cell_id = face_cell_list[0]; + size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0); + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim]; + + Rdxd MatriceProj(identity); + MatriceProj -= tensorProduct(normal, normal); + vectorSym = MatriceProj * vectorSym; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); + + const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); + + const auto& edge_list = bc.edgeList(); + + for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { + const EdgeId edge_id = edge_list[i_edge]; + + const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; + // Assert(face_cell_list.size() == 1); + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim]; + + Rdxd MatriceProj(identity); + MatriceProj -= tensorProduct(normal, normal); + vectorSym = MatriceProj * vectorSym; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + } + } + } + }, + boundary_condition); + } + } + + void + _applyNeumannflatBoundaryCondition(const BoundaryConditionList& bc_list, + const MeshType& mesh, + NodeValuePerCell<Rp>& stateNode, + EdgeValuePerCell<Rp>& stateEdge, + FaceValuePerCell<Rp>& stateFace) const + { + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<NeumannflatBoundaryCondition, T>) { + // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + std::cout << " Traitement WALL \n"; + const Rd& normal = bc.outgoingNormal(); + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); + // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& face_list = bc.faceList(); + const auto& node_list = bc.nodeList(); + + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + const NodeId node_id = node_list[i_node]; + + const auto& node_cell_list = node_to_cell_matrix[node_id]; + // Assert(face_cell_list.size() == 1); + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id); + + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim]; + // stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim]; + } + } + + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; + + const auto& face_cell_list = face_to_cell_matrix[face_id]; + Assert(face_cell_list.size() == 1); + + CellId face_cell_id = face_cell_list[0]; + size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0); + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); + + const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); + + const auto& edge_list = bc.edgeList(); + + for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { + const EdgeId edge_id = edge_list[i_edge]; + + const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; + // Assert(face_cell_list.size() == 1); + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + + Rd vectorSym(zero); + for (size_t dim = 0; dim < Dimension; ++dim) + vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim]; + + vectorSym -= dot(vectorSym, normal) * normal; + + for (size_t dim = 0; dim < Dimension; ++dim) + stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim]; + } + } + } + } + }, + boundary_condition); + } + } + + void + _applyInflowBoundaryCondition(const BoundaryConditionList& bc_list, + const MeshType& mesh, + NodeValuePerCell<Rp>& stateNode, + EdgeValuePerCell<Rp>& stateEdge, + FaceValuePerCell<Rp>& stateFace) const + { + for (const auto& boundary_condition : bc_list) { + std::visit( + [&](auto&& bc) { + using T = std::decay_t<decltype(bc)>; + if constexpr (std::is_same_v<InflowListBoundaryCondition, T>) { + // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); + std::cout << " Traitement INFLOW \n"; + + const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); + const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); + + const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); + const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); + // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& face_list = bc.faceList(); + const auto& node_list = bc.nodeList(); + + const auto& face_array_list = bc.faceArrayList(); + const auto& node_array_list = bc.nodeArrayList(); + + for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { + const NodeId node_id = node_list[i_node]; + + const auto& node_cell_list = node_to_cell_matrix[node_id]; + // Assert(face_cell_list.size() == 1); + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id); + + for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { + CellId node_cell_id = node_cell_list[i_cell]; + size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; + + for (size_t dim = 0; dim < Dimension + 2; ++dim) + stateNode[node_cell_id][node_local_number_in_cell][dim] = node_array_list[i_node][dim]; + } + } + + for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { + const FaceId face_id = face_list[i_face]; + + const auto& face_cell_list = face_to_cell_matrix[face_id]; + Assert(face_cell_list.size() == 1); + + CellId face_cell_id = face_cell_list[0]; + size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0); + + for (size_t dim = 0; dim < Dimension + 2; ++dim) + stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][dim]; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); + + const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); + // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); + + const auto& edge_list = bc.edgeList(); + + const auto& edge_array_list = bc.edgeArrayList(); + + for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { + const EdgeId edge_id = edge_list[i_edge]; + + const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); + + // Assert(face_cell_list.size() == 1); + + for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { + CellId edge_cell_id = edge_cell_list[i_cell]; + size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; + + for (size_t dim = 0; dim < Dimension + 2; ++dim) + stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim]; + } + } + } + } + }, + boundary_condition); + } + } + + public: + inline double + pression(const double rho, const double epsilon, const double gam) const + { + return (gam - 1) * rho * epsilon; + } + + // void + // computeLimitorVolumicScalarQuantityMinModDukowicz(const DiscreteFunctionDPk<Dimension, double>& q_bar, + // CellValue<double>& Limitor_q) const + + std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>> + solve(const std::shared_ptr<const MeshType>& p_mesh, + const DiscreteFunctionP0<const double>& rho_n, + const DiscreteFunctionP0<const Rd>& u_n, + const DiscreteFunctionP0<const double>& E_n, + const double& gamma, + const DiscreteFunctionP0<const double>& c_n, + const DiscreteFunctionP0<const double>& p_n, + const size_t& degree, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, + const double& dt, + const bool checkLocalConservation) const + { + auto rho = copy(rho_n); + auto u = copy(u_n); + auto E = copy(E_n); + // auto c = copy(c_n); + // auto p = copy(p_n); + + std::cout << " degre " << degree << "\n"; + + std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list; + + for (auto&& bc_descriptor : bc_descriptor_list) { + if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) { + symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared()); + } + } + + PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::boundary, + std::max(size_t(1), degree), + symmetry_boundary_descriptor_list); + + std::cout << " Apres reconstuction descriptor " + << "\n"; + + auto _epsilon_limiter = [=, this](const MeshType& mesh, const DiscreteFunctionP0<const double>& epsilon, + auto epsilon_R, CellValue<double>& lambda_epsilon) { + 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() + .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 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); + } + 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 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 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, std::min(epsilon_x1, epsilon_x2))); + epsilon_R_max = std::max(epsilon_R_max, std::max(epsilon_x0, std::max(epsilon_x1, epsilon_x2))); + } + } else { + { + } + } + 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))); + }); + }; + + auto epsilon = E - 0.5 * dot(u, u); + + // assemblage direct + // auto reconstructions = PolynomialReconstruction{reconstruction_descriptor}.build(rho, rho * u, rho * E); + + // Pour assemblage Leibniz + // on passe par la densite et les variables specifiques + auto reconstructions = PolynomialReconstruction{reconstruction_descriptor}.build(rho, u, epsilon); + std::cout << " Apres reconstruction build " + << "\n"; + + // Fonction HP et D + auto remove_mean = [=, this]<typename DataType>(DiscreteFunctionDPk<Dimension, DataType>& HP_S) { + /* + const QuadratureFormula<Dimension> qf = + QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor(m_quadrature_degree)); + + auto cell_type = p_mesh->connectivity().cellType(); + + auto Vj = MeshDataManager::instance().getMeshData(*p_mesh).Vj(); + auto xj = MeshDataManager::instance().getMeshData(*p_mesh).xj(); + + auto xr = p_mesh->xr(); + auto xl = MeshDataManager::instance().getMeshData(*p_mesh).xl(); // p_mesh->xl(); + + auto cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix(); + auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix(); + + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + const auto HP_S_coefs = HP_S.coefficients(cell_id); + DataType value; + if constexpr (std::is_arithmetic_v<DataType>) { + value = 0; + } else { + value = zero; + } + HP_S_coefs[0] = value; + + switch (cell_type[cell_id]) { + case CellType::Quadrangle: { + const auto cell_to_node = cell_to_node_matrix[cell_id]; + const auto cell_to_face = cell_to_face_matrix[cell_id]; + + const Rd& a00 = xr[cell_to_node[0]]; + const Rd& a01 = xl[cell_to_face[3]][0]; + const Rd& a02 = xr[cell_to_node[3]]; + + const Rd& a10 = xl[cell_to_face[0]][0]; + const Rd& a11 = xj[cell_id]; + const Rd& a20 = xr[cell_to_node[1]]; + + const Rd& a21 = xl[cell_to_face[1]][0]; + const Rd& a12 = xl[cell_to_face[2]][0]; + const Rd& a22 = xr[cell_to_node[2]]; + + const SquareTransformation<2> T(a00, a01, a02, // + a10, a11, a12, // + a20, a21, a22); + + for (size_t i_point = 0; i_point < qf.numberOfPoints(); ++i_point) { + const auto xi = qf.point(i_point); + const auto& x = T(xi); + value += qf.weight(i_point) * T.jacobianDeterminant(xi) * HP_S[cell_id](x); + } + + break; + } + default: { + throw NotImplementedError("invalid cell type"); + } + } + + HP_S_coefs[0] = -1. / Vj[cell_id] * value; + }); + */ + return HP_S; + }; + + auto compute_HP = [=]<typename DataType>(const DiscreteFunctionDPk<Dimension, const double>& rho_L, + const DiscreteFunctionDPk<Dimension, DataType>& P_S) + -> DiscreteFunctionDPk<Dimension, std::remove_const_t<DataType>> { + DiscreteFunctionDPk<Dimension, std::remove_const_t<DataType>> HP_S{p_mesh, degree}; + + if constexpr (Dimension == 2) { + switch (degree) { + case 0: + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + auto HP_S_coefs = HP_S.coefficients(cell_id); + + if constexpr (std::is_arithmetic_v<DataType>) { + HP_S_coefs[0] = 0; + } else { + HP_S_coefs[0] = zero; + } + }); + break; + case 1: + // 1 x + // y + // + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + const auto rho_coefs = rho_L.coefficients(cell_id); + + const auto& rho_j = rho[cell_id]; + // const auto& rho_x = rho_coefs[1]; + // const auto& rho_y = rho_coefs[2]; + + const auto P_S_coefs = P_S.coefficients(cell_id); + + const auto& S_x = P_S_coefs[1]; + const auto& S_y = P_S_coefs[2]; + + auto HP_S_coefs = HP_S.coefficients(cell_id); + + if constexpr (std::is_arithmetic_v<DataType>) { + HP_S_coefs[0] = 0; + } else { + HP_S_coefs[0] = zero; + } + HP_S_coefs[1] = rho_j * S_x; + HP_S_coefs[2] = rho_j * S_y; + }); + break; + case 2: + // 1, x , x^2 + // y xy + // y^2 + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + const auto rho_coefs = rho_L.coefficients(cell_id); + + const auto& rho_j = rho[cell_id]; + const auto& rho_x = rho_coefs[1]; + const auto& rho_y = rho_coefs[3]; + + const auto P_S_coefs = P_S.coefficients(cell_id); + + const auto& S_x = P_S_coefs[1]; + const auto& S_xx = P_S_coefs[2]; + const auto& S_y = P_S_coefs[3]; + const auto& S_xy = P_S_coefs[4]; + const auto& S_yy = P_S_coefs[5]; + + auto HP_S_coefs = HP_S.coefficients(cell_id); + + if constexpr (std::is_arithmetic_v<DataType>) { + HP_S_coefs[0] = 0; + } else { + HP_S_coefs[0] = zero; + } + HP_S_coefs[1] = rho_j * S_x; + HP_S_coefs[2] = .5 * (rho_j * S_xx + 2 * rho_x * S_x); + HP_S_coefs[3] = rho_j * S_y; + HP_S_coefs[4] = .5 * (rho_j * S_xy + rho_x * S_y + rho_y * S_x); + HP_S_coefs[5] = .5 * (rho_j * S_yy + 2 * rho_y * S_y); + // HP_S_coefs *= .5; + }); + break; + case 3: + // 1, + // x , x^2, x^3 + // y, xy, x^2 y + // y^2, x y^2 + // y^3 + + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + const auto rho_coefs = rho_L.coefficients(cell_id); + + const auto& rho_j = rho[cell_id]; + + const auto& rho_x = rho_coefs[1]; // ro_xx coef2, ro_xxx coef3 + const auto& rho_xx = rho_coefs[2]; // ro_xx coef2, ro_xxx coef3 + const auto& rho_xxx = rho_coefs[3]; // ro_xx coef2, ro_xxx coef3 + + const auto& rho_y = rho_coefs[4]; // + const auto& rho_xy = rho_coefs[5]; // + const auto& rho_xxy = rho_coefs[6]; // + + const auto& rho_yy = rho_coefs[7]; // + const auto& rho_xyy = rho_coefs[8]; // + const auto& rho_yyy = rho_coefs[9]; // + + const auto P_S_coefs = P_S.coefficients(cell_id); + + const auto& S_x = P_S_coefs[1]; + const auto& S_xx = P_S_coefs[2]; + const auto& S_xxx = P_S_coefs[3]; + + const auto& S_y = P_S_coefs[4]; + const auto& S_xy = P_S_coefs[5]; + const auto& S_xxy = P_S_coefs[6]; + + const auto& S_yy = P_S_coefs[7]; + const auto& S_xyy = P_S_coefs[8]; + const auto& S_yyy = P_S_coefs[9]; + + auto HP_S_coefs = HP_S.coefficients(cell_id); + + if constexpr (std::is_arithmetic_v<DataType>) { + HP_S_coefs[0] = 0; + } else { + HP_S_coefs[0] = zero; + } + + // Cf formule.. (44) du CRAS + HP_S_coefs[1] = rho_j * S_x; + HP_S_coefs[2] = rho_j * S_xx + 0.5 * rho_x * S_x; + HP_S_coefs[3] = rho_j * S_y; + HP_S_coefs[4] = rho_j * S_xy + rho_x * S_y + rho_y * S_x; + HP_S_coefs[5] = rho_j * S_yy + 0.5 * rho_y * S_y; + HP_S_coefs[6] = HP_S_coefs[0]; + HP_S_coefs[7] = HP_S_coefs[0]; + HP_S_coefs[8] = HP_S_coefs[0]; + HP_S_coefs[9] = HP_S_coefs[0]; + }); + break; + default: { + throw NotImplementedError("NotImplementedError degree > 3"); + } + } + } else { + throw NotImplementedError("NotImplementedError in 3D"); + } + return HP_S; + }; + + auto compute_D = + [=](const DiscreteFunctionDPk<Dimension, const Rd>& P_u) -> DiscreteFunctionDPk<Dimension, double> { + DiscreteFunctionDPk<Dimension, double> D{p_mesh, degree}; + + switch (degree) { + case 0: + case 1: + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + const auto D_coefs = D.coefficients(cell_id); + + for (size_t i = 0; i < D_coefs.size(); ++i) + D_coefs[i] = 0; + }); + + break; + + case 2: + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { + const auto& rho_j = rho[cell_id]; + + const auto P_u_coefs = P_u.coefficients(cell_id); + + const auto& u_x = P_u_coefs[1]; // (u1,u2)_x + const auto& u_y = P_u_coefs[3]; // (u1,u2)_y + + const auto& u1x = u_x[0]; + const auto& u2x = u_x[1]; + + const auto& u1y = u_y[0]; + const auto& u2y = u_y[1]; + + const auto D_coefs = D.coefficients(cell_id); + + D_coefs[0] = 0; + D_coefs[1] = 0; + D_coefs[2] = .5 * rho_j * dot(u_x, u_x); + D_coefs[3] = 0; + D_coefs[4] = rho_j * (u1x * u1y + u2x * u2y); + D_coefs[5] = .5 * rho_j * dot(u_y, u_y); + }); + + break; + default: { + throw NotImplementedError("NotImplementedError degree > 3"); + } + } + return D; + }; + std::cout << " Apres HP, D " + << " \n"; + DiscreteFunctionDPk rho_bar = reconstructions[0]->template get<DiscreteFunctionDPk<Dimension, const double>>(); + + DiscreteFunctionDPk P_u = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>(); + DiscreteFunctionDPk P_epsilon = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>(); + + DiscreteFunctionDPk rho_L = copy(rho_bar); + Limitor.density_limiter(*p_mesh, symmetry_boundary_descriptor_list, rho, rho_L); + std::cout << " Apres Limitor density " + << " \n"; + + DiscreteFunctionDPk HP_epsilon = compute_HP(rho_L, P_epsilon); + std::cout << " Apres Hp eps " + << " \n"; + + remove_mean(HP_epsilon); + DiscreteFunctionDPk HP_u = compute_HP(rho_L, P_u); + std::cout << " Apres Hp u " + << " \n"; + + remove_mean(HP_u); + + std::cout << " Avant D " + << " \n"; + DiscreteFunctionDPk D = compute_D(P_u); + std::cout << " Apres D " + << " \n"; + + remove_mean(D); + std::cout << " Apres Comput HP_e, Hpu et D " + << " \n"; + + CellValue<double> lambda_epsilon(p_mesh->connectivity()); + lambda_epsilon.fill(1); + + auto epsilon_bar = [=](const CellId cell_id, const Rd& x) { + const double tau = 1. / rho_L[cell_id](x); + const Rd tau_HP_u = tau * HP_u[cell_id](x); + + return epsilon[cell_id] + lambda_epsilon[cell_id] * + (tau * HP_epsilon[cell_id](x) + tau * D[cell_id](x) - 0.5 * dot(tau_HP_u, tau_HP_u)); + }; + // specific_internal_nrj_limiter(const MeshType& mesh, const DiscreteFunctionP0<const double>& rho, + // const DiscreteFunctionDPk<Dimension, double>& rho_L, + // const DiscreteFunctionP0<const double>& epsilon, + // const DiscreteFunctionDPk<Dimension, double>& epsilon_R, + // CellValue<const double>& lambda_epsilon) + // Limitor.specific_internal_nrj_limiter(*p_mesh, epsilon, epsilon_bar, lambda_epsilon); + + std::cout << " Avant epsilon limit " + << "\n"; + _epsilon_limiter(*p_mesh, epsilon, epsilon_bar, lambda_epsilon); + // DiscreteFunctionDPk rho_u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>(); + // DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const + // double>>(); + + CellValue<double> lambda_u{p_mesh->connectivity()}; + parallel_for( + p_mesh->numberOfCells(), + PUGS_LAMBDA(const CellId cell_id) { lambda_u[cell_id] = std::sqrt(lambda_epsilon[cell_id]); }); + + // Du coup on a la vitesse + auto u_bar = [=](const CellId cell_id, const Rd& x) { + return u[cell_id] + lambda_u[cell_id] * (1. / rho_L[cell_id](x)) * HP_u[cell_id](x); + }; + + auto E_bar = // epsilon_bar + .5 * dot(u_bar, u_bar); + [=](const CellId cell_id, const Rd& x) { + return epsilon_bar(cell_id, x) + .5 * dot(u_bar(cell_id, x), u_bar(cell_id, x)); + }; + // Les var conservatives assemblees : rho U et rho E + auto rho_u_bar = // rho_L * u_bar; + [=](const CellId cell_id, const Rd& x) { + return rho_L[cell_id](x) * u[cell_id] + lambda_u[cell_id] * HP_u[cell_id](x); + }; + + auto rho_E_bar = // rho_L * E_bar; + [=](const CellId cell_id, const Rd& x) { return rho_L[cell_id](x) * E_bar(cell_id, x); }; + + // Pour les flux .. + + // auto p_bar = [&rho_L, &epsilon_bar](const CellId cell_id, const Rd& x) { + // const double rho_epsilon = (rho_L[cell_id](x) * epsilon_bar(cell_id, x)); + + // constexpr double gam = 1.4; + // return (gam - 1) * rho_epsilon; // pression(rho_L[cell_id](x), epsilon_bar(cell_id, x), gam); + // }; + + // + // Creation des variables conservatives limitées + // + + auto bc_list = this->_getBCList(*p_mesh, bc_descriptor_list); + + auto rhoE = rho * E; + auto rhoU = rho * u; + + // Construction du vecteur des variables conservatives + // + // Ci dessous juste ordre 1 + // + CellValue<Rp> State{p_mesh->connectivity()}; + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + State[j][0] = rho[j]; + for (size_t d = 0; d < Dimension; ++d) + State[j][1 + d] = rhoU[j][d]; + State[j][1 + Dimension] = rhoE[j]; + }); + + // CellValue<Rp> State{p_mesh->connectivity()}; + // + // Remplir ici les reconstructions d'ordre élevé + + // + const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix(); + const auto& cell_to_edge_matrix = p_mesh->connectivity().cellToEdgeMatrix(); + const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix(); + + const auto xr = p_mesh->xr(); + auto xl = MeshDataManager::instance().getMeshData(*p_mesh).xl(); + auto xe = MeshDataManager::instance().getMeshData(*p_mesh).xe(); + + NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()}; + StateAtNode.fill(zero); + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + // StateAtNode[j].fill(State[j]); + const auto& cell_to_node = cell_to_node_matrix[j]; + + for (size_t l = 0; l < cell_to_node.size(); ++l) { + const NodeId& node = cell_to_node[l]; + const double rhoj_node = rho_L[j](xr[node]); + StateAtNode[j][l][0] = rhoj_node; + for (size_t dim = 0; dim < Dimension; ++dim) + StateAtNode[j][l][1 + dim] = rho_u_bar(j, xr[node])[dim]; + StateAtNode[j][l][1 + Dimension] = rho_E_bar(j, xr[node]); + } + }); + + EdgeValuePerCell<Rp> StateAtEdge{p_mesh->connectivity()}; + StateAtEdge.fill(zero); + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + // eStateAtEdge[j].fill(State[j]); + const auto& cell_to_edge = cell_to_edge_matrix[j]; + + for (size_t l = 0; l < cell_to_edge.size(); ++l) { + const EdgeId& edge = cell_to_edge[l]; + + StateAtEdge[j][l][0] = rho_L[j](xe[edge]); + for (size_t dim = 0; dim < Dimension; ++dim) + StateAtEdge[j][l][1 + dim] = rho_u_bar(j, xe[edge])[dim]; + StateAtEdge[j][l][1 + Dimension] = rho_E_bar(j, xe[edge]); + } + }); + FaceValuePerCell<Rp> StateAtFace{p_mesh->connectivity()}; + StateAtFace.fill(zero); + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + // StateAtFace[j].fill(State[j]); + const auto& cell_to_face = cell_to_face_matrix[j]; + + for (size_t l = 0; l < cell_to_face.size(); ++l) { + const FaceId& face = cell_to_face[l]; + + StateAtFace[j][l][0] = rho_L[j](xl[face]); + for (size_t dim = 0; dim < Dimension; ++dim) + StateAtFace[j][l][1 + dim] = rho_u_bar(j, xl[face])[dim]; + StateAtFace[j][l][1 + Dimension] = rho_E_bar(j, xl[face]); + } + }); + + // Conditions aux limites des etats conservatifs + + _applyInflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + //_applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + _applySymmetricBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + _applyNeumannflatBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + + // + // Construction du flux .. ok pour l'ordre 1 + // + CellValue<Rdxd> rhoUtensU{p_mesh->connectivity()}; + CellValue<Rdxd> Pid(p_mesh->connectivity()); + Pid.fill(identity); + CellValue<Rdxd> rhoUtensUPlusPid{p_mesh->connectivity()}; + rhoUtensUPlusPid.fill(zero); + + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + rhoUtensU[j] = tensorProduct(rhoU[j], u[j]); + Pid[j] *= p_n[j]; + rhoUtensUPlusPid[j] = rhoUtensU[j] + Pid[j]; + }); + + auto Flux_rho = rhoU; + auto Flux_qtmvt = rhoUtensUPlusPid; // rhoUtensU + Pid; + auto Flux_totnrj = (rhoE + p_n) * u; + + // Flux avec prise en compte des states at Node/Edge/Face + + NodeValuePerCell<Rd> Flux_rhoAtCellNode{p_mesh->connectivity()}; + EdgeValuePerCell<Rd> Flux_rhoAtCellEdge{p_mesh->connectivity()}; + FaceValuePerCell<Rd> Flux_rhoAtCellFace{p_mesh->connectivity()}; + /* + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_to_node = cell_to_node_matrix[j]; + + for (size_t l = 0; l < cell_to_node.size(); ++l) { + for (size_t dim = 0; dim < Dimension; ++dim) + Flux_rhoAtCellNode[j][l][dim] = StateAtNode[j][l][0] * StateAtNode[j][l][dim]; + } + + const auto& cell_to_face = cell_to_face_matrix[j]; + + for (size_t l = 0; l < cell_to_face.size(); ++l) { + for (size_t dim = 0; dim < Dimension; ++dim) + Flux_rhoAtCellFace[j][l][dim] = StateAtFace[j][l][0] * StateAtFace[j][l][dim]; + } + + const auto& cell_to_edge = cell_to_edge_matrix[j]; + + for (size_t l = 0; l < cell_to_edge.size(); ++l) { + for (size_t dim = 0; dim < Dimension; ++dim) + Flux_rhoAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim]; + } + }); + */ + NodeValuePerCell<Rdxd> Flux_qtmvtAtCellNode{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid; + EdgeValuePerCell<Rdxd> Flux_qtmvtAtCellEdge{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid; + FaceValuePerCell<Rdxd> Flux_qtmvtAtCellFace{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid; + /* + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_to_node = cell_to_node_matrix[j]; + for (size_t l = 0; l < cell_to_node.size(); ++l) { + for (size_t dim = 0; dim < Dimension; ++dim) + Flux_qtmvtAtCellNode[j][l][dim] = StateAtNode[j][l][0] * StateAtNode[j][l][dim]; + } + + const auto& cell_to_face = cell_to_face_matrix[j]; + + for (size_t l = 0; l < cell_to_face.size(); ++l) { + for (size_t dim = 0; dim < Dimension; ++dim) + Flux_qtmvtAtCellFace[j][l][dim] = StateAtFace[j][l][0] * StateAtFace[j][l][dim]; + } + + const auto& cell_to_edge = cell_to_edge_matrix[j]; + + for (size_t l = 0; l < cell_to_edge.size(); ++l) { + for (size_t dim = 0; dim < Dimension; ++dim) + Flux_qtmvtAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim]; + } + }); + */ + // parallel_for( + // p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + // Flux_qtmvtAtCellNode[j] = Flux_qtmvtAtCellEdge[j] = Flux_qtmvtAtCellFace[j] = Flux_qtmvt[j]; + // }); + + NodeValuePerCell<Rd> Flux_totnrjAtCellNode{p_mesh->connectivity()}; + EdgeValuePerCell<Rd> Flux_totnrjAtCellEdge{p_mesh->connectivity()}; + FaceValuePerCell<Rd> Flux_totnrjAtCellFace{p_mesh->connectivity()}; + + Flux_rhoAtCellEdge.fill(zero); + Flux_totnrjAtCellEdge.fill(zero); + Flux_qtmvtAtCellEdge.fill(zero); + + // Les flux aux nodes/edge/faces + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_to_node = cell_to_node_matrix[j]; + + for (size_t l = 0; l < cell_to_node.size(); ++l) { + // Etats conservatifs aux noeuds + const double rhonode = StateAtNode[j][l][0]; + Rd Unode; + for (size_t dim = 0; dim < Dimension; ++dim) + Unode[dim] = StateAtNode[j][l][dim + 1] / rhonode; + const double rhoEnode = StateAtNode[j][l][Dimension + 1]; + // + const double Enode = rhoEnode / rhonode; + const double epsilonnode = Enode - .5 * dot(Unode, Unode); + const Rd rhoUnode = rhonode * Unode; + const Rdxd rhoUtensUnode = tensorProduct(rhoUnode, Unode); + + const double Pressionnode = pression(rhonode, epsilonnode, gamma); + + const double rhoEnodePlusP = rhoEnode + Pressionnode; + + Rdxd rhoUtensUPlusPidnode(identity); + rhoUtensUPlusPidnode *= Pressionnode; + rhoUtensUPlusPidnode += rhoUtensUnode; + + Flux_rhoAtCellNode[j][l] = rhoUnode; + Flux_qtmvtAtCellNode[j][l] = rhoUtensUPlusPidnode; + Flux_totnrjAtCellNode[j][l] = rhoEnodePlusP * Unode; + } + + const auto& cell_to_face = cell_to_face_matrix[j]; + + for (size_t l = 0; l < cell_to_face.size(); ++l) { + const double rhoface = StateAtFace[j][l][0]; + Rd Uface; + for (size_t dim = 0; dim < Dimension; ++dim) + Uface[dim] = StateAtFace[j][l][dim + 1] / rhoface; + const double rhoEface = StateAtFace[j][l][Dimension + 1]; + // + const double Eface = rhoEface / rhoface; + const double epsilonface = Eface - .5 * dot(Uface, Uface); + const Rd rhoUface = rhoface * Uface; + const Rdxd rhoUtensUface = tensorProduct(rhoUface, Uface); + + const double Pressionface = pression(rhoface, epsilonface, gamma); + + const double rhoEfacePlusP = rhoEface + Pressionface; + + Rdxd rhoUtensUPlusPidface(identity); + rhoUtensUPlusPidface *= Pressionface; + rhoUtensUPlusPidface += rhoUtensUface; + + Flux_rhoAtCellFace[j][l] = rhoUface; + Flux_qtmvtAtCellFace[j][l] = rhoUtensUPlusPidface; + Flux_totnrjAtCellFace[j][l] = rhoEfacePlusP * Uface; + } + + if constexpr (Dimension == 3) { + const auto& cell_to_edge = cell_to_edge_matrix[j]; + + for (size_t l = 0; l < cell_to_edge.size(); ++l) { + const double rhoedge = StateAtEdge[j][l][0]; + Rd Uedge; + for (size_t dim = 0; dim < Dimension; ++dim) + Uedge[dim] = StateAtEdge[j][l][dim + 1] / rhoedge; + const double rhoEedge = StateAtEdge[j][l][Dimension + 1]; + // + const double Eedge = rhoEedge / rhoedge; + const double epsilonedge = Eedge - .5 * dot(Uedge, Uedge); + const Rd rhoUedge = rhoedge * Uedge; + const Rdxd rhoUtensUedge = tensorProduct(rhoUedge, Uedge); + + const double Pressionedge = pression(rhoedge, epsilonedge, gamma); + + const double rhoEedgePlusP = rhoEedge + Pressionedge; + + Rdxd rhoUtensUPlusPidedge(identity); + rhoUtensUPlusPidedge *= Pressionedge; + rhoUtensUPlusPidedge += rhoUtensUedge; + + Flux_rhoAtCellEdge[j][l] = rhoUedge; + Flux_qtmvtAtCellEdge[j][l] = rhoUtensUPlusPidedge; + Flux_totnrjAtCellEdge[j][l] = rhoEedgePlusP * Uedge; + } + } + }); + + // parallel_for( + // p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + // Flux_totnrjAtCellNode[j] = Flux_totnrjAtCellEdge[j] = Flux_totnrjAtCellFace[j] = Flux_totnrj[j]; + // }); + + MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh); + + auto volumes_cell = mesh_data.Vj(); + + // Calcul des matrices de viscosité aux node/edge/face + + const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr(); + const NodeValuePerCell<const Rd> njr = mesh_data.njr(); + + const FaceValuePerCell<const Rd> Cjf = mesh_data.Cjf(); + const FaceValuePerCell<const Rd> njf = mesh_data.njf(); + + const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix(); + const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells(); + + const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix(); + const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells(); + + // We compute Gjr, Gjf + + NodeValuePerCell<Rp> Gjr{p_mesh->connectivity()}; + Gjr.fill(zero); + EdgeValuePerCell<Rp> Gje{p_mesh->connectivity()}; + Gje.fill(zero); + FaceValuePerCell<Rp> Gjf{p_mesh->connectivity()}; + Gjf.fill(zero); + + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_to_node = cell_to_node_matrix[j]; + + for (size_t l = 0; l < cell_to_node.size(); ++l) { + const NodeId& node = cell_to_node[l]; + const auto& node_to_cell = node_to_cell_matrix[node]; + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node); + + const Rd& Cjr_loc = Cjr(j, l); + + for (size_t k = 0; k < node_to_cell.size(); ++k) { + const CellId K = node_to_cell[k]; + const size_t R = node_local_number_in_its_cells[k]; + const Rd& Ckr_loc = Cjr(K, R); + + // Une moyenne entre les etats jk + + Rd uNode = .5 * (u_n[j] + u_n[K]); + double cNode = .5 * (c_n[j] + c_n[K]); + + // Viscosity j k + Rpxp ViscosityMatrixJK(identity); + const double MaxmaxabsVpNormjk = + std::max(toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uNode, cNode, + Cjr_loc), + toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uNode, cNode, + Ckr_loc)); + + ViscosityMatrixJK *= MaxmaxabsVpNormjk; + const Rd& u_Cjr = Flux_qtmvtAtCellNode[K][R] * Cjr_loc; // Flux_qtmvt[K] * Cjr_loc; + + const Rp& statediff = StateAtNode[j][l] - StateAtNode[K][R]; // State[j] - State[K]; + const Rp& diff = ViscosityMatrixJK * statediff; + + Gjr[j][l][0] += dot(Flux_rhoAtCellNode[K][R], Cjr_loc); // dot(Flux_rho[K], Cjr_loc); + for (size_t d = 0; d < Dimension; ++d) + Gjr[j][l][1 + d] += u_Cjr[d]; + Gjr[j][l][1 + Dimension] += dot(Flux_totnrjAtCellNode[K][R], Cjr_loc); // dot(Flux_totnrj[K], Cjr_loc); + + Gjr[j][l] += diff; + } + + Gjr[j][l] *= 1. / node_to_cell.size(); + } + }); + synchronize(Gjr); + if (checkLocalConservation) { + auto is_boundary_node = p_mesh->connectivity().isBoundaryNode(); + + NodeValue<double> MaxErrorConservationNode(p_mesh->connectivity()); + MaxErrorConservationNode.fill(0.); + // double MaxErrorConservationNode = 0; + parallel_for( + p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId l) { + const auto& node_to_cell = node_to_cell_matrix[l]; + const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(l); + + if (not is_boundary_node[l]) { + Rp SumGjr(zero); + for (size_t k = 0; k < node_to_cell.size(); ++k) { + const CellId K = node_to_cell[k]; + const unsigned int R = node_local_number_in_its_cells[k]; + SumGjr += Gjr[K][R]; + } + // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr)); + MaxErrorConservationNode[l] = l2Norm(SumGjr); + } + }); + std::cout << " Max Error Node " << max(MaxErrorConservationNode) << "\n"; + } + // + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + // Edge + const auto& cell_to_face = cell_to_face_matrix[j]; + + for (size_t l = 0; l < cell_to_face.size(); ++l) { + const FaceId& face = cell_to_face[l]; + const auto& face_to_cell = face_to_cell_matrix[face]; + const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(face); + + const Rd& Cjf_loc = Cjf(j, l); + + for (size_t k = 0; k < face_to_cell.size(); ++k) { + const CellId K = face_to_cell[k]; + const unsigned int R = face_local_number_in_its_cells[k]; + const Rd& Ckf_loc = Cjf(K, R); + // Une moyenne entre les etats jk + + Rd uFace = .5 * (u_n[j] + u_n[K]); + double cFace = .5 * (c_n[j] + c_n[K]); + + // Viscosity j k + Rpxp ViscosityMatrixJK(identity); + const double MaxmaxabsVpNormjk = + std::max(toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uFace, cFace, + Cjf_loc), + toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uFace, cFace, + Ckf_loc)); + + ViscosityMatrixJK *= MaxmaxabsVpNormjk; + + const Rd& u_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc; // Flux_qtmvt[K] * Cjf_loc; + + const Rp& statediff = StateAtFace[j][l] - StateAtFace[K][R]; // State[j] - State[K]; + const Rp& diff = ViscosityMatrixJK * statediff; + + Gjf[j][l][0] += dot(Flux_rhoAtCellFace[K][R], Cjf_loc); // dot(Flux_rho[K], Cjf_loc); + for (size_t d = 0; d < Dimension; ++d) + Gjf[j][l][1 + d] += u_Cjf[d]; + Gjf[j][l][1 + Dimension] += dot(Flux_totnrjAtCellFace[K][R], Cjf_loc); // dot(Flux_totnrj[K], Cjf_loc); + + Gjf[j][l] += diff; + } + + Gjf[j][l] *= 1. / face_to_cell.size(); + } + }); + synchronize(Gjf); + if (checkLocalConservation) { + auto is_boundary_face = p_mesh->connectivity().isBoundaryFace(); + + FaceValue<double> MaxErrorConservationFace(p_mesh->connectivity()); + MaxErrorConservationFace.fill(0.); + + parallel_for( + p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) { + const auto& face_to_cell = face_to_cell_matrix[l]; + const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(l); + + if (not is_boundary_face[l]) { + Rp SumGjf(zero); + for (size_t k = 0; k < face_to_cell.size(); ++k) { + const CellId K = face_to_cell[k]; + const unsigned int R = face_local_number_in_its_cells[k]; + SumGjf += Gjf[K][R]; + } + MaxErrorConservationFace[l] = l2Norm(SumGjf); + // MaxErrorConservationFace = std::max(MaxErrorConservationFace, l2Norm(SumGjf)); + } + }); + std::cout << " Max Error Face " << max(MaxErrorConservationFace) << "\n"; + } + + if constexpr (Dimension == 3) { + const auto& edge_to_cell_matrix = p_mesh->connectivity().edgeToCellMatrix(); + const auto& edge_local_numbers_in_their_cells = p_mesh->connectivity().edgeLocalNumbersInTheirCells(); + + const EdgeValuePerCell<const Rd> Cje = mesh_data.Cje(); + const EdgeValuePerCell<const Rd> nje = mesh_data.nje(); + + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + // Edge + const auto& cell_to_edge = cell_to_edge_matrix[j]; + + for (size_t l = 0; l < cell_to_edge.size(); ++l) { + const EdgeId& edge = cell_to_edge[l]; + const auto& edge_to_cell = edge_to_cell_matrix[edge]; + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge); + + const Rd& Cje_loc = Cje(j, l); + + for (size_t k = 0; k < edge_to_cell.size(); ++k) { + const CellId K = edge_to_cell[k]; + const size_t R = edge_local_number_in_its_cells[k]; + + const Rd& Cke_loc = Cje(K, R); + + // Une moyenne entre les etats jk + + Rd uEdge = .5 * (u_n[j] + u_n[K]); + double cEdge = .5 * (c_n[j] + c_n[K]); + + // Viscosity j k + Rpxp ViscosityMatrixJK(identity); + const double MaxmaxabsVpNormjk = + std::max(toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uEdge, cEdge, + Cje_loc), + toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uEdge, cEdge, + Cke_loc)); + + ViscosityMatrixJK *= MaxmaxabsVpNormjk; + + const Rd& u_Cje = Flux_qtmvtAtCellEdge[K][R] * Cje_loc; // Flux_qtmvt[K] * Cje_loc; + + const Rp& statediff = StateAtEdge[j][l] - StateAtEdge[K][R]; // State[j] - State[K]; + const Rp& diff = ViscosityMatrixJK * statediff; + + Gje[j][l][0] += dot(Flux_rhoAtCellEdge[K][R], Cje_loc); // dot(Flux_rho[K], Cje_loc); + for (size_t d = 0; d < Dimension; ++d) + Gje[j][l][1 + d] += u_Cje[d]; + Gje[j][l][1 + Dimension] += dot(Flux_totnrjAtCellEdge[K][R], Cje_loc); // dot(Flux_totnrj[K], Cje_loc); + + Gje[j][l] += diff; + } + + Gje[j][l] *= 1. / edge_to_cell.size(); + } + }); + synchronize(Gje); + + if (checkLocalConservation) { + auto is_boundary_edge = p_mesh->connectivity().isBoundaryEdge(); + + EdgeValue<double> MaxErrorConservationEdge(p_mesh->connectivity()); + MaxErrorConservationEdge.fill(0.); + // double MaxErrorConservationEdge = 0; + parallel_for( + p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId l) { + const auto& edge_to_cell = edge_to_cell_matrix[l]; + const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(l); + + if (not is_boundary_edge[l]) { + Rp SumGje(zero); + for (size_t k = 0; k < edge_to_cell.size(); ++k) { + const CellId K = edge_to_cell[k]; + const unsigned int R = edge_local_number_in_its_cells[k]; + SumGje += Gje[K][R]; + } + // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2Norm(SumGje)); + MaxErrorConservationEdge[l] = l2Norm(SumGje); + } + }); + std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << "\n"; + } + } // dim 3 + + // Pour les assemblages + double theta = 2. / 3.; //.5; + double eta = 1. / 6.; //.2; + if constexpr (Dimension == 2) { + eta = 0; + } + // else{ + // theta = 1. / 3.; + // eta = 1. / 3.; + // theta = .5; + // eta = 0; + //} + // + parallel_for( + p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { + const auto& cell_to_node = cell_to_node_matrix[j]; + + Rp SumFluxesNode(zero); + + for (size_t l = 0; l < cell_to_node.size(); ++l) { + SumFluxesNode += Gjr[j][l]; + } + // SumFluxesNode *= (1 - theta); + + const auto& cell_to_edge = cell_to_edge_matrix[j]; + Rp SumFluxesEdge(zero); + + for (size_t l = 0; l < cell_to_edge.size(); ++l) { + SumFluxesEdge += Gje[j][l]; + } + + const auto& cell_to_face = cell_to_face_matrix[j]; + Rp SumFluxesFace(zero); + + for (size_t l = 0; l < cell_to_face.size(); ++l) { + SumFluxesFace += Gjf[j][l]; + } + // SumFluxesEdge *= (theta); + + const Rp SumFluxes = (1 - theta - eta) * SumFluxesNode + eta * SumFluxesEdge + theta * SumFluxesFace; + + State[j] -= dt / volumes_cell[j] * SumFluxes; + + rho[j] = State[j][0]; + for (size_t d = 0; d < Dimension; ++d) + u[j][d] = State[j][1 + d] / rho[j]; + E[j] = State[j][1 + Dimension] / rho[j]; + }); + + return std::make_tuple(std::make_shared<const DiscreteFunctionVariant>(rho), + std::make_shared<const DiscreteFunctionVariant>(u), + std::make_shared<const DiscreteFunctionVariant>(E)); + } + RusanovEulerianCompositeSolver_v2_order_n() : m_quadrature_degree(8) {} + // RusanovEulerianCompositeSolver_v2_order_n() = default; + ~RusanovEulerianCompositeSolver_v2_order_n() = default; +}; + +template <MeshConcept MeshType> +class RusanovEulerianCompositeSolver_v2_order_n<MeshType>::WallBoundaryCondition +{ +}; + +template <> +class RusanovEulerianCompositeSolver_v2_order_n<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 RusanovEulerianCompositeSolver_v2_order_n<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 RusanovEulerianCompositeSolver_v2_order_n<MeshType>::NeumannflatBoundaryCondition +{ +}; +template <> +class RusanovEulerianCompositeSolver_v2_order_n<Mesh<2>>::NeumannflatBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary; + const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary; + + public: + const Rd& + outgoingNormal() const + { + return m_mesh_flat_node_boundary.outgoingNormal(); + } + + size_t + numberOfNodes() const + { + return m_mesh_flat_node_boundary.nodeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_flat_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_flat_node_boundary.nodeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_flat_face_boundary.faceList(); + } + + NeumannflatBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary, + const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary) + : m_mesh_flat_node_boundary(mesh_flat_node_boundary), m_mesh_flat_face_boundary(mesh_flat_face_boundary) + { + ; + } + + ~NeumannflatBoundaryCondition() = default; +}; + +template <> +class RusanovEulerianCompositeSolver_v2_order_n<Mesh<3>>::NeumannflatBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary; + const MeshFlatEdgeBoundary<MeshType> m_mesh_flat_edge_boundary; + const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary; + + public: + const Rd& + outgoingNormal() const + { + return m_mesh_flat_node_boundary.outgoingNormal(); + } + + size_t + numberOfNodes() const + { + return m_mesh_flat_node_boundary.nodeList().size(); + } + + size_t + numberOfEdges() const + { + return m_mesh_flat_edge_boundary.edgeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_flat_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_flat_node_boundary.nodeList(); + } + + const Array<const EdgeId>& + edgeList() const + { + return m_mesh_flat_edge_boundary.edgeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_flat_face_boundary.faceList(); + } + + NeumannflatBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary, + const MeshFlatEdgeBoundary<MeshType>& mesh_flat_edge_boundary, + const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary) + : m_mesh_flat_node_boundary(mesh_flat_node_boundary), + m_mesh_flat_edge_boundary(mesh_flat_edge_boundary), + m_mesh_flat_face_boundary(mesh_flat_face_boundary) + { + ; + } + + ~NeumannflatBoundaryCondition() = default; +}; + +template <MeshConcept MeshType> +class RusanovEulerianCompositeSolver_v2_order_n<MeshType>::SymmetryBoundaryCondition +{ +}; + +template <> +class RusanovEulerianCompositeSolver_v2_order_n<Mesh<2>>::SymmetryBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary; + const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary; + + public: + const Rd& + outgoingNormal() const + { + return m_mesh_flat_node_boundary.outgoingNormal(); + } + + size_t + numberOfNodes() const + { + return m_mesh_flat_node_boundary.nodeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_flat_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_flat_node_boundary.nodeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_flat_face_boundary.faceList(); + } + + SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary, + const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary) + : m_mesh_flat_node_boundary(mesh_flat_node_boundary), m_mesh_flat_face_boundary(mesh_flat_face_boundary) + { + ; + } + + ~SymmetryBoundaryCondition() = default; +}; + +template <> +class RusanovEulerianCompositeSolver_v2_order_n<Mesh<3>>::SymmetryBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary; + const MeshFlatEdgeBoundary<MeshType> m_mesh_flat_edge_boundary; + const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary; + + public: + const Rd& + outgoingNormal() const + { + return m_mesh_flat_node_boundary.outgoingNormal(); + } + + size_t + numberOfNodes() const + { + return m_mesh_flat_node_boundary.nodeList().size(); + } + + size_t + numberOfEdges() const + { + return m_mesh_flat_edge_boundary.edgeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_flat_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_flat_node_boundary.nodeList(); + } + + const Array<const EdgeId>& + edgeList() const + { + return m_mesh_flat_edge_boundary.edgeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_flat_face_boundary.faceList(); + } + + SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary, + const MeshFlatEdgeBoundary<MeshType>& mesh_flat_edge_boundary, + const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary) + : m_mesh_flat_node_boundary(mesh_flat_node_boundary), + m_mesh_flat_edge_boundary(mesh_flat_edge_boundary), + m_mesh_flat_face_boundary(mesh_flat_face_boundary) + { + ; + } + + ~SymmetryBoundaryCondition() = default; +}; + +template <MeshConcept MeshType> +class RusanovEulerianCompositeSolver_v2_order_n<MeshType>::InflowListBoundaryCondition +{ +}; + +template <> +class RusanovEulerianCompositeSolver_v2_order_n<Mesh<2>>::InflowListBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshNodeBoundary m_mesh_node_boundary; + const MeshFaceBoundary m_mesh_face_boundary; + const Table<const double> m_node_array_list; + const Table<const double> m_face_array_list; + + public: + size_t + numberOfNodes() const + { + return m_mesh_node_boundary.nodeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_node_boundary.nodeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_face_boundary.faceList(); + } + + const Table<const double>& + nodeArrayList() const + { + return m_node_array_list; + } + + const Table<const double>& + faceArrayList() const + { + return m_face_array_list; + } + + InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, + const MeshFaceBoundary& mesh_face_boundary, + const Table<const double>& node_array_list, + const Table<const double>& face_array_list) + : m_mesh_node_boundary(mesh_node_boundary), + m_mesh_face_boundary(mesh_face_boundary), + m_node_array_list(node_array_list), + m_face_array_list(face_array_list) + { + ; + } + + ~InflowListBoundaryCondition() = default; +}; + +template <> +class RusanovEulerianCompositeSolver_v2_order_n<Mesh<3>>::InflowListBoundaryCondition +{ + public: + using Rd = TinyVector<Dimension, double>; + + private: + const MeshNodeBoundary m_mesh_node_boundary; + const MeshEdgeBoundary m_mesh_edge_boundary; + const MeshFaceBoundary m_mesh_face_boundary; + const Table<const double> m_node_array_list; + const Table<const double> m_edge_array_list; + const Table<const double> m_face_array_list; + + public: + size_t + numberOfNodes() const + { + return m_mesh_node_boundary.nodeList().size(); + } + + size_t + numberOfEdges() const + { + return m_mesh_edge_boundary.edgeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_node_boundary.nodeList(); + } + + const Array<const EdgeId>& + edgeList() const + { + return m_mesh_edge_boundary.edgeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_face_boundary.faceList(); + } + + const Table<const double>& + nodeArrayList() const + { + return m_node_array_list; + } + + const Table<const double>& + edgeArrayList() const + { + return m_edge_array_list; + } + + const Table<const double>& + faceArrayList() const + { + return m_face_array_list; + } + + InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, + const MeshEdgeBoundary& mesh_edge_boundary, + const MeshFaceBoundary& mesh_face_boundary, + const Table<const double>& node_array_list, + const Table<const double>& edge_array_list, + const Table<const double>& face_array_list) + : m_mesh_node_boundary(mesh_node_boundary), + m_mesh_edge_boundary(mesh_edge_boundary), + m_mesh_face_boundary(mesh_face_boundary), + m_node_array_list(node_array_list), + m_edge_array_list(edge_array_list), + m_face_array_list(face_array_list) + { + ; + } + + ~InflowListBoundaryCondition() = default; +}; + +template <MeshConcept MeshType> +class RusanovEulerianCompositeSolver_v2_order_n<MeshType>::OutflowBoundaryCondition +{ +}; + +template <> +class RusanovEulerianCompositeSolver_v2_order_n<Mesh<2>>::OutflowBoundaryCondition +{ + using Rd = TinyVector<Dimension, double>; + + private: + const MeshNodeBoundary m_mesh_node_boundary; + const MeshFaceBoundary m_mesh_face_boundary; + + public: + size_t + numberOfNodes() const + { + return m_mesh_node_boundary.nodeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_node_boundary.nodeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_face_boundary.faceList(); + } + + OutflowBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) + : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary) + { + ; + } +}; + +template <> +class RusanovEulerianCompositeSolver_v2_order_n<Mesh<3>>::OutflowBoundaryCondition +{ + using Rd = TinyVector<Dimension, double>; + + private: + const MeshNodeBoundary m_mesh_node_boundary; + const MeshEdgeBoundary m_mesh_edge_boundary; + const MeshFaceBoundary m_mesh_face_boundary; + + public: + size_t + numberOfNodes() const + { + return m_mesh_node_boundary.nodeList().size(); + } + size_t + numberOfEdges() const + { + return m_mesh_edge_boundary.edgeList().size(); + } + + size_t + numberOfFaces() const + { + return m_mesh_face_boundary.faceList().size(); + } + + const Array<const NodeId>& + nodeList() const + { + return m_mesh_node_boundary.nodeList(); + } + + const Array<const EdgeId>& + edgeList() const + { + return m_mesh_edge_boundary.edgeList(); + } + + const Array<const FaceId>& + faceList() const + { + return m_mesh_face_boundary.faceList(); + } + + OutflowBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, + const MeshEdgeBoundary& mesh_edge_boundary, + const MeshFaceBoundary& mesh_face_boundary) + : m_mesh_node_boundary(mesh_node_boundary), + + m_mesh_edge_boundary(mesh_edge_boundary), + + m_mesh_face_boundary(mesh_face_boundary) + { + ; + } +}; + +std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>> +rusanovEulerianCompositeSolver_v2_order_n( + const std::shared_ptr<const DiscreteFunctionVariant>& rho_v, + const std::shared_ptr<const DiscreteFunctionVariant>& u_v, + const std::shared_ptr<const DiscreteFunctionVariant>& E_v, + const double& gamma, + const std::shared_ptr<const DiscreteFunctionVariant>& c_v, + const std::shared_ptr<const DiscreteFunctionVariant>& p_v, + const size_t& degree, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, + const double& dt, + const bool check) +{ + std::shared_ptr mesh_v = getCommonMesh({rho_v, u_v, E_v, c_v, p_v}); + if (not mesh_v) { + throw NormalError("discrete functions are not defined on the same mesh"); + } + + if (not checkDiscretizationType({rho_v, u_v, E_v}, DiscreteFunctionType::P0)) { + throw NormalError("acoustic solver expects P0 functions"); + } + + return std::visit( + PUGS_LAMBDA(auto&& p_mesh) + ->std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>> { + using MeshType = mesh_type_t<decltype(p_mesh)>; + static constexpr size_t Dimension = MeshType::Dimension; + using Rd = TinyVector<Dimension>; + + if constexpr (Dimension == 1) { + throw NormalError("RusanovEulerianCompositeSolver v2 order n is not available in 1D"); + } else { + if constexpr (is_polygonal_mesh_v<MeshType>) { + return RusanovEulerianCompositeSolver_v2_order_n<MeshType>{} + .solve(p_mesh, rho_v->get<DiscreteFunctionP0<const double>>(), u_v->get<DiscreteFunctionP0<const Rd>>(), + E_v->get<DiscreteFunctionP0<const double>>(), gamma, c_v->get<DiscreteFunctionP0<const double>>(), + p_v->get<DiscreteFunctionP0<const double>>(), degree, bc_descriptor_list, dt, check); + } else { + throw NormalError("RusanovEulerianCompositeSolver v2 order n is only defined on polygonal meshes"); + } + } + }, + mesh_v->variant()); +} diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2_order_n.hpp b/src/scheme/RusanovEulerianCompositeSolver_v2_order_n.hpp new file mode 100644 index 0000000000000000000000000000000000000000..f5e57539fd3807640e56c34f166f8e677bda9447 --- /dev/null +++ b/src/scheme/RusanovEulerianCompositeSolver_v2_order_n.hpp @@ -0,0 +1,30 @@ +#ifndef RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_ORDER_N_HPP +#define RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_ORDER_N_HPP + +#include <memory> +#include <mesh/MeshVariant.hpp> +#include <scheme/CellbyCellLimitation.hpp> +#include <scheme/DiscreteFunctionVariant.hpp> +#include <scheme/IBoundaryConditionDescriptor.hpp> +#include <scheme/RusanovEulerianCompositeSolverTools.hpp> +#include <vector> + +// double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v, +// const std::shared_ptr<const DiscreteFunctionVariant>& c_v); + +std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>, + std::shared_ptr<const DiscreteFunctionVariant>> +rusanovEulerianCompositeSolver_v2_order_n( + const std::shared_ptr<const DiscreteFunctionVariant>& rho, + const std::shared_ptr<const DiscreteFunctionVariant>& u, + const std::shared_ptr<const DiscreteFunctionVariant>& E, + const double& gamma, + const std::shared_ptr<const DiscreteFunctionVariant>& c, + const std::shared_ptr<const DiscreteFunctionVariant>& p, + const size_t& degree, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, + const double& dt, + const bool check = false); + +#endif // RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_ORDER_N_HPP