diff --git a/src/scheme/GKSNavier.cpp b/src/scheme/GKSNavier.cpp
index f85ae84ce007453b12f5ffe12565599c51517fc5..48b98b66b1375148e80b011389d195d0bb1c5777 100644
--- a/src/scheme/GKSNavier.cpp
+++ b/src/scheme/GKSNavier.cpp
@@ -31,10 +31,8 @@ gks_inv_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
CellValue<double> local_inv_dt{mesh.connectivity()};
parallel_for(
- mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
- // local_inv_dt[cell_id] = (c[cell_id] + abs(U[cell_id])) / (Vj[cell_id]);
- local_inv_dt[cell_id] = (c[cell_id] + abs(U[cell_id])) / (Vj[cell_id] * Vj[cell_id]);
- });
+ mesh.numberOfCells(),
+ PUGS_LAMBDA(CellId cell_id) { local_inv_dt[cell_id] = (c[cell_id] + abs(U[cell_id])) / (Vj[cell_id]); });
return max(local_inv_dt);
} else {
@@ -75,9 +73,9 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
double p = rho / (2 * lambda);
M(0, 0) = rho;
- M(1, 1) = rhoU[0] * U[0] + p;
+ M(1, 1) = rho * U_2 + p;
M(2, 2) =
- 0.25 * rho * U_2 * U_2 + (delta + 4) * 0.5 * U_2 * p + 0.125 * (delta * delta + 6 * delta + 8) * p / lambda;
+ 0.25 * rho * U_2 * U_2 + 0.5 * (delta + 3) * U_2 * p + 0.125 * (delta * delta + 4 * delta + 3) * p / lambda;
M(1, 0) = rhoU[0];
M(0, 1) = M(1, 0);
@@ -85,11 +83,10 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
M(2, 0) = rhoE;
M(0, 2) = M(2, 0);
- M(2, 1) = 0.5 * rhoU[0] * U_2 + (delta + 4) * 0.5 * U[0] * p;
+ M(2, 1) = 0.5 * (rhoU[0] * U_2 + (delta + 3) * U[0] * p);
M(1, 2) = M(2, 1);
- const TinyMatrix<SIZEproblem> M_bis = M;
- const TinyMatrix<SIZEproblem> inverseM = inverse(M_bis);
+ const TinyMatrix<SIZEproblem> inverseM = inverse(M);
return inverseM;
}
@@ -104,12 +101,11 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
TinyMatrix<SIZEproblem> invMatrixM(zero);
TinyMatrix<SIZEproblem> MatrixP(zero);
- TinyMatrix<SIZEproblem> invMatrixP(zero);
TinyMatrix<SIZEproblem> invMatrixMtilde(zero);
invMatrixMtilde(0, 0) = 1;
invMatrixMtilde(1, 1) = 1;
- invMatrixMtilde(2, 2) = 2 * sqrt2 * lambda / (std::sqrt(delta + 1) * sqrt_rho);
+ invMatrixMtilde(2, 2) = (8 * lambda * lambda) / (rho * (delta + 1));
MatrixP(0, 0) = 1 / sqrt_rho;
MatrixP(1, 1) = sqrt2 * sqrt_lambda / sqrt_rho;
@@ -118,16 +114,22 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
MatrixP(2, 0) = 0.5 * U_2 - 0.25 * (delta + 1) / lambda;
MatrixP(2, 1) = -U[0];
- invMatrixP(0, 0) = sqrt_rho;
- invMatrixP(1, 1) = sqrt_rho / (sqrt2 * sqrt_lambda);
- invMatrixP(2, 2) = 1;
- invMatrixP(1, 0) = sqrt_rho * U[0];
- invMatrixP(2, 0) = sqrt_rho * (0.5 * U_2 + 0.25 * (delta + 1) / lambda);
- invMatrixP(2, 1) = U[0] * sqrt_rho / (sqrt2 * sqrt_lambda);
+ // TinyMatrix<SIZEproblem> MatrixM(zero);
+ // TinyMatrix<SIZEproblem> invMatrixP(zero);
+ // TinyMatrix<SIZEproblem> MatrixMtilde(zero);
- invMatrixM = MatrixP * invMatrixMtilde * invMatrixP;
+ // MatrixMtilde(0, 0) = 1;
+ // MatrixMtilde(1, 1) = 1;
+ // MatrixMtilde(2, 2) = (rho * (delta + 1)) / (8 * lambda * lambda);
- // std::cout << std ::endl << invMatrixP << std ::endl << std ::endl;
+ // invMatrixP(0, 0) = sqrt_rho;
+ // invMatrixP(1, 1) = sqrt_rho / (sqrt2 * sqrt_lambda);
+ // invMatrixP(2, 2) = 1;
+ // invMatrixP(1, 0) = sqrt_rho * U[0];
+ // invMatrixP(2, 0) = sqrt_rho * (0.5 * U_2 + 0.25 * (delta + 1) / lambda);
+ // invMatrixP(2, 1) = U[0] * sqrt_rho / (sqrt2 * sqrt_lambda);
+
+ invMatrixM = transpose(MatrixP) * invMatrixMtilde * (MatrixP);
return invMatrixM;
}
@@ -160,7 +162,7 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
return xi2_moments;
}
- std::tuple<CellValue<const TinyVector<SIZEproblem>>, CellValue<const TinyVector<SIZEproblem>>>
+ CellValue<const TinyVector<SIZEproblem>>
_computeGradConservVariables_for_a(const MeshType& mesh,
const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& rhoU,
@@ -174,8 +176,6 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
CellValue<TinyVector<SIZEproblem>> Gradj{mesh.connectivity()};
Gradj.fill(zero);
- CellValue<TinyVector<SIZEproblem>> Gradjpu{mesh.connectivity()};
- Gradjpu.fill(zero);
NodeValue<double> rho_node_mean{mesh.connectivity()};
NodeValue<Rd> rhoU_node_mean{mesh.connectivity()};
@@ -201,23 +201,11 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
rhoU_node_mean[node_id] *= 0.5;
rhoE_node_mean[node_id] *= 0.5;
- for (size_t l = 0; l < node_cells.size(); l++) {
- if (node_cells.size() == 1) {
- continue;
- }
- if (l == 0) {
- Gradj[node_cells[l]][0] = (rho_node_mean[node_id] - rho[node_cells[l]]) / Vj[node_cells[l]];
- Gradj[node_cells[l]][1] = (rhoU_node_mean[node_id][0] - rhoU[node_cells[l]][0]) / Vj[node_cells[l]];
- Gradj[node_cells[l]][2] = (rhoE_node_mean[node_id] - rhoE[node_cells[l]]) / Vj[node_cells[l]];
- } else {
- Gradjpu[node_cells[l]][0] = (rho[node_cells[l]] - rho_node_mean[node_id]) / Vj[node_cells[l]];
- Gradjpu[node_cells[l]][1] = (rhoU[node_cells[l]][0] - rhoU_node_mean[node_id][0]) / Vj[node_cells[l]];
- Gradjpu[node_cells[l]][2] = (rhoE[node_cells[l]] - rhoE_node_mean[node_id]) / Vj[node_cells[l]];
- ;
- }
- }
+ Gradj[node_cells[0]][0] = (rho_node_mean[node_id] - rho[node_cells[0]]) / Vj[node_cells[0]];
+ Gradj[node_cells[0]][1] = (rhoU_node_mean[node_id][0] - rhoU[node_cells[0]][0]) / Vj[node_cells[0]];
+ Gradj[node_cells[0]][2] = (rhoE_node_mean[node_id] - rhoE[node_cells[0]]) / Vj[node_cells[0]];
});
- return {Gradj, Gradjpu};
+ return Gradj;
}
TinyVector<SIZEproblem>
@@ -247,7 +235,7 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
return B;
}
- std::tuple<CellValue<TinyVector<SIZEproblem>>, CellValue<TinyVector<SIZEproblem>>>
+ CellValue<TinyVector<SIZEproblem>>
_computeGradConservVariables_for_b(const MeshType& mesh,
const DiscreteScalarFunction& rho,
const DiscreteVectorFunction& rhoU,
@@ -264,31 +252,17 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
CellValue<TinyVector<SIZEproblem>> Gradj{mesh.connectivity()};
Gradj.fill(zero);
- CellValue<TinyVector<SIZEproblem>> Gradjpu{mesh.connectivity()};
- Gradjpu.fill(zero);
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_cells = node_to_cell_matrix[node_id];
- for (size_t l = 0; l < node_cells.size(); l++) {
- if (node_cells.size() == 1) {
- continue;
- }
- if (l == 0) {
- Gradj[node_cells[l]][0] = (rho_node[node_id] - rho[node_cells[l]]) / Vj[node_cells[l]];
- Gradj[node_cells[l]][1] = (rhoU_node[node_id][0] - rhoU[node_cells[l]][0]) / Vj[node_cells[l]];
- Gradj[node_cells[l]][2] = (rhoE_node[node_id] - rhoE[node_cells[l]]) / Vj[node_cells[l]];
- } else {
- Gradjpu[node_cells[l]][0] = (rho[node_cells[l]] - rho_node[node_id]) / Vj[node_cells[l]];
- Gradjpu[node_cells[l]][1] = (rhoU[node_cells[l]][0] - rhoU_node[node_id][0]) / Vj[node_cells[l]];
- Gradjpu[node_cells[l]][2] = (rhoE[node_cells[l]] - rhoE_node[node_id]) / Vj[node_cells[l]];
- ;
- }
- }
+ Gradj[node_cells[0]][0] = (rho_node[node_id] - rho[node_cells[0]]) / Vj[node_cells[0]];
+ Gradj[node_cells[0]][1] = (rhoU_node[node_id][0] - rhoU[node_cells[0]][0]) / Vj[node_cells[0]];
+ Gradj[node_cells[0]][2] = (rhoE_node[node_id] - rhoE[node_cells[0]]) / Vj[node_cells[0]];
});
- return {Gradj, Gradjpu};
+ return Gradj;
}
std::tuple<TinyMatrix<SIZEproblem>, TinyMatrix<SIZEproblem>, TinyMatrix<SIZEproblem>>
@@ -444,41 +418,11 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
- double U_2 = U[cell_id][0] * U[cell_id][0];
err_function_term[cell_id] = std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0]);
- exp_term[cell_id] = std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
+ exp_term[cell_id] =
+ std::exp(-lambda[cell_id] * U[cell_id][0] * U[cell_id][0]) / std::sqrt(pi * lambda[cell_id]);
});
- // NodeValue<double> rho_node_mean{mesh.connectivity()};
- // NodeValue<Rd> rhoU_node_mean{mesh.connectivity()};
- // NodeValue<double> rhoE_node_mean{mesh.connectivity()};
- // NodeValue<Rd> U_node_mean{mesh.connectivity()};
- // NodeValue<double> lambda_node_mean{mesh.connectivity()};
-
- // rho_node_mean.fill(1);
- // rhoU_node_mean.fill(zero);
- // rhoE_node_mean.fill(1);
- // U_node_mean.fill(zero);
- // lambda_node_mean.fill(1);
-
- // parallel_for(
- // mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
- // const auto& node_cells = node_to_cell_matrix[node_id];
-
- // for (size_t l = 0; l < node_cells.size(); l++) {
- // if (node_cells.size() == 1)
- // continue;
-
- // rho_node_mean[node_id] += rho[node_cells[l]];
- // rhoU_node_mean[node_id] += rhoU[node_cells[l]];
- // rhoE_node_mean[node_id] += rhoE[node_cells[l]];
- // }
-
- // rho_node_mean[node_id] *= 0.5;
- // rhoU_node_mean[node_id] *= 0.5;
- // rhoE_node_mean[node_id] *= 0.5;
- // });
-
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_cells = node_to_cell_matrix[node_id];
@@ -541,9 +485,8 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
}
}
- auto [Gradj_for_a, Gradjpu_for_a] = this->_computeGradConservVariables_for_a(mesh, rho, rhoU, rhoE);
- auto [Gradj_for_b, Gradjpu_for_b] =
- this->_computeGradConservVariables_for_b(mesh, rho, rhoU, rhoE, rho_node, rhoU_node, rhoE_node);
+ auto Gradj_for_a = this->_computeGradConservVariables_for_a(mesh, rho, rhoU, rhoE);
+ auto Gradj_for_b = this->_computeGradConservVariables_for_b(mesh, rho, rhoU, rhoE, rho_node, rhoU_node, rhoE_node);
NodeValue<TinyVector<SIZEproblem>> Fluxes_term{mesh.connectivity()};
Fluxes_term.fill(zero);
@@ -567,8 +510,6 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
rhoU_right_node[0] = rhoU[node_cells[1]][0] - 0.5 * Vj[node_cells[1]] * Gradj_for_a[node_cells[1]][1];
const double rhoE_right_node = rhoE[node_cells[1]] - 0.5 * Vj[node_cells[1]] * Gradj_for_a[node_cells[1]][2];
- // peut etre des + 0.5* ici (au dessus)
-
Rd U_right_node;
U_right_node[0] = rhoU_right_node[0] / rho_right_node;
const double lambda_right_node =
@@ -580,45 +521,44 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
TinyVector<7> u_moments_right = this->_compute_u_moments(rho_right_node, U_right_node, lambda_right_node);
TinyVector<2> xi2_moments_right = this->_compute_xi2_moments(lambda_right_node, delta);
- // TinyMatrix<SIZEproblem> InvM_left =
- // this->_computeInvMatrixM(rho_left_node, U_left_node, lambda_left_node, delta);
- // TinyMatrix<SIZEproblem> InvM_right =
- // this->_computeInvMatrixM(rho_right_node, U_right_node, lambda_right_node, delta);
- // TinyMatrix<SIZEproblem> InvMInterface =
- // this->_computeInvMatrixM(rho_node[node_id], U_node[node_id], lambda_node[node_id], delta);
-
TinyMatrix<SIZEproblem> InvM_left =
- this->_computeInvM_bis(rho_left_node, rhoU_left_node, rhoE_left_node, U_left_node, lambda_left_node, delta);
- TinyMatrix<SIZEproblem> InvM_right = this->_computeInvM_bis(rho_right_node, rhoU_right_node, rhoE_right_node,
- U_right_node, lambda_right_node, delta);
+ this->_computeInvMatrixM(rho_left_node, U_left_node, lambda_left_node, delta);
+ TinyMatrix<SIZEproblem> InvM_right =
+ this->_computeInvMatrixM(rho_right_node, U_right_node, lambda_right_node, delta);
TinyMatrix<SIZEproblem> InvMInterface =
- this->_computeInvM_bis(rho_node[node_id], rhoU_node[node_id], rhoE_node[node_id], U_node[node_id],
- lambda_node[node_id], delta);
+ this->_computeInvMatrixM(rho_node[node_id], U_node[node_id], lambda_node[node_id], delta);
- // std::cout << std::endl << "InvM : " << InvM << "\t and \t InvM_bis" << InvM_bis << std::endl;
+ // TinyMatrix<SIZEproblem> InvM_left =
+ // this->_computeInvM_bis(rho_left_node, rhoU_left_node, rhoE_left_node, U_left_node, lambda_left_node,
+ // delta);
+ // TinyMatrix<SIZEproblem> InvM_right = this->_computeInvM_bis(rho_right_node, rhoU_right_node, rhoE_right_node,
+ // U_right_node, lambda_right_node, delta);
+ // TinyMatrix<SIZEproblem> InvMInterface =
+ // this->_computeInvM_bis(rho_node[node_id], rhoU_node[node_id], rhoE_node[node_id], U_node[node_id],
+ // lambda_node[node_id], delta);
TinyVector<SIZEproblem> al;
TinyVector<SIZEproblem> ar;
TinyVector<SIZEproblem> bl;
TinyVector<SIZEproblem> br;
- for (size_t l = 0; l < node_cells.size(); l++) {
- if (node_cells.size() == 1) {
- continue;
- }
- if (l == 0) {
- al = this->_compute_a_and_b(InvM_left, Gradj_for_a[node_cells[l]]);
- bl = this->_compute_a_and_b(InvMInterface, Gradj_for_b[node_cells[l]]);
- } else {
- ar = this->_compute_a_and_b(InvM_right, Gradjpu_for_a[node_cells[l]]);
- br = this->_compute_a_and_b(InvMInterface, Gradjpu_for_b[node_cells[l]]);
- }
- }
-
- // TinyMatrix<SIZEproblem> Matrix_for_A_left =
- // this->_computeMatricesC1_full_moments(u_moments_left, xi2_moments_left);
- // TinyMatrix<SIZEproblem> Matrix_for_A_right =
- // this->_computeMatricesC1_full_moments(u_moments_right, xi2_moments_right);
+ al = this->_compute_a_and_b(InvM_left, Gradj_for_a[node_cells[0]]);
+ bl = this->_compute_a_and_b(InvMInterface, Gradj_for_b[node_cells[0]]);
+ ar = this->_compute_a_and_b(InvM_right, Gradj_for_a[node_cells[0]]);
+ br = this->_compute_a_and_b(InvMInterface, Gradj_for_b[node_cells[0]]);
+
+ // for (size_t l = 0; l < node_cells.size(); l++) {
+ // if (node_cells.size() == 1) {
+ // continue;
+ // }
+ // if (l == 0) {
+ // al = this->_compute_a_and_b(InvM_left, Gradj_for_a[node_cells[l]]);
+ // bl = this->_compute_a_and_b(InvMInterface, Gradj_for_b[node_cells[l]]);
+ // } else {
+ // ar = this->_compute_a_and_b(InvM_right, -Gradj_for_a[node_cells[l]]);
+ // br = this->_compute_a_and_b(InvMInterface, -Gradj_for_b[node_cells[l]]);
+ // }
+ // }
auto [C0_semi_pos_left, C1_semi_pos_left, C2_semi_pos_left] =
this->_computeMatricesC0C1C2_semi_moments(rho_left_node, U_left_node, lambda_left_node, u_moments_left,
@@ -668,8 +608,8 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
F_fn[1] = C1_semi_pos_left(1, 0) + C1_semi_neg_right(1, 0);
F_fn[2] = C1_semi_pos_left(2, 0) + C1_semi_neg_right(2, 0);
- Fluxes_term[node_id] = (1 - eta[node_id]) * G + eta[node_id] * F_fn + Fluxes_from_a + Fluxes_from_b +
- Fluxes_from_B + Fluxes_for_Navier;
+ Fluxes_term[node_id] = (1 - eta[node_id]) * G + eta[node_id] * F_fn + Fluxes_for_Navier + Fluxes_from_a +
+ Fluxes_from_b + Fluxes_from_B;
});
NodeValuePerCell<TinyVector<SIZEproblem>> Fluxes_term_apply{mesh.connectivity()};