Select Git revision
-
Stéphane Del Pino authored
This function allow to interrupt the script execution. The parameter exit code is returned by pugs when leaving.
Stéphane Del Pino authoredThis function allow to interrupt the script execution. The parameter exit code is returned by pugs when leaving.
HybridHLLcRusanovEulerianCompositeSolver_v2.cpp 80.69 KiB
#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 double rhoj = StateAtFace[j][l][0];
const double rhoK = StateAtFace[K][R][0];
Rd Uj;
Rd UK;
for (size_t dim = 0; dim < Dimension; ++dim){
Uj[dim] = StateAtFace[j][l][dim + 1] / rhoj;
UK[dim] = StateAtFace[K][R][dim + 1] / rhoK;
}
const double uL = Uj[0];
const double uR = UK[0];
const double rhoEj = StateAtFace[j][l][Dimension + 1];
const double rhoEK = StateAtFace[K][R][Dimension + 1];
const double Ej = rhoEj / rhoj;
const double EK = rhoEK / rhoK;
const double epsilonj = Ej - .5 * dot(Uj, Uj);
const double epsilonK = EK - .5 * dot(UK, UK);
const double Pressionj = pression(rhoj, epsilonj, gamma);
const double PressionK = pression(rhoK, epsilonK, gamma);
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) / l2Norm(Cjf_loc); // SL
const double MaxVpNormjk = std::max(MinMaxVpNormj.second, MinMaxVpNormk.second) / l2Norm(Cjf_loc); // SR
const double diffP = PressionK - Pressionj;
const double diffVelocityj = MinVpNormjk - uL;
const double diffVelocityK = MaxVpNormjk - uR;
const double SC = (diffP + rhoj * uL * diffVelocityj - rhoK * uR * diffVelocityK)/(rhoj * diffVelocityj - rhoK * diffVelocityK); // SC
if (MinVpNormjk >= 0 ){ // SL ≥ 0
const Rd& uj_Cjf = Flux_qtmvtAtCellFace[j][l] * Cjf_loc; // Flux_qtmvt[j] * Cjf_loc;
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
const Rd& uk_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc; // Flux_qtmvt[K] * Cjf_loc;
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{
if (SC >= 0){ // SL ≤ 0 ≤ SC
const Rd& uj_Cjf = Flux_qtmvtAtCellFace[j][l] * Cjf_loc; // Flux_qtmvt[j] * Cjf_loc;
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)
Rp DL;
DL[0] = 1;
DL[1] = SC;
for (size_t d = 2; d < (Dimension + 1); ++d){
DL[d] = Uj[d-1];
}
DL[Dimension + 1] = Ej + rhoj * (SC - uL) * (SC + (Pressionj / (rhoj * diffVelocityj)));
const Rp UCL = rhoj * (diffVelocityj / (MinVpNormjk - SC)) * DL;
const Rp diffStates = MinVpNormjk * l2Norm(Cjf_loc) * (UCL - StateAtFace[j][l]);
Gjf[j][l] += diffStates;
}
else{
// SC ≤ 0 ≤ SR
const Rd& uk_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc; // Flux_qtmvt[K] * Cjf_loc;
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)
Rp DR;
DR[0] = 1;
DR[1] = SC;
for (size_t d = 2; d < (Dimension + 1); ++d){
DR[d] = UK[d-1];
}
DR[Dimension + 1] = EK + rhoK * (SC - uR) * (SC + (PressionK / (rhoK * diffVelocityK)));
const Rp UCR = rhoK * (diffVelocityK / (MaxVpNormjk - SC)) * DR;
const Rp diffStates = MaxVpNormjk * l2Norm(Cjf_loc) * (UCR - StateAtFace[K][R]);
Gjf[j][l] += diffStates;
}
}
} }
});
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());
}