diff --git a/src/scheme/EulerKineticSolverMultiD.cpp b/src/scheme/EulerKineticSolverMultiD.cpp
index f299da04d71838d25f34446caadb581e2e7f8ce3..20077749b71094bed6ebf343795b264fdd3e42c6 100644
--- a/src/scheme/EulerKineticSolverMultiD.cpp
+++ b/src/scheme/EulerKineticSolverMultiD.cpp
@@ -46,8 +46,8 @@ getLambdaVector2D(const std::shared_ptr<const MeshVariant>& mesh,
for (size_t m = 1; m < 4 * M + 1; ++m) {
size_t i_wave = (l - 1) * 4 * M + m - 1;
double ll = l;
- lambda_vector[i_wave] = TinyVector<2>{(ll / L) * lambda * std::cos((m * pi) / (2 * M) + 0.25 * pi),
- (ll / L) * lambda * std::sin((m * pi) / (2 * M) + 0.25 * pi)};
+ lambda_vector[i_wave] = TinyVector<2>{(ll / L) * lambda * std::cos((m * pi) / (2 * M)), // + 0.2 * pi),
+ (ll / L) * lambda * std::sin((m * pi) / (2 * M))}; // + 0.2 * pi)};
}
}
// lambda_vector[4 * L * M] = TinyVector<2>{0, 0}; //!!!!!!!!!!!!!!!!!!!!
@@ -597,7 +597,7 @@ class EulerKineticSolverMultiD
double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
- if (li_nlr > 0) {
+ if (li_nlr > 1e-14) {
Fr[node_id][i_wave] += li_nlr * Fn[cell_id][i_wave];
sum += li_nlr;
}
@@ -647,7 +647,7 @@ class EulerKineticSolverMultiD
double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
- if (li_nl > 0) {
+ if (li_nl > 1e-14) {
Fl[face_id][i_wave] += Fn[cell_id][i_wave];
}
}
@@ -697,24 +697,20 @@ class EulerKineticSolverMultiD
const FaceId face_id = face_list[i_face];
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];
+ const auto& face_to_cell = face_to_cell_matrix[face_id];
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];
+ const CellId cell_id = face_to_cell[i_cell];
- TinyVector<Dimension> n_face = nl[face_id];
- const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
+ double li_nl = dot(m_lambda_vector[i_wave], n);
- double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
- if (li_nl < 0) {
+ if (li_nl < -1e-14) {
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]) * dot(rho_u[cell_id], rho_u[cell_id]);
+ double rho_e = rho_Ej_prime - 0.5 * (1. / rhoj_prime) * dot(rho_uj, rho_uj);
double p = (m_gamma - 1) * rho_e;
TinyVector<Dimension> A_rho_prime = rho_uj_prime;
@@ -762,7 +758,7 @@ class EulerKineticSolverMultiD
auto txt = I - nxn;
double li_nl = dot(m_lambda_vector[i_wave], n);
- if (li_nl < 0) {
+ if (li_nl < -1e-14) {
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;
@@ -813,7 +809,7 @@ class EulerKineticSolverMultiD
auto n = sign * nl[face_id];
double li_nl = dot(m_lambda_vector[i_wave], n);
- if (li_nl < 0) {
+ if (li_nl < -1e-14) {
Fl_rho[face_id][i_wave] += Fn_rho[cell_id][i_wave];
Fl_rho_u[face_id][i_wave] += Fn_rho_u[cell_id][i_wave];
Fl_rho_E[face_id][i_wave] += Fn_rho_E[cell_id][i_wave];
@@ -891,7 +887,7 @@ class EulerKineticSolverMultiD
double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell));
- if (li_njr > 0) {
+ if (li_njr > 1e-14) {
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;
@@ -902,7 +898,7 @@ class EulerKineticSolverMultiD
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) {
+ if (li_njr_prime > 1e-14) {
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];
@@ -976,7 +972,7 @@ class EulerKineticSolverMultiD
double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell));
- if (li_njr > 0) {
+ if (li_njr > 1e-14) {
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;
@@ -987,7 +983,7 @@ class EulerKineticSolverMultiD
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) {
+ if (li_njr_prime > 1e-14) {
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];
@@ -1060,7 +1056,7 @@ class EulerKineticSolverMultiD
double li_njr = dot(m_lambda_vector[i_wave], njr(cell_id, i_node_cell));
- if (li_njr > 0) {
+ if (li_njr > 1e-14) {
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;
@@ -1071,7 +1067,7 @@ class EulerKineticSolverMultiD
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) {
+ if (li_njr_prime > 1e-14) {
Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_njr_prime;
Fr_rho_u[node_id][i_wave] += li_njr_prime * Fn_rho_u[cell_id][i_wave];
Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_njr_prime;
@@ -1178,7 +1174,7 @@ class EulerKineticSolverMultiD
double li_nlr = dot(m_lambda_vector[i_wave], n_face);
- if (li_nlr > 0) {
+ if (li_nlr > 1e-14) {
Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_nlr;
Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u[cell_id][i_wave];
Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_nlr;
@@ -1189,7 +1185,7 @@ class EulerKineticSolverMultiD
TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
double li_nlr_prime = dot(m_lambda_vector[i_wave], nlr_prime);
- if (li_nlr_prime > 0) {
+ if (li_nlr_prime > 1e-14) {
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];
@@ -1276,7 +1272,7 @@ class EulerKineticSolverMultiD
double li_nlr = dot(m_lambda_vector[i_wave], n_face);
- if (li_nlr > 0) {
+ if (li_nlr > 1e-14) {
Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_nlr;
Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u[cell_id][i_wave];
Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_nlr;
@@ -1287,7 +1283,7 @@ class EulerKineticSolverMultiD
TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
double li_nlr_prime = dot(m_lambda_vector[i_wave], nlr_prime);
- if (li_nlr_prime > 0) {
+ if (li_nlr_prime > 1e-14) {
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];
@@ -1374,7 +1370,7 @@ class EulerKineticSolverMultiD
double li_nlr = dot(m_lambda_vector[i_wave], n_face);
- if (li_nlr > 0) {
+ if (li_nlr > 1e-14) {
Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_nlr;
Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u[cell_id][i_wave];
Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_nlr;
@@ -1385,7 +1381,7 @@ class EulerKineticSolverMultiD
TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
double li_nlr_prime = dot(m_lambda_vector[i_wave], nlr_prime);
- if (li_nlr_prime > 0) {
+ if (li_nlr_prime > 1e-14) {
Fr_rho[node_id][i_wave] += Fn_rho[cell_id][i_wave] * li_nlr_prime;
Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u[cell_id][i_wave];
Fr_rho_E[node_id][i_wave] += Fn_rho_E[cell_id][i_wave] * li_nlr_prime;
@@ -1583,20 +1579,18 @@ class EulerKineticSolverMultiD
DiscreteFunctionP0Vector<TinyVector<Dimension>> Fnp1_rho_u(p_mesh, nb_waves);
DiscreteFunctionP0Vector<double> Fnp1_rho_E(p_mesh, nb_waves);
- // Compute first order F
-
- 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);
+ const CellValue<const double> rho = compute_PFn(Fn_rho);
+ const CellValue<const TinyVector<Dimension>> rho_u = compute_PFn(Fn_rho_u);
+ const CellValue<const double> rho_E = compute_PFn(Fn_rho_E);
CellArray<double> deltaLFn_rho{p_mesh->connectivity(), nb_waves};
CellArray<TinyVector<Dimension>> deltaLFn_rho_u{p_mesh->connectivity(), nb_waves};
CellArray<double> deltaLFn_rho_E{p_mesh->connectivity(), nb_waves};
if (m_scheme == 1) {
- 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);
+ FaceArray<double> Fl_rho = compute_Flux1_face(Fn_rho);
+ FaceArray<TinyVector<Dimension>> Fl_rho_u = compute_Flux1_face(Fn_rho_u);
+ FaceArray<double> Fl_rho_E = compute_Flux1_face(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);
@@ -1604,9 +1598,9 @@ class EulerKineticSolverMultiD
synchronize(Fl_rho_u);
synchronize(Fl_rho_E);
- deltaLFn_rho = compute_deltaLFn_face<double>(Fl_rho);
- deltaLFn_rho_u = compute_deltaLFn_face<TinyVector<Dimension>>(Fl_rho_u);
- deltaLFn_rho_E = compute_deltaLFn_face<double>(Fl_rho_E);
+ deltaLFn_rho = compute_deltaLFn_face(Fl_rho);
+ deltaLFn_rho_u = compute_deltaLFn_face(Fl_rho_u);
+ deltaLFn_rho_E = compute_deltaLFn_face(Fl_rho_E);
} else if (m_scheme == 2) {
NodeArray<double> Fr_rho = compute_Flux1_glace<double>(Fn_rho, is_boundary_node);
@@ -1615,6 +1609,7 @@ class EulerKineticSolverMultiD
NodeArray<double> Fr_rho_E = compute_Flux1_glace<double>(Fn_rho_E, is_boundary_node);
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);
+
synchronize(Fr_rho);
synchronize(Fr_rho_u);
synchronize(Fr_rho_E);
diff --git a/src/scheme/EulerKineticSolverMultiDOrder2.cpp b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
index 1df0243085bd8ec08b0c007b7669662ef4a1d694..acefc2734a3475e14ca6a3269a6e9c22fe795a8f 100644
--- a/src/scheme/EulerKineticSolverMultiDOrder2.cpp
+++ b/src/scheme/EulerKineticSolverMultiDOrder2.cpp
@@ -432,7 +432,7 @@ class EulerKineticSolverMultiDOrder2
double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
- if (li_nlr > 0) {
+ if (li_nlr > 1e-14) {
Flr[node_id][i_face][i_wave] += Fn[cell_id][i_wave] * li_nlr;
sum += li_nlr;
}
@@ -483,7 +483,7 @@ class EulerKineticSolverMultiDOrder2
double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
- if (li_nl > 0) {
+ if (li_nl > 1e-14) {
if ((not is_wall_boundary_face[face_id])) { //} and (not is_tiny_vector_v<T>)) {
Fl[face_id][i_wave] += DPk_Fn(cell_id, i_wave)(m_xl[face_id]);
} else {
@@ -589,7 +589,7 @@ class EulerKineticSolverMultiDOrder2
double li_nlr = sign * dot(m_lambda_vector[i_wave], n_face);
- if (li_nlr > 0) {
+ if (li_nlr > 1e-14) {
Fr[node_id][i_wave] += li_nlr * DPk_Fn(cell_id, i_wave)(m_xl[face_id]);
sum += li_nlr;
}
@@ -669,38 +669,45 @@ class EulerKineticSolverMultiDOrder2
const double sign = face_cell_is_reversed(cell_id, i_face_cell) ? -1 : 1;
double li_nl = sign * dot(m_lambda_vector[i_wave], n_face);
- if (li_nl < 0) {
- double rhol_prime = rhol;
- TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul;
- double rho_El_prime = rho_El;
-
- double rho_el = rho_El - 0.5 * (1. / rhol) * dot(rho_ul, rho_ul);
- double pl = (m_gamma - 1) * rho_el;
-
- TinyVector<Dimension> A_rho_prime = rho_ul_prime;
- TinyMatrix<Dimension> A_rho_u_prime =
- 1. / rhol_prime * tensorProduct(rho_ul_prime, rho_ul_prime) + pl * I;
- TinyVector<Dimension> A_rho_E_prime = (rho_El_prime + pl) / rhol_prime * rho_ul_prime;
-
- double Fn_rho_prime = rhol_prime / nb_waves;
- TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ul_prime;
- double Fn_rho_E_prime = rho_El_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];
+ // if (li_nl < -1e-14) {
+ double rhol_prime = rhol;
+ TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul;
+ double rho_El_prime = rho_El;
+
+ double rho_el = rho_El - 0.5 * (1. / rhol) * dot(rho_ul, rho_ul);
+ double pl = (m_gamma - 1) * rho_el;
+
+ TinyVector<Dimension> A_rho_prime = rho_ul_prime;
+ TinyMatrix<Dimension> A_rho_u_prime =
+ 1. / rhol_prime * tensorProduct(rho_ul_prime, rho_ul_prime) + pl * I;
+ TinyVector<Dimension> A_rho_E_prime = (rho_El_prime + pl) / rhol_prime * rho_ul_prime;
+
+ double Fn_rho_prime = rhol_prime / nb_waves;
+ TinyVector<Dimension> Fn_rho_u_prime = 1. / nb_waves * rho_ul_prime;
+ double Fn_rho_E_prime = rho_El_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];
+ }
+ if (li_nl < -1e-14) {
Fl_rho[face_id][i_wave] += Fn_rho_prime;
Fl_rho_u[face_id][i_wave] += Fn_rho_u_prime;
Fl_rho_E[face_id][i_wave] += Fn_rho_E_prime;
}
+ // if (cell_id >= 500 and cell_id < 502) {
+ // std::cout << "rho'(" << cell_id << ") = " << rhol_prime << "\n";
+ // std::cout << "rho_u'(" << cell_id << ") = " << rho_ul_prime << "\n";
+ // std::cout << "rho_E'(" << cell_id << ") = " << rho_El_prime << "\n";
+ // std::cout << "F_rho_u'(" << i_wave << ") = " << Fn_rho_u_prime << "\n";
+ // }
}
}
}
+ // std::cout << "\n";
} else if constexpr (std::is_same_v<BCType, WallBoundaryCondition>) {
auto face_list = bc.faceList();
@@ -732,7 +739,7 @@ class EulerKineticSolverMultiDOrder2
auto txt = I - nxn;
double li_nl = dot(m_lambda_vector[i_wave], n);
- if (li_nl < 0) {
+ if (li_nl < -1e-14) {
double rhol_prime = rhol;
TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul;
double rho_El_prime = rho_El;
@@ -782,7 +789,7 @@ class EulerKineticSolverMultiDOrder2
auto n = sign * nl[face_id];
double li_nl = dot(m_lambda_vector[i_wave], n);
- if (li_nl < 0) {
+ if (li_nl < -1e-14) {
Fl_rho[face_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]);
Fl_rho_u[face_id][i_wave] += Fn_rho_u(cell_id, i_wave)(m_xl[face_id]);
Fl_rho_E[face_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]);
@@ -1172,7 +1179,7 @@ class EulerKineticSolverMultiDOrder2
double li_nlr = dot(m_lambda_vector[i_wave], n_face);
- if (li_nlr > 0) {
+ if (li_nlr > 1e-14) {
Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]);
Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
@@ -1183,7 +1190,7 @@ class EulerKineticSolverMultiDOrder2
TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
double li_nlr_prime = dot(m_lambda_vector[i_wave], nlr_prime);
- if (li_nlr_prime > 0) {
+ if (li_nlr_prime > 1e-14) {
double rhol_prime = rhol;
TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul;
double rho_El_prime = rho_El;
@@ -1280,7 +1287,7 @@ class EulerKineticSolverMultiDOrder2
double li_nlr = dot(m_lambda_vector[i_wave], n_face);
- if (li_nlr > 0) {
+ if (li_nlr > 1e-14) {
Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]);
Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
@@ -1291,7 +1298,7 @@ class EulerKineticSolverMultiDOrder2
TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
double li_nlr_prime = dot(m_lambda_vector[i_wave], nlr_prime);
- if (li_nlr_prime > 0) {
+ if (li_nlr_prime > 1e-14) {
double rhol_prime = rhol;
TinyVector<Dimension> rho_ul_prime = txt * rho_ul - nxn * rho_ul;
double rho_El_prime = rho_El;
@@ -1378,7 +1385,7 @@ class EulerKineticSolverMultiDOrder2
double li_nlr = dot(m_lambda_vector[i_wave], n_face);
- if (li_nlr > 0) {
+ if (li_nlr > 1e-14) {
Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
Fr_rho_u[node_id][i_wave] += li_nlr * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]);
Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * li_nlr;
@@ -1389,7 +1396,7 @@ class EulerKineticSolverMultiDOrder2
TinyVector<Dimension> nlr_prime = txt * n_face - nxn * n_face;
double li_nlr_prime = dot(m_lambda_vector[i_wave], nlr_prime);
- if (li_nlr_prime > 0) {
+ if (li_nlr_prime > 1e-14) {
Fr_rho[node_id][i_wave] += Fn_rho(cell_id, i_wave)(m_xl[face_id]) * li_nlr_prime;
Fr_rho_u[node_id][i_wave] += li_nlr_prime * Fn_rho_u(cell_id, i_wave)(m_xl[face_id]);
Fr_rho_E[node_id][i_wave] += Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) * li_nlr_prime;
@@ -1735,7 +1742,7 @@ class EulerKineticSolverMultiDOrder2
symmetry_boundary_descriptor_list);
PolynomialReconstructionDescriptor reconstruction_descriptor_vector(IntegrationMethodType::cell_center, 1,
- boundary_descriptor_list);
+ symmetry_boundary_descriptor_list);
auto reconstruction_scalar =
PolynomialReconstruction{reconstruction_descriptor_scalar}.build(DiscreteFunctionP0Vector(p_mesh, Fn_rho),
@@ -1751,31 +1758,34 @@ class EulerKineticSolverMultiDOrder2
auto DPk_Fn_rho_u =
reconstruction_vector[0]->template get<DiscreteFunctionDPkVector<Dimension, const TinyVector<Dimension>>>();
- const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
- const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
-
- for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
- const auto cell_to_face = cell_to_face_matrix[cell_id];
- if ((cell_id >= 500) and (cell_id < 510)) {
- size_t i_face = 2; // for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
-
- FaceId face_id = cell_to_face[i_face];
-
- const double sign = face_cell_is_reversed(cell_id, i_face) ? -1 : 1;
- const TinyVector<Dimension> n = sign * nl[face_id];
-
- // for (size_t i_wave = 0; i_wave < m_lambda_vector.size(); ++i_wave) {
- size_t i_wave = 0;
- std::cout << "n = " << n << ", x = " << m_xj[cell_id] << ", lambda_1 = " << m_lambda_vector[i_wave] << "\n";
- std::cout << "F^rho_" << i_wave << "(" << cell_id << "," << face_id
- << ") = " << DPk_Fn_rho(cell_id, i_wave)(m_xl[face_id]) << "; F^rho_u_" << i_wave << "(" << cell_id
- << "," << face_id << ") = " << DPk_Fn_rho_u(cell_id, i_wave)(m_xl[face_id]) << "; F^rho_E_" << i_wave
- << "(" << cell_id << "," << face_id << ") = " << DPk_Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) << "\n";
- // }
- //}
- }
- }
- exit(0);
+ // const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
+ // // const auto& face_cell_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
+
+ // for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
+ // const auto cell_to_face = cell_to_face_matrix[cell_id];
+ // if ((cell_id >= 500) and (cell_id < 510)) {
+ // size_t i_face = 2; // for (size_t i_face = 0; i_face < cell_to_face.size(); ++i_face) {
+
+ // FaceId face_id = cell_to_face[i_face];
+
+ // // const double sign = face_cell_is_reversed(cell_id, i_face) ? -1 : 1;
+ // // const TinyVector<Dimension> n = sign * nl[face_id];
+
+ // // for (size_t i_wave = 0; i_wave < m_lambda_vector.size(); ++i_wave) {
+ // size_t i_wave = 0;
+ // // std::cout << "n = " << n << ", x = " << m_xj[cell_id] << ", lambda_1 = " << m_lambda_vector[i_wave] <<
+ // "\n"; std::cout << "F^rho_" << i_wave << "(" << cell_id << "," << face_id
+ // << ") = " << DPk_Fn_rho(cell_id, i_wave)(m_xl[face_id]) << "; F^rho_u_" << i_wave << "(" << cell_id
+ // << "," << face_id << ") = " << DPk_Fn_rho_u(cell_id, i_wave)(m_xl[face_id]) << "; F^rho_E_" <<
+ // i_wave
+ // << "(" << cell_id << "," << face_id << ") = " << DPk_Fn_rho_E(cell_id, i_wave)(m_xl[face_id]) <<
+ // "\n";
+ // // }
+ // //}
+ // }
+ // }
+ // std::cout << "\n";
+ // // exit(0);
DiscreteFunctionDPkVector<Dimension, double> Fn_rho_lim = copy(DPk_Fn_rho);
DiscreteFunctionDPkVector<Dimension, double> Fn_rho_E_lim = copy(DPk_Fn_rho_E);
@@ -1787,6 +1797,19 @@ class EulerKineticSolverMultiDOrder2
vector_limiter(Fn_rho_u, Fn_rho_u_lim);
+ // const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
+ // for (CellId cell_id = 500; cell_id < 502; ++cell_id) {
+ // std::cout << "rho(" << cell_id << ") = " << rho[cell_id] << "\n";
+ // std::cout << "rho_u(" << cell_id << ") = " << rho_u[cell_id] << "\n";
+ // std::cout << "rho_E(" << cell_id << ") = " << rho_E[cell_id] << "\n";
+ // const auto cell_to_face = cell_to_face_matrix[cell_id];
+ // FaceId face_id = cell_to_face[0];
+ // for (size_t i_wave = 0; i_wave < nb_waves; ++i_wave) {
+ // std::cout << "F_rho_u(" << i_wave << ") = " << Fn_rho_u_lim(cell_id, i_wave)(m_xl[face_id]) << "\n";
+ // }
+ // }
+ // std::cout << "\n";
+
CellArray<double> deltaLFn_rho{p_mesh->connectivity(), nb_waves};
CellArray<TinyVector<Dimension>> deltaLFn_rho_u{p_mesh->connectivity(), nb_waves};
CellArray<double> deltaLFn_rho_E{p_mesh->connectivity(), nb_waves};
@@ -1866,7 +1889,6 @@ class EulerKineticSolverMultiDOrder2
DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_lim = copy(DPk_Fnp1_rho);
DiscreteFunctionDPkVector<Dimension, double> Fnp1_rho_E_lim = copy(DPk_Fnp1_rho_E);
- // const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
// for (CellId cell_id = 0; cell_id < m_mesh.numberOfCells(); ++cell_id) {
// const auto cell_to_face = cell_to_face_matrix[cell_id];
// if ((cell_id >= 500) and (cell_id < 510)) {