From af488bd097d543ae2a37327847544fe1bdf2e4a9 Mon Sep 17 00:00:00 2001
From: Clovis <clovis.schoeck@etudiant.univ-rennes.fr>
Date: Fri, 20 Sep 2024 16:52:44 +0200
Subject: [PATCH] Navier-Stokes GKS solver programmed
---
src/language/modules/SchemeModule.cpp | 6 +-
src/scheme/GKS2.cpp | 76 ++-
src/scheme/GKSNavier.cpp | 697 +++++++++++++++++---------
src/scheme/GKSNavier.hpp | 11 +-
4 files changed, 499 insertions(+), 291 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 342a440db..1361721b1 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -459,11 +459,11 @@ SchemeModule::SchemeModule()
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoU,
- const std::shared_ptr<const DiscreteFunctionVariant>& rhoE,
- const std::shared_ptr<const DiscreteFunctionVariant>& tau, const double& delta, const double& dt)
+ const std::shared_ptr<const DiscreteFunctionVariant>& rhoE, const double& viscosity, const double& delta,
+ const double& dt)
-> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
- return GKSHandler{getCommonMesh({rho, rhoU, rhoE})}.solver().gksNavier(rho, rhoU, rhoE, tau, delta, dt);
+ return GKSHandler{getCommonMesh({rho, rhoU, rhoE})}.solver().gksNavier(rho, rhoU, rhoE, viscosity, delta, dt);
}
));
diff --git a/src/scheme/GKS2.cpp b/src/scheme/GKS2.cpp
index 6ba5583e0..cd5a89ff5 100644
--- a/src/scheme/GKS2.cpp
+++ b/src/scheme/GKS2.cpp
@@ -33,7 +33,7 @@ class GKS2
if (tau_n[cell_id] == 0)
eta[cell_id] = 0;
else
- eta[cell_id] = tau_n[cell_id]; //(tau_n[cell_id] / dt) * (1 - std::exp(-dt / tau_n[cell_id]));
+ eta[cell_id] = (tau_n[cell_id] / dt) * (1 - std::exp(-dt / tau_n[cell_id]));
}
// std::cout << "eta = " << eta << std::endl;
@@ -83,6 +83,10 @@ class GKS2
CellValue<Rd> U{p_mesh->connectivity()};
// U.fill(Rd(0));
+ CellValue<double> erf_term_pos(mesh.connectivity());
+ CellValue<double> erf_term_neg(mesh.connectivity());
+ CellValue<double> exp_term(mesh.connectivity());
+
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
U[cell_id][0] = rho_U_n[cell_id][0] / rho_n[cell_id];
@@ -90,6 +94,14 @@ class GKS2
lambda[cell_id] = 0.5 * (1. + delta) * rho_n[cell_id] / (2 * rho_E_n[cell_id] - rho_U_2);
});
+ parallel_for(
+ mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+ double U_2 = U[cell_id][0] * U[cell_id][0];
+ erf_term_pos[cell_id] = 1. + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0]);
+ erf_term_neg[cell_id] = 1. - 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]);
+ });
+
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_cells = node_to_cell_matrix[node_id];
@@ -106,35 +118,22 @@ class GKS2
continue;
if (l == 0) {
- double U_2_left = U[node_cells[l]][0] * U[node_cells[l]][0];
-
- rho_cell_left =
- rho_n[node_cells[l]] * (1 + std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0]));
+ rho_cell_left = rho_n[node_cells[l]] * erf_term_pos[node_cells[l]];
rho_U_cell_left[0] =
- rho_U_n[node_cells[l]][0] * (1. + std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) +
- rho_n[node_cells[l]] * std::exp(-lambda[node_cells[l]] * U_2_left) /
- std::sqrt(pi * lambda[node_cells[l]]);
-
- rho_E_cell_left =
- rho_E_n[node_cells[l]] * (1. + std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) +
- 0.5 * rho_U_n[node_cells[l]][0] * std::exp(-lambda[node_cells[l]] * U_2_left) /
- std::sqrt(pi * lambda[node_cells[l]]);
- } else {
- double U_2_right = U[node_cells[l]][0] * U[node_cells[l]][0];
+ rho_U_n[node_cells[l]][0] * erf_term_pos[node_cells[l]] + rho_n[node_cells[l]] * exp_term[node_cells[l]];
+ rho_E_cell_left = rho_E_n[node_cells[l]] * erf_term_pos[node_cells[l]] +
+ 0.5 * rho_U_n[node_cells[l]][0] * exp_term[node_cells[l]];
+ } else {
rho_cell_right =
- rho_n[node_cells[l]] * (1 - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0]));
+ rho_n[node_cells[l]] * (1. - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0]));
rho_U_cell_right[0] =
- rho_U_n[node_cells[l]][0] * (1. - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) -
- rho_n[node_cells[l]] * std::exp(-lambda[node_cells[l]] * U_2_right) /
- std::sqrt(pi * lambda[node_cells[l]]);
-
- rho_E_cell_right =
- rho_E_n[node_cells[l]] * (1. - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) -
- 0.5 * rho_U_n[node_cells[l]][0] * std::exp(-lambda[node_cells[l]] * U_2_right) /
- std::sqrt(pi * lambda[node_cells[l]]);
+ rho_U_n[node_cells[l]][0] * erf_term_neg[node_cells[l]] - rho_n[node_cells[l]] * exp_term[node_cells[l]];
+
+ rho_E_cell_right = rho_E_n[node_cells[l]] * erf_term_neg[node_cells[l]] -
+ 0.5 * rho_U_n[node_cells[l]][0] * exp_term[node_cells[l]];
}
}
rho_node[node_id] = 0.5 * (rho_cell_left + rho_cell_right);
@@ -159,28 +158,26 @@ class GKS2
double U_2_left = U[node_cells[l]][0] * U[node_cells[l]][0];
double rho_U_2_left = rho_U_n[node_cells[l]][0] * U[node_cells[l]][0];
- F2_fn_left[0] = (rho_U_2_left + 0.5 * rho_n[node_cells[l]] / lambda[node_cells[l]]) *
- (1. + std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) +
- rho_U_n[node_cells[l]][0] * std::exp(-lambda[node_cells[l]] * U_2_left) /
- std::sqrt(pi * lambda[node_cells[l]]);
+ F2_fn_left[0] =
+ (rho_U_2_left + 0.5 * rho_n[node_cells[l]] / lambda[node_cells[l]]) * erf_term_pos[node_cells[l]] +
+ rho_U_n[node_cells[l]][0] * exp_term[node_cells[l]];
F3_fn_left[0] = 0.5 * rho_U_n[node_cells[l]][0] * (U_2_left + 0.5 * (delta + 3) / lambda[node_cells[l]]) *
- (1. + std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) +
+ erf_term_pos[node_cells[l]] +
0.5 * rho_n[node_cells[l]] * (U_2_left + 0.5 * (delta + 2) / lambda[node_cells[l]]) *
- std::exp(-lambda[node_cells[l]] * U_2_left) / std::sqrt(pi * lambda[node_cells[l]]);
+ exp_term[node_cells[l]];
} else {
double U_2_right = U[node_cells[l]][0] * U[node_cells[l]][0];
double rho_U_2_right = rho_U_n[node_cells[l]][0] * U[node_cells[l]][0];
- F2_fn_right[0] = (rho_U_2_right + 0.5 * rho_n[node_cells[l]] / lambda[node_cells[l]]) *
- (1. - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) -
- rho_U_n[node_cells[l]][0] * std::exp(-lambda[node_cells[l]] * U_2_right) /
- std::sqrt(pi * lambda[node_cells[l]]);
+ F2_fn_right[0] =
+ (rho_U_2_right + 0.5 * rho_n[node_cells[l]] / lambda[node_cells[l]]) * erf_term_neg[node_cells[l]] -
+ rho_U_n[node_cells[l]][0] * exp_term[node_cells[l]];
F3_fn_right[0] = 0.5 * rho_U_n[node_cells[l]][0] * (U_2_right + 0.5 * (delta + 3) / lambda[node_cells[l]]) *
- (1. - std::erf(std::sqrt(lambda[node_cells[l]]) * U[node_cells[l]][0])) -
+ erf_term_neg[node_cells[l]] -
0.5 * rho_n[node_cells[l]] * (U_2_right + 0.5 * (delta + 2) / lambda[node_cells[l]]) *
- std::exp(-lambda[node_cells[l]] * U_2_right) / std::sqrt(pi * lambda[node_cells[l]]);
+ exp_term[node_cells[l]];
}
}
@@ -215,13 +212,6 @@ class GKS2
const size_t& nb_nodes = rho_flux_G.numberOfSubValues(cell_id);
for (size_t r = 0; r < nb_nodes; r++) {
rho[cell_id] -= dt / Vj[cell_id] * rho_flux_G(cell_id, r)[0];
- // rho_U[cell_id] -= dt / Vj[cell_id] *
- // ((1 - eta[cell_id]) * rho_U_flux_G(cell_id, r) + eta[cell_id] *
- // (rho_U_flux_Ffn(cell_id,
- // r)));
- // rho_E[cell_id] -=
- // dt / Vj[cell_id] *
- // ((1 - eta[cell_id]) * rho_E_flux_G(cell_id, r)[0] + eta[cell_id] * (rho_E_flux_Ffn(cell_id, r)[0]));
rho_U[cell_id] -=
dt / Vj[cell_id] *
(rho_U_flux_G(cell_id, r) + eta[cell_id] * (rho_U_flux_Ffn(cell_id, r) - rho_U_flux_G(cell_id, r)));
diff --git a/src/scheme/GKSNavier.cpp b/src/scheme/GKSNavier.cpp
index 29c4d8cf0..f85ae84ce 100644
--- a/src/scheme/GKSNavier.cpp
+++ b/src/scheme/GKSNavier.cpp
@@ -31,8 +31,10 @@ 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]; });
+ 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]);
+ });
return max(local_inv_dt);
} else {
@@ -46,7 +48,8 @@ template <MeshConcept MeshType>
class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
{
private:
- static constexpr size_t Dimension = MeshType::Dimension;
+ static constexpr size_t Dimension = MeshType::Dimension;
+ static constexpr size_t SIZEproblem = 2 + Dimension;
using Rd = TinyVector<Dimension>;
using Rdxd = TinyMatrix<Dimension>;
@@ -59,67 +62,78 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
const double pi = std::acos(-1);
const double sqrt2 = std::sqrt(2);
- TinyMatrix<2 + Dimension>
- _computeInvMatrixM(const double& rho, const double& U, const double& lambda, const double& delta) const
+ TinyMatrix<SIZEproblem>
+ _computeInvM_bis(const double& rho,
+ const Rd& rhoU,
+ const double& rhoE,
+ const Rd& U,
+ const double& lambda,
+ const double& delta) const
{
- const double sqrt_rho = std::sqrt(rho);
- const double sqrt_lambda = std::sqrt(lambda);
- const double U_2 = U * U;
+ TinyMatrix<SIZEproblem> M;
+ double U_2 = U[0] * U[0];
+ double p = rho / (2 * lambda);
- TinyMatrix<2 + Dimension> invMatrixM(zero);
- TinyMatrix<2 + Dimension> MatrixP(zero);
- TinyMatrix<2 + Dimension> transposeMatrixP(zero);
- TinyMatrix<2 + Dimension> invMatrixMtilde(zero);
+ M(0, 0) = rho;
+ M(1, 1) = rhoU[0] * U[0] + 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;
- invMatrixMtilde(1, 1) = 1;
- invMatrixMtilde(2, 2) = 1;
- invMatrixMtilde(3, 3) = 2 * sqrt2 * lambda / (std::sqrt(delta + 1) * sqrt_rho);
+ M(1, 0) = rhoU[0];
+ M(0, 1) = M(1, 0);
- MatrixP(1, 1) = 1 / sqrt_rho;
- MatrixP(2, 2) = sqrt2 * sqrt_lambda / sqrt_rho;
- MatrixP(3, 3) = 1;
- MatrixP(2, 1) = -sqrt2 * sqrt_lambda * U / sqrt_rho;
- MatrixP(3, 1) = 0.5 * U_2 - 0.25 * (delta + 1) / lambda;
- MatrixP(3, 2) = -U;
+ M(2, 0) = rhoE;
+ M(0, 2) = M(2, 0);
- invMatrixM = MatrixP * invMatrixMtilde * transpose(MatrixP);
+ M(2, 1) = 0.5 * rhoU[0] * U_2 + (delta + 4) * 0.5 * U[0] * p;
+ M(1, 2) = M(2, 1);
- return invMatrixM;
+ const TinyMatrix<SIZEproblem> M_bis = M;
+ const TinyMatrix<SIZEproblem> inverseM = inverse(M_bis);
+
+ return inverseM;
}
- TinyVector<2 + Dimension>
- _compute_al(const TinyMatrix<2 + Dimension>& invMatrixM,
- const DiscreteScalarFunction& rho,
- const DiscreteVectorFunction& rhoU,
- const DiscreteScalarFunction& rhoE)
- {}
+ TinyMatrix<SIZEproblem>
+ _computeInvMatrixM(const double& rho, const Rd& U, const double& lambda, const double& delta) const
+ {
+ const double sqrt_rho = std::sqrt(rho);
+
+ const double sqrt_lambda = std::sqrt(lambda);
+ const double U_2 = U[0] * U[0];
+
+ TinyMatrix<SIZEproblem> invMatrixM(zero);
+ TinyMatrix<SIZEproblem> MatrixP(zero);
+ TinyMatrix<SIZEproblem> invMatrixP(zero);
+ TinyMatrix<SIZEproblem> invMatrixMtilde(zero);
- TinyVector<2 + Dimension>
- _compute_ar(const TinyMatrix<2 + Dimension>& invMatrixM)
- {}
+ invMatrixMtilde(0, 0) = 1;
+ invMatrixMtilde(1, 1) = 1;
+ invMatrixMtilde(2, 2) = 2 * sqrt2 * lambda / (std::sqrt(delta + 1) * sqrt_rho);
- TinyVector<2 + Dimension>
- _compute_Al(const TinyMatrix<2 + Dimension>& invMatrixM)
- {}
+ MatrixP(0, 0) = 1 / sqrt_rho;
+ MatrixP(1, 1) = sqrt2 * sqrt_lambda / sqrt_rho;
+ MatrixP(2, 2) = 1;
+ MatrixP(1, 0) = -sqrt2 * sqrt_lambda * U[0] / sqrt_rho;
+ MatrixP(2, 0) = 0.5 * U_2 - 0.25 * (delta + 1) / lambda;
+ MatrixP(2, 1) = -U[0];
- TinyVector<2 + Dimension>
- _compute_Ar(const TinyMatrix<2 + Dimension>& invMatrixM)
- {}
+ 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);
- TinyVector<2 + Dimension>
- _compute_bl(const TinyMatrix<2 + Dimension>& invMatrixM)
- {}
+ invMatrixM = MatrixP * invMatrixMtilde * invMatrixP;
- TinyVector<2 + Dimension>
- _compute_br(const TinyMatrix<2 + Dimension>& invMatrixM)
- {}
+ // std::cout << std ::endl << invMatrixP << std ::endl << std ::endl;
- TinyVector<2 + Dimension>
- _compute_B(const TinyMatrix<2 + Dimension>& invMatrixM)
- {}
+ return invMatrixM;
+ }
TinyVector<7>
- _compute_u_moments(const double& rho, const Rd& U, const double& lambda)
+ _compute_u_moments(const double& rho, const Rd& U, const double& lambda) const
{
const double U_2 = U[0] * U[0];
const double U_4 = U_2 * U_2;
@@ -133,11 +147,11 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
u_moments[4] = rho * (U_4 + 3 * U_2 / lambda + 0.75 / lambda_2);
u_moments[5] = rho * U[0] * (U_4 + 5 * U_2 / lambda + 3.75 / lambda_2);
u_moments[6] = rho * (U_4 * U_2 + 7.5 * U_4 / lambda + 11.25 * U_2 / lambda_2 + 1.875 / (lambda_2 * lambda));
-
return u_moments;
}
+
TinyVector<2>
- _compute_xi2_moments(const double& lambda, const double& delta)
+ _compute_xi2_moments(const double& lambda, const double& delta) const
{
TinyVector<2> xi2_moments;
xi2_moments[0] = delta * 0.5 / lambda;
@@ -146,126 +160,222 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
return xi2_moments;
}
- TinyMatrix<2 + Dimension>
- _computeMatrixC1(const double& rho,
- const Rd& U,
- const double& lambda,
- TinyVector<7> u_moments,
- TinyVector<2> xi2_moments)
+ std::tuple<CellValue<const TinyVector<SIZEproblem>>, CellValue<const TinyVector<SIZEproblem>>>
+ _computeGradConservVariables_for_a(const MeshType& mesh,
+ const DiscreteScalarFunction& rho,
+ const DiscreteVectorFunction& rhoU,
+ const DiscreteScalarFunction& rhoE) const
{
- const double U_2 = U[0] * U[0];
- const double U_4 = U_2 * U_2;
- const double lambda_2 = lambda * lambda;
+ auto& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- const double erf_term = (1 + std::erf(std::sqrt(lambda) * U[0]));
- const double exp_term = std::exp(-lambda * U_2) / std::sqrt(pi * lambda);
+ auto Vj = mesh_data.Vj();
- TinyMatrix<2 + Dimension> MatrixC1(zero);
+ auto node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
- MatrixC1(1, 1) = u_moments[2] * erf_term + rho * U * exp_term;
- MatrixC1(2, 2) = u_moments[4] * erf_term + rho * U * (U_2 + 2.5 / lambda) * exp_term;
- MatrixC1(3, 3) =
- (u_moments[6] + 2 * xi2_moments[0] * u_moments[4] + xi2_moments[1] * u_moments[2]) * erf_term +
- rho * (U * (U_4 + 7 * U_2 / lambda + 8.25 / lambda_2) + 2 * xi2_moments[0] * U * (U_2 + 2.5 / lambda)) * exp_term;
- MatrixC1(3, 3) *= 0.25;
+ CellValue<TinyVector<SIZEproblem>> Gradj{mesh.connectivity()};
+ Gradj.fill(zero);
+ CellValue<TinyVector<SIZEproblem>> Gradjpu{mesh.connectivity()};
+ Gradjpu.fill(zero);
- MatrixC1(1, 1) *= 0.25;
- MatrixC1(2, 2) *= 0.25;
- MatrixC1(3, 3) *= 0.25;
+ NodeValue<double> rho_node_mean{mesh.connectivity()};
+ NodeValue<Rd> rhoU_node_mean{mesh.connectivity()};
+ NodeValue<double> rhoE_node_mean{mesh.connectivity()};
+ rho_node_mean.fill(0);
+ rhoU_node_mean.fill(zero);
+ rhoE_node_mean.fill(0);
- MatrixC1(2, 1) = u_moments[3] + rho * (U_2 + 1 / lambda) * exp_term;
- MatrixC1(2, 1) *= 0.5;
- MatrixC1(3, 1) = 0.5 * (u_moments[4] + u_moments[2] * xi2_moments[1]) * erf_term +
- rho * U[0] * ((U_2 + 2.5 / lambda) + xi2_moments[0]) * exp_term;
- MatrixC1(3, 1) *= 0.5;
- MatrixC1(3, 2) = (u_moments[5] + u_moments[3] * xi2_moments[0]) * erf_term +
- rho * (U_4 + 4.5 * U_2 / lambda + 2 / lambda_2 + (U_2 + 1 / lambda) * xi2_moments[0]) * exp_term;
- MatrixC1(3, 2) *= 0.5;
-
- return (MatrixC1 + transpose(MatrixC1));
+ 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;
+
+ 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]];
+ ;
+ }
+ }
+ });
+ return {Gradj, Gradjpu};
}
- TinyMatrix<2 + Dimension>
- _computeMatrixC2(const double& rho,
- const Rd& U,
- const double& lambda,
- TinyVector<7> u_moments,
- TinyVector<2> xi2_moments)
+ TinyVector<SIZEproblem>
+ _compute_a_and_b(const TinyMatrix<SIZEproblem>& invMatrixM, const TinyVector<SIZEproblem>& Grad) const
+ {
+ TinyVector<SIZEproblem> a = 2 * invMatrixM * Grad;
+ return a;
+ }
+
+ TinyVector<SIZEproblem>
+ _compute_A(const TinyMatrix<SIZEproblem>& invMatrixM,
+ const TinyMatrix<SIZEproblem>& C1_full,
+ const TinyVector<SIZEproblem>& a) const
+ {
+ TinyVector<SIZEproblem> A = (invMatrixM * C1_full) * a;
+ return A;
+ }
+
+ TinyVector<SIZEproblem>
+ _compute_B(const TinyMatrix<SIZEproblem>& invMatrixM_interface,
+ const TinyMatrix<SIZEproblem>& C0_semi_pos,
+ const TinyMatrix<SIZEproblem>& C0_semi_neg,
+ const TinyVector<SIZEproblem>& al,
+ const TinyVector<SIZEproblem>& ar) const
+ {
+ TinyVector<SIZEproblem> B = invMatrixM_interface * (C0_semi_pos * al + C0_semi_neg * ar);
+ return B;
+ }
+
+ std::tuple<CellValue<TinyVector<SIZEproblem>>, CellValue<TinyVector<SIZEproblem>>>
+ _computeGradConservVariables_for_b(const MeshType& mesh,
+ const DiscreteScalarFunction& rho,
+ const DiscreteVectorFunction& rhoU,
+ const DiscreteScalarFunction& rhoE,
+ const NodeValue<double>& rho_node,
+ const NodeValue<Rd>& rhoU_node,
+ const NodeValue<double>& rhoE_node) const
+ {
+ auto& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ auto Vj = mesh_data.Vj();
+
+ auto node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+
+ 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]];
+ ;
+ }
+ }
+ });
+
+ return {Gradj, Gradjpu};
+ }
+
+ std::tuple<TinyMatrix<SIZEproblem>, TinyMatrix<SIZEproblem>, TinyMatrix<SIZEproblem>>
+ _computeMatricesC0C1C2_semi_moments(const double& rho,
+ const Rd& U,
+ const double& lambda,
+ const TinyVector<7>& u_moments,
+ const TinyVector<2>& xi2_moments,
+ const int& sign_u) const
{
const double U_2 = U[0] * U[0];
const double U_4 = U_2 * U_2;
const double lambda_2 = lambda * lambda;
- const double erf_term = (1 - std::erf(std::sqrt(lambda) * U[0]));
- const double exp_term = std::exp(-lambda * U_2) / std::sqrt(pi * lambda);
-
- TinyMatrix<2 + Dimension> MatrixC2(zero);
+ const double erf_term = (1 + sign_u * std::erf(std::sqrt(lambda) * U[0]));
+ const double exp_term = rho * sign_u * std::exp(-lambda * U_2) / std::sqrt(pi * lambda);
+
+ TinyMatrix<SIZEproblem> MatrixC0(zero);
+ MatrixC0(0, 0) = 0.5 * u_moments[0] * erf_term;
+ MatrixC0(1, 1) = 0.5 * (u_moments[2] * erf_term + U[0] * exp_term);
+ MatrixC0(2, 2) = 0.5 * (u_moments[4] * erf_term + U[0] * (U_2 + 2.5 / lambda) * exp_term) +
+ MatrixC0(0, 0) * xi2_moments[1] + 2 * MatrixC0(1, 1) * xi2_moments[0];
+ MatrixC0(2, 2) *= 0.25;
+
+ MatrixC0(1, 0) = 0.5 * (u_moments[1] * erf_term + exp_term);
+ MatrixC0(0, 1) = MatrixC0(1, 0);
+ MatrixC0(2, 0) = MatrixC0(1, 1) + MatrixC0(0, 0) * xi2_moments[0];
+ MatrixC0(2, 0) *= 0.5;
+ MatrixC0(0, 2) = MatrixC0(2, 0);
+ MatrixC0(2, 1) = 0.5 * (u_moments[3] * erf_term + (U_2 + 1 / lambda) * exp_term) + MatrixC0(1, 0) * xi2_moments[0];
+ MatrixC0(2, 1) *= 0.5;
+ MatrixC0(1, 2) = MatrixC0(2, 1);
+
+ TinyMatrix<SIZEproblem> MatrixC1(zero);
+ MatrixC1(0, 0) = MatrixC0(1, 0);
+ MatrixC1(1, 1) = 0.5 * (u_moments[3] * erf_term + (U_2 + 1 / lambda) * exp_term);
+ MatrixC1(2, 2) = 0.5 * (u_moments[5] * erf_term + (U_4 + 4.5 * U_2 / lambda + 2 / lambda_2) * exp_term) +
+ MatrixC1(0, 0) * xi2_moments[1] + 2 * MatrixC1(1, 1) * xi2_moments[0];
+ MatrixC1(2, 2) *= 0.25;
- MatrixC2(1, 1) = u_moments[2] * erf_term - rho * U * exp_term;
- MatrixC2(2, 2) = u_moments[4] * erf_term - rho * U * (U_2 + 2.5 / lambda) * exp_term;
- MatrixC2(3, 3) =
- (u_moments[6] + 2 * xi2_moments[0] * u_moments[4] + xi2_moments[1] * u_moments[2]) * erf_term -
- rho * (U * (U_4 + 7 * U_2 / lambda + 8.25 / lambda_2) + 2 * xi2_moments[0] * U * (U_2 + 2.5 / lambda)) * exp_term;
- MatrixC2(3, 3) *= 0.25;
+ MatrixC1(1, 0) = MatrixC0(1, 1);
+ MatrixC1(0, 1) = MatrixC1(1, 0);
+ MatrixC1(2, 0) = MatrixC1(1, 1) + MatrixC1(0, 0) * xi2_moments[0];
+ MatrixC1(2, 0) *= 0.5;
+ MatrixC1(0, 2) = MatrixC1(2, 0);
+ MatrixC1(2, 1) =
+ 0.5 * (u_moments[4] * erf_term + U[0] * (U_2 + 2.5 / lambda) * exp_term) + MatrixC1(1, 0) * xi2_moments[0];
+ MatrixC1(2, 1) *= 0.5;
+ MatrixC1(1, 2) = MatrixC1(2, 1);
- MatrixC2(1, 1) *= 0.25;
+ TinyMatrix<SIZEproblem> MatrixC2(zero);
+ MatrixC2(0, 0) = MatrixC1(1, 0);
+ MatrixC2(1, 1) = 0.5 * (u_moments[4] * erf_term + U[0] * (U_2 + 2.5 / lambda) * exp_term);
+ MatrixC2(2, 2) = 0.5 * (u_moments[6] * erf_term + U[0] * (U_4 + 7 * U_2 / lambda + 8.25 / lambda_2) * exp_term) +
+ MatrixC2(0, 0) * xi2_moments[1] + 2 * MatrixC2(1, 1) * xi2_moments[0];
MatrixC2(2, 2) *= 0.25;
- MatrixC2(3, 3) *= 0.25;
- MatrixC2(2, 1) = u_moments[3] * erf_term - rho * (U_2 + 1 / lambda) * exp_term;
+ MatrixC2(1, 0) = MatrixC1(1, 1);
+ MatrixC2(0, 1) = MatrixC2(1, 0);
+ MatrixC2(2, 0) = MatrixC2(1, 1) + MatrixC2(0, 0) * xi2_moments[0];
+ MatrixC2(2, 0) *= 0.5;
+ MatrixC2(0, 2) = MatrixC2(2, 0);
+ MatrixC2(2, 1) = 0.5 * (u_moments[5] * erf_term + (U_4 + 4.5 * U_2 / lambda + 2 / lambda_2) * exp_term) +
+ MatrixC2(1, 0) * xi2_moments[0];
MatrixC2(2, 1) *= 0.5;
- MatrixC2(3, 1) = 0.5 * (u_moments[4] + u_moments[2] * xi2_moments[1]) * erf_term -
- rho * U[0] * ((U_2 + 2.5 / lambda) + xi2_moments[0]) * exp_term;
- MatrixC2(3, 1) *= 0.5;
- MatrixC2(3, 2) = (u_moments[5] + u_moments[3] * xi2_moments[0]) * erf_term -
- rho * (U_4 + 4.5 * U_2 / lambda + 2 / lambda_2 + (U_2 + 1 / lambda) * xi2_moments[0]) * exp_term;
- MatrixC2(3, 2) *= 0.5;
-
- return (MatrixC2 + transpose(MatrixC2));
+ MatrixC2(1, 2) = MatrixC2(2, 1);
+
+ return {MatrixC0, MatrixC1, MatrixC2};
}
- TinyMatrix<2 + Dimension>
- _computeMatrixC3(const double& rho,
- const Rd& U,
- const double& lambda,
- TinyVector<7> u_moments,
- TinyVector<2> xi2_moments)
+ TinyMatrix<SIZEproblem>
+ _computeMatricesC1_full_moments(const TinyVector<7>& u_moments, const TinyVector<2>& xi2_moments) const
{
- const double U_2 = U[0] * U[0];
- const double U_4 = U_2 * U_2;
- const double lambda_2 = lambda * lambda;
+ TinyMatrix<SIZEproblem> MatrixC1(zero);
+
+ MatrixC1(0, 0) = u_moments[1];
+ MatrixC1(1, 1) = u_moments[3];
+ MatrixC1(2, 2) = 0.25 * (u_moments[5] + 2 * xi2_moments[0] * u_moments[3] + xi2_moments[1] * u_moments[1]);
+
+ MatrixC1(1, 0) = u_moments[2];
+ MatrixC1(0, 1) = MatrixC1(1, 0);
+ MatrixC1(2, 0) = 0.5 * (u_moments[3] + u_moments[1] * xi2_moments[0]);
+ MatrixC1(0, 2) = MatrixC1(2, 0);
+ MatrixC1(2, 1) = (u_moments[4] + u_moments[2] * xi2_moments[0]);
+ MatrixC1(1, 2) = MatrixC1(2, 1);
- const double erf_term = (1 - std::erf(std::sqrt(lambda) * U[0]));
- const double exp_term = std::exp(-lambda * U_2) / std::sqrt(pi * lambda);
-
- TinyMatrix<2 + Dimension> MatrixC3(zero);
-
- MatrixC3(1, 1) = u_moments[2] * erf_term - rho * U * exp_term;
- MatrixC3(2, 2) = u_moments[4] * erf_term - rho * U * (U_2 + 2.5 / lambda) * exp_term;
- MatrixC3(3, 3) =
- (u_moments[6] + 2 * xi2_moments[0] * u_moments[4] + xi2_moments[1] * u_moments[2]) * erf_term -
- rho * (U * (U_4 + 7 * U_2 / lambda + 8.25 / lambda_2) + 2 * xi2_moments[0] * U * (U_2 + 2.5 / lambda)) * exp_term;
- MatrixC3(3, 3) *= 0.25;
-
- MatrixC3(1, 1) *= 0.5;
- MatrixC3(2, 2) *= 0.5;
- MatrixC3(3, 3) *= 0.5;
-
- MatrixC3(2, 1) = u_moments[3] * erf_term - rho * (U_2 + 1 / lambda) * exp_term;
- MatrixC3(2, 1) *= 0.5;
- MatrixC3(3, 1) = 0.5 * (u_moments[4] + u_moments[2] * xi2_moments[1]) * erf_term -
- rho * U[0] * ((U_2 + 2.5 / lambda) + xi2_moments[0]) * exp_term;
- MatrixC3(3, 1) *= 0.5;
- MatrixC3(3, 2) = (u_moments[5] + u_moments[3] * xi2_moments[0]) * erf_term -
- rho * (U_4 + 4.5 * U_2 / lambda + 2 / lambda_2 + (U_2 + 1 / lambda) * xi2_moments[0]) * exp_term;
- MatrixC3(3, 2) *= 0.5;
-
- // Divide by 2 the diagonal because there are double terms on it after the addition of the matrix and its transpose
- MatrixC3(1, 1) *= 0.5;
- MatrixC3(2, 2) *= 0.5;
- MatrixC3(3, 3) *= 0.5;
-
- return (MatrixC3 + transpose(MatrixC3));
+ return {MatrixC1};
}
public:
@@ -275,51 +385,35 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
compute_fluxes(const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoU_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoE_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& tau_v,
+ const double& viscosity,
const double& delta,
const double& dt) const
{
auto mesh_v = getCommonMesh({rho_v, rhoU_v, rhoE_v});
const MeshType& mesh = *mesh_v->get<MeshType>();
- DiscreteFunctionP0<const double> tau_n = tau_v->get<DiscreteScalarFunction>();
- CellValue<double> eta(mesh.connectivity());
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- if (tau_n[cell_id] == 0)
- eta[cell_id] = 0;
- else
- eta[cell_id] = (tau_n[cell_id] / dt) * (1 - std::exp(-dt / tau_n[cell_id]));
- }
-
DiscreteScalarFunction rho_n = rho_v->get<DiscreteScalarFunction>();
DiscreteVectorFunction rhoU_n = rhoU_v->get<DiscreteVectorFunction>();
DiscreteScalarFunction rhoE_n = rhoE_v->get<DiscreteScalarFunction>();
- DiscreteFunctionP0<double> rho = copy(rho_n);
- DiscreteFunctionP0<Rd> rhoU = copy(rhoU_n);
- DiscreteFunctionP0<double> rhoE = copy(rhoE_n);
+ DiscreteScalarFunction rho = copy(rho_n);
+ DiscreteVectorFunction rhoU = copy(rhoU_n);
+ DiscreteScalarFunction rhoE = copy(rhoE_n);
+
+ CellValue<double> lambda{mesh.connectivity()};
+ lambda.fill(1);
+ CellValue<Rd> U{mesh.connectivity()};
+ U.fill(zero);
auto& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+ auto Vj = mesh_data.Vj();
+
auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
auto node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const NodeValuePerCell<const Rd> njr = mesh_data.njr();
- NodeValuePerCell<Rd> rho_flux_G(mesh.connectivity());
- NodeValuePerCell<Rd> rhoU_flux_G(mesh.connectivity());
- NodeValuePerCell<Rd> rhoE_flux_G(mesh.connectivity());
- rho_flux_G.fill(zero);
- rhoU_flux_G.fill(zero);
- rhoE_flux_G.fill(zero);
-
- NodeValuePerCell<Rd> rho_flux_Ffn(mesh.connectivity());
- NodeValuePerCell<Rd> rhoU_flux_Ffn(mesh.connectivity());
- NodeValuePerCell<Rd> rhoE_flux_Ffn(mesh.connectivity());
- rho_flux_Ffn.fill(zero);
- rhoU_flux_Ffn.fill(zero);
- rhoE_flux_Ffn.fill(zero);
-
CellValue<double> rho_fluxes(mesh.connectivity());
CellValue<Rd> rhoU_fluxes(mesh.connectivity());
CellValue<double> rhoE_fluxes(mesh.connectivity());
@@ -327,21 +421,16 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
rhoU_fluxes.fill(zero);
rhoE_fluxes.fill(0);
- NodeValue<Rd> F2fn(mesh.connectivity());
- NodeValue<Rd> F3fn(mesh.connectivity());
- F2fn.fill(zero);
- F3fn.fill(zero);
-
NodeValue<double> rho_node(mesh.connectivity());
NodeValue<Rd> rhoU_node(mesh.connectivity());
NodeValue<double> rhoE_node(mesh.connectivity());
- rho_node.fill(0);
+ NodeValue<Rd> U_node(mesh.connectivity());
+ NodeValue<double> lambda_node(mesh.connectivity());
+ rho_node.fill(1);
rhoU_node.fill(zero);
- rhoE_node.fill(0);
- CellValue<double> lambda{mesh.connectivity()};
- lambda.fill(1);
- CellValue<Rd> U{mesh.connectivity()};
- U.fill(zero);
+ rhoE_node.fill(1);
+ U_node.fill(zero);
+ lambda_node.fill(1);
CellValue<double> err_function_term(mesh.connectivity());
CellValue<double> exp_term(mesh.connectivity());
@@ -360,6 +449,36 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
exp_term[cell_id] = std::exp(-lambda[cell_id] * U_2) / 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];
@@ -398,77 +517,179 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
rhoE_node[node_id] = 0.5 * (rhoE_cell_left + rhoE_cell_right);
});
+ parallel_for(
+ mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ U_node[node_id][0] = rhoU_node[node_id][0] / rho_node[node_id];
+ double rhoU2_node = rhoU_node[node_id][0] * rhoU_node[node_id][0] / rho_node[node_id];
+ lambda_node[node_id] = (1 + delta) * rho_node[node_id] / (2 * (2 * rhoE_node[node_id] - rhoU2_node));
+ });
+
+ NodeValue<double> tau{mesh.connectivity()};
+ tau.fill(1);
+ NodeValue<double> eta{mesh.connectivity()};
+ eta.fill(1);
+ NodeValue<double> zeta{mesh.connectivity()};
+ zeta.fill(1);
+ for (NodeId node_id = 0; node_id < mesh.numberOfNodes(); ++node_id) {
+ tau[node_id] = viscosity * rho_node[node_id] / (2 * lambda_node[node_id]);
+ if (viscosity == 0) {
+ eta[node_id] = 0;
+ zeta[node_id] = 0;
+ } else {
+ zeta[node_id] = std::exp(-dt / tau[node_id]);
+ eta[node_id] = (tau[node_id] / dt) * (1 - zeta[node_id]);
+ }
+ }
+
+ 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);
+
+ NodeValue<TinyVector<SIZEproblem>> Fluxes_term{mesh.connectivity()};
+ Fluxes_term.fill(zero);
+
parallel_for(
mesh.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
const auto& node_cells = node_to_cell_matrix[node_id];
- double F2_fn_left = 0;
- double F3_fn_left = 0;
- double F2_fn_right = 0;
- double F3_fn_right = 0;
+ const double rho_left_node = rho[node_cells[0]] + 0.5 * Vj[node_cells[0]] * Gradj_for_a[node_cells[0]][0];
+ Rd rhoU_left_node;
+ rhoU_left_node[0] = rhoU[node_cells[0]][0] + 0.5 * Vj[node_cells[0]] * Gradj_for_a[node_cells[0]][1];
+ const double rhoE_left_node = rhoE[node_cells[0]] + 0.5 * Vj[node_cells[0]] * Gradj_for_a[node_cells[0]][2];
- for (size_t l = 0; l < node_cells.size(); l++) {
- if (node_cells.size() == 1)
- continue;
+ Rd U_left_node;
+ U_left_node[0] = rhoU_left_node[0] / rho_left_node;
+ const double lambda_left_node =
+ (1 + delta) * rho_left_node / (2 * (2 * rhoE_left_node - rhoU_left_node[0] * U_left_node[0]));
- if (l == 0) {
- double U_2_left = U[node_cells[l]][0] * U[node_cells[l]][0];
- double rhoU_2_left = rhoU_n[node_cells[l]][0] * U[node_cells[l]][0];
+ const double rho_right_node = rho[node_cells[1]] - 0.5 * Vj[node_cells[1]] * Gradj_for_a[node_cells[1]][0];
+ Rd rhoU_right_node;
+ 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];
- F2_fn_left = (rhoU_2_left + 0.5 * rho_n[node_cells[l]] / lambda[node_cells[l]]) *
- (1. + err_function_term[node_cells[l]]) +
- rhoU_n[node_cells[l]][0] * exp_term[node_cells[l]];
+ // peut etre des + 0.5* ici (au dessus)
- F3_fn_left = 0.5 * rhoU_n[node_cells[l]][0] * (U_2_left + 0.5 * (delta + 3) / lambda[node_cells[l]]) *
- (1. + err_function_term[node_cells[l]]) +
- 0.5 * rho_n[node_cells[l]] * (U_2_left + 0.5 * (delta + 2) / lambda[node_cells[l]]) *
- exp_term[node_cells[l]];
- } else {
- double U_2_right = U[node_cells[l]][0] * U[node_cells[l]][0];
- double rhoU_2_right = rhoU_n[node_cells[l]][0] * U[node_cells[l]][0];
+ Rd U_right_node;
+ U_right_node[0] = rhoU_right_node[0] / rho_right_node;
+ const double lambda_right_node =
+ (1 + delta) * rho_right_node / (2 * (2 * rhoE_right_node - rhoU_right_node[0] * U_right_node[0]));
+
+ TinyVector<7> u_moments_left = this->_compute_u_moments(rho_left_node, U_left_node, lambda_left_node);
+ TinyVector<2> xi2_moments_left = this->_compute_xi2_moments(lambda_left_node, delta);
+
+ 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);
- F2_fn_right = (rhoU_2_right + 0.5 * rho_n[node_cells[l]] / lambda[node_cells[l]]) *
- (1. - err_function_term[node_cells[l]]) -
- rhoU_n[node_cells[l]][0] * exp_term[node_cells[l]];
+ 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);
- F3_fn_right = 0.5 * rhoU_n[node_cells[l]][0] * (U_2_right + 0.5 * (delta + 3) / lambda[node_cells[l]]) *
- (1. - err_function_term[node_cells[l]]) -
- 0.5 * rho_n[node_cells[l]] * (U_2_right + 0.5 * (delta + 2) / lambda[node_cells[l]]) *
- exp_term[node_cells[l]];
+ // std::cout << std::endl << "InvM : " << InvM << "\t and \t InvM_bis" << InvM_bis << std::endl;
+
+ 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]]);
}
}
- F2fn[node_id][0] = 0.5 * (F2_fn_left + F2_fn_right);
- F3fn[node_id][0] = 0.5 * (F3_fn_left + F3_fn_right);
+ // 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);
+
+ 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,
+ xi2_moments_left, 1);
+
+ auto [C0_semi_neg_right, C1_semi_neg_right, C2_semi_neg_right] =
+ this->_computeMatricesC0C1C2_semi_moments(rho_right_node, U_right_node, lambda_right_node, u_moments_right,
+ xi2_moments_left, -1);
+
+ TinyVector<SIZEproblem> Al = this->_compute_A(InvM_left, C1_semi_pos_left, -al);
+ TinyVector<SIZEproblem> Ar = this->_compute_A(InvM_right, C1_semi_neg_right, -ar);
+
+ TinyVector<SIZEproblem> B = this->_compute_B(InvMInterface, C0_semi_pos_left, C0_semi_neg_right, al, ar);
+
+ TinyVector<SIZEproblem> Fluxes_from_a = (tau[node_id] * zeta[node_id] - tau[node_id] * eta[node_id]) *
+ (C2_semi_pos_left * al + C2_semi_neg_right * ar);
+ TinyVector<SIZEproblem> Fluxes_for_Navier =
+ -tau[node_id] * eta[node_id] *
+ ((C1_semi_pos_left * Al + C1_semi_neg_right * Ar) + (C2_semi_pos_left * al + C2_semi_neg_right * ar));
+
+ TinyVector<7> u_moments_interface =
+ this->_compute_u_moments(rho_node[node_id], U_node[node_id], lambda_node[node_id]);
+ TinyVector<2> xi2_moments_interface = this->_compute_xi2_moments(lambda_node[node_id], delta);
+
+ auto [C0_semi_pos_interface, C1_semi_pos_interface, C2_semi_pos_interface] =
+ this->_computeMatricesC0C1C2_semi_moments(rho_node[node_id], U_node[node_id], lambda_node[node_id],
+ u_moments_interface, xi2_moments_interface, 1);
+ auto [C0_semi_neg_interface, C1_semi_neg_interface, C2_semi_neg_interface] =
+ this->_computeMatricesC0C1C2_semi_moments(rho_node[node_id], U_node[node_id], lambda_node[node_id],
+ u_moments_interface, xi2_moments_interface, -1);
+
+ TinyMatrix<SIZEproblem> Matrix_for_B =
+ this->_computeMatricesC1_full_moments(u_moments_interface, xi2_moments_interface);
+
+ TinyVector<SIZEproblem> Fluxes_from_b = (2 * tau[node_id] * eta[node_id] - tau[node_id] * (1 + zeta[node_id])) *
+ (C2_semi_pos_interface * bl + C2_semi_neg_interface * br);
+ TinyVector<SIZEproblem> Fluxes_from_B =
+ (dt / 2 - tau[node_id] + tau[node_id] * eta[node_id]) * Matrix_for_B * B;
+
+ TinyVector<SIZEproblem> G;
+ G[0] = u_moments_interface[1];
+ G[1] = u_moments_interface[2];
+ G[2] = 0.5 * (u_moments_interface[3] + u_moments_interface[1] * xi2_moments_interface[0]);
+
+ TinyVector<SIZEproblem> F_fn;
+ F_fn[0] = C1_semi_pos_left(0, 0) + C1_semi_neg_right(0, 0);
+ 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;
});
+ NodeValuePerCell<TinyVector<SIZEproblem>> Fluxes_term_apply{mesh.connectivity()};
+ Fluxes_term_apply.fill(zero);
+
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
const auto& cell_nodes = cell_to_node_matrix[cell_id];
for (size_t r = 0; r < cell_nodes.size(); r++) {
- double rhoU_2_node = rhoU_node[cell_nodes[r]][0] * rhoU_node[cell_nodes[r]][0] / rho_node[cell_nodes[r]];
-
- rho_flux_G(cell_id, r) = rhoU_node[cell_nodes[r]][0] * njr(cell_id, r);
- rhoU_flux_G(cell_id, r) =
- (delta * rhoU_2_node / (1. + delta) + 2 * rhoE_node[cell_nodes[r]] / (1. + delta)) * njr(cell_id, r);
- rhoE_flux_G(cell_id, r) =
- rhoU_node[cell_nodes[r]][0] / rho_node[cell_nodes[r]] *
- ((3. + delta) * rhoE_node[cell_nodes[r]] / (1. + delta) - rhoU_2_node / (1. + delta)) * njr(cell_id, r);
-
- rhoU_flux_Ffn(cell_id, r) = F2fn[cell_nodes[r]][0] * njr(cell_id, r);
- rhoE_flux_Ffn(cell_id, r) = F3fn[cell_nodes[r]][0] * njr(cell_id, r);
+ Fluxes_term_apply(cell_id, r) = njr(cell_id, r)[0] * Fluxes_term[cell_nodes[r]];
}
});
- for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
- const size_t& nb_nodes = rho_flux_G.numberOfSubValues(cell_id);
+ for (CellId cell_id = 2; cell_id < mesh.numberOfCells() - 2; ++cell_id) {
+ const size_t& nb_nodes = Fluxes_term_apply.numberOfSubValues(cell_id);
for (size_t r = 0; r < nb_nodes; r++) {
- rho_fluxes[cell_id] += rho_flux_G(cell_id, r)[0];
- rhoU_fluxes[cell_id] +=
- ((1 - eta[cell_id]) * rhoU_flux_G(cell_id, r) + eta[cell_id] * (rhoU_flux_Ffn(cell_id, r)));
- rhoE_fluxes[cell_id] +=
- ((1 - eta[cell_id]) * rhoE_flux_G(cell_id, r)[0] + eta[cell_id] * (rhoE_flux_Ffn(cell_id, r)[0]));
+ rho_fluxes[cell_id] += Fluxes_term_apply(cell_id, r)[0];
+ rhoU_fluxes[cell_id][0] += Fluxes_term_apply(cell_id, r)[1];
+ rhoE_fluxes[cell_id] += Fluxes_term_apply(cell_id, r)[2];
}
}
@@ -547,13 +768,13 @@ class GKSHandler::GKSNAVIER final : public GKSHandler::IGKSNAVIER
gksNavier(const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoU_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoE_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& tau_v,
+ const double& viscosity,
const double& delta,
const double& dt) const
{
std::shared_ptr mesh_v = getCommonMesh({rho_v, rhoU_v, rhoE_v});
- auto [rho_fluxes, rhoU_fluxes, rhoE_fluxes] = compute_fluxes(rho_v, rhoU_v, rhoE_v, tau_v, delta, dt);
+ auto [rho_fluxes, rhoU_fluxes, rhoE_fluxes] = compute_fluxes(rho_v, rhoU_v, rhoE_v, viscosity, delta, dt);
return apply_fluxes(rho_v, rhoU_v, rhoE_v, rho_fluxes, rhoU_fluxes, rhoE_fluxes, dt);
}
@@ -572,7 +793,11 @@ GKSHandler::GKSHandler(const std::shared_ptr<const MeshVariant>& mesh_v)
[&](auto&& mesh) {
using MeshType = mesh_type_t<decltype(mesh)>;
if constexpr (is_polygonal_mesh_v<MeshType>) {
- m_gks_navier = std::make_unique<GKSNAVIER<MeshType>>();
+ if constexpr (MeshType::Dimension == 1) {
+ m_gks_navier = std::make_unique<GKSNAVIER<MeshType>>();
+ } else {
+ throw NormalError("unexpected mesh type");
+ }
} else {
throw NormalError("unexpected mesh type");
}
diff --git a/src/scheme/GKSNavier.hpp b/src/scheme/GKSNavier.hpp
index 571f87d83..2410d32d3 100644
--- a/src/scheme/GKSNavier.hpp
+++ b/src/scheme/GKSNavier.hpp
@@ -19,13 +19,6 @@ double gks_inv_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c,
class GKSHandler
{
- public:
- enum class SolverType
- {
- Glace,
- Eucclhyd
- };
-
private:
struct IGKSNAVIER
{
@@ -35,7 +28,7 @@ class GKSHandler
compute_fluxes(const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoU_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoE_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& tau_v,
+ const double& viscosity,
const double& delta,
const double& dt) const = 0;
@@ -56,7 +49,7 @@ class GKSHandler
gksNavier(const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoU_v,
const std::shared_ptr<const DiscreteFunctionVariant>& rhoE_v,
- const std::shared_ptr<const DiscreteFunctionVariant>& tau_v,
+ const double& viscosity,
const double& delta,
const double& dt) const = 0;
--
GitLab