From 25b37797740829bcc28e86b41ced9de46276584a Mon Sep 17 00:00:00 2001
From: HOCH PHILIPPE <philippe.hoch@gmail.com>
Date: Fri, 23 Aug 2024 05:56:29 +0200
Subject: [PATCH] Adding second order in space for 3 current solver with
induced limitation
---
src/language/modules/SchemeModule.cpp | 111 +++
src/scheme/CMakeLists.txt | 3 +
.../RusanovEulerianCompositeSolver_v2.cpp | 684 +-----------------
3 files changed, 120 insertions(+), 678 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 32f689003..13ce019b4 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -44,8 +44,11 @@
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
#include <scheme/RoeViscousFormEulerianCompositeSolver_v2.hpp>
+#include <scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.hpp>
#include <scheme/RusanovEulerianCompositeSolver.hpp>
+#include <scheme/RusanovEulerianCompositeSolver_o2.hpp>
#include <scheme/RusanovEulerianCompositeSolver_v2.hpp>
+#include <scheme/RusanovEulerianCompositeSolver_v2_o2.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
#include <scheme/VariableBCDescriptor.hpp>
#include <scheme/test_reconstruction.hpp>
@@ -543,6 +546,44 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("rusanov_eulerian_composite_solver_version1_o2",
+ std::function(
+
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return rusanovEulerianCompositeSolver_o2(rho, u, E, gamma, c, p, degree,
+ bc_descriptor_list, dt);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("rusanov_eulerian_composite_solver_version1_o2_with_checks",
+ std::function(
+
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return rusanovEulerianCompositeSolver_o2(rho, u, E, gamma, c, p, degree,
+ bc_descriptor_list, dt, true);
+ }
+
+ ));
+
this->_addBuiltinFunction("rusanov_eulerian_composite_solver_version2",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
@@ -612,6 +653,76 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("rusanov_eulerian_composite_solver_version2_o2",
+ std::function(
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return rusanovEulerianCompositeSolver_v2_o2(rho, u, E, gamma, c, p, degree,
+ bc_descriptor_list, dt);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("rusanov_eulerian_composite_solver_version2_o2_with_checks",
+ std::function(
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return rusanovEulerianCompositeSolver_v2_o2(rho, u, E, gamma, c, p, degree,
+ bc_descriptor_list, dt, true);
+ }));
+
+ this->_addBuiltinFunction("roe_viscosityform_eulerian_composite_solver_version2_o2",
+ std::function(
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return roeViscousFormEulerianCompositeSolver_v2_o2(rho, u, E, gamma, c, p, degree,
+ bc_descriptor_list, dt);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("roe_viscosityform_eulerian_composite_solver_version2_o2_with_checks",
+ std::function(
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return roeViscousFormEulerianCompositeSolver_v2_o2(rho, u, E, gamma, c, p, degree,
+ bc_descriptor_list, dt, true);
+ }
+
+ ));
+
this->_addBuiltinFunction("eucclhyd_fluxes",
std::function(
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index bae8c3faa..14739b777 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -26,6 +26,9 @@ add_library(
RusanovEulerianCompositeSolver.cpp
RusanovEulerianCompositeSolver_v2.cpp
RoeViscousFormEulerianCompositeSolver_v2.cpp
+ RusanovEulerianCompositeSolver_o2.cpp
+ RusanovEulerianCompositeSolver_v2_o2.cpp
+ RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
test_reconstruction.cpp
)
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
index 5e8bc7ffe..c00f2311b 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
@@ -12,13 +12,8 @@
#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>
@@ -602,517 +597,6 @@ 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>>
@@ -1131,126 +615,6 @@ 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>>();
-
- 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
- // Pour rho (plus ou moins classique)
- computeLimitorVolumicScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, Limitor_rho);
-
- // Pour les variables d'interet (physique)
- // Pour eps ici G.P. (avec Leibniz et toute l artillerie)
- computeLimitorInternalNrjScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, valeps, grad_eps_cell,
- // valu,
- grad_u_cell, Limitor_eps);
- // Pour u : relation deduite du principe d'induction
- parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Limitor_u[j] = sqrt(Limitor_eps[j]); });
-
- //
- // Creation des variables conservatives limitées
- //
- auto rho_bar_lim = copy(rho_bar);
- auto rho_u_bar_lim = copy(u_bar); // a modifier .. du a la densite
- auto rho_E_bar_lim = copy(eps_bar); // a modifier et a completer : densité et energie cinetique
-
- // Assemblage des ctes et ordres élevés
- parallel_for(
- p_mesh->numberOfCells(),
- PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j](x) - rho[j]);
- // rho_u_bar_lim.coefficients(j)[0] = valrho[j] * valu[j]; // * valu[j];
- // rho_E_bar_lim.coefficients(j)[0] = valrho[j] * (valeps[j] + .5 * dot(valu[j], valu[j]));
-
- // rho_bar_lim.coefficients(j)[0] // cte inchangée
- rho_u_bar_lim.coefficients(j)[0] *= valrho[j]; // du coup (rho u)
- rho_E_bar_lim.coefficients(j)[0] *= valrho[j]; // du coup (rho*epsilon)
- rho_E_bar_lim.coefficients(j)[0] += valrho[j] * .5 * dot(valu[j], valu[j]);
- // Rdxd grad_rhou_lim=
- // valrho[j]*Limitor_u[j]*grad_u_cell[j]+Limitor_rho[j]*tensorProduct(valu[j],grad_rho_cell[j]);
- //
- Rdxd gradU_lim_transpose = Limitor_u[j] * transpose(grad_u_cell[j]);
- const Rd gradUtu_lim = gradU_lim_transpose * valu[j];
-
- for (size_t dim = 0; dim < Dimension; ++dim) {
- rho_bar_lim.coefficients(j)[1 + dim] *= Limitor_rho[j];
- rho_u_bar_lim.coefficients(j)[1 + dim] *= Limitor_u[j] * valrho[j]; // c'est la partie gradu limite
- rho_u_bar_lim.coefficients(j)[1 + dim] += Limitor_rho[j] * grad_rho_cell[j][dim] * valu[j];
-
- rho_E_bar_lim.coefficients(j)[1 + dim] *= Limitor_eps[j] * valrho[j];
- rho_E_bar_lim.coefficients(j)[1 + dim] += gradUtu_lim[dim] * valrho[j];
-
- rho_E_bar_lim.coefficients(j)[1 + dim] +=
- (valeps[j] + .5 * dot(valu[j], valu[j])) * Limitor_rho[j] * grad_rho_cell[j][dim];
- }
- });
-
// auto c = copy(c_n);
// auto p = copy(p_n);
@@ -1281,59 +645,23 @@ class RusanovEulerianCompositeSolver_v2
const auto& cell_to_edge_matrix = p_mesh->connectivity().cellToEdgeMatrix();
const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
- const auto xr = p_mesh->xr();
- auto xl = MeshDataManager::instance().getMeshData(*p_mesh).xl();
- auto xe = MeshDataManager::instance().getMeshData(*p_mesh).xe();
+ // const auto xr = p_mesh->xr();
+ // auto xl = MeshDataManager::instance().getMeshData(*p_mesh).xl();
+ // auto xe = MeshDataManager::instance().getMeshData(*p_mesh).xe();
NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()};
StateAtNode.fill(zero);
parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
- // StateAtNode[j].fill(State[j]);
- const auto& cell_to_node = cell_to_node_matrix[j];
-
- for (size_t l = 0; l < cell_to_node.size(); ++l) {
- const NodeId& node = cell_to_node[l];
-
- StateAtNode[j][l][0] = rho_bar_lim[j](xr[node]);
- for (size_t dim = 0; dim < Dimension; ++dim)
- StateAtNode[j][l][1 + dim] = rho_u_bar_lim[j](xr[node])[dim];
- StateAtNode[j][l][1 + Dimension] = rho_E_bar_lim[j](xr[node]);
- }
- });
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtNode[j].fill(State[j]); });
EdgeValuePerCell<Rp> StateAtEdge{p_mesh->connectivity()};
StateAtEdge.fill(zero);
parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
- // eStateAtEdge[j].fill(State[j]);
- const auto& cell_to_edge = cell_to_edge_matrix[j];
-
- for (size_t l = 0; l < cell_to_edge.size(); ++l) {
- const EdgeId& edge = cell_to_edge[l];
-
- StateAtEdge[j][l][0] = rho_bar_lim[j](xe[edge]);
- for (size_t dim = 0; dim < Dimension; ++dim)
- StateAtEdge[j][l][1 + dim] = rho_u_bar_lim[j](xe[edge])[dim];
- StateAtEdge[j][l][1 + Dimension] = rho_E_bar_lim[j](xe[edge]);
- }
- });
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtEdge[j].fill(State[j]); });
FaceValuePerCell<Rp> StateAtFace{p_mesh->connectivity()};
StateAtFace.fill(zero);
parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
- // StateAtFace[j].fill(State[j]);
- const auto& cell_to_face = cell_to_face_matrix[j];
-
- for (size_t l = 0; l < cell_to_face.size(); ++l) {
- const FaceId& face = cell_to_face[l];
-
- StateAtFace[j][l][0] = rho_bar_lim[j](xl[face]);
- for (size_t dim = 0; dim < Dimension; ++dim)
- StateAtFace[j][l][1 + dim] = rho_u_bar_lim[j](xl[face])[dim];
- StateAtFace[j][l][1 + Dimension] = rho_E_bar_lim[j](xl[face]);
- }
- });
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtFace[j].fill(State[j]); });
// Conditions aux limites des etats conservatifs
--
GitLab