diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.cpp b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
index 02c995fdd75728747366ba4c1afb2563af28340d..8651c43999c61cb380933657cb1db271a9055ef1 100644
--- a/src/scheme/EulerKineticSolverMultiDOrder2.cpp
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
@@ -452,7 +452,9 @@ class EulerKineticSolverMultiDOrder2
template <typename T>
FaceArray<T>
- compute_Flux_face(const DiscreteFunctionDPkVector<Dimension, const T> DPk_Fn)
+ compute_Flux_face(const DiscreteFunctionDPkVector<Dimension, const T>& DPk_Fn,
+ const CellArray<const T>& Fn,
+ CellValue<bool>& is_wall_boundary_cell)
{
const size_t nb_waves = m_lambda_vector.size();
FaceArray<T> Fl(m_mesh.connectivity(), nb_waves);
@@ -482,7 +484,11 @@ class EulerKineticSolverMultiDOrder2
double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
if (li_nl > 0) {
- Fl[face_id][i_wave] += DPk_Fn(cell_id, i_wave)(m_xl[face_id]);
+ if (not is_wall_boundary_cell[cell_id]) {
+ Fl[face_id][i_wave] += DPk_Fn(cell_id, i_wave)(m_xl[face_id]);
+ } else {
+ Fl[face_id][i_wave] += Fn[cell_id][i_wave];
+ }
}
}
}
@@ -609,6 +615,9 @@ class EulerKineticSolverMultiDOrder2
DiscreteFunctionDPkVector<Dimension, double> Fn_rho,
DiscreteFunctionDPkVector<Dimension, TinyVector<Dimension>> Fn_rho_u,
DiscreteFunctionDPkVector<Dimension, double> Fn_rho_E,
+ const CellArray<const double> Fj_rho,
+ const CellArray<const TinyVector<Dimension>> Fj_rho_u,
+ const CellArray<const double> Fj_rho_E,
const BoundaryConditionList& bc_list)
{
const size_t nb_waves = m_lambda_vector.size();
@@ -711,9 +720,9 @@ class EulerKineticSolverMultiDOrder2
double rho_El = 0;
for (size_t k_wave = 0; k_wave < nb_waves; ++k_wave) {
- rhol += Fn_rho(cell_id, k_wave)(m_xl[face_id]);
- rho_ul += Fn_rho_u(cell_id, k_wave)(m_xl[face_id]);
- rho_El += Fn_rho_E(cell_id, k_wave)(m_xl[face_id]);
+ rhol += Fj_rho[cell_id][k_wave];
+ rho_ul += Fj_rho_u[cell_id][k_wave];
+ rho_El += Fj_rho_E[cell_id][k_wave];
}
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
@@ -1667,6 +1676,31 @@ class EulerKineticSolverMultiDOrder2
bc_v);
}
+ CellValue<bool> is_wall_boundary_cell{p_mesh->connectivity()};
+ is_wall_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, WallBoundaryCondition>) {
+ auto node_list = bc.nodeList();
+
+ is_wall_boundary_cell.fill(0);
+
+ auto node_to_cell_matrix = m_mesh.connectivity().nodeToCellMatrix();
+
+ 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_wall_boundary_cell[cell_id] = 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);
@@ -1728,11 +1762,13 @@ class EulerKineticSolverMultiDOrder2
CellArray<double> deltaLFn_rho_E{p_mesh->connectivity(), nb_waves};
if (m_scheme == 1) {
- FaceArray<double> Fl_rho = compute_Flux_face<double>(Fn_rho_lim);
- FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux_face<TinyVector<Dimension>>(Fn_rho_u_lim);
- FaceArray<double> Fl_rho_E = compute_Flux_face<double>(Fn_rho_E_lim);
+ FaceArray<double> Fl_rho = compute_Flux_face<double>(Fn_rho_lim, Fn_rho, is_wall_boundary_cell);
+ FaceArray<TinyVector<Dimension>> Fl_rho_u =
+ compute_Flux_face<TinyVector<Dimension>>(Fn_rho_u_lim, Fn_rho_u, is_wall_boundary_cell);
+ FaceArray<double> Fl_rho_E = compute_Flux_face<double>(Fn_rho_E_lim, Fn_rho_E, is_wall_boundary_cell);
- apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, bc_list);
+ apply_face_bc(Fl_rho, Fl_rho_u, Fl_rho_E, Fn_rho_lim, Fn_rho_u_lim, Fn_rho_E_lim, Fn_rho, Fn_rho_u, Fn_rho_E,
+ bc_list);
synchronize(Fl_rho);
synchronize(Fl_rho_u);
@@ -1808,11 +1844,13 @@ class EulerKineticSolverMultiDOrder2
CellArray<double> deltaLFnp1_rho_E{p_mesh->connectivity(), nb_waves};
if (m_scheme == 1) {
- FaceArray<double> Fl_rho_np1 = compute_Flux_face<double>(Fnp1_rho_lim);
- FaceArray<TinyVector<Dimension>> Fl_rho_u_np1 = compute_Flux_face<TinyVector<Dimension>>(Fnp1_rho_u_lim);
- FaceArray<double> Fl_rho_E_np1 = compute_Flux_face<double>(Fnp1_rho_E_lim);
+ FaceArray<double> Fl_rho_np1 = compute_Flux_face<double>(Fnp1_rho_lim, Fn_rho, is_wall_boundary_cell);
+ FaceArray<TinyVector<Dimension>> Fl_rho_u_np1 =
+ compute_Flux_face<TinyVector<Dimension>>(Fnp1_rho_u_lim, Fn_rho_u, is_wall_boundary_cell);
+ FaceArray<double> Fl_rho_E_np1 = compute_Flux_face<double>(Fnp1_rho_E_lim, Fn_rho_E, is_wall_boundary_cell);
- apply_face_bc(Fl_rho_np1, Fl_rho_u_np1, Fl_rho_E_np1, Fnp1_rho_lim, Fnp1_rho_u_lim, Fnp1_rho_E_lim, bc_list);
+ apply_face_bc(Fl_rho_np1, Fl_rho_u_np1, Fl_rho_E_np1, Fnp1_rho_lim, Fnp1_rho_u_lim, Fnp1_rho_E_lim, Fn_rho,
+ Fn_rho_u, Fn_rho_E, bc_list);
synchronize(Fl_rho_np1);
synchronize(Fl_rho_u_np1);