Select Git revision
EulerKineticSolverMultiD.cpp
EulerKineticSolverMultiD.cpp 71.61 KiB
#include <scheme/EulerKineticSolverMultiD.hpp>
#include <language/utils/EvaluateAtPoints.hpp>
#include <language/utils/InterpolateItemArray.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshFaceBoundary.hpp>
#include <mesh/MeshFlatFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshVariant.hpp>
#include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionDPkVector.hpp>
#include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
#include <scheme/DiscreteFunctionP0Vector.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/IBoundaryConditionDescriptor.hpp>
#include <scheme/IDiscreteFunctionDescriptor.hpp>
#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <scheme/PolynomialReconstructionDescriptor.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <scheme/WallBoundaryConditionDescriptor.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/LineTransformation.hpp>
#include <tuple>
std::vector<TinyVector<2>>
getLambdaVector2D(const std::shared_ptr<const MeshVariant>& mesh,
const double& lambda,
const size_t& L,
const size_t& M)
{
return std::visit(
[&](auto&& p_mesh) -> std::vector<TinyVector<2>> {
using MeshType = mesh_type_t<decltype(p_mesh)>;
if constexpr (is_polygonal_mesh_v<MeshType> and (MeshType::Dimension == 2)) {
const double pi = std::acos(-1);
std::vector<TinyVector<2>> lambda_vector(4 * L * M); //!!!!!!!!!!!!!!!!!!!!
for (size_t l = 1; l < L + 1; ++l) {
for (size_t m = 1; m < 4 * M + 1; ++m) {
size_t i_wave = (l - 1) * 4 * M + m - 1;
double ll = l;
lambda_vector[i_wave] = TinyVector<2>{(ll / L) * lambda * std::cos((m * pi) / (2 * M)),
(ll / L) * lambda * std::sin((m * pi) / (2 * M))};
}
}
// lambda_vector[4 * L * M] = TinyVector<2>{0, 0}; //!!!!!!!!!!!!!!!!!!!!
return lambda_vector;
} else {
throw NotImplementedError("Invalid mesh type for MultiD Euler Kinetic solver");
}
},
mesh->variant());
}
std::vector<TinyVector<3>>
getLambdaVector3D(const std::shared_ptr<const MeshVariant>& mesh, const double& lambda)
{
return std::visit(
[&](auto&& p_mesh) -> std::vector<TinyVector<3>> {
using MeshType = mesh_type_t<decltype(p_mesh)>;
if constexpr (is_polygonal_mesh_v<MeshType> and (MeshType::Dimension == 3)) {
std::vector<TinyVector<3>> lambda_vector{6};
lambda_vector[0] = TinyVector<3>{lambda, 0., 0.};
lambda_vector[1] = TinyVector<3>{-lambda, 0., 0.};
lambda_vector[2] = TinyVector<3>{0., lambda, 0.};
lambda_vector[3] = TinyVector<3>{0., -lambda, 0.};
lambda_vector[4] = TinyVector<3>{0., 0., lambda};
lambda_vector[5] = TinyVector<3>{0., 0., -lambda};
return lambda_vector;
} else {
throw NotImplementedError("Invalid mesh type for MultiD Euler Kinetic solver");
}
},
mesh->variant());
}
template <MeshConcept MeshType>
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
GetEulerKineticWavesMultiD(const std::shared_ptr<const MeshType> p_mesh,
const std::vector<TinyVector<MeshType::Dimension>> lambda_vector,
DiscreteFunctionP0<const double> rho,
DiscreteFunctionP0<const TinyVector<MeshType::Dimension>> rho_u,
DiscreteFunctionP0<const double> rho_E,
const double& gamma)
{
if (lambda_vector.size() < 3) {
throw NormalError("lambda vector must be of dimension greater than 3");
}
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
const NodeValue<const TinyVector<MeshType::Dimension>> m_xr = p_mesh->xr();
const CellValue<const TinyVector<MeshType::Dimension>> m_xj = mesh_data.xj();
const CellValue<const double> m_Vj = mesh_data.Vj();
DiscreteFunctionP0Vector<double> Fn_rho(p_mesh, lambda_vector.size());
DiscreteFunctionP0Vector<TinyVector<MeshType::Dimension>> Fn_rho_u(p_mesh, lambda_vector.size());
DiscreteFunctionP0Vector<double> Fn_rho_E(p_mesh, lambda_vector.size());
// const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
CellValue<TinyVector<MeshType::Dimension>> A_rho{p_mesh->connectivity()};
CellValue<TinyMatrix<MeshType::Dimension>> A_rho_u{p_mesh->connectivity()};
CellValue<TinyVector<MeshType::Dimension>> A_rho_E{p_mesh->connectivity()};
const TinyMatrix<MeshType::Dimension> I = identity;
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
double p = (gamma - 1) * rho_e;
A_rho[cell_id] = rho_u[cell_id];
A_rho_u[cell_id] = 1. / rho[cell_id] * tensorProduct(rho_u[cell_id], rho_u[cell_id]) + p * I;
A_rho_E[cell_id] = (rho_E[cell_id] + p) / rho[cell_id] * rho_u[cell_id];
});
const size_t nb_waves = lambda_vector.size();
TinyVector<MeshType::Dimension> inv_S = zero;
for (size_t d = 0; d < MeshType::Dimension; ++d) {
for (size_t i = 0; i < nb_waves; ++i) {
inv_S[d] += lambda_vector[i][d] * lambda_vector[i][d];
}
inv_S[d] = 1. / inv_S[d];
}
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
for (size_t i = 0; i < nb_waves; ++i) { Fn_rho[cell_id][i] = (1. / nb_waves) * rho[cell_id];
Fn_rho_u[cell_id][i] = (1. / nb_waves) * rho_u[cell_id];
Fn_rho_E[cell_id][i] = (1. / nb_waves) * rho_E[cell_id];
for (size_t d1 = 0; d1 < MeshType::Dimension; ++d1) {
Fn_rho[cell_id][i] += inv_S[d1] * lambda_vector[i][d1] * A_rho[cell_id][d1];
for (size_t d2 = 0; d2 < MeshType::Dimension; ++d2) {
Fn_rho_u[cell_id][i][d1] += inv_S[d2] * lambda_vector[i][d2] * A_rho_u[cell_id](d2, d1);
}
Fn_rho_E[cell_id][i] += inv_S[d1] * lambda_vector[i][d1] * A_rho_E[cell_id][d1];
}
}
});
return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(Fn_rho),
std::make_shared<DiscreteFunctionVariant>(Fn_rho_u),
std::make_shared<DiscreteFunctionVariant>(Fn_rho_E));
}
template <size_t Dimension>
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
getEulerKineticWavesMultiD(const std::vector<TinyVector<Dimension>>& lambda_vector,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_u,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_E,
const double& gamma)
{
const std::shared_ptr<const MeshVariant> mesh_v = getCommonMesh({rho, rho_u, rho_E});
if (mesh_v.use_count() == 0) {
throw NormalError("rho rho_u and rho_E have to be defined on the same mesh");
}
return std::visit(
[&](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)>;
if constexpr (is_polygonal_mesh_v<MeshType> and (MeshType::Dimension == Dimension)) {
auto rho_h = rho->get<DiscreteFunctionP0<const double>>();
auto rho_u_h = rho_u->get<DiscreteFunctionP0<const TinyVector<MeshType::Dimension>>>();
auto rho_E_h = rho_E->get<DiscreteFunctionP0<const double>>();
return GetEulerKineticWavesMultiD(p_mesh, lambda_vector, rho_h, rho_u_h, rho_E_h, gamma);
} else {
throw NotImplementedError("Invalid mesh type for MultiD Euler Kinetic solver");
}
},
mesh_v->variant());
}
template std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
getEulerKineticWavesMultiD(const std::vector<TinyVector<2>>& lambda_vector,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_u,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_E,
const double& gamma);
template std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
getEulerKineticWavesMultiD(const std::vector<TinyVector<3>>& lambda_vector,
const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_u,
const std::shared_ptr<const DiscreteFunctionVariant>& rho_E,
const double& gamma);
template <MeshConcept MeshType>
class EulerKineticSolverMultiD
{
private:
constexpr static size_t Dimension = MeshType::Dimension;
using Rd = TinyVector<Dimension>;
class SymmetryBoundaryCondition;
class WallBoundaryCondition;
class InflowListBoundaryCondition;
class OutflowBoundaryCondition;
using BoundaryCondition = std::
variant<SymmetryBoundaryCondition, WallBoundaryCondition, InflowListBoundaryCondition, OutflowBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
const double m_dt;
const double m_eps;
const std::vector<TinyVector<Dimension>>& m_lambda_vector;
const size_t m_SpaceOrder;
const CellArray<const double> m_Fn_rho;
const CellArray<const TinyVector<Dimension>> m_Fn_rho_u;
const CellArray<const double> m_Fn_rho_E;
const double m_gamma;
const std::shared_ptr<const MeshType> p_mesh;
const MeshType& m_mesh;
const CellValue<const double> m_Vj = MeshDataManager::instance().getMeshData(m_mesh).Vj();
const FaceValue<const TinyVector<Dimension>> nl = MeshDataManager::instance().getMeshData(m_mesh).nl();
CellValue<const double> m_inv_Vj;
const CellValue<const TinyVector<Dimension>> m_xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
const NodeValue<const TinyVector<Dimension>> m_xr = m_mesh.xr();
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::symmetry: {
bc_list.emplace_back(
SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
break;
}
case IBoundaryConditionDescriptor::Type::wall: {
bc_list.emplace_back(WallBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFaceBoundary(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());
}
auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
CellValue<size_t> is_boundary_cell{p_mesh->connectivity()};
is_boundary_cell.fill(0);
auto node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
auto node_list = node_boundary.nodeList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
auto cell_list = node_to_cell_matrix[node_list[i_node]];
for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
const CellId cell_id = cell_list[i_cell];
is_boundary_cell[cell_id] = 1;
}
}
size_t nb_boundary_cells = sum(is_boundary_cell);
Array<CellId> cell_list{nb_boundary_cells};
size_t index = 0;
for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
if (is_boundary_cell[cell_id]) {
cell_list[index++] = cell_id;
}
}
auto xj = MeshDataManager::instance().getMeshData(m_mesh).xj();
Table<const double> cell_array_list =
InterpolateItemArray<double(Rd)>::template interpolate<ItemType::cell>(inflow_list_bc_descriptor
.functionSymbolIdList(),
xj, cell_list);
bc_list.emplace_back(InflowListBoundaryCondition(cell_array_list, cell_list));
break;
}
case IBoundaryConditionDescriptor::Type::outflow: {
bc_list.emplace_back(OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
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 acoustic solver";
throw NormalError(error_msg.str());
}
}
return bc_list;
}
public:
CellArray<TinyVector<MeshType::Dimension>>
getA(const DiscreteFunctionP0<const double>& rho,
const DiscreteFunctionP0<const double>& rho_u1,
const DiscreteFunctionP0<const double>& rho_u2,
const DiscreteFunctionP0<const double>& rho_E) const
{
// const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
CellArray<TinyVector<Dimension>> vec_A{m_mesh.connectivity(), 4};
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) *
(rho_u1[cell_id] * rho_u1[cell_id] + rho_u2[cell_id] * rho_u2[cell_id]);
double p = (m_gamma - 1) * rho_e;
vec_A[cell_id][0][0] = rho_u1[cell_id];
vec_A[cell_id][1][0] = rho_u1[cell_id] * rho_u1[cell_id] / rho[cell_id] + p;
vec_A[cell_id][2][0] = rho_u1[cell_id] * rho_u2[cell_id] / rho[cell_id]; vec_A[cell_id][3][0] = (rho_E[cell_id] + p) * rho_u1[cell_id] / rho[cell_id];
vec_A[cell_id][0][1] = rho_u2[cell_id];
vec_A[cell_id][1][1] = rho_u1[cell_id] * rho_u2[cell_id] / rho[cell_id];
vec_A[cell_id][2][1] = rho_u2[cell_id] * rho_u2[cell_id] / rho[cell_id] + p;
vec_A[cell_id][3][1] = (rho_E[cell_id] + p) * rho_u2[cell_id] / rho[cell_id];
});
return vec_A;
}
const CellValue<const TinyVector<Dimension>>
getA_rho(const CellValue<const TinyVector<Dimension>>& rho_u) const
{
return rho_u;
}
const CellValue<const TinyMatrix<Dimension>>
getA_rho_u(const CellValue<const double>& rho,
const CellValue<const TinyVector<Dimension>>& rho_u,
const CellValue<const double>& rho_E) const
{
CellValue<TinyMatrix<Dimension>> vec_A{m_mesh.connectivity()};
const TinyMatrix<Dimension> I = identity;
CellValue<double> rho_e{m_mesh.connectivity()};
CellValue<double> p{m_mesh.connectivity()};
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
rho_e[cell_id] = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
p[cell_id] = (m_gamma - 1) * rho_e[cell_id];
vec_A[cell_id] = 1. / rho[cell_id] * tensorProduct(rho_u[cell_id], rho_u[cell_id]) + p[cell_id] * I;
});
return vec_A;
}
const CellValue<const TinyVector<Dimension>>
getA_rho_E(const CellValue<const double>& rho,
const CellValue<const TinyVector<Dimension>>& rho_u,
const CellValue<const double>& rho_E) const
{
CellValue<TinyVector<Dimension>> vec_A{m_mesh.connectivity()};
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
double p = (m_gamma - 1) * rho_e;
vec_A[cell_id] = (rho_E[cell_id] + p) / rho[cell_id] * rho_u[cell_id];
});
return vec_A;
}
const CellArray<const double>
compute_scalar_M(const CellValue<const double>& u, const CellValue<const TinyVector<Dimension>>& A_u)
{
const size_t nb_waves = m_lambda_vector.size();
CellArray<double> M{p_mesh->connectivity(), nb_waves};
TinyVector<Dimension> inv_S = zero;
for (size_t d = 0; d < MeshType::Dimension; ++d) {
for (size_t i = 0; i < nb_waves; ++i) {
inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
}
inv_S[d] = 1. / inv_S[d];
}
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { for (size_t i = 0; i < nb_waves; ++i) {
M[cell_id][i] = (1. / nb_waves) * u[cell_id];
for (size_t d = 0; d < Dimension; ++d) {
M[cell_id][i] += inv_S[d] * m_lambda_vector[i][d] * A_u[cell_id][d];
}
}
});
return M;
}
const CellArray<const TinyVector<Dimension>>
compute_vector_M(const CellValue<const TinyVector<Dimension>>& u, const CellValue<const TinyMatrix<Dimension>>& A_u)
{
const size_t nb_waves = m_lambda_vector.size();
CellArray<TinyVector<Dimension>> M{p_mesh->connectivity(), nb_waves};
TinyVector<MeshType::Dimension> inv_S = zero;
for (size_t d = 0; d < MeshType::Dimension; ++d) {
for (size_t i = 0; i < nb_waves; ++i) {
inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
}
inv_S[d] = 1. / inv_S[d];
}
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
M[cell_id][i_wave] = 1. / nb_waves * u[cell_id];
for (size_t d1 = 0; d1 < Dimension; ++d1) {
for (size_t d2 = 0; d2 < Dimension; ++d2) {
M[cell_id][i_wave][d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_u[cell_id](d2, d1);
}
}
}
});
return M;
}
template <typename T>
const CellValue<const T>
compute_PFnp1(const CellArray<const T>& Fn, const CellArray<const T>& deltaLFn)
{
CellValue<T> PFnp1{p_mesh->connectivity()};
if constexpr (is_tiny_vector_v<T>) {
PFnp1.fill(zero);
} else {
PFnp1.fill(0.);
}
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
const double inv_Vj = m_inv_Vj[cell_id];
const double dt_over_Vj = m_dt * inv_Vj;
for (size_t i_wave = 0; i_wave < m_lambda_vector.size(); ++i_wave) {
PFnp1[cell_id] += Fn[cell_id][i_wave] - dt_over_Vj * deltaLFn[cell_id][i_wave];
}
});
return PFnp1;
}
template <typename T>
const CellValue<const T>
compute_PFn(const CellArray<const T>& F)
{
const size_t nb_waves = m_lambda_vector.size(); CellValue<T> PFn{p_mesh->connectivity()};
if constexpr (is_tiny_vector_v<T>) {
PFn.fill(zero);
} else {
PFn.fill(0.);
}
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
for (size_t i = 0; i < nb_waves; ++i) {
PFn[cell_id] += F[cell_id][i];
}
});
return PFn;
}
template <typename T>
NodeArray<T>
compute_Flux1_glace(const CellArray<const T>& Fn, NodeValue<bool> is_symmetry_boundary_node)
{
if constexpr (Dimension > 1) {
const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
const size_t nb_waves = m_lambda_vector.size();
NodeArray<T> Fr(m_mesh.connectivity(), nb_waves);
const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
if constexpr (is_tiny_vector_v<T>) {
Fr.fill(zero);
} else {
Fr.fill(0.);
}
parallel_for(
p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
const auto& node_to_cell = node_to_cell_matrix[node_id];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
double sum_li_njr = 0;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
const unsigned int i_node = node_local_number_in_its_cell[i_cell];
double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node));
if (li_njr > 0) {
Fr[node_id][i_wave] += li_njr * Fn[cell_id][i_wave];
sum_li_njr += li_njr;
}
}
if ((sum_li_njr != 0) and (not is_symmetry_boundary_node[node_id])) {
Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave];
}
}
});
return Fr;
} else {
throw NormalError("Glace not defined in 1d");
}
}
template <typename T>
FaceArrayPerNode<T>
compute_Flux1_eucclhyd(const CellArray<const T>& Fn, NodeValue<bool> is_symmetry_boundary_node)
{
if constexpr (Dimension > 1) {
const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
const size_t nb_waves = m_lambda_vector.size();
FaceArrayPerNode<T> Flr(m_mesh.connectivity(), nb_waves); const auto& node_local_numbers_in_their_faces = p_mesh->connectivity().nodeLocalNumbersInTheirFaces();
const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
const auto& node_to_face_matrix = p_mesh->connectivity().nodeToFaceMatrix();
const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
if constexpr (is_tiny_vector_v<T>) {
Flr.fill(zero);
} else {
Flr.fill(0.);
}
parallel_for(
p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_local_number_in_its_face = node_local_numbers_in_their_faces.itemArray(node_id);
const auto& node_to_face = node_to_face_matrix[node_id];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
const FaceId face_id = node_to_face[i_face];
const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
const auto& face_to_cell = face_to_cell_matrix[face_id];
const size_t i_node_face = node_local_number_in_its_face[i_face];
double sum = 0;
for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
const CellId cell_id = face_to_cell[i_cell];
const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
TinyVector<Dimension> n_face = nlr(face_id, i_node_face);
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
if (li_nlr > 0) {
Flr[node_id][i_face][i_wave] += li_nlr * Fn[cell_id][i_wave];
sum += li_nlr;
}
}
if ((sum != 0) and (not is_symmetry_boundary_node[node_id])) {
Flr[node_id][i_face][i_wave] = 1. / sum * Flr[node_id][i_face][i_wave];
}
}
}
});
return Flr;
} else {
throw NormalError("Eucclhyd not defined in 1d");
}
}
template <typename T>
FaceArray<T>
compute_Flux1_face(const CellArray<const T> Fn)
{
const size_t nb_waves = m_lambda_vector.size();
FaceArray<T> Fl(m_mesh.connectivity(), nb_waves);
const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
if constexpr (is_tiny_vector_v<T>) {
Fl.fill(zero);
} else {
Fl.fill(0.);
}
parallel_for( p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
const auto& face_to_cell = face_to_cell_matrix[face_id];
for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
const CellId cell_id = face_to_cell[i_cell];
const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
TinyVector<Dimension> n_face = nl[face_id];
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
if (li_nl > 0) {
Fl[face_id][i_wave] += Fn[cell_id][i_wave];
}
}
}
});
return Fl;
}
void
apply_face_bc(FaceArray<double> Fl_rho,
FaceArray<TinyVector<Dimension>> Fl_rho_u,
FaceArray<double> Fl_rho_E,
const CellValue<const double> rho,
const CellValue<const TinyVector<Dimension>> rho_u,
const CellValue<const double> rho_E,
const CellArray<const double> Fn_rho,
const CellArray<const TinyVector<Dimension>> Fn_rho_u,
const CellArray<const double> Fn_rho_E,
const BoundaryConditionList& bc_list)
{
const size_t nb_waves = m_lambda_vector.size();
const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
TinyVector<MeshType::Dimension> inv_S = zero;
for (size_t d = 0; d < MeshType::Dimension; ++d) {
for (size_t i = 0; i < nb_waves; ++i) {
inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
}
inv_S[d] = 1. / inv_S[d];
}
for (auto&& bc_v : bc_list) {
std::visit(
[=, this](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
auto face_list = bc.faceList();
auto n = bc.outgoingNormal();
auto nxn = tensorProduct(n, n);
TinyMatrix<Dimension> I = identity;
auto txt = I - nxn;
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
const auto& face_to_cell = face_to_cell_matrix[face_id];
for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
const CellId cell_id = face_to_cell[i_cell];
const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
TinyVector<Dimension> n_face = nl[face_id];
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
if (li_nl < 0) {
double rhoj_prime = rho[cell_id];
TinyVector<Dimension> rho_uj = rho_u[cell_id];
TinyVector<Dimension> rho_uj_prime = txt * rho_uj - nxn * rho_uj;
double rho_Ej_prime = rho_E[cell_id];
double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
double p = (m_gamma - 1) * rho_e;
TinyVector<Dimension> A_rho_prime = rho_uj_prime;
TinyMatrix<Dimension> A_rho_u_prime =
1. / rhoj_prime * tensorProduct(rho_uj_prime, rho_uj_prime) + p * I;
TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime;
double Fn_rho_prime = rhoj_prime / nb_waves;
TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_uj_prime;
double Fn_rho_E_prime = rho_Ej_prime / nb_waves;
for (size_t d1 = 0; d1 < Dimension; ++d1) {
Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1];
for (size_t d2 = 0; d2 < Dimension; ++d2) {
Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1);
}
Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1];
}
Fl_rho[face_id][i_wave] += Fn_rho_prime;
Fl_rho_u[face_id][i_wave] += Fn_rho_u_prime;
Fl_rho_E[face_id][i_wave] += Fn_rho_E_prime;
}
}
}
}
} else if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
auto face_list = bc.faceList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
const auto& face_to_cell = face_to_cell_matrix[face_id];
for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
const CellId cell_id = face_to_cell[i_cell];
const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
auto n = sign * nl[face_id];
auto nxn = tensorProduct(n, n);
TinyMatrix<Dimension> I = identity;
auto txt = I - nxn;
double li_nl = dot(m_lambda_vector[i_wave], n);
if (li_nl < 0) {
double rhoj_prime = rho[cell_id];
TinyVector<Dimension> rho_uj = rho_u[cell_id];
TinyVector<Dimension> rho_uj_prime = txt * rho_uj - nxn * rho_uj;
double rho_Ej_prime = rho_E[cell_id];
double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
double p = (m_gamma - 1) * rho_e;
TinyVector<Dimension> A_rho_prime = rho_uj_prime;
TinyMatrix<Dimension> A_rho_u_prime = 1. / rhoj_prime * tensorProduct(rho_uj_prime, rho_uj_prime) + p * I;
TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime;
double Fn_rho_prime = rhoj_prime / nb_waves;
TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_uj_prime;
double Fn_rho_E_prime = rho_Ej_prime / nb_waves;
for (size_t d1 = 0; d1 < Dimension; ++d1) {
Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1];
for (size_t d2 = 0; d2 < Dimension; ++d2) {
Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1);
}
Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1];
}
Fl_rho[face_id][i_wave] += Fn_rho_prime;
Fl_rho_u[face_id][i_wave] += Fn_rho_u_prime;
Fl_rho_E[face_id][i_wave] += Fn_rho_E_prime;
}
}
}
}
} else if constexpr (std::is_same_v<BCType, OutflowBoundaryCondition>) {
auto face_list = bc.faceList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
const auto& face_to_cell = face_to_cell_matrix[face_id];
for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
const CellId cell_id = face_to_cell[i_cell];
const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
auto n = sign * nl[face_id];
double li_nl = dot(m_lambda_vector[i_wave], n);
if (li_nl < 0) {
Fl_rho[face_id][i_wave] += Fn_rho[cell_id][i_wave];
Fl_rho_u[face_id][i_wave] += Fn_rho_u[cell_id][i_wave];
Fl_rho_E[face_id][i_wave] += Fn_rho_E[cell_id][i_wave];
}
}
}
}
}
},
bc_v);
}
}
void
apply_glace_bc(NodeArray<double> Fr_rho,
NodeArray<TinyVector<Dimension>> Fr_rho_u,
NodeArray<double> Fr_rho_E,
const CellValue<const double> rho,
const CellValue<const TinyVector<Dimension>> rho_u,
const CellValue<const double> rho_E,
const BoundaryConditionList& bc_list)
{
const size_t nb_waves = m_lambda_vector.size();
const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
TinyVector<MeshType::Dimension> inv_S = zero;
for (size_t d = 0; d < MeshType::Dimension; ++d) { for (size_t i = 0; i < nb_waves; ++i) {
inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
}
inv_S[d] = 1. / inv_S[d];
}
NodeValue<bool> node_is_done{p_mesh->connectivity()};
node_is_done.fill(false);
for (auto&& bc_v : bc_list) {
std::visit(
[=, this](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
auto node_list = bc.nodeList();
auto n = bc.outgoingNormal();
auto nxn = tensorProduct(n, n);
TinyMatrix<Dimension> I = identity;
auto txt = I - nxn;
// for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
// const CellId cell_id = node_to_cell[i_cell];
// const unsigned int i_node = node_local_number_in_its_cell[i_cell];
// double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node));
// if (li_njr > 0) {
// Fr[node_id][i_wave] += Fn[cell_id][i_wave] * li_njr;
// sum_li_njr += li_njr;
// }
// }
// if (sum_li_njr != 0) {
// Fr[node_id][i_wave] /= sum_li_njr;
// }
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
if (not node_is_done[node_id]) {
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
const auto& node_to_cell = node_to_cell_matrix[node_id];
double sum_li_njr = 0;
for (size_t i_cell = 0; i_cell < node_to_cell.size(); ++i_cell) {
const CellId cell_id = node_to_cell[i_cell];
const size_t i_node_cell = node_local_number_in_its_cell[i_cell];
double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell));
if (li_njr > 0) {
sum_li_njr += li_njr;
}
TinyVector<Dimension> njr_prime = txt * njr(cell_id, i_node_cell) - nxn * njr(cell_id, i_node_cell);
double li_njr_prime = dot(m_lambda_vector[i_wave], njr_prime);
if (li_njr_prime > 0) {
double rhoj_prime = rho[cell_id];
TinyVector<Dimension> rho_uj_prime = txt * rho_u[cell_id] - nxn * rho_u[cell_id];
double rho_Ej_prime = rho_E[cell_id];
double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
double p = (m_gamma - 1) * rho_e;
TinyVector<Dimension> A_rho_prime = rho_uj_prime;
TinyMatrix<Dimension> A_rho_u_prime =
1. / rhoj_prime * tensorProduct(rho_uj_prime, rho_uj_prime) + p * I;
TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime;
double Fn_rho_prime = rhoj_prime / nb_waves; TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_uj_prime;
double Fn_rho_E_prime = rho_Ej_prime / nb_waves;
for (size_t d1 = 0; d1 < Dimension; ++d1) {
Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1];
for (size_t d2 = 0; d2 < Dimension; ++d2) {
Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1);
}
Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1];
}
Fr_rho[node_id][i_wave] += Fn_rho_prime * li_njr_prime;
Fr_rho_u[node_id][i_wave] += 1. / li_njr_prime * Fn_rho_u_prime;
Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_njr_prime;
sum_li_njr += li_njr_prime;
}
}
if (sum_li_njr != 0) {
Fr_rho[node_id][i_wave] /= sum_li_njr;
Fr_rho_u[node_id][i_wave] = 1. / sum_li_njr * Fr_rho_u[node_id][i_wave];
Fr_rho_E[node_id][i_wave] /= sum_li_njr;
}
}
node_is_done[node_id] = true;
}
}
}
},
bc_v);
}
}
void
apply_eucclhyd_bc(FaceArrayPerNode<double> Flr_rho,
FaceArrayPerNode<TinyVector<Dimension>> Flr_rho_u,
FaceArrayPerNode<double> Flr_rho_E,
const CellValue<const double> rho,
const CellValue<const TinyVector<Dimension>> rho_u,
const CellValue<const double> rho_E,
const BoundaryConditionList& bc_list)
{
const size_t nb_waves = m_lambda_vector.size();
// const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
const auto& face_local_numbers_in_their_nodes = p_mesh->connectivity().faceLocalNumbersInTheirNodes();
const auto& face_to_node_matrix = p_mesh->connectivity().faceToNodeMatrix();
const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
TinyVector<MeshType::Dimension> inv_S = zero;
for (size_t d = 0; d < MeshType::Dimension; ++d) {
for (size_t i = 0; i < nb_waves; ++i) {
inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
}
inv_S[d] = 1. / inv_S[d];
}
FaceArrayPerNode<double> sum_li_nlr(m_mesh.connectivity(), nb_waves);
sum_li_nlr.fill(0);
for (auto&& bc_v : bc_list) {
std::visit(
[=, this](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
auto face_list = bc.faceList();
auto n = bc.outgoingNormal(); auto nxn = tensorProduct(n, n);
TinyMatrix<Dimension> I = identity;
auto txt = I - nxn;
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_to_node = face_to_node_matrix[face_id];
const auto& face_local_number_in_its_cell = face_local_numbers_in_their_cells.itemArray(face_id);
const auto& face_to_cell = face_to_cell_matrix[face_id];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
for (size_t i_node = 0; i_node < face_to_node.size(); ++i_node) {
const NodeId node_id = face_to_node[i_node];
const auto& face_local_number_in_its_node = face_local_numbers_in_their_nodes.itemArray(face_id);
const size_t i_face_node = face_local_number_in_its_node[i_node];
for (size_t i_cell = 0; i_cell < face_to_cell.size(); ++i_cell) {
const CellId cell_id = face_to_cell[i_cell];
const size_t i_face_cell = face_local_number_in_its_cell[i_cell];
TinyVector<Dimension> n_face = nlr(face_id, i_node);
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
double li_nlr_prime = sign * dot(m_lambda_vector[i_wave], nlr_prime);
double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
if (li_nlr > 0) {
sum_li_nlr[node_id][i_face_node][i_wave] += li_nlr;
}
if (li_nlr_prime > 0) {
double rhoj_prime = rho[cell_id];
TinyVector<Dimension> rho_uj_prime = txt * rho_u[cell_id] - nxn * rho_u[cell_id];
double rho_Ej_prime = rho_E[cell_id];
double rho_e = rho_E[cell_id] - 0.5 * (1. / rho[cell_id]) * dot(rho_u[cell_id], rho_u[cell_id]);
double p = (m_gamma - 1) * rho_e;
TinyVector<Dimension> A_rho_prime = rho_uj_prime;
TinyMatrix<Dimension> A_rho_u_prime =
1. / rhoj_prime * tensorProduct(rho_uj_prime, rho_uj_prime) + p * I;
TinyVector<Dimension> A_rho_E_prime = (rho_Ej_prime + p) / rhoj_prime * rho_uj_prime;
double Fn_rho_prime = rhoj_prime / nb_waves;
TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_uj_prime;
double Fn_rho_E_prime = rho_Ej_prime / nb_waves;
for (size_t d1 = 0; d1 < Dimension; ++d1) {
Fn_rho_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_prime[d1];
for (size_t d2 = 0; d2 < Dimension; ++d2) {
Fn_rho_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_rho_u_prime(d2, d1);
}
Fn_rho_E_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_rho_E_prime[d1];
}
Flr_rho[node_id][i_face_node][i_wave] += Fn_rho_prime * li_nlr_prime;
Flr_rho_u[node_id][i_face_node][i_wave] += li_nlr_prime * Fn_rho_u_prime;
Flr_rho_E[node_id][i_face_node][i_wave] += Fn_rho_E_prime * li_nlr_prime;
sum_li_nlr[node_id][i_face_node][i_wave] += li_nlr_prime;
}
}
}
}
}
}
},
bc_v); }
for (auto&& bc_v : bc_list) {
std::visit(
[=](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
auto face_list = bc.faceList();
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_to_node = face_to_node_matrix[face_id];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
for (size_t i_node = 0; i_node < face_to_node.size(); ++i_node) {
const NodeId node_id = face_to_node[i_node];
const auto& face_local_number_in_its_node = face_local_numbers_in_their_nodes.itemArray(face_id);
const size_t i_face_node = face_local_number_in_its_node[i_node];
if (sum_li_nlr[node_id][i_face_node][i_wave] != 0) {
Flr_rho[node_id][i_face_node][i_wave] =
1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho[node_id][i_face_node][i_wave];
Flr_rho_u[node_id][i_face_node][i_wave] =
1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho_u[node_id][i_face_node][i_wave];
Flr_rho_E[node_id][i_face_node][i_wave] =
1. / sum_li_nlr[node_id][i_face_node][i_wave] * Flr_rho[node_id][i_face_node][i_wave];
}
}
}
}
}
},
bc_v);
}
}
template <typename T>
CellArray<T>
compute_deltaLFn_glace(NodeArray<T> Fr)
{
if constexpr (Dimension > 1) {
const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
const size_t nb_waves = m_lambda_vector.size();
CellArray<T> deltaLFn(p_mesh->connectivity(), nb_waves);
const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
if constexpr (is_tiny_vector_v<T>) {
deltaLFn.fill(zero);
} else {
deltaLFn.fill(0.);
}
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
const auto& cell_to_node = cell_to_node_matrix[cell_id];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
const NodeId node_id = cell_to_node[i_node];
double li_Cjr = dot(m_lambda_vector[i_wave], Cjr(cell_id, i_node));
deltaLFn[cell_id][i_wave] += li_Cjr * Fr[node_id][i_wave];
}
}
});
return deltaLFn;
} else { throw NormalError("Glace not defined in 1d");
}
}
template <typename T>
CellArray<T>
compute_deltaLFn_eucclhyd(FaceArrayPerNode<T> Flr)
{
if constexpr (Dimension > 1) {
const size_t nb_waves = m_lambda_vector.size();
CellArray<T> deltaLFn(p_mesh->connectivity(), nb_waves);
if constexpr (is_tiny_vector_v<T>) {
deltaLFn.fill(zero);
} else {
deltaLFn.fill(0.);
}
const NodeValuePerFace<const TinyVector<Dimension>> Nlr = MeshDataManager::instance().getMeshData(m_mesh).Nlr();
const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
const auto& face_to_node_matrix = p_mesh->connectivity().faceToNodeMatrix();
const auto& node_to_face_matrix = p_mesh->connectivity().nodeToFaceMatrix();
const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
const auto& cell_to_face = cell_to_face_matrix[cell_id];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
for (size_t i_face_cell = 0; i_face_cell < cell_to_face.size(); ++i_face_cell) {
const FaceId face_id = cell_to_face[i_face_cell];
const auto& face_to_node = face_to_node_matrix[face_id];
for (size_t i_node_face = 0; i_node_face < face_to_node.size(); ++i_node_face) {
const NodeId node_id = face_to_node[i_node_face];
const auto& node_to_face = node_to_face_matrix[node_id];
auto local_face_number_in_node = [&](FaceId face_number) {
for (size_t i_face = 0; i_face < node_to_face.size(); ++i_face) {
if (face_number == node_to_face[i_face]) {
return i_face;
}
}
return std::numeric_limits<size_t>::max();
};
const size_t i_face_node = local_face_number_in_node(face_id);
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
const double li_Nlr = sign * dot(m_lambda_vector[i_wave], Nlr(face_id, i_node_face));
deltaLFn[cell_id][i_wave] += li_Nlr * Flr[node_id][i_face_node][i_wave];
}
}
}
});
return deltaLFn;
} else {
throw NormalError("Eucclhyd not defined in 1d");
}
}
template <typename T>
CellArray<T>
compute_deltaLFn_face(FaceArray<T> Fl)
{
if constexpr (Dimension > 1) {
const FaceValue<const TinyVector<Dimension>> Nl = MeshDataManager::instance().getMeshData(m_mesh).Nl(); const size_t nb_waves = m_lambda_vector.size();
CellArray<T> deltaLFn(p_mesh->connectivity(), nb_waves);
const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
if constexpr (is_tiny_vector_v<T>) {
deltaLFn.fill(zero);
} else {
deltaLFn.fill(0.);
}
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
const auto& cell_to_face = cell_to_face_matrix[cell_id];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
for (size_t i_face_cell = 0; i_face_cell < cell_to_face.size(); ++i_face_cell) {
const FaceId face_id = cell_to_face[i_face_cell];
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
const double li_Nl = sign * dot(m_lambda_vector[i_wave], Nl[face_id]);
deltaLFn[cell_id][i_wave] += li_Nl * Fl[face_id][i_wave];
}
}
});
return deltaLFn;
} else {
throw NormalError("Nl in meaningless in 1d");
}
}
CellValue<bool>
getInflowBoundaryCells(const BoundaryConditionList& bc_list)
{
CellValue<bool> is_boundary_cell{p_mesh->connectivity()};
is_boundary_cell.fill(false);
for (auto&& bc_v : bc_list) {
std::visit(
[=](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<BCType, InflowListBoundaryCondition>) {
auto cell_list = bc.cellList();
for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
const CellId cell_id = cell_list[i_cell];
is_boundary_cell[cell_id] = true;
}
}
},
bc_v);
}
return is_boundary_cell;
}
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
apply(const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
BoundaryConditionList bc_list = this->_getBCList(m_mesh, bc_descriptor_list);
const size_t nb_waves = m_lambda_vector.size();
CellValue<bool> is_inflow_boundary_cell = getInflowBoundaryCells(bc_list);
const CellValue<double> rho_ex{p_mesh->connectivity()};
const CellValue<TinyVector<Dimension>> rho_u_ex{p_mesh->connectivity()};
const CellValue<double> rho_E_ex{p_mesh->connectivity()};
rho_ex.fill(1); // !! A modifier en ne parcourant que les bords
rho_u_ex.fill(zero);
rho_E_ex.fill(0);
for (auto&& bc_v : bc_list) {
std::visit(
[=](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<BCType, InflowListBoundaryCondition>) {
auto cell_list = bc.cellList();
auto cell_array_list = bc.cellArrayList();
for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
const CellId cell_id = cell_list[i_cell];
rho_ex[cell_id] = cell_array_list[i_cell][0];
rho_u_ex[cell_id][0] = cell_array_list[i_cell][1];
rho_u_ex[cell_id][1] = cell_array_list[i_cell][2];
rho_E_ex[cell_id] = cell_array_list[i_cell][3];
}
}
},
bc_v);
}
NodeValue<bool> is_symmetry_boundary_node{p_mesh->connectivity()};
is_symmetry_boundary_node.fill(false);
for (auto&& bc_v : bc_list) {
std::visit(
[=](auto&& bc) {
using BCType = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
auto node_list = bc.nodeList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
is_symmetry_boundary_node[node_list[i_node]] = 1;
}
}
},
bc_v);
}
const CellValue<const TinyVector<Dimension>> A_rho_ex = getA_rho(rho_u_ex);
const CellValue<const TinyMatrix<Dimension>> A_rho_u_ex = getA_rho_u(rho_ex, rho_u_ex, rho_E_ex);
const CellValue<const TinyVector<Dimension>> A_rho_E_ex = getA_rho_E(rho_ex, rho_u_ex, rho_E_ex);
const CellArray<const double> Fn_rho_ex = compute_scalar_M(rho_ex, A_rho_ex);
const CellArray<const TinyVector<Dimension>> Fn_rho_u_ex = compute_vector_M(rho_u_ex, A_rho_u_ex);
const CellArray<const double> Fn_rho_E_ex = compute_scalar_M(rho_E_ex, A_rho_E_ex);
const CellArray<const double> Fn_rho = copy(m_Fn_rho);
const CellArray<const TinyVector<Dimension>> Fn_rho_u = copy(m_Fn_rho_u);
const CellArray<const double> Fn_rho_E = copy(m_Fn_rho_E);
DiscreteFunctionP0Vector<double> Fnp1_rho(p_mesh, nb_waves);
DiscreteFunctionP0Vector<TinyVector<Dimension>> Fnp1_rho_u(p_mesh, nb_waves);
DiscreteFunctionP0Vector<double> Fnp1_rho_E(p_mesh, nb_waves);
// Compute first order F
const CellValue<const double> rho = compute_PFn<double>(Fn_rho);
const CellValue<const TinyVector<Dimension>> rho_u = compute_PFn<TinyVector<Dimension>>(Fn_rho_u);
const CellValue<const double> rho_E = compute_PFn<double>(Fn_rho_E);
// NodeArray<double> Fr_rho = compute_Flux1_glace<double>(Fn_rho, is_symmetry_boundary_node);
// NodeArray<TinyVector<Dimension>> Fr_rho_u =
// compute_Flux1_glace<TinyVector<Dimension>>(Fn_rho_u, is_symmetry_boundary_node);
// NodeArray<double> Fr_rho_E = compute_Flux1_glace<double>(Fn_rho_E, is_symmetry_boundary_node);
FaceArrayPerNode<double> Flr_rho = compute_Flux1_eucclhyd<double>(Fn_rho, is_symmetry_boundary_node);
FaceArrayPerNode<TinyVector<Dimension>> Flr_rho_u =
compute_Flux1_eucclhyd<TinyVector<Dimension>>(Fn_rho_u, is_symmetry_boundary_node);
FaceArrayPerNode<double> Flr_rho_E = compute_Flux1_eucclhyd<double>(Fn_rho_E, is_symmetry_boundary_node);
// FaceArray<double> Fl_rho = compute_Flux1_face<double>(Fn_rho);
// FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux1_face<TinyVector<Dimension>>(Fn_rho_u);
// FaceArray<double> Fl_rho_E = compute_Flux1_face<double>(Fn_rho_E);
// apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
// apply_glace_bc(Fr_rho, Fr_rho_u, Fr_rho_E, rho, rho_u, rho_E, Fn_rho, Fn_rho_u, Fn_rho_E, bc_list);
apply_eucclhyd_bc(Flr_rho, Flr_rho_u, Flr_rho_E, rho, rho_u, rho_E, bc_list);
// synchronize(Flr_rho);
// synchronize(Flr_rho_u);
// synchronize(Flr_rho_E);
// const CellArray<const double> deltaLFn_rho = compute_deltaLFn_glace(Fr_rho);
// const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = compute_deltaLFn_glace(Fr_rho_u);
// const CellArray<const double> deltaLFn_rho_E = compute_deltaLFn_glace(Fr_rho_E);
const CellArray<const double> deltaLFn_rho = compute_deltaLFn_eucclhyd(Flr_rho);
const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u = compute_deltaLFn_eucclhyd(Flr_rho_u);
const CellArray<const double> deltaLFn_rho_E = compute_deltaLFn_eucclhyd(Flr_rho_E);
// const CellArray<const double> deltaLFn_rho = compute_deltaLFn_face<double>(Fl_rho);
// const CellArray<const TinyVector<Dimension>> deltaLFn_rho_u =
// compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u);
// const CellArray<const double> deltaLFn_rho_E = compute_deltaLFn_face<double>(Fl_rho_E);
const CellValue<const TinyVector<Dimension>> A_rho = getA_rho(rho_u);
const CellValue<const TinyMatrix<Dimension>> A_rho_u = getA_rho_u(rho, rho_u, rho_E);
const CellValue<const TinyVector<Dimension>> A_rho_E = getA_rho_E(rho, rho_u, rho_E);
const CellArray<const double> M_rho = compute_scalar_M(rho, A_rho);
const CellArray<const TinyVector<Dimension>> M_rho_u = compute_vector_M(rho_u, A_rho_u);
const CellArray<const double> M_rho_E = compute_scalar_M(rho_E, A_rho_E);
const CellValue<const double> rho_np1 = compute_PFnp1<double>(Fn_rho, deltaLFn_rho);
const CellValue<const TinyVector<Dimension>> rho_u_np1 =
compute_PFnp1<TinyVector<Dimension>>(Fn_rho_u, deltaLFn_rho_u);
const CellValue<const double> rho_E_np1 = compute_PFnp1<double>(Fn_rho_E, deltaLFn_rho_E);
const CellValue<const TinyVector<Dimension>> A_rho_np1 = getA_rho(rho_u_np1);
const CellValue<const TinyMatrix<Dimension>> A_rho_u_np1 = getA_rho_u(rho_np1, rho_u_np1, rho_E_np1);
const CellValue<const TinyVector<Dimension>> A_rho_E_np1 = getA_rho_E(rho_np1, rho_u_np1, rho_E_np1);
const CellArray<const double> M_rho_np1 = compute_scalar_M(rho_np1, A_rho_np1);
const CellArray<const TinyVector<Dimension>> M_rho_u_np1 = compute_vector_M(rho_u_np1, A_rho_u_np1);
const CellArray<const double> M_rho_E_np1 = compute_scalar_M(rho_E_np1, A_rho_E_np1);
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
if (not is_inflow_boundary_cell[cell_id]) {
const double c1 = 1. / (m_eps + m_dt);
const double c2 = m_eps * m_dt * m_inv_Vj[cell_id];
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
Fnp1_rho[cell_id][i_wave] = c1 * (m_eps * Fn_rho[cell_id][i_wave] - c2 * deltaLFn_rho[cell_id][i_wave] +
m_dt * M_rho_np1[cell_id][i_wave]);
Fnp1_rho_u[cell_id][i_wave] =
c1 * (m_eps * Fn_rho_u[cell_id][i_wave] - c2 * deltaLFn_rho_u[cell_id][i_wave] +
m_dt * M_rho_u_np1[cell_id][i_wave]);
Fnp1_rho_E[cell_id][i_wave] =
c1 * (m_eps * Fn_rho_E[cell_id][i_wave] - c2 * deltaLFn_rho_E[cell_id][i_wave] +
m_dt * M_rho_E_np1[cell_id][i_wave]);
}
} else {
for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
Fnp1_rho[cell_id][i_wave] = Fn_rho_ex[cell_id][i_wave];
Fnp1_rho_u[cell_id][i_wave] = Fn_rho_u_ex[cell_id][i_wave]; Fnp1_rho_E[cell_id][i_wave] = Fn_rho_E_ex[cell_id][i_wave];
}
}
});
return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(Fnp1_rho),
std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_u),
std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_E));
}
EulerKineticSolverMultiD(std::shared_ptr<const MeshType> mesh,
const double& dt,
const double& eps,
const double& gamma,
const std::vector<TinyVector<MeshType::Dimension>>& lambda_vector,
const size_t& SpaceOrder,
const CellArray<const double>& Fn_rho,
const CellArray<const TinyVector<MeshType::Dimension>>& Fn_rho_u,
const CellArray<const double>& Fn_rho_E)
: m_dt{dt},
m_eps{eps},
m_lambda_vector{lambda_vector},
m_SpaceOrder{SpaceOrder},
m_Fn_rho{Fn_rho},
m_Fn_rho_u{Fn_rho_u},
m_Fn_rho_E{Fn_rho_E},
m_gamma{gamma},
p_mesh{mesh},
m_mesh{*mesh}
{
if ((lambda_vector.size() != Fn_rho[CellId(0)].size()) or (lambda_vector.size() != Fn_rho_u[CellId(0)].size()) or
(lambda_vector.size() != Fn_rho_u[CellId(0)].size()) or
((lambda_vector.size() != Fn_rho_E[CellId(0)].size()))) {
throw NormalError("Incompatible dimensions of lambda vector and Fn");
}
m_inv_Vj = [&] {
CellValue<double> inv_Vj{m_mesh.connectivity()};
parallel_for(
m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { inv_Vj[cell_id] = 1. / m_Vj[cell_id]; });
return inv_Vj;
}();
}
~EulerKineticSolverMultiD() = default;
};
template <MeshConcept MeshType>
class EulerKineticSolverMultiD<MeshType>::SymmetryBoundaryCondition
{
public:
static constexpr const size_t Dimension = MeshType::Dimension;
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(); }
const Array<const NodeId>&
nodeList() const
{
return m_mesh_flat_node_boundary.nodeList();
}
size_t
numberOfFaces() const
{
return m_mesh_flat_face_boundary.faceList().size();
}
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 <MeshConcept MeshType>
class EulerKineticSolverMultiD<MeshType>::WallBoundaryCondition
{
public:
static constexpr const size_t Dimension = MeshType::Dimension;
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();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
size_t
numberOfFaces() const
{
return m_mesh_face_boundary.faceList().size();
}
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)
{ ;
}
~WallBoundaryCondition() = default;
};
template <MeshConcept MeshType>
class EulerKineticSolverMultiD<MeshType>::InflowListBoundaryCondition
{
public:
static constexpr const size_t Dimension = MeshType::Dimension;
using Rd = TinyVector<Dimension, double>;
private:
const Table<const double> m_cell_array_list;
const Array<const CellId> m_cell_list;
public:
size_t
numberOfCells() const
{
return m_cell_list.size();
}
const Array<const CellId>&
cellList() const
{
return m_cell_list;
}
const Table<const double>&
cellArrayList() const
{
return m_cell_array_list;
}
InflowListBoundaryCondition(const Table<const double>& cell_array_list, const Array<const CellId> cell_list)
: m_cell_array_list(cell_array_list), m_cell_list(cell_list)
{
;
}
~InflowListBoundaryCondition() = default;
};
template <MeshConcept MeshType>
class EulerKineticSolverMultiD<MeshType>::OutflowBoundaryCondition
{
public:
static constexpr const size_t Dimension = MeshType::Dimension;
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();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
size_t
numberOfFaces() const
{
return m_mesh_face_boundary.faceList().size();
}
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)
{
;
}
~OutflowBoundaryCondition() = default;
};
template <size_t Dimension>
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
eulerKineticSolverMultiD(const double& dt,
const double& gamma,
const std::vector<TinyVector<Dimension>>& lambda_vector,
const double& eps,
const size_t& SpaceOrder,
const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
{
const std::shared_ptr<const MeshVariant> mesh_v = getCommonMesh({F_rho, F_rho_u, F_rho_E});
if (mesh_v.use_count() == 0) {
throw NormalError("F_rho, F_rho_u1, F_rho_u2 and F_rho_E have to be defined on the same mesh");
}
if (SpaceOrder > 1) {
throw NotImplementedError("Euler kinetic solver in Multi D not implemented at order > 1");
}
return std::visit(
[&](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)>;
if constexpr (is_polygonal_mesh_v<MeshType> and (MeshType::Dimension == Dimension)) {
auto Fn_rho = F_rho->get<DiscreteFunctionP0Vector<const double>>();
auto Fn_rho_u = F_rho_u->get<DiscreteFunctionP0Vector<const TinyVector<MeshType::Dimension>>>();
auto Fn_rho_E = F_rho_E->get<DiscreteFunctionP0Vector<const double>>();
EulerKineticSolverMultiD solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, Fn_rho.cellArrays(),
Fn_rho_u.cellArrays(), Fn_rho_E.cellArrays());
return solver.apply(bc_descriptor_list);
} else {
throw NotImplementedError("Invalid mesh type for Multi D Euler Kinetic solver");
}
},
mesh_v->variant());
}
template std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
eulerKineticSolverMultiD(const double& dt,
const double& gamma, const std::vector<TinyVector<2>>& lambda_vector,
const double& eps,
const size_t& SpaceOrder,
const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
template std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
eulerKineticSolverMultiD(const double& dt,
const double& gamma,
const std::vector<TinyVector<3>>& lambda_vector,
const double& eps,
const size_t& SpaceOrder,
const std::shared_ptr<const DiscreteFunctionVariant>& F_rho,
const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);