diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt
index 1ecbe1db54c85a4282ed1a55a043115f09bc7525..b94d4249c0bc2d013d62384d7329e9914feccf84 100644
--- a/src/language/modules/CMakeLists.txt
+++ b/src/language/modules/CMakeLists.txt
@@ -20,6 +20,7 @@ add_library(
target_link_libraries(
PugsLanguageModules
${HIGHFIVE_TARGET}
+ PugsScheme
)
add_dependencies(
@@ -27,5 +28,6 @@ add_dependencies(
PugsLanguageModules
PugsLanguageAlgorithms
PugsUtils
+ PugsScheme
PugsMesh
)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index d9c6f0d141d0d3ed5b259e673c70dbad56f4061e..ef3e5dbc9dc6e501e185aface5817480b94f4ac5 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -36,7 +36,6 @@
#include <scheme/DiscreteFunctionVectorIntegrator.hpp>
#include <scheme/DiscreteFunctionVectorInterpoler.hpp>
#include <scheme/EulerKineticSolver.hpp>
-#include <scheme/EulerKineticSolver2D.hpp>
#include <scheme/EulerKineticSolverAcoustic2LagrangeFVOrder2.hpp>
#include <scheme/EulerKineticSolverAcousticLagrangeFV.hpp>
#include <scheme/EulerKineticSolverFirstOrderFV.hpp>
@@ -44,6 +43,7 @@
#include <scheme/EulerKineticSolverMeanFluxMood.hpp>
#include <scheme/EulerKineticSolverMoodFD.hpp>
#include <scheme/EulerKineticSolverMoodFV.hpp>
+#include <scheme/EulerKineticSolverMultiD.hpp>
#include <scheme/EulerKineticSolverOneFluxMood.hpp>
#include <scheme/EulerKineticSolverThirdOrderFD.hpp>
#include <scheme/EulerKineticSolverThirdOrderFV.hpp>
@@ -956,47 +956,62 @@ SchemeModule::SchemeModule()
));
- this->_addBuiltinFunction("euler_kinetic_2D",
+ this->_addBuiltinFunction("euler_kinetic_MultiD",
std::function(
[](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_u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u2,
+ const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u,
const std::shared_ptr<const DiscreteFunctionVariant>& F_E,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho,
+ F_rho_u, F_E, bc_descriptor_list);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("euler_kinetic_MultiD",
+ std::function(
+
+ [](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_E,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
- return eulerKineticSolver2D(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho, F_rho_u1,
- F_rho_u2, F_E, bc_descriptor_list);
+ return eulerKineticSolverMultiD(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho,
+ F_rho_u, F_E, bc_descriptor_list);
}
));
- this->_addBuiltinFunction("get_lambda_vector",
+ this->_addBuiltinFunction("get_lambda_vector_2d",
std::function(
[](const std::shared_ptr<const MeshVariant>& mesh, const double& lambda, const size_t& L,
const size_t& M) -> std::vector<TinyVector<2>> {
- return getLambdaVector(mesh, lambda, L, M);
+ return getLambdaVector2D(mesh, lambda, L, M);
}
));
- this->_addBuiltinFunction("get_velocity", std::function(
+ this->_addBuiltinFunction("get_lambda_vector_3d", std::function(
- [](const std::shared_ptr<const DiscreteFunctionVariant>& Fn_rho_u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& Fn_rho_u2)
- -> std::shared_ptr<const DiscreteFunctionVariant> {
- return getVelocity(Fn_rho_u1, Fn_rho_u2);
- }
+ [](const std::shared_ptr<const MeshVariant>& mesh,
+ const double& lambda) -> std::vector<TinyVector<3>> {
+ return getLambdaVector3D(mesh, lambda);
+ }
- ));
+ ));
- this->_addBuiltinFunction("get_euler_kinetic_waves_2D",
+ this->_addBuiltinFunction("get_euler_kinetic_waves_MultiD",
std::function(
[](const std::vector<TinyVector<2>>& lambda_vector,
@@ -1005,9 +1020,23 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& rho_E,
const double& gamma) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return getEulerKineticWavesMultiD(lambda_vector, rho, rho_u, rho_E, gamma);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("get_euler_kinetic_waves_MultiD",
+ std::function(
+
+ [](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) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
- return getEulerKineticWaves2D(lambda_vector, rho, rho_u, rho_E, gamma);
+ return getEulerKineticWavesMultiD(lambda_vector, rho, rho_u, rho_E, gamma);
}
));
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index e26e157ec3c145d4d53397859a4ff6f4de1e0c3f..e243b52496e81225803fe7ae28bd320ff44748be 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -21,7 +21,7 @@ add_library(
EulerKineticSolverOneFluxMood.cpp
EulerKineticSolverMoodFD.cpp
EulerKineticSolverMoodFV.cpp
- EulerKineticSolver2D.cpp
+ EulerKineticSolverMultiD.cpp
EulerKineticSolverThirdOrderMoodFV.cpp
EulerKineticSolverThirdOrderMoodFD.cpp
EulerKineticSolverThirdOrderFV.cpp
diff --git a/src/scheme/EulerKineticSolver2D.hpp b/src/scheme/EulerKineticSolver2D.hpp
deleted file mode 100644
index a7937ee7ae067059b671d10b05674e109b7d9b88..0000000000000000000000000000000000000000
--- a/src/scheme/EulerKineticSolver2D.hpp
+++ /dev/null
@@ -1,43 +0,0 @@
-#ifndef EULER_KINETIC_SOLVER_2D_HPP
-#define EULER_KINETIC_SOLVER_2D_HPP
-
-#include <language/utils/FunctionSymbolId.hpp>
-#include <scheme/DiscreteFunctionVariant.hpp>
-
-class IBoundaryConditionDescriptor;
-
-std::vector<TinyVector<2>> getLambdaVector(const std::shared_ptr<const MeshVariant>& mesh,
- const double& lambda,
- const size_t& L,
- const size_t& M);
-
-std::shared_ptr<const DiscreteFunctionVariant> getVelocity(
- const std::shared_ptr<const DiscreteFunctionVariant>& Fn_rho_u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& Fn_rho_u2);
-
-std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
-getEulerKineticWaves2D(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);
-
-std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>>
-eulerKineticSolver2D(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_u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u2,
- const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
-
-#endif // EULER_KINETIC_SOLVER_2D_HPP
diff --git a/src/scheme/EulerKineticSolver2D.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
similarity index 69%
rename from src/scheme/EulerKineticSolver2D.cpp
rename to src/scheme/EulerKineticSolverMultiD.cpp
index 8f7f0af61e0d92441db063036041bcecfb9810e8..3ea60272a74641a139f9acfb34174554e770c4bd 100644
--- a/src/scheme/EulerKineticSolver2D.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -1,4 +1,4 @@
-#include <scheme/EulerKineticSolver2D.hpp>
+#include <scheme/EulerKineticSolverMultiD.hpp>
#include <language/utils/EvaluateAtPoints.hpp>
#include <language/utils/InterpolateItemArray.hpp>
@@ -29,7 +29,10 @@
#include <tuple>
std::vector<TinyVector<2>>
-getLambdaVector(const std::shared_ptr<const MeshVariant>& mesh, const double& lambda, const size_t& L, const size_t& M)
+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>> {
@@ -51,59 +54,46 @@ getLambdaVector(const std::shared_ptr<const MeshVariant>& mesh, const double& la
return lambda_vector;
} else {
- throw NotImplementedError("Invalid mesh type for 2D Euler Kinetic solver");
+ throw NotImplementedError("Invalid mesh type for MultiD Euler Kinetic solver");
}
},
mesh->variant());
}
-std::shared_ptr<const DiscreteFunctionVariant>
-getVelocity(const std::shared_ptr<const DiscreteFunctionVariant>& Fn_rho_u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& Fn_rho_u2)
+std::vector<TinyVector<3>>
+getLambdaVector3D(const std::shared_ptr<const MeshVariant>& mesh, const double& lambda)
{
- const std::shared_ptr<const MeshVariant> mesh_v = getCommonMesh({Fn_rho_u1, Fn_rho_u2});
- if (mesh_v.use_count() == 0) {
- throw NormalError("Fn_rho_u1 and Fn_rho_u2 have to be defined on the same mesh");
- }
return std::visit(
- [&](auto&& p_mesh) -> std::shared_ptr<const DiscreteFunctionVariant> {
+ [&](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 == 2)) {
- auto F_rho_u1 = Fn_rho_u1->get<DiscreteFunctionP0Vector<const double>>();
- auto F_rho_u2 = Fn_rho_u2->get<DiscreteFunctionP0Vector<const double>>();
-
- DiscreteFunctionP0<TinyVector<MeshType::Dimension>> rho_u{p_mesh};
-
- parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- rho_u[cell_id][0] = 0;
- rho_u[cell_id][1] = 0;
- size_t nb_waves = F_rho_u1[cell_id].size();
- for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
- rho_u[cell_id][0] += F_rho_u1[cell_id][i_wave];
- rho_u[cell_id][1] += F_rho_u2[cell_id][i_wave];
- }
- });
+ if constexpr (is_polygonal_mesh_v<MeshType> and (MeshType::Dimension == 3)) {
+ std::vector<TinyVector<3>> lambda_vector{6};
- return std::make_shared<DiscreteFunctionVariant>(rho_u);
+ 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 2D Euler Kinetic solver");
+ throw NotImplementedError("Invalid mesh type for MultiD Euler Kinetic solver");
}
},
- mesh_v->variant());
+ mesh->variant());
}
template <MeshConcept MeshType>
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
-getEulerKineticWaves2D(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)
+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");
@@ -115,71 +105,67 @@ getEulerKineticWaves2D(const std::shared_ptr<const MeshType> p_mesh,
const CellValue<const double> m_Vj = mesh_data.Vj();
DiscreteFunctionP0Vector<double> Fn_rho(p_mesh, lambda_vector.size());
- DiscreteFunctionP0Vector<double> Fn_rho_u1(p_mesh, lambda_vector.size());
- DiscreteFunctionP0Vector<double> Fn_rho_u2(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();
- CellArray<TinyVector<MeshType::Dimension>> vec_A{p_mesh->connectivity(), 4};
+ 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;
- vec_A[cell_id][0][0] = rho_u[cell_id][0];
- vec_A[cell_id][1][0] = rho_u[cell_id][0] * rho_u[cell_id][0] / rho[cell_id] + p;
- vec_A[cell_id][2][0] = rho_u[cell_id][0] * rho_u[cell_id][1] / rho[cell_id];
- vec_A[cell_id][3][0] = (rho_E[cell_id] + p) * rho_u[cell_id][0] / rho[cell_id];
-
- vec_A[cell_id][0][1] = rho_u[cell_id][1];
- vec_A[cell_id][1][1] = rho_u[cell_id][0] * rho_u[cell_id][1] / rho[cell_id];
- vec_A[cell_id][2][1] = rho_u[cell_id][1] * rho_u[cell_id][1] / rho[cell_id] + p;
- vec_A[cell_id][3][1] = (rho_E[cell_id] + p) * rho_u[cell_id][1] / rho[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 i = 0; i < nb_waves; ++i) {
- for (size_t d = 0; d < MeshType::Dimension; ++d) {
- inv_S[d] += std::pow(lambda_vector[i][d], 2);
- }
- }
+
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] + inv_S[0] * vec_A[cell_id][0][0] * lambda_vector[i][0] +
- inv_S[1] * vec_A[cell_id][0][1] * lambda_vector[i][1];
- Fn_rho_u1[cell_id][i] = (1. / nb_waves) * rho_u[cell_id][0] +
- inv_S[0] * vec_A[cell_id][1][0] * lambda_vector[i][0] +
- inv_S[1] * vec_A[cell_id][1][1] * lambda_vector[i][1];
- Fn_rho_u2[cell_id][i] = (1. / nb_waves) * rho_u[cell_id][1] +
- inv_S[0] * vec_A[cell_id][2][0] * lambda_vector[i][0] +
- inv_S[1] * vec_A[cell_id][2][1] * lambda_vector[i][1];
- Fn_rho_E[cell_id][i] = (1. / nb_waves) * rho_E[cell_id] +
- inv_S[0] * vec_A[cell_id][3][0] * lambda_vector[i][0] +
- inv_S[1] * vec_A[cell_id][3][1] * lambda_vector[i][1];
+ 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_u1),
- std::make_shared<DiscreteFunctionVariant>(Fn_rho_u2),
+ 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>,
std::shared_ptr<const DiscreteFunctionVariant>>
-getEulerKineticWaves2D(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)
+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) {
@@ -189,23 +175,41 @@ getEulerKineticWaves2D(const std::vector<TinyVector<2>>& lambda_vector,
return std::visit(
[&](auto&& p_mesh)
-> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
- 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 == 2)) {
+ 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 getEulerKineticWaves2D(p_mesh, lambda_vector, rho_h, rho_u_h, rho_E_h, gamma);
+ return GetEulerKineticWavesMultiD(p_mesh, lambda_vector, rho_h, rho_u_h, rho_E_h, gamma);
} else {
- throw NotImplementedError("Invalid mesh type for 2D Euler Kinetic solver");
+ 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 EulerKineticSolver2D
+class EulerKineticSolverMultiD
{
private:
constexpr static size_t Dimension = MeshType::Dimension;
@@ -225,18 +229,14 @@ class EulerKineticSolver2D
const std::vector<TinyVector<Dimension>>& m_lambda_vector;
const size_t m_SpaceOrder;
const CellArray<const double> m_Fn_rho;
- const CellArray<TinyVector<Dimension>> m_Fn_rho_u;
+ 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 NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
- const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
- const NodeValuePerFace<const TinyVector<Dimension>> Nlr = MeshDataManager::instance().getMeshData(m_mesh).Nlr();
- const NodeValuePerFace<const TinyVector<Dimension>> nlr = MeshDataManager::instance().getMeshData(m_mesh).nlr();
- const FaceValue<const TinyVector<Dimension>> nl = MeshDataManager::instance().getMeshData(m_mesh).nl();
- const FaceValue<const TinyVector<Dimension>> Nl = MeshDataManager::instance().getMeshData(m_mesh).Nl();
+ 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();
@@ -508,97 +508,107 @@ class EulerKineticSolver2D
NodeArray<T>
compute_Flux1(CellArray<T>& Fn)
{
- 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.);
- }
+ 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;
+ 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];
+ 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;
+ 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;
}
}
- if (sum_li_njr != 0) {
- Fr[node_id][i_wave] /= sum_li_njr;
- }
- }
- });
+ });
- return Fr;
+ return Fr;
+ } else {
+ throw NormalError("Glace not defined in 1d");
+ }
}
template <typename T>
FaceArrayPerNode<T>
compute_Flux1_eucchlyd(CellArray<T> Fn)
{
- 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.);
- }
+ 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];
+ 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];
+ 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;
+ 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];
+ 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;
+ 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);
+ double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
- if (li_nlr > 0) {
- Flr[node_id][i_face][i_wave] += Fn[cell_id][i_wave] * li_nlr;
- sum += li_nlr;
+ if (li_nlr > 0) {
+ Flr[node_id][i_face][i_wave] += Fn[cell_id][i_wave] * li_nlr;
+ sum += li_nlr;
+ }
+ }
+ if (sum != 0) {
+ Flr[node_id][i_face][i_wave] /= sum;
}
- }
- if (sum != 0) {
- Flr[node_id][i_face][i_wave] /= sum;
}
}
- }
- });
+ });
- return Flr;
+ return Flr;
+ } else {
+ throw NormalError("Eucclhyd not defined in 1d");
+ }
}
template <typename T>
@@ -815,110 +825,126 @@ class EulerKineticSolver2D
DiscreteFunctionP0Vector<double>
compute_deltaLFn(NodeArray<double> Fr)
{
- const size_t nb_waves = m_lambda_vector.size();
- DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
+ 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();
+ DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
- const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+ const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
- parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
- const auto& cell_to_node = cell_to_node_matrix[cell_id];
+ 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) {
- deltaLFn[cell_id][i_wave] = 0;
- for (size_t i_node = 0; i_node < cell_to_node.size(); ++i_node) {
- const NodeId node_id = cell_to_node[i_node];
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ deltaLFn[cell_id][i_wave] = 0;
+ 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] += Fr[node_id][i_wave] * li_Cjr;
+ double li_Cjr = dot(m_lambda_vector[i_wave], Cjr(cell_id, i_node));
+ deltaLFn[cell_id][i_wave] += Fr[node_id][i_wave] * li_Cjr;
+ }
}
- }
- });
- return deltaLFn;
+ });
+ return deltaLFn;
+ } else {
+ throw NormalError("Glace not defined in 1d");
+ }
}
DiscreteFunctionP0Vector<double>
compute_deltaLFn_eucclhyd(FaceArrayPerNode<double> Fr)
{
- const size_t nb_waves = m_lambda_vector.size();
- DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
+ if constexpr (Dimension > 1) {
+ const size_t nb_waves = m_lambda_vector.size();
+ DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
- 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();
+ const NodeValuePerFace<const TinyVector<Dimension>> Nlr = MeshDataManager::instance().getMeshData(m_mesh).Nlr();
- parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
- const auto& cell_to_face = cell_to_face_matrix[cell_id];
+ 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();
- for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
- deltaLFn[cell_id][i_wave] = 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) {
+ deltaLFn[cell_id][i_wave] = 0;
- 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_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];
+ 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;
+ 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();
- };
+ return std::numeric_limits<size_t>::max();
+ };
- const size_t i_face_node = local_face_number_in_node(face_id);
+ 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 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));
+ const double li_Nlr = sign * dot(m_lambda_vector[i_wave], Nlr(face_id, i_node_face));
- deltaLFn[cell_id][i_wave] += Fr[node_id][i_face_node][i_wave] * li_Nlr;
+ deltaLFn[cell_id][i_wave] += Fr[node_id][i_face_node][i_wave] * li_Nlr;
+ }
}
}
- }
- });
- return deltaLFn;
+ });
+ return deltaLFn;
+ } else {
+ throw NormalError("Eucclhyd not defined in 1d");
+ }
}
template <typename T>
CellArray<T>
compute_deltaLFn_face(FaceArray<T> Fl)
{
- const size_t nb_waves = m_lambda_vector.size();
- CellArray<T> deltaLFn(p_mesh->connectivity(), nb_waves);
+ 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();
+ 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.);
- }
+ 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];
+ 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];
+ 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 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]);
+ 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];
+ deltaLFn[cell_id][i_wave] += li_Nl * Fl[face_id][i_wave];
+ }
}
- }
- });
- return deltaLFn;
+ });
+ return deltaLFn;
+ } else {
+ throw NormalError("Nl in meaningless in 1d");
+ }
}
CellValue<bool>
@@ -947,7 +973,6 @@ class EulerKineticSolver2D
}
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
- 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)
@@ -998,8 +1023,7 @@ class EulerKineticSolver2D
const CellArray<const double> Fn_rho_E = copy(m_Fn_rho_E);
DiscreteFunctionP0Vector<double> Fnp1_rho(p_mesh, nb_waves);
- DiscreteFunctionP0Vector<double> Fnp1_rho_u1(p_mesh, nb_waves);
- DiscreteFunctionP0Vector<double> Fnp1_rho_u2(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
@@ -1091,41 +1115,36 @@ class EulerKineticSolver2D
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_u1[cell_id][i_wave] =
- c1 * (m_eps * Fn_rho_u[cell_id][i_wave][0] - c2 * deltaLFn_rho_u[cell_id][i_wave][0] +
- m_dt * M_rho_u_np1[cell_id][i_wave][0]);
- Fnp1_rho_u2[cell_id][i_wave] =
- c1 * (m_eps * Fn_rho_u[cell_id][i_wave][1] - c2 * deltaLFn_rho_u[cell_id][i_wave][1] +
- m_dt * M_rho_u_np1[cell_id][i_wave][1]);
+ 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_u1[cell_id][i_wave] = Fn_rho_u_ex[cell_id][i_wave][0];
- Fnp1_rho_u2[cell_id][i_wave] = Fn_rho_u_ex[cell_id][i_wave][1];
- Fnp1_rho_E[cell_id][i_wave] = Fn_rho_E_ex[cell_id][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_u1),
- std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_u2),
+ std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_u),
std::make_shared<DiscreteFunctionVariant>(Fnp1_rho_E));
}
- EulerKineticSolver2D(std::shared_ptr<const MeshType> mesh,
- const double& dt,
- const double& eps,
- const double& gamma,
- const std::vector<TinyVector<2>>& lambda_vector,
- const size_t& SpaceOrder,
- const CellArray<const double>& Fn_rho,
- const CellArray<TinyVector<2>>& Fn_rho_u,
- const CellArray<const double>& Fn_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},
@@ -1151,11 +1170,11 @@ class EulerKineticSolver2D
}();
}
- ~EulerKineticSolver2D() = default;
+ ~EulerKineticSolverMultiD() = default;
};
template <MeshConcept MeshType>
-class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition
+class EulerKineticSolverMultiD<MeshType>::SymmetryBoundaryCondition
{
public:
static constexpr const size_t Dimension = MeshType::Dimension;
@@ -1208,7 +1227,7 @@ class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition
};
template <MeshConcept MeshType>
-class EulerKineticSolver2D<MeshType>::WallBoundaryCondition
+class EulerKineticSolverMultiD<MeshType>::WallBoundaryCondition
{
public:
static constexpr const size_t Dimension = MeshType::Dimension;
@@ -1254,7 +1273,7 @@ class EulerKineticSolver2D<MeshType>::WallBoundaryCondition
};
template <MeshConcept MeshType>
-class EulerKineticSolver2D<MeshType>::InflowListBoundaryCondition
+class EulerKineticSolverMultiD<MeshType>::InflowListBoundaryCondition
{
public:
static constexpr const size_t Dimension = MeshType::Dimension;
@@ -1294,7 +1313,7 @@ class EulerKineticSolver2D<MeshType>::InflowListBoundaryCondition
};
template <MeshConcept MeshType>
-class EulerKineticSolver2D<MeshType>::OutflowBoundaryCondition
+class EulerKineticSolverMultiD<MeshType>::OutflowBoundaryCondition
{
public:
static constexpr const size_t Dimension = MeshType::Dimension;
@@ -1339,60 +1358,71 @@ class EulerKineticSolver2D<MeshType>::OutflowBoundaryCondition
~OutflowBoundaryCondition() = default;
};
+template <size_t Dimension>
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
- std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
-eulerKineticSolver2D(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_u1,
- const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_u2,
- const std::shared_ptr<const DiscreteFunctionVariant>& F_rho_E,
- const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+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_u1, F_rho_u2, F_rho_E});
+ 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 2D not implemented at order > 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>, 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 == 2)) {
- auto Fn_rho = F_rho->get<DiscreteFunctionP0Vector<const double>>();
- auto Fn_rho_u1 = F_rho_u1->get<DiscreteFunctionP0Vector<const double>>();
- auto Fn_rho_u2 = F_rho_u2->get<DiscreteFunctionP0Vector<const double>>();
- auto Fn_rho_E = F_rho_E->get<DiscreteFunctionP0Vector<const double>>();
-
- const double nb_waves = lambda_vector.size();
- CellArray<const double> Fn_rho_array = Fn_rho.cellArrays();
- CellArray<TinyVector<MeshType::Dimension>> Fn_rho_u_array(p_mesh->connectivity(), nb_waves);
- CellArray<const double> Fn_rho_E_array = Fn_rho_E.cellArrays();
-
- parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
- Fn_rho_u_array[cell_id][i_wave][0] = Fn_rho_u1[cell_id][i_wave];
- Fn_rho_u_array[cell_id][i_wave][1] = Fn_rho_u2[cell_id][i_wave];
- }
- });
+ 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>>();
- EulerKineticSolver2D solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, Fn_rho_array, Fn_rho_u_array,
- Fn_rho_E_array);
+ 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 2D Euler Kinetic solver");
+ 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);
diff --git a/src/scheme/EulerKineticSolverMultiD.hpp b/src/scheme/EulerKineticSolverMultiD.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..987f10f8a97189ca0db782cf0545e92b0d5bb933
--- /dev/null
+++ b/src/scheme/EulerKineticSolverMultiD.hpp
@@ -0,0 +1,40 @@
+#ifndef EULER_KINETIC_SOLVER_MULTID_HPP
+#define EULER_KINETIC_SOLVER_MULTID_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+
+class IBoundaryConditionDescriptor;
+
+std::vector<TinyVector<2>> getLambdaVector2D(const std::shared_ptr<const MeshVariant>& mesh,
+ const double& lambda,
+ const size_t& L,
+ const size_t& M);
+
+std::vector<TinyVector<3>> getLambdaVector3D(const std::shared_ptr<const MeshVariant>& mesh, const double& lambda);
+
+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);
+
+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);
+
+#endif // EULER_KINETIC_SOLVER_MULTID_HPP