diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index e1af0c003f55eff3855b1a16ae9a0c452d6ce8d4..7d144f7e981313141f16ad40f5da0baa4211fb6a 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -941,17 +941,13 @@ SchemeModule::SchemeModule()
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_E,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex,
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>,
std::shared_ptr<const DiscreteFunctionVariant>> {
return eulerKineticSolver2D(dt, gamma, lambda_vector, eps, SpaceOrder, F_rho, F_rho_u1,
- F_rho_u2, F_E, rho_ex, rho_u_ex, rho_E_ex,
- bc_descriptor_list);
+ F_rho_u2, F_E, bc_descriptor_list);
}
));
diff --git a/src/scheme/EulerKineticSolver2D.cpp b/src/scheme/EulerKineticSolver2D.cpp
index f2130599b3b744015108d7514ad2fc1e8d137844..281db07dc6b593de71798d906a764043e3136bf6 100644
--- a/src/scheme/EulerKineticSolver2D.cpp
+++ b/src/scheme/EulerKineticSolver2D.cpp
@@ -217,13 +217,9 @@ class EulerKineticSolver2D
const double m_eps;
const std::vector<TinyVector<Dimension>>& m_lambda_vector;
const size_t m_SpaceOrder;
- const DiscreteFunctionP0Vector<const double> m_Fn_rho;
- const DiscreteFunctionP0Vector<const double> m_Fn_rho_u1;
- const DiscreteFunctionP0Vector<const double> m_Fn_rho_u2;
- const DiscreteFunctionP0Vector<const double> m_Fn_rho_E;
- const DiscreteFunctionP0<const double> m_rho_ex;
- const DiscreteFunctionP0<const TinyVector<Dimension>> m_rho_u_ex;
- const DiscreteFunctionP0<const double> m_rho_E_ex;
+ const CellArray<const double> m_Fn_rho;
+ const CellArray<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;
@@ -344,43 +340,123 @@ class EulerKineticSolver2D
return vec_A;
}
- DiscreteFunctionP0Vector<double>
- compute_M(const DiscreteFunctionP0<const double>& u, const DiscreteFunctionP0<TinyVector<MeshType::Dimension>> A_u)
+ const CellValue<const TinyVector<Dimension>>
+ getA_rho(const CellValue<const TinyVector<Dimension>>& rho_u) const
{
- const size_t nb_waves = m_lambda_vector.size();
- DiscreteFunctionP0Vector<double> M{p_mesh, nb_waves};
+ return rho_u;
+ }
- 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(m_lambda_vector[i][d], 2);
+ 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 = 0; i < nb_waves; ++i) {
- M[cell_id][i] = (1. / nb_waves) * u[cell_id] + inv_S[0] * A_u[cell_id][0] * m_lambda_vector[i][0] +
- inv_S[1] * A_u[cell_id][1] * m_lambda_vector[i][1];
+ 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;
}
- DiscreteFunctionP0<double>
- compute_PFnp1(const DiscreteFunctionP0Vector<const double>& Fn, DiscreteFunctionP0Vector<double>& deltaLFn)
+ template <typename T>
+ const CellValue<const T>
+ compute_PFnp1(const CellArray<const T>& Fn, const CellArray<const T>& deltaLFn)
{
- DiscreteFunctionP0<double> PFnp1{p_mesh};
+ 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;
- PFnp1[cell_id] = 0;
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];
}
@@ -389,37 +465,49 @@ class EulerKineticSolver2D
return PFnp1;
}
- DiscreteFunctionP0<double>
- compute_PFn(const DiscreteFunctionP0Vector<const double> F)
+ template <typename T>
+ const CellValue<const T>
+ compute_PFn(const CellArray<const T>& F)
{
- DiscreteFunctionP0<double> PFn{p_mesh};
+ 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) {
- double PF = 0;
- for (size_t i = 0; i < m_lambda_vector.size(); ++i) {
- PF += F[cell_id][i];
+ for (size_t i = 0; i < nb_waves; ++i) {
+ PFn[cell_id] += F[cell_id][i];
}
- PFn[cell_id] = PF;
});
return PFn;
}
- NodeArray<double>
- compute_Flux1(DiscreteFunctionP0Vector<double> Fn)
+ template <typename T>
+ NodeArray<T>
+ compute_Flux1(CellArray<T>& Fn)
{
const size_t nb_waves = m_lambda_vector.size();
- NodeArray<double> Fr(m_mesh.connectivity(), nb_waves);
+ 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;
- Fr[node_id][i_wave] = 0;
+ 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];
@@ -439,17 +527,24 @@ class EulerKineticSolver2D
return Fr;
}
- FaceArrayPerNode<double>
- compute_Flux1_eucchlyd(DiscreteFunctionP0Vector<double> Fn)
+ template <typename T>
+ FaceArrayPerNode<T>
+ compute_Flux1_eucchlyd(CellArray<T> Fn)
{
const size_t nb_waves = m_lambda_vector.size();
- FaceArrayPerNode<double> Fr(m_mesh.connectivity(), nb_waves);
+ 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);
@@ -461,8 +556,8 @@ class EulerKineticSolver2D
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];
- Fr[node_id][i_face][i_wave] = 0;
- 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];
@@ -474,37 +569,42 @@ class EulerKineticSolver2D
double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
if (li_nlr > 0) {
- Fr[node_id][i_face][i_wave] += Fn[cell_id][i_wave] * li_nlr;
+ Flr[node_id][i_face][i_wave] += Fn[cell_id][i_wave] * li_nlr;
sum += li_nlr;
}
}
if (sum != 0) {
- Fr[node_id][i_face][i_wave] /= sum;
+ Flr[node_id][i_face][i_wave] /= sum;
}
}
}
});
- return Fr;
+ return Flr;
}
- FaceArray<double>
- compute_Flux1_face(const DiscreteFunctionP0Vector<const double> Fn)
+ template <typename T>
+ FaceArray<T>
+ compute_Flux1_face(const CellArray<const T> Fn)
{
const size_t nb_waves = m_lambda_vector.size();
- FaceArray<double> Fl(m_mesh.connectivity(), nb_waves);
+ 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];
- Fl[face_id][i_wave] = 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];
@@ -525,15 +625,13 @@ class EulerKineticSolver2D
}
void
- apply_scalar_bc(FaceArray<double> Fl_rho,
- FaceArray<double> Fl_rho_u1,
- FaceArray<double> Fl_rho_u2,
- FaceArray<double> Fl_rho_E,
- const DiscreteFunctionP0<const double> rho,
- const DiscreteFunctionP0<const double> rho_u1,
- const DiscreteFunctionP0<const double> rho_u2,
- const DiscreteFunctionP0<const double> rho_E,
- const BoundaryConditionList& bc_list)
+ apply_symmetry_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 BoundaryConditionList& bc_list)
{
const size_t nb_waves = m_lambda_vector.size();
const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
@@ -541,12 +639,10 @@ class EulerKineticSolver2D
const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
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(m_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] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
+ }
inv_S[d] = 1. / inv_S[d];
}
@@ -577,39 +673,33 @@ class EulerKineticSolver2D
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_u1[cell_id], rho_u2[cell_id]};
+ 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]) *
- (rho_u1[cell_id] * rho_u1[cell_id] + rho_u2[cell_id] * rho_u2[cell_id]);
- double p = (m_gamma - 1) * rho_e;
+ 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;
- TinyVector<Dimension> A_rho_u1_prime{rho_uj_prime[0] * rho_uj_prime[0] / rhoj_prime + p,
- rho_uj_prime[0] * rho_uj_prime[1] / rhoj_prime};
- TinyVector<Dimension> A_rho_u2_prime{rho_uj_prime[0] * rho_uj_prime[1] / rhoj_prime,
- rho_uj_prime[1] * rho_uj_prime[1] / rhoj_prime + p};
+ 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 +
- inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_prime[0] +
- inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_prime[1];
- double Fn_rho_u1_prime = rho_uj_prime[0] / nb_waves +
- inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_u1_prime[0] +
- inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_u1_prime[1];
- double Fn_rho_u2_prime = rho_uj_prime[1] / nb_waves +
- inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_u2_prime[0] +
- inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_u2_prime[1];
- double Fn_rho_E_prime = rho_Ej_prime / nb_waves +
- inv_S[0] * m_lambda_vector[i_wave][0] * A_rho_E_prime[0] +
- inv_S[1] * m_lambda_vector[i_wave][1] * A_rho_E_prime[1];
+ 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_u1[face_id][i_wave] += Fn_rho_u1_prime;
- Fl_rho_u2[face_id][i_wave] += Fn_rho_u2_prime;
+ Fl_rho_u[face_id][i_wave] += Fn_rho_u_prime;
Fl_rho_E[face_id][i_wave] += Fn_rho_E_prime;
}
}
@@ -695,22 +785,27 @@ class EulerKineticSolver2D
return deltaLFn;
}
- DiscreteFunctionP0Vector<double>
- compute_deltaLFn_face(FaceArray<double> Fl)
+ template <typename T>
+ CellArray<T>
+ compute_deltaLFn_face(FaceArray<T> Fl)
{
const size_t nb_waves = m_lambda_vector.size();
- DiscreteFunctionP0Vector<double> deltaLFn(p_mesh, nb_waves);
+ 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) {
- 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];
@@ -718,7 +813,7 @@ class EulerKineticSolver2D
const double li_Nl = sign * dot(m_lambda_vector[i_wave], Nl[face_id]);
- deltaLFn[cell_id][i_wave] += Fl[face_id][i_wave] * li_Nl;
+ deltaLFn[cell_id][i_wave] += li_Nl * Fl[face_id][i_wave];
}
}
});
@@ -733,8 +828,6 @@ class EulerKineticSolver2D
is_boundary_cell.fill(false);
- auto node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
-
for (auto&& bc_v : bc_list) {
std::visit(
[=](auto&& bc) {
@@ -765,14 +858,12 @@ class EulerKineticSolver2D
CellValue<bool> is_boundary_cell = getBoundaryCells(bc_list);
- DiscreteFunctionP0<double> rho_ex{p_mesh};
- DiscreteFunctionP0<double> rho_u1_ex{p_mesh};
- DiscreteFunctionP0<double> rho_u2_ex{p_mesh};
- DiscreteFunctionP0<double> rho_E_ex{p_mesh};
+ 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);
- rho_u1_ex.fill(0);
- rho_u2_ex.fill(0);
+ 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) {
@@ -785,8 +876,8 @@ class EulerKineticSolver2D
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_u1_ex[cell_id] = cell_array_list[i_cell][1];
- rho_u2_ex[cell_id] = cell_array_list[i_cell][2];
+ 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];
}
}
@@ -794,29 +885,17 @@ class EulerKineticSolver2D
bc_v);
}
- const CellArray<const TinyVector<Dimension>> APFn_ex = getA(rho_ex, rho_u1_ex, rho_u2_ex, rho_E_ex);
- const DiscreteFunctionP0<TinyVector<Dimension>> A_rho_ex{p_mesh};
- const DiscreteFunctionP0<TinyVector<Dimension>> A_rho_u1_ex{p_mesh};
- const DiscreteFunctionP0<TinyVector<Dimension>> A_rho_u2_ex{p_mesh};
- const DiscreteFunctionP0<TinyVector<Dimension>> A_rho_E_ex{p_mesh};
-
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- A_rho_ex[cell_id] = APFn_ex[cell_id][0];
- A_rho_u1_ex[cell_id] = APFn_ex[cell_id][1];
- A_rho_u2_ex[cell_id] = APFn_ex[cell_id][2];
- A_rho_E_ex[cell_id] = APFn_ex[cell_id][3];
- });
+ 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);
- DiscreteFunctionP0Vector<double> Fn_rho_ex = compute_M(rho_ex, A_rho_ex);
- DiscreteFunctionP0Vector<double> Fn_rho_u1_ex = compute_M(rho_u1_ex, A_rho_u1_ex);
- DiscreteFunctionP0Vector<double> Fn_rho_u2_ex = compute_M(rho_u2_ex, A_rho_u2_ex);
- DiscreteFunctionP0Vector<double> Fn_rho_E_ex = compute_M(rho_E_ex, A_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);
- DiscreteFunctionP0Vector<double> Fn_rho = copy(m_Fn_rho);
- DiscreteFunctionP0Vector<double> Fn_rho_u1 = copy(m_Fn_rho_u1);
- DiscreteFunctionP0Vector<double> Fn_rho_u2 = copy(m_Fn_rho_u2);
- DiscreteFunctionP0Vector<double> Fn_rho_E = copy(m_Fn_rho_E);
+ 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<double> Fnp1_rho_u1(p_mesh, nb_waves);
@@ -825,105 +904,82 @@ class EulerKineticSolver2D
// Compute first order F
- const DiscreteFunctionP0<const double> rho = compute_PFn(Fn_rho);
- const DiscreteFunctionP0<const double> rho_u1 = compute_PFn(Fn_rho_u1);
- const DiscreteFunctionP0<const double> rho_u2 = compute_PFn(Fn_rho_u2);
- const DiscreteFunctionP0<const double> rho_E = compute_PFn(Fn_rho_E);
+ 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);
std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
- // Reconstruction uniquement à l'intérieur
+ // Reconstruction uniquement à l'intérieur donc pas utile
// for (auto&& bc_descriptor : bc_descriptor_list) {
// if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
// symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
// }
// }
- PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
- symmetry_boundary_descriptor_list);
+ // PolynomialReconstructionDescriptor reconstruction_descriptor(IntegrationMethodType::cell_center, 1,
+ // symmetry_boundary_descriptor_list);
- auto reconstruction =
- PolynomialReconstruction{reconstruction_descriptor}.build(Fn_rho, Fn_rho_u1, Fn_rho_u2, Fn_rho_E);
- auto DPk_Fn_rho = reconstruction[0]->get<DiscreteFunctionDPkVector<Dimension, const double>>();
+ // auto reconstruction = PolynomialReconstruction{reconstruction_descriptor}.build(Fn_rho, Fn_rho_u, Fn_rho_E);
+ // auto DPk_Fn_rho = reconstruction[0]->get<DiscreteFunctionDPkVector<Dimension, const double>>();
+ // auto DPk_Fn_rho_u = reconstruction[1]->get<DiscreteFunctionDPkVector<Dimension, const
+ // TinyVector<Dimension>>>(); auto DPk_Fn_rho_E = reconstruction[2]->get<DiscreteFunctionDPkVector<Dimension,
+ // const double>>();
// NodeArray<double> Fr_rho = compute_Flux1(Fn_rho);
// NodeArray<double> Fr_rho_u1 = compute_Flux1(Fn_rho_u1);
// NodeArray<double> Fr_rho_u2 = compute_Flux1(Fn_rho_u2);
// NodeArray<double> Fr_rho_E = compute_Flux1(Fn_rho_E);
- // FaceArrayPerNode<double> Fr_rho = compute_Flux1_eucchlyd(Fn_rho);
- // FaceArrayPerNode<double> Fr_rho_u1 = compute_Flux1_eucchlyd(Fn_rho_u1);
- // FaceArrayPerNode<double> Fr_rho_u2 = compute_Flux1_eucchlyd(Fn_rho_u2);
- // FaceArrayPerNode<double> Fr_rho_E = compute_Flux1_eucchlyd(Fn_rho_E);
+ // FaceArrayPerNode<double> Flr_rho = compute_Flux1_eucchlyd(Fn_rho);
+ // FaceArrayPerNode<double> Flr_rho_u1 = compute_Flux1_eucchlyd(Fn_rho_u1);
+ // FaceArrayPerNode<double> Flr_rho_u2 = compute_Flux1_eucchlyd(Fn_rho_u2);
+ // FaceArrayPerNode<double> Flr_rho_E = compute_Flux1_eucchlyd(Fn_rho_E);
- FaceArray<double> Fr_rho = compute_Flux1_face(Fn_rho);
- FaceArray<double> Fr_rho_u1 = compute_Flux1_face(Fn_rho_u1);
- FaceArray<double> Fr_rho_u2 = compute_Flux1_face(Fn_rho_u2);
- FaceArray<double> Fr_rho_E = compute_Flux1_face(Fn_rho_E);
+ 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_scalar_bc(Fr_rho, Fr_rho_u1, Fr_rho_u2, Fr_rho_E, rho, rho_u1, rho_u2, rho_E, bc_list);
+ apply_symmetry_bc(Fl_rho, Fl_rho_u, Fl_rho_E, rho, rho_u, rho_E, bc_list);
- synchronize(Fr_rho);
- synchronize(Fr_rho_u1);
- synchronize(Fr_rho_u2);
- synchronize(Fr_rho_E);
+ synchronize(Fl_rho);
+ synchronize(Fl_rho_u);
+ synchronize(Fl_rho_E);
// DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn(Fr_rho);
// DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn(Fr_rho_u1);
// DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn(Fr_rho_u2);
// DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn(Fr_rho_E);
- // DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn_eucclhyd(Fr_rho);
- // DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn_eucclhyd(Fr_rho_u1);
- // DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn_eucclhyd(Fr_rho_u2);
- // DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn_eucclhyd(Fr_rho_E);
+ // DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn_eucclhyd(Flr_rho);
+ // DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn_eucclhyd(Flr_rho_u1);
+ // DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn_eucclhyd(Flr_rho_u2);
+ // DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn_eucclhyd(Flr_rho_E);
- DiscreteFunctionP0Vector<double> deltaLFn_rho = compute_deltaLFn_face(Fr_rho);
- DiscreteFunctionP0Vector<double> deltaLFn_rho_u1 = compute_deltaLFn_face(Fr_rho_u1);
- DiscreteFunctionP0Vector<double> deltaLFn_rho_u2 = compute_deltaLFn_face(Fr_rho_u2);
- DiscreteFunctionP0Vector<double> deltaLFn_rho_E = compute_deltaLFn_face(Fr_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 CellArray<const TinyVector<Dimension>> APFn = getA(rho, rho_u1, rho_u2, rho_E);
- const DiscreteFunctionP0<TinyVector<Dimension>> APFn_rho{p_mesh};
- const DiscreteFunctionP0<TinyVector<Dimension>> APFn_rho_u1{p_mesh};
- const DiscreteFunctionP0<TinyVector<Dimension>> APFn_rho_u2{p_mesh};
- const DiscreteFunctionP0<TinyVector<Dimension>> APFn_rho_E{p_mesh};
+ 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);
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- APFn_rho[cell_id] = APFn[cell_id][0];
- APFn_rho_u1[cell_id] = APFn[cell_id][1];
- APFn_rho_u2[cell_id] = APFn[cell_id][2];
- APFn_rho_E[cell_id] = APFn[cell_id][3];
- });
+ 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 DiscreteFunctionP0Vector<const double> MPFn_rho = compute_M(rho, APFn_rho);
- const DiscreteFunctionP0Vector<const double> MPFn_rho_u1 = compute_M(rho_u1, APFn_rho_u1);
- const DiscreteFunctionP0Vector<const double> MPFn_rho_u2 = compute_M(rho_u2, APFn_rho_u2);
- const DiscreteFunctionP0Vector<const double> MPFn_rho_E = compute_M(rho_E, APFn_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 DiscreteFunctionP0<const double> rho_np1 = compute_PFnp1(Fn_rho, deltaLFn_rho);
- const DiscreteFunctionP0<const double> rho_u1_np1 = compute_PFnp1(Fn_rho_u1, deltaLFn_rho_u1);
- const DiscreteFunctionP0<const double> rho_u2_np1 = compute_PFnp1(Fn_rho_u2, deltaLFn_rho_u2);
- const DiscreteFunctionP0<const double> rho_E_np1 = compute_PFnp1(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 TinyVector<Dimension>> APFnp1 = getA(rho_np1, rho_u1_np1, rho_u2_np1, rho_E_np1);
- const DiscreteFunctionP0<TinyVector<Dimension>> APFnp1_rho{p_mesh};
- const DiscreteFunctionP0<TinyVector<Dimension>> APFnp1_rho_u1{p_mesh};
- const DiscreteFunctionP0<TinyVector<Dimension>> APFnp1_rho_u2{p_mesh};
- const DiscreteFunctionP0<TinyVector<Dimension>> APFnp1_rho_E{p_mesh};
-
- parallel_for(
- m_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
- APFnp1_rho[cell_id] = APFnp1[cell_id][0];
- APFnp1_rho_u1[cell_id] = APFnp1[cell_id][1];
- APFnp1_rho_u2[cell_id] = APFnp1[cell_id][2];
- APFnp1_rho_E[cell_id] = APFnp1[cell_id][3];
- });
-
- const DiscreteFunctionP0Vector<const double> MPFnp1_rho = compute_M(rho_np1, APFnp1_rho);
- const DiscreteFunctionP0Vector<const double> MPFnp1_rho_u1 = compute_M(rho_u1_np1, APFnp1_rho_u1);
- const DiscreteFunctionP0Vector<const double> MPFnp1_rho_u2 = compute_M(rho_u2_np1, APFnp1_rho_u2);
- const DiscreteFunctionP0Vector<const double> MPFnp1_rho_E = compute_M(rho_E_np1, APFnp1_rho_E);
+ 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) {
@@ -932,22 +988,22 @@ class EulerKineticSolver2D
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 * MPFnp1_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_u1[cell_id][i_wave] - c2 * deltaLFn_rho_u1[cell_id][i_wave] +
- m_dt * MPFnp1_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_u2[cell_id][i_wave] - c2 * deltaLFn_rho_u2[cell_id][i_wave] +
- m_dt * MPFnp1_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_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 * MPFnp1_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_u1_ex[cell_id][i_wave];
- Fnp1_rho_u2[cell_id][i_wave] = Fn_rho_u2_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];
}
}
@@ -965,30 +1021,23 @@ class EulerKineticSolver2D
const double& gamma,
const std::vector<TinyVector<2>>& lambda_vector,
const size_t& SpaceOrder,
- const DiscreteFunctionP0Vector<const double>& Fn_rho,
- const DiscreteFunctionP0Vector<const double>& Fn_rho_u1,
- const DiscreteFunctionP0Vector<const double>& Fn_rho_u2,
- const DiscreteFunctionP0Vector<const double>& Fn_rho_E,
- const DiscreteFunctionP0<const double>& rho_n_ex,
- const DiscreteFunctionP0<const TinyVector<2>>& rho_u_n_ex,
- const DiscreteFunctionP0<const double>& rho_E_n_ex)
+ const CellArray<const double>& Fn_rho,
+ const CellArray<TinyVector<2>>& 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_u1{Fn_rho_u1},
- m_Fn_rho_u2{Fn_rho_u2},
+ m_Fn_rho_u{Fn_rho_u},
m_Fn_rho_E{Fn_rho_E},
- m_rho_ex{rho_n_ex},
- m_rho_u_ex{rho_u_n_ex},
- m_rho_E_ex{rho_E_n_ex},
m_gamma{gamma},
p_mesh{mesh},
m_mesh{*mesh}
{
- if ((lambda_vector.size() != Fn_rho.size()) or (lambda_vector.size() != Fn_rho_u1.size()) or
- (lambda_vector.size() != Fn_rho_u2.size()) or ((lambda_vector.size() != Fn_rho_E.size()))) {
+ 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");
}
@@ -1003,58 +1052,6 @@ class EulerKineticSolver2D
~EulerKineticSolver2D() = default;
};
-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::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex,
- 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, rho_ex, rho_u_ex, rho_E_ex});
- 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");
- }
-
- 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>> {
- 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>>();
- auto rho_n_ex = rho_ex->get<DiscreteFunctionP0<const double>>();
- auto rho_u_n_ex = rho_u_ex->get<DiscreteFunctionP0<const TinyVector<MeshType::Dimension>>>();
- auto rho_E_n_ex = rho_E_ex->get<DiscreteFunctionP0<const double>>();
-
- EulerKineticSolver2D solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, Fn_rho, Fn_rho_u1, Fn_rho_u2,
- Fn_rho_E, rho_n_ex, rho_u_n_ex, rho_E_n_ex);
-
- return solver.apply(bc_descriptor_list);
- } else {
- throw NotImplementedError("Invalid mesh type for 2D Euler Kinetic solver");
- }
- },
- mesh_v->variant());
-}
-
template <MeshConcept MeshType>
class EulerKineticSolver2D<MeshType>::SymmetryBoundaryCondition
{
@@ -1147,3 +1144,61 @@ class EulerKineticSolver2D<MeshType>::InflowListBoundaryCondition
~InflowListBoundaryCondition() = default;
};
+
+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)
+{
+ const std::shared_ptr<const MeshVariant> mesh_v = getCommonMesh({F_rho, F_rho_u1, F_rho_u2, 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");
+ }
+
+ 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>> {
+ 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];
+ }
+ });
+
+ EulerKineticSolver2D solver(p_mesh, dt, eps, gamma, lambda_vector, SpaceOrder, Fn_rho_array, Fn_rho_u_array,
+ Fn_rho_E_array);
+
+ return solver.apply(bc_descriptor_list);
+ } else {
+ throw NotImplementedError("Invalid mesh type for 2D Euler Kinetic solver");
+ }
+ },
+ mesh_v->variant());
+}
diff --git a/src/scheme/EulerKineticSolver2D.hpp b/src/scheme/EulerKineticSolver2D.hpp
index 76d024c82658527ee2fe895ada184929fabbb6ef..a7937ee7ae067059b671d10b05674e109b7d9b88 100644
--- a/src/scheme/EulerKineticSolver2D.hpp
+++ b/src/scheme/EulerKineticSolver2D.hpp
@@ -38,9 +38,6 @@ eulerKineticSolver2D(const double& dt,
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::shared_ptr<const DiscreteFunctionVariant>& rho_ex,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho_u_ex,
- const std::shared_ptr<const DiscreteFunctionVariant>& rho_E_ex,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
#endif // EULER_KINETIC_SOLVER_2D_HPP