From 61a3f496d5b231b32d88d0ab677b794d444a8ed8 Mon Sep 17 00:00:00 2001
From: HOCH PHILIPPE <philippe.hoch@gmail.com>
Date: Mon, 19 Aug 2024 23:40:24 +0200
Subject: [PATCH] debut mise en place limitation induite o2
---
.../RusanovEulerianCompositeSolver_v2.cpp | 606 +++++++++++++++++-
1 file changed, 603 insertions(+), 3 deletions(-)
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
index bd7a407b0..5111335a2 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
@@ -12,9 +12,13 @@
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshTraits.hpp>
#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionDPk.hpp>
+#include <scheme/DiscreteFunctionDPkVariant.hpp>
+#include <scheme/DiscreteFunctionDPkVector.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
-
+#include <scheme/PolynomialReconstruction.hpp>
+#include <utils/PugsTraits.hpp>
#include <variant>
template <MeshConcept MeshTypeT>
@@ -598,6 +602,517 @@ class RusanovEulerianCompositeSolver_v2
return (gam - 1) * rho * epsilon;
}
+ // void
+ // computeLimitorVolumicScalarQuantityMinModDukowicz(const DiscreteFunctionDPk<Dimension, double>& q_bar,
+ // CellValue<double>& Limitor_q) const
+ void
+ computeLimitorVolumicScalarQuantityMinModDukowiczGradient(const MeshType& mesh,
+ const CellValue<double>& q,
+ const CellValue<Rd>& gradq,
+ CellValue<double>& Limitor_q,
+ const bool enableWeakBoundPositivityOnly = false) const
+ {
+ Assert(Dimension == 2 or Dimension == 3);
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+
+ const auto& xr = mesh.xr();
+
+ MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ // auto volumes_cell = mesh_data.Vj();
+
+ const auto& xcentroid = mesh_data.xj();
+ const auto& xface = mesh_data.xl();
+
+ if (Dimension == 2) {
+ parallel_for(
+ mesh.numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+ const double State = q[j];
+ const Rd& Cent = xcentroid[j]; // mesh.xi(j);
+ const Rd& Gradl = gradq[j];
+
+ // Min et Max des valeurs voisines
+
+ double minvalvois(State);
+ double maxvalvois(State);
+
+ // Calcul des extremas de la fonction
+
+ double minval(std::numeric_limits<double>::max());
+ double maxval(-std::numeric_limits<double>::max());
+
+ // At nodes
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ // At faces
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ const NodeId& node = cell_to_node[l];
+ const auto& node_to_cell = node_to_cell_matrix[node];
+
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const CellId& cellvois = node_to_cell[k];
+ const double statevoiscell = q[cellvois];
+ const Rd Centvois = xcentroid[cellvois]; // mesh.element()[mayvois].centroid().x();
+ const Rd Mvois(Centvois - Cent);
+
+ const double valVois = State + dot(Gradl, Mvois); // Partie MinMod
+
+ minval = std::min(minval, valVois);
+ maxval = std::max(maxval, valVois);
+
+ minvalvois = std::min(minvalvois, statevoiscell);
+ maxvalvois = std::max(maxvalvois, statevoiscell);
+ }
+
+ const Rd M(xr[node] - Cent);
+
+ const double valLineO2node = State + dot(Gradl, M); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2node);
+ maxval = std::max(maxval, valLineO2node);
+ }
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const FaceId& face = cell_to_face[l];
+
+ const Rd Mb(xface[face] - Cent);
+
+ const double valLineO2bras = State + dot(Gradl, Mb); // Idem aux demi faces
+
+ minval = std::min(minval, valLineO2bras);
+ maxval = std::max(maxval, valLineO2bras);
+ }
+
+ static const double epsZer0 = 1e-15;
+ // on applique le traitement idem de l'ordre 2..
+ double coef1 = 0, coef2 = 0;
+ if (enableWeakBoundPositivityOnly) {
+ coef1 = 1;
+ const double minGlobal = 1e-12;
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minGlobal - State) / ((minval - State));
+
+ } else {
+ if (std::fabs(maxval - State) <= epsZer0) // epsIsEqualAbs())
+ coef1 = 1.;
+ else
+ coef1 = (maxvalvois - State) / ((maxval - State));
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minvalvois - State) / ((minval - State));
+ }
+ const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+ Limitor_q[j] = alfa;
+ // grad(i) *= alfa;
+ });
+ } else {
+ // throw NormalError(" Limiteur Scal Grad non encore implemente 3D");
+ const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
+ const auto& xedge = mesh_data.xe();
+
+ parallel_for(
+ mesh.numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+ const double State = q[j];
+ const Rd& Cent = xcentroid[j]; // mesh.xi(j);
+ const Rd& Gradl = gradq[j];
+
+ // Min et Max des valeurs voisines
+
+ double minvalvois(State);
+ double maxvalvois(State);
+
+ // Calcul des extremas de la fonction
+
+ double minval(std::numeric_limits<double>::max());
+ double maxval(-std::numeric_limits<double>::max());
+
+ // At nodes
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ // At faces
+ const auto& cell_to_face = cell_to_face_matrix[j];
+ // At faces
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ const NodeId& node = cell_to_node[l];
+ const auto& node_to_cell = node_to_cell_matrix[node];
+
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const CellId& cellvois = node_to_cell[k];
+ const double statevoiscell = q[cellvois];
+ const Rd Centvois = xcentroid[cellvois]; // mesh.element()[mayvois].centroid().x();
+ const Rd Mvois(Centvois - Cent);
+
+ const double valVois = State + dot(Gradl, Mvois); // Partie MinMod
+
+ minval = std::min(minval, valVois);
+ maxval = std::max(maxval, valVois);
+
+ minvalvois = std::min(minvalvois, statevoiscell);
+ maxvalvois = std::max(maxvalvois, statevoiscell);
+ }
+
+ const Rd M(xr[node] - Cent);
+
+ const double valLineO2node = State + dot(Gradl, M); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2node);
+ maxval = std::max(maxval, valLineO2node);
+ }
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const FaceId& face = cell_to_face[l];
+
+ const Rd Mb(xface[face] - Cent);
+
+ const double valLineO2face = State + dot(Gradl, Mb); // Idem aux demi faces
+
+ minval = std::min(minval, valLineO2face);
+ maxval = std::max(maxval, valLineO2face);
+ }
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ const EdgeId& edge = cell_to_edge[l];
+
+ const Rd Me(xedge[edge] - Cent);
+
+ const double valLineO2edge = State + dot(Gradl, Me); // Idem aux demi faces
+
+ minval = std::min(minval, valLineO2edge);
+ maxval = std::max(maxval, valLineO2edge);
+ }
+
+ static const double epsZer0 = 1e-15;
+ // on applique le traitement idem de l'ordre 2..
+ double coef1 = 0, coef2 = 0;
+ if (enableWeakBoundPositivityOnly) {
+ coef1 = 1;
+ const double minGlobal = 1e-12;
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minGlobal - State) / ((minval - State));
+
+ } else {
+ if (std::fabs(maxval - State) <= epsZer0) // epsIsEqualAbs())
+ coef1 = 1.;
+ else
+ coef1 = (maxvalvois - State) / ((maxval - State));
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minvalvois - State) / ((minval - State));
+ }
+ const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+ Limitor_q[j] = alfa;
+ // grad(i) *= alfa;
+ });
+ }
+ }
+
+ void
+ computeLimitorInternalNrjScalarQuantityMinModDukowiczGradient(const MeshType& mesh,
+ const CellValue<double>& rho,
+ const CellValue<Rd>& gradrho,
+ const CellValue<double>& internalnrj,
+ const CellValue<Rd>& gradinternalnrj,
+ // const CellValue<Rd>& u,
+ const CellValue<Rdxd>& gradu,
+ CellValue<double>& Limitor_eps,
+ const bool enableWeakBoundPositivityOnly = false) const
+ {
+ Assert(Dimension == 2 or Dimension == 3);
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& xr = mesh.xr();
+
+ MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+ const auto& xcentroid = mesh_data.xj();
+ const auto& xface = mesh_data.xl();
+
+ if (Dimension == 2) {
+ parallel_for(
+ mesh.numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+ const double rhol = rho[j];
+ const Rd& Cent = xcentroid[j];
+ const Rd& grad_rhol = gradrho[j];
+
+ // const Rd& ul = u[j];
+ const Rdxd& grad_ul = gradu[j];
+
+ const double internal_nrj = internalnrj[j];
+ // const double kinetic_nrj = .5 * dot(ul, ul);
+
+ const Rd& grad_internal_nrj = gradinternalnrj[j];
+
+ // const Rd& grad_kinetic_internal_nrj = gradeps[j];
+
+ // const double qtt_limitor_density = rhol / (rhol + dot(grad,(x-xc));
+ // At nodes
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ // At faces
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ // Min et Max des valeurs voisines
+ const double State = internal_nrj;
+
+ double minvalvois(State);
+ double maxvalvois(State);
+
+ // Calcul des extremas de la fonction
+
+ double minval(std::numeric_limits<double>::max());
+ double maxval(-std::numeric_limits<double>::max());
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ const NodeId& node = cell_to_node[l];
+ const auto& node_to_cell = node_to_cell_matrix[node];
+
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const CellId& cellvois = node_to_cell[k];
+ const Rd Centvois = xcentroid[cellvois];
+ const Rd Mvois(Centvois - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mvois));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mvois);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mvois);
+
+ const double valVois = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie MinMod
+
+ minval = std::min(minval, valVois);
+ maxval = std::max(maxval, valVois);
+
+ const double statevoiscell = internalnrj[cellvois];
+
+ minvalvois = std::min(minvalvois, statevoiscell);
+ maxvalvois = std::max(maxvalvois, statevoiscell);
+ }
+
+ const Rd Mr(xr[node] - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mr));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mr);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mr);
+
+ const double valLineO2node = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2node);
+ maxval = std::max(maxval, valLineO2node);
+ }
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const FaceId& face = cell_to_face[l];
+
+ const Rd Mb(xface[face] - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mb));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mb);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mb);
+
+ const double valLineO2face = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2face);
+ maxval = std::max(maxval, valLineO2face);
+ }
+
+ static const double epsZer0 = 1e-15;
+
+ // on applique le traitement idem de l'ordre 2..
+ double coef1 = 0, coef2 = 0;
+ if (enableWeakBoundPositivityOnly) {
+ coef1 = 1;
+ const double minGlobal = 1e-12;
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minGlobal - State) / ((minval - State));
+
+ } else {
+ if (std::fabs(maxval - State) <= epsZer0) // epsIsEqualAbs())
+ coef1 = 1.;
+ else
+ coef1 = (maxvalvois - State) / ((maxval - State));
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minvalvois - State) / ((minval - State));
+ }
+ const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+ Limitor_eps[j] = alfa;
+ });
+ } else {
+ // throw NormalError(" Limiteur Scal Grad internal nrj non encore implemente 3D");
+ const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
+
+ const auto& xedge = mesh_data.xe();
+
+ parallel_for(
+ mesh.numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j] - rho[j]); });
+ const double rhol = rho[j];
+ const Rd& Cent = xcentroid[j];
+ const Rd& grad_rhol = gradrho[j];
+
+ // const Rd& ul = u[j];
+ const Rdxd& grad_ul = gradu[j];
+
+ const double internal_nrj = internalnrj[j];
+ // const double kinetic_nrj = .5 * dot(ul, ul);
+
+ const Rd& grad_internal_nrj = gradinternalnrj[j];
+
+ // const Rd& grad_kinetic_internal_nrj = gradeps[j];
+
+ // const double qtt_limitor_density = rhol / (rhol + dot(grad,(x-xc));
+ // At nodes
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ // At faces
+ const auto& cell_to_face = cell_to_face_matrix[j];
+ // At edges
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+
+ // Min et Max des valeurs voisines
+ const double State = internal_nrj;
+
+ double minvalvois(State);
+ double maxvalvois(State);
+
+ // Calcul des extremas de la fonction
+
+ double minval(std::numeric_limits<double>::max());
+ double maxval(-std::numeric_limits<double>::max());
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ const NodeId& node = cell_to_node[l];
+ const auto& node_to_cell = node_to_cell_matrix[node];
+
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const CellId& cellvois = node_to_cell[k];
+ const Rd Centvois = xcentroid[cellvois];
+ const Rd Mvois(Centvois - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mvois));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mvois);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mvois);
+
+ const double valVois = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie MinMod
+
+ minval = std::min(minval, valVois);
+ maxval = std::max(maxval, valVois);
+
+ const double statevoiscell = internalnrj[cellvois];
+
+ minvalvois = std::min(minvalvois, statevoiscell);
+ maxvalvois = std::max(maxvalvois, statevoiscell);
+ }
+
+ const Rd Mr(xr[node] - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mr));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mr);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mr);
+
+ const double valLineO2node = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2node);
+ maxval = std::max(maxval, valLineO2node);
+ }
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const FaceId& face = cell_to_face[l];
+
+ const Rd Mb(xface[face] - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Mb));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Mb);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Mb);
+
+ const double valLineO2face = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2face);
+ maxval = std::max(maxval, valLineO2face);
+ }
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ const EdgeId& edge = cell_to_edge[l];
+
+ const Rd Me(xedge[edge] - Cent);
+
+ const double qtt_limitor_density = rhol / (rhol + dot(grad_rhol, Me));
+
+ const Rd& High_ordervelocity = qtt_limitor_density * (grad_ul * Me);
+ const double High_orderinternalnrj = qtt_limitor_density * dot(grad_internal_nrj, Me);
+
+ const double valLineO2edge = internal_nrj + High_orderinternalnrj -
+ .5 * dot(High_ordervelocity, High_ordervelocity); // Partie Dukowicz
+
+ minval = std::min(minval, valLineO2edge);
+ maxval = std::max(maxval, valLineO2edge);
+ }
+
+ static const double epsZer0 = 1e-15;
+
+ // on applique le traitement idem de l'ordre 2..
+ double coef1 = 0, coef2 = 0;
+ if (enableWeakBoundPositivityOnly) {
+ coef1 = 1;
+ const double minGlobal = 1e-12;
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minGlobal - State) / ((minval - State));
+
+ } else {
+ if (std::fabs(maxval - State) <= epsZer0) // epsIsEqualAbs())
+ coef1 = 1.;
+ else
+ coef1 = (maxvalvois - State) / ((maxval - State));
+
+ if (std::fabs(minval - State) <= epsZer0) // epsIsEqualAbs())
+ coef2 = 1.;
+ else
+ coef2 = (minvalvois - State) / ((minval - State));
+ }
+ const double alfa = std::max(0., std::min(1., std::min(coef1, coef2)));
+
+ Limitor_eps[j] = alfa;
+ });
+ }
+ }
+
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
@@ -615,6 +1130,88 @@ class RusanovEulerianCompositeSolver_v2
auto rho = copy(rho_n);
auto u = copy(u_n);
auto E = copy(E_n);
+ // using DiscreteScalarFunction = DiscreteFunctionP0<const double>;
+ // using DiscreteVectorFunction = DiscreteFunctionP0<const Rd>;
+ // DiscreteScalarFunction rhoD = rho_n.get<DiscreteScalarFunction>();
+ // DiscreteVectorFunction uD = u_n.get<DiscreteVectorFunction>();
+ // DiscreteScalarFunction ED = E_n.get<DiscreteScalarFunction>();
+
+ std::vector<std::shared_ptr<const IBoundaryDescriptor>> symmetry_boundary_descriptor_list;
+
+ for (auto&& bc_descriptor : bc_descriptor_list) {
+ if (bc_descriptor->type() == IBoundaryConditionDescriptor::Type::symmetry) {
+ symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
+ }
+ }
+ static const int degree = 1;
+ PolynomialReconstructionDescriptor
+ reconstruction_descriptor(IntegrationMethodType::cell_center, // IntegrationMethodType::boundary,
+ degree, symmetry_boundary_descriptor_list);
+ //
+ //
+ //
+ CellValue<double> valrho({p_mesh->connectivity()});
+ CellValue<Rd> valu({p_mesh->connectivity()});
+ CellValue<double> valeps({p_mesh->connectivity()});
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ valrho[j] = rho[j];
+ valu[j] = u[j];
+ valeps[j] = E[j] - .5 * dot(u[j], u[j]);
+ });
+
+ DiscreteFunctionP0<double> eps(p_mesh, valeps);
+
+ // assemblage direct
+ // auto reconstructions = PolynomialReconstruction{reconstruction_descriptor}.build(rho, rho * u, rho * E);
+
+ // Pour assemblage Leibniz
+ // ca calcul au moins les gradients
+ auto reconstructions = PolynomialReconstruction{reconstruction_descriptor}.build(rho, u, eps);
+
+ DiscreteFunctionDPk rho_bar = reconstructions[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ DiscreteFunctionDPk u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ DiscreteFunctionDPk eps_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+ // DiscreteFunctionDPk rho_u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
+ // DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+ auto rho_bar_lim = rho_bar;
+ auto rho_u_bar_lim = u_bar;
+ auto rho_E_bar_lim = eps_bar;
+
+ CellValue<double> Limitor_rho(p_mesh->connectivity());
+ CellValue<double> Limitor_u(p_mesh->connectivity());
+ CellValue<double> Limitor_eps(p_mesh->connectivity());
+
+ CellValue<Rd> grad_rho_cell(p_mesh->connectivity());
+ CellValue<Rdxd> grad_u_cell(p_mesh->connectivity());
+ CellValue<Rd> grad_eps_cell(p_mesh->connectivity());
+
+ // DiscreteFunctionDPk<Dimension, const double> rho_bar_lim;
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j](x) - rho[j]);
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ grad_rho_cell[j][dim] = rho_bar.coefficients(j)[1 + dim]; // si base : x, puis y , puis z
+ grad_eps_cell[j][dim] = eps_bar.coefficients(j)[1 + dim]; // idem ci dessus
+ const Rd gradu = u_bar.coefficients(j)[1 + dim];
+ for (size_t dim2 = 0; dim2 < Dimension; ++dim2) {
+ grad_u_cell[j](dim, dim2) = gradu[dim2];
+ }
+ }
+ grad_u_cell[j] = transpose(grad_u_cell[j]);
+ });
+
+ // Appels des limiteurs cell by cell
+ computeLimitorVolumicScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, Limitor_rho);
+
+ computeLimitorInternalNrjScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, valeps, grad_eps_cell,
+ // valu,
+ grad_u_cell, Limitor_eps);
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Limitor_u[j] = sqrt(Limitor_eps[j]); });
+
// auto c = copy(c_n);
// auto p = copy(p_n);
@@ -637,6 +1234,9 @@ class RusanovEulerianCompositeSolver_v2
});
// CellValue<Rp> State{p_mesh->connectivity()};
+ //
+ // Remplir ici les reconstructions d'ordre élevé
+
NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()};
StateAtNode.fill(zero);
parallel_for(
@@ -706,7 +1306,7 @@ class RusanovEulerianCompositeSolver_v2
Flux_rhoAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim];
}
});
-*/
+ */
NodeValuePerCell<Rdxd> Flux_qtmvtAtCellNode{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
EdgeValuePerCell<Rdxd> Flux_qtmvtAtCellEdge{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
FaceValuePerCell<Rdxd> Flux_qtmvtAtCellFace{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
@@ -734,7 +1334,7 @@ class RusanovEulerianCompositeSolver_v2
Flux_qtmvtAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim];
}
});
-*/
+ */
// parallel_for(
// p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
// Flux_qtmvtAtCellNode[j] = Flux_qtmvtAtCellEdge[j] = Flux_qtmvtAtCellFace[j] = Flux_qtmvt[j];
--
GitLab