Skip to content
Snippets Groups Projects
Commit 5f33da9e authored by Philippe Hoch's avatar Philippe Hoch
Browse files

adding (not yet fully implemented) induced limitation (CellbyCell)

parent a85bc537
No related branches found
No related tags found
No related merge requests found
...@@ -17,12 +17,15 @@ ...@@ -17,12 +17,15 @@
#include <mesh/MeshNodeBoundary.hpp> #include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshTraits.hpp> #include <mesh/MeshTraits.hpp>
#include <mesh/MeshVariant.hpp> #include <mesh/MeshVariant.hpp>
#include <mesh/StencilManager.hpp>
#include <scheme/DiscreteFunctionDPk.hpp> #include <scheme/DiscreteFunctionDPk.hpp>
#include <scheme/DiscreteFunctionDPkVariant.hpp> #include <scheme/DiscreteFunctionDPkVariant.hpp>
#include <scheme/DiscreteFunctionDPkVector.hpp> #include <scheme/DiscreteFunctionDPkVector.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp> #include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/InflowListBoundaryConditionDescriptor.hpp> #include <scheme/InflowListBoundaryConditionDescriptor.hpp>
#include <scheme/PolynomialReconstruction.hpp> #include <scheme/PolynomialReconstruction.hpp>
#include <scheme/PolynomialReconstructionDescriptor.hpp>
#include <utils/PugsTraits.hpp> #include <utils/PugsTraits.hpp>
#include <variant> #include <variant>
...@@ -39,6 +42,7 @@ class CellByCellLimitation ...@@ -39,6 +42,7 @@ class CellByCellLimitation
public: public:
void void
density_limiter(const MeshType& mesh, density_limiter(const MeshType& mesh,
const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
const DiscreteFunctionP0<const double>& rho, const DiscreteFunctionP0<const double>& rho,
DiscreteFunctionDPk<Dimension, double>& rho_L) const DiscreteFunctionDPk<Dimension, double>& rho_L) const
{ {
...@@ -46,16 +50,20 @@ class CellByCellLimitation ...@@ -46,16 +50,20 @@ class CellByCellLimitation
const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix(); const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
const auto& xr = mesh.xr(); const auto& xr = mesh.xr();
const auto& xl = mesh.xl(); // const auto& xl = mesh.xl();
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
auto stencil = StencilManager::instance().getStencil(mesh.connectivity(), 1); // auto stencil = StencilManager::instance()._getStencilArray(mesh.connectivity(), 1);
auto stencil = StencilManager::instance()
.getCellToCellStencilArray(mesh.connectivity(),
StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
symmetry_boundary_descriptor_list);
auto xj = mesh_data.xj(); auto xj = mesh_data.xj();
const auto& xl = mesh_data.xl();
const QuadratureFormula<1> qf = // const QuadratureFormula<1> qf =
QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree)); // QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
...@@ -81,21 +89,31 @@ class CellByCellLimitation ...@@ -81,21 +89,31 @@ class CellByCellLimitation
rho_bar_max = std::max(rho_bar_max, rho_xk); rho_bar_max = std::max(rho_bar_max, rho_xk);
} }
if constexpr (Dimension == 2) {
auto face_list = cell_to_face_matrix[cell_id]; auto face_list = cell_to_face_matrix[cell_id];
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) { for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face]; const FaceId face_id = face_list[i_face];
const Rd& x0 = xr[face_to_node_matrix[face_id][0]]; const Rd& x0 = xr[face_to_node_matrix[face_id][0]];
const Rd& x1 = xl[face_id][0];
const Rd& x2 = xr[face_to_node_matrix[face_id][1]]; const Rd& x2 = xr[face_to_node_matrix[face_id][1]];
const Rd& x1 = .5 * (x0 + x2);
const LineParabolicTransformation<Dimension> t(x0, x1, x2); const double rho_x0 = rho_L[cell_id](x0);
for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) { const double rho_x1 = rho_L[cell_id](x1);
const double rho_xk = rho_L[cell_id](t(qf.point(iq))); const double rho_x2 = rho_L[cell_id](x2);
rho_bar_min = std::min(rho_bar_min, rho_xk); rho_bar_min = std::min(rho_bar_min, std::min(rho_x0, std::min(rho_x1, rho_x2)));
rho_bar_max = std::max(rho_bar_max, rho_xk); rho_bar_max = std::max(rho_bar_max, std::max(rho_x0, std::max(rho_x1, rho_x2)));
} }
// const LineParabolicTransformation<Dimension> t(x0, x1, x2);
// for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) {
// const double rho_xk = rho_L[cell_id](t(qf.point(iq)));
// rho_bar_min = std::min(rho_bar_min, rho_xk);
// rho_bar_max = std::max(rho_bar_max, rho_xk);
//}
} else {
} }
const double eps = 1E-14; const double eps = 1E-14;
...@@ -122,11 +140,17 @@ class CellByCellLimitation ...@@ -122,11 +140,17 @@ class CellByCellLimitation
} }
void void
specific_internal_nrj_limiter(const MeshType& mesh, specific_internal_nrj_limiter(
const DiscreteFunctionP0<const double>& rho, const MeshType& mesh,
const DiscreteFunctionDPk<Dimension, double>& rho_L const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
DiscreteFunctionP0<const double>& epsilon,
DiscreteFunctionDPk<Dimension, double>& epsilon_R) const // const DiscreteFunctionP0<const double>& rho,
// const DiscreteFunctionDPk<Dimension, double>& rho_L,
const DiscreteFunctionP0<const double>& epsilon,
// const DiscreteFunctionDPk<Dimension, double>&
auto epsilon_R(const CellId cell_id, const Rd& x),
CellValue<const double>& lambda_epsilon) // const
// DiscreteFunctionDPk<Dimension, double>& epsilon_R) const
{ {
const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix(); const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix(); const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
...@@ -136,12 +160,15 @@ class CellByCellLimitation ...@@ -136,12 +160,15 @@ class CellByCellLimitation
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh); MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
auto stencil = StencilManager::instance().getStencil(mesh.connectivity(), 1); // auto stencil = StencilManager::instance()._getStencilArray(mesh.connectivity(), 1);
auto stencil = StencilManager::instance()
.getCellToCellStencilArray(mesh.connectivity(),
StencilDescriptor{1, StencilDescriptor::ConnectionType::by_nodes},
symmetry_boundary_descriptor_list);
auto xj = mesh_data.xj(); auto xj = mesh_data.xj();
const QuadratureFormula<1> qf = // const QuadratureFormula<1> qf =
QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree)); // QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
...@@ -174,14 +201,20 @@ class CellByCellLimitation ...@@ -174,14 +201,20 @@ class CellByCellLimitation
const Rd& x0 = xr[face_to_node_matrix[face_id][0]]; const Rd& x0 = xr[face_to_node_matrix[face_id][0]];
const Rd& x1 = xl[face_id][0]; const Rd& x1 = xl[face_id][0];
const Rd& x2 = xr[face_to_node_matrix[face_id][1]]; const Rd& x2 = xr[face_to_node_matrix[face_id][1]];
const double epsilon_x0 = epsilon_R(cell_id, x0);
const double epsilon_x1 = epsilon_R(cell_id, x1);
const double epsilon_x2 = epsilon_R(cell_id, x2);
const LineParabolicTransformation<Dimension> t(x0, x1, x2); epsilon_R_min = std::min(epsilon_R_min, std::min(epsilon_x0, std::min(epsilon_x1, epsilon_x2)));
for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) { epsilon_R_max = std::max(epsilon_R_max, std::max(epsilon_x0, std::max(epsilon_x1, epsilon_x2)));
const double epsilon_xk = epsilon_R(cell_id, t(qf.point(iq)));
epsilon_R_min = std::min(epsilon_R_min, epsilon_xk); // const LineParabolicTransformation<Dimension> t(x0, x1, x2);
epsilon_R_max = std::max(epsilon_R_max, epsilon_xk); // for (size_t iq = 0; iq < qf.numberOfPoints(); ++iq) {
} // const double epsilon_xk = epsilon_R(cell_id, t(qf.point(iq)));
// epsilon_R_min = std::min(epsilon_R_min, epsilon_xk);
// epsilon_R_max = std::max(epsilon_R_max, epsilon_xk);
//}
} }
const double eps = 1E-14; const double eps = 1E-14;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment