diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
index 4022bfda3422af3d18f9281b254d108b9ede1a51..9aeb37209da9d5f340c1f19f4245eaa7aaf085c5 100644
--- a/src/scheme/EulerKineticSolverMultiD.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -523,23 +523,25 @@ class EulerKineticSolverMultiD
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 (not is_symmetry_boundary_node[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) {
+ Fr[node_id][i_wave] = 1. / sum_li_njr * Fr[node_id][i_wave];
}
- }
- 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];
}
}
});
@@ -829,6 +831,9 @@ class EulerKineticSolverMultiD
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();
@@ -844,95 +849,178 @@ class EulerKineticSolverMultiD
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;
+ auto node_list = bc.nodeList();
- // 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];
+ TinyMatrix<Dimension> I = identity;
- // 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;
- // }
+ NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()};
+ nr.fill(zero);
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];
+ const NodeId node_id = node_list[i_node];
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(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];
- 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];
+ nr[node_id] += njr(cell_id, i_node_cell);
+ }
+ nr[node_id] = 1. / sqrt(dot(nr[node_id], nr[node_id])) * nr[node_id];
+ auto nxn = tensorProduct(nr[node_id], nr[node_id]);
+ auto txt = I - nxn;
- double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell));
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ double sum = 0;
+ Fr_rho[node_id][i_wave] = 0;
+ Fr_rho_u[node_id][i_wave] = zero;
+ Fr_rho_E[node_id][i_wave] = 0;
- if (li_njr > 0) {
- sum_li_njr += li_njr;
- }
+ 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];
- 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);
+ double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell));
- 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];
+ if (li_njr > 0) {
+ Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_njr;
+ Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u[cell_id][i_wave];
+ Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_njr;
- 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;
+ sum += li_njr;
+ }
- 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;
+ 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);
- 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;
+ 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];
- 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];
- }
+ 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;
- 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;
+ 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;
- sum_li_njr += li_njr_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] += li_njr_prime * Fn_rho_u_prime;
+ Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_njr_prime;
+
+ sum += li_njr_prime;
+ }
+ }
+ if (sum != 0) {
+ Fr_rho[node_id][i_wave] /= sum;
+ Fr_rho_u[node_id][i_wave] = 1. / sum * Fr_rho_u[node_id][i_wave];
+ Fr_rho_E[node_id][i_wave] /= sum;
+ }
+ }
+ }
+ } else if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
+ auto node_list = bc.nodeList();
+
+ TinyMatrix<Dimension> I = identity;
+
+ NodeValue<TinyVector<Dimension>> nr{p_mesh->connectivity()};
+ nr.fill(zero);
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ const NodeId node_id = node_list[i_node];
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cell = node_local_numbers_in_their_cells.itemArray(node_id);
+
+ 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];
+
+ nr[node_id] += njr(cell_id, i_node_cell);
+ }
+ nr[node_id] = 1. / sqrt(dot(nr[node_id], nr[node_id])) * nr[node_id];
+ auto nxn = tensorProduct(nr[node_id], nr[node_id]);
+ auto txt = I - nxn;
+
+ for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ double sum = 0;
+ Fr_rho[node_id][i_wave] = 0;
+ Fr_rho_u[node_id][i_wave] = zero;
+ Fr_rho_E[node_id][i_wave] = 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) {
+ Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_njr;
+ Fr_rho_u[node_id][i_wave] += li_njr * Fn_rho_u[cell_id][i_wave];
+ Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_njr;
+
+ sum += li_njr;
}
- 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;
+
+ 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] += li_njr_prime * Fn_rho_u_prime;
+ Fr_rho_E[node_id][i_wave] += Fn_rho_E_prime * li_njr_prime;
+
+ sum += li_njr_prime;
}
}
- node_is_done[node_id] = true;
+ if (sum != 0) {
+ Fr_rho[node_id][i_wave] /= sum;
+ Fr_rho_u[node_id][i_wave] = 1. / sum * Fr_rho_u[node_id][i_wave];
+ Fr_rho_E[node_id][i_wave] /= sum;
+ }
}
}
}
@@ -1323,35 +1411,35 @@ class EulerKineticSolverMultiD
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);
+ 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);
+ // 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);
+ 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);
+ synchronize(Fr_rho);
+ synchronize(Fr_rho_u);
+ synchronize(Fr_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_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_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 =