#include <scheme/HybridHLLcRusanovEulerianCompositeSolver_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 HybridHLLcRusanovEulerianCompositeSolver_v2 { private: using MeshType = MeshTypeT; static constexpr size_t Dimension = MeshType::Dimension; using Rdxd = TinyMatrix<Dimension>; using Rd = TinyVector<Dimension>; using Rpxp = TinyMatrix<Dimension + 2>; using Rp = TinyVector<Dimension + 2>; class SymmetryBoundaryCondition; class InflowListBoundaryCondition; class OutflowBoundaryCondition; class WallBoundaryCondition; class NeumannflatBoundaryCondition; using BoundaryCondition = std::variant<SymmetryBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition, NeumannflatBoundaryCondition, WallBoundaryCondition>; using BoundaryConditionList = std::vector<BoundaryCondition>; BoundaryConditionList _getBCList(const MeshType& mesh, const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const { BoundaryConditionList bc_list; for (const auto& bc_descriptor : bc_descriptor_list) { bool is_valid_boundary_condition = true; switch (bc_descriptor->type()) { case IBoundaryConditionDescriptor::Type::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 Eulerian Composite solver"; throw NormalError(error_msg.str()); } } return bc_list; } 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& 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 + 1] = 0; } } 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]; // stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = 0; } 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]; // stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = 0; } } } } }, 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] = 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 _applyInflowBoundaryCondition(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValuePerCell<Rp>& stateNode, EdgeValuePerCell<Rp>& stateEdge, FaceValuePerCell<Rp>& stateFace) const { for (const auto& boundary_condition : bc_list) { std::visit( [&](auto&& bc) { using T = std::decay_t<decltype(bc)>; if constexpr (std::is_same_v<InflowListBoundaryCondition, T>) { // MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); std::cout << " Traitement INFLOW \n"; const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix(); const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix(); const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells(); const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells(); // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); const auto& face_list = bc.faceList(); const auto& node_list = bc.nodeList(); const auto& face_array_list = bc.faceArrayList(); const auto& node_array_list = bc.nodeArrayList(); for (size_t i_node = 0; i_node < node_list.size(); ++i_node) { const NodeId node_id = node_list[i_node]; const auto& node_cell_list = node_to_cell_matrix[node_id]; // Assert(face_cell_list.size() == 1); const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id); for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) { CellId node_cell_id = node_cell_list[i_cell]; size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; for (size_t dim = 0; dim < Dimension + 2; ++dim) stateNode[node_cell_id][node_local_number_in_cell][dim] = node_array_list[i_node][dim]; } } for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { const FaceId face_id = face_list[i_face]; const auto& face_cell_list = face_to_cell_matrix[face_id]; Assert(face_cell_list.size() == 1); CellId face_cell_id = face_cell_list[0]; size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0); for (size_t dim = 0; dim < Dimension + 2; ++dim) stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][dim]; } if constexpr (Dimension == 3) { const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix(); const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells(); // const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed(); const auto& edge_list = bc.edgeList(); const auto& edge_array_list = bc.edgeArrayList(); for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) { const EdgeId edge_id = edge_list[i_edge]; const auto& edge_cell_list = edge_to_cell_matrix[edge_id]; const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id); // Assert(face_cell_list.size() == 1); for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) { CellId edge_cell_id = edge_cell_list[i_cell]; size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell]; for (size_t dim = 0; dim < Dimension + 2; ++dim) stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim]; } } } } }, boundary_condition); } } public: inline double pression(const double rho, const double epsilon, const double gam) const { return (gam - 1) * rho * epsilon; } std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>> solve(const std::shared_ptr<const MeshType>& p_mesh, const DiscreteFunctionP0<const double>& rho_n, const DiscreteFunctionP0<const Rd>& u_n, const DiscreteFunctionP0<const double>& E_n, const double& gamma, const DiscreteFunctionP0<const double>& c_n, const DiscreteFunctionP0<const double>& p_n, // const size_t degree, const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, const double& dt, const bool checkLocalConservation) const { auto rho = copy(rho_n); auto u = copy(u_n); auto E = copy(E_n); // auto c = copy(c_n); // auto p = copy(p_n); 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]); }); 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()}; 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); 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 += Gjr[K][R]; maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjr[K][R])); } // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr)); MaxErrorConservationNode[l] = l2Norm(SumGjr) / maxGjrAbs; } }); 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); 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]; } } const std::pair<double, double> MinMaxVpNormj = toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(u_n[j], c_n[j], Cjf_loc); const std::pair<double, double> MinMaxVpNormk = toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(u_n[K], c_n[K], Cjf_loc); const double MinVpNormjk = std::min(MinMaxVpNormj.first, MinMaxVpNormk.first); // SL const double MaxVpNormjk = std::max(MinMaxVpNormj.second, MinMaxVpNormk.second); // SR Rp SC = zero; const Rd& uk_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc; // Flux_qtmvt[K] * Cjf_loc; const Rd& uj_Cjf = Flux_qtmvtAtCellFace[j][l] * Cjf_loc; // Flux_qtmvt[j] * Cjf_loc; const Rp& statediff = StateAtFace[K][R] - StateAtFace[j][l]; // State[K] - State[j]; const Rp& diff = MinVpNormjk * statediff; SC[0] = (dot(Flux_rhoAtCellFace[K][R], Cjf_loc) - dot(Flux_rhoAtCellFace[j][l], Cjf_loc) + diff[1 + Dimension]) / diff[0]; for (size_t d = 0; d < Dimension; ++d) SC[1 + d] = (uk_Cjf[d] - uj_Cjf[d] + diff[d]) / diff[d]; SC[1 + Dimension] = (dot(Flux_totnrjAtCellFace[K][R], Cjf_loc) - dot(Flux_totnrjAtCellFace[j][l], Cjf_loc) + diff[1 + Dimension]) / diff[1 + Dimension]; if (MinVpNormjk >= 0 ){ // SL >=0 Gjf[j][l][0] = dot(Flux_rhoAtCellFace[j][l], Cjf_loc); // dot(Flux_roh[j] , Cjf_loc) for (size_t d = 0; d < Dimension; ++d) Gjf[j][l][1 + d] = uj_Cjf[d]; Gjf[j][l][1 + Dimension] = dot(Flux_totnrjAtCellFace[j][l], Cjf_loc); // dot(Flux_totnrj[K] , Cjf_loc) } else{ if (MaxVpNormjk <= 0 ){ // SR <= 0 Gjf[j][l][0] = dot(Flux_rhoAtCellFace[K][R], Cjf_loc); // dot(Flux_roh[K] , Cjf_loc) for (size_t d = 0; d < Dimension; ++d) Gjf[j][l][1 + d] = uk_Cjf[d]; Gjf[j][l][1 + Dimension] = dot(Flux_totnrjAtCellFace[K][R], Cjf_loc); // dot(Flux_totnrj[K] , Cjf_loc) } else{ Rp diffL = zero; Rp diffR = zero; if (SC[0] >= 0 ){// SL <= 0 <= SC diffL[0] = MinVpNormjk * SC[0] * statediff[0]; Gjf[j][l][0] = (SC[0] * dot(Flux_rhoAtCellFace[j][l], Cjf_loc) - MinVpNormjk * dot(Flux_rhoAtCellFace[K][R], Cjf_loc) + diffL[0]) / (SC[0] - MinVpNormjk); } else{ // SC <= 0 <= SR diffR[0] = SC[0] * MaxVpNormjk * statediff[0]; Gjf[j][l][0] = (MaxVpNormjk * dot(Flux_rhoAtCellFace[j][l], Cjf_loc) - SC[0] * dot(Flux_rhoAtCellFace[K][R], Cjf_loc) + diffR[0]) / (MaxVpNormjk - SC[0]); } for (size_t d = 0; d < Dimension; ++d){ if (SC[1 + d] >= 0 ){// SL <= 0 <= SC diffL[1 + d] = MinVpNormjk * SC[1 + d] * statediff[1 + d]; Gjf[j][l][1 + d] = (SC[1 + d] * uj_Cjf[d] - MinVpNormjk * uk_Cjf[d] + diffL[d]) / (SC[1 + d] - MinVpNormjk); } else{ // SC <= 0 <= SR diffR[1 + d] = SC[1 + d] * MaxVpNormjk * statediff[1 + d]; Gjf[j][l][1 + d] = (MaxVpNormjk * uj_Cjf[d] - SC[1 + d] * uk_Cjf[d] + diffR[d]) / (MaxVpNormjk - SC[1 + d]); } } if (SC[1 + Dimension] >= 0 ){// SL <= 0 <= SC diffL[1 + Dimension] = MinVpNormjk * SC[1 + Dimension] * statediff[1 + Dimension]; Gjf[j][l][1 + Dimension] = (SC[1 + Dimension] * dot(Flux_totnrjAtCellFace[j][l], Cjf_loc) - MinVpNormjk * dot(Flux_totnrjAtCellFace[K][R], Cjf_loc) + diffL[1 + Dimension]) / (SC[1 + Dimension] - MinVpNormjk); } else{ // SC <= 0 <= SR diffR[1 + Dimension] = SC[1 + Dimension] * MaxVpNormjk * statediff[1 + Dimension]; Gjf[j][l][1 + Dimension] = (MaxVpNormjk * dot(Flux_totnrjAtCellFace[j][l], Cjf_loc) - SC[1 + Dimension] * dot(Flux_totnrjAtCellFace[K][R], Cjf_loc) + diffR[1 + Dimension]) / (MaxVpNormjk - SC[1 + Dimension]); } } } } }); 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); 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 += Gjf[K][R]; maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjf[K][R])); } MaxErrorConservationFace[l] = l2Norm(SumGjf) / maxGjrAbs; // 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); 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 += Gje[K][R]; maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gje[K][R])); } // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2Norm(SumGje)); MaxErrorConservationEdge[l] = l2Norm(SumGje) / maxGjrAbs; } }); std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << "\n"; } } // dim 3 // Pour les assemblages double theta = 2. / 3.; //.5; 2. / 3. double eta = 1. / 6.; //.2; 1. / 6. 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)); } HybridHLLcRusanovEulerianCompositeSolver_v2() = default; ~HybridHLLcRusanovEulerianCompositeSolver_v2() = default; }; template <MeshConcept MeshType> class HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>::WallBoundaryCondition { }; template <> class HybridHLLcRusanovEulerianCompositeSolver_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 HybridHLLcRusanovEulerianCompositeSolver_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 HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>::NeumannflatBoundaryCondition { }; template <> class HybridHLLcRusanovEulerianCompositeSolver_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 HybridHLLcRusanovEulerianCompositeSolver_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 HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>::SymmetryBoundaryCondition { }; template <> class HybridHLLcRusanovEulerianCompositeSolver_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 HybridHLLcRusanovEulerianCompositeSolver_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 HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>::InflowListBoundaryCondition { }; template <> class HybridHLLcRusanovEulerianCompositeSolver_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 HybridHLLcRusanovEulerianCompositeSolver_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 HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>::OutflowBoundaryCondition { }; template <> class HybridHLLcRusanovEulerianCompositeSolver_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 HybridHLLcRusanovEulerianCompositeSolver_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>> hybridHLLcRusanovEulerianCompositeSolver_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("Hybrid HLLc/Rusanov EulerianCompositeSolver v2 is not available in 1D"); } else { if constexpr (is_polygonal_mesh_v<MeshType>) { return HybridHLLcRusanovEulerianCompositeSolver_v2<MeshType>{} .solve(p_mesh, rho_v->get<DiscreteFunctionP0<const double>>(), u_v->get<DiscreteFunctionP0<const Rd>>(), E_v->get<DiscreteFunctionP0<const double>>(), gamma, c_v->get<DiscreteFunctionP0<const double>>(), p_v->get<DiscreteFunctionP0<const double>>(), bc_descriptor_list, dt, check); } else { throw NormalError("Hybrid HLLc/Rusanov EulerianCompositeSolver v2 is only defined on polygonal meshes"); } } }, mesh_v->variant()); }