diff --git a/src/scheme/EulerKineticSolverLagrangeMultiD.cpp b/src/scheme/EulerKineticSolverLagrangeMultiD.cpp
index 70423dfc07916d6f20cec371604b3bd53e1b37f4..3e8eb0b8a170b63565ea9a377745e65f6c255fd0 100644
--- a/src/scheme/EulerKineticSolverLagrangeMultiD.cpp
+++ b/src/scheme/EulerKineticSolverLagrangeMultiD.cpp
@@ -251,6 +251,241 @@ class EulerKineticSolverLagrangeMultiD
return Fr;
}
+ void
+ apply_bc(NodeArray<double>& Fr_p,
+ NodeArray<TinyVector<Dimension>>& Fr_u,
+ const CellValue<const double>& p,
+ const CellValue<const TinyVector<Dimension>>& u,
+ const CellArray<const double>& Fn_p,
+ const CellArray<const TinyVector<Dimension>>& Fn_u,
+ const double& lambda,
+ const BoundaryConditionList& bc_list)
+ {
+ const size_t nb_waves = m_lambda_vector.size();
+ const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
+ const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
+ const NodeValuePerCell<const TinyVector<Dimension>> njr = MeshDataManager::instance().getMeshData(m_mesh).njr();
+ const NodeValuePerCell<const TinyVector<Dimension>> Cjr = MeshDataManager::instance().getMeshData(m_mesh).Cjr();
+
+ TinyVector<MeshType::Dimension> inv_S = zero;
+ for (size_t d = 0; d < MeshType::Dimension; ++d) {
+ for (size_t i = 0; i < nb_waves; ++i) {
+ inv_S[d] += m_lambda_vector[i][d] * m_lambda_vector[i][d];
+ }
+ inv_S[d] = 1. / inv_S[d];
+ }
+
+ for (auto&& bc_v : bc_list) {
+ std::visit(
+ [=, this](auto&& bc) {
+ using BCType = std::decay_t<decltype(bc)>;
+ if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+ auto 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] += Cjr(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_p[node_id][i_wave] = 0;
+ Fr_u[node_id][i_wave] = zero;
+
+ 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 > 1e-14) {
+ Fr_p[node_id][i_wave] += Fn_p[cell_id][i_wave] * li_njr;
+ Fr_u[node_id][i_wave] += li_njr * Fn_u[cell_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 > 1e-14) {
+ double p_prime = p[cell_id];
+ TinyVector<Dimension> u_prime = txt * u[cell_id] - nxn * u[cell_id];
+
+ TinyVector<Dimension> A_p_prime = lambda * lambda * u_prime;
+ TinyMatrix<Dimension> A_u_prime = p_prime * I;
+
+ double Fn_p_prime = p_prime / nb_waves;
+ TinyVector<Dimension> Fn_u_prime = 1. / nb_waves * u_prime;
+
+ for (size_t d1 = 0; d1 < Dimension; ++d1) {
+ Fn_p_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_p_prime[d1];
+ for (size_t d2 = 0; d2 < Dimension; ++d2) {
+ Fn_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_u_prime(d2, d1);
+ }
+ }
+
+ Fr_p[node_id][i_wave] += Fn_p_prime * li_njr_prime;
+ Fr_u[node_id][i_wave] += li_njr_prime * Fn_u_prime;
+
+ sum += li_njr_prime;
+ }
+ }
+ if (sum != 0) {
+ Fr_p[node_id][i_wave] /= sum;
+ Fr_u[node_id][i_wave] = 1. / sum * Fr_u[node_id][i_wave];
+ }
+ }
+ }
+ } 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] += Cjr(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_p[node_id][i_wave] = 0;
+ Fr_u[node_id][i_wave] = zero;
+
+ 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 > 1e-14) {
+ Fr_p[node_id][i_wave] += Fn_p[cell_id][i_wave] * li_njr;
+ Fr_u[node_id][i_wave] += li_njr * Fn_u[cell_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 > 1e-14) {
+ double p_prime = p[cell_id];
+ TinyVector<Dimension> u_prime = txt * u[cell_id] - nxn * u[cell_id];
+
+ TinyVector<Dimension> A_p_prime = lambda * lambda * u_prime;
+ TinyMatrix<Dimension> A_u_prime = p_prime * I;
+
+ double Fn_p_prime = p_prime / nb_waves;
+ TinyVector<Dimension> Fn_u_prime = 1. / nb_waves * u_prime;
+
+ for (size_t d1 = 0; d1 < Dimension; ++d1) {
+ Fn_p_prime += inv_S[d1] * m_lambda_vector[i_wave][d1] * A_p_prime[d1];
+ for (size_t d2 = 0; d2 < Dimension; ++d2) {
+ Fn_u_prime[d1] += inv_S[d2] * m_lambda_vector[i_wave][d2] * A_u_prime(d2, d1);
+ }
+ }
+
+ Fr_p[node_id][i_wave] += Fn_p_prime * li_njr_prime;
+ Fr_u[node_id][i_wave] += li_njr_prime * Fn_u_prime;
+
+ sum += li_njr_prime;
+ }
+ }
+ if (sum != 0) {
+ Fr_p[node_id][i_wave] /= sum;
+ Fr_u[node_id][i_wave] = 1. / sum * Fr_u[node_id][i_wave];
+ }
+ }
+ }
+ } else if constexpr (std::is_same_v<BCType, OutflowBoundaryCondition>) {
+ 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] += Cjr(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_p[node_id][i_wave] = 0;
+ Fr_u[node_id][i_wave] = zero;
+
+ 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 > 1e-14) {
+ Fr_p[node_id][i_wave] += Fn_p[cell_id][i_wave] * li_njr;
+ Fr_u[node_id][i_wave] += li_njr * Fn_u[cell_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 > 1e-14) {
+ Fr_p[node_id][i_wave] += Fn_p[cell_id][i_wave] * li_njr_prime;
+ Fr_u[node_id][i_wave] += li_njr_prime * Fn_u[cell_id][i_wave];
+
+ sum += li_njr_prime;
+ }
+ }
+ if (sum != 0) {
+ Fr_p[node_id][i_wave] /= sum;
+ Fr_u[node_id][i_wave] = 1. / sum * Fr_u[node_id][i_wave];
+ }
+ }
+ }
+ }
+ },
+ bc_v);
+ }
+ }
+
const NodeValue<TinyVector<Dimension>>
compute_velocity_flux(NodeArray<double>& F, const double& lambda)
{
@@ -409,6 +644,8 @@ class EulerKineticSolverLagrangeMultiD
NodeArray<double> Fr_p = compute_Flux<double>(F_p, is_boundary_node);
NodeArray<TinyVector<Dimension>> Fr_u = compute_Flux<TinyVector<Dimension>>(F_u, is_boundary_node);
+ apply_bc(Fr_p, Fr_u, p, u, F_p, F_u, lambda, bc_list);
+
const NodeValue<TinyVector<Dimension>> ur = compute_velocity_flux(Fr_p, lambda);
const NodeValue<TinyVector<Dimension>> pr = compute_pressure_flux(Fr_u);