Adding Roe composite under Flux Form 2D and 3D

src/scheme/CMakeLists.txt

@@ -28,6 +28,7 @@ add_library(
   RusanovEulerianCompositeSolver.cpp
   RusanovEulerianCompositeSolver_v2.cpp
   RoeViscousFormEulerianCompositeSolver_v2.cpp
+  RoeFluxFormEulerianCompositeSolver_v2.cpp
   RusanovEulerianCompositeSolver_o2.cpp
   RusanovEulerianCompositeSolver_v2_o2.cpp
   RusanovEulerianCompositeSolver_v2_order_n.cpp

src/scheme/CellbyCellLimitation.hpp

 #ifndef CELL_BY_CELL_LIMITATION_HPP
 #define CELL_BY_CELL_LIMITATION_HPP
 
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
-
 #include <language/utils/InterpolateItemArray.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
@@ -24,6 +21,7 @@
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/InflowListBoundaryConditionDescriptor.hpp>
+#include <scheme/LimitationTools.hpp>
 #include <scheme/PolynomialReconstruction.hpp>
 #include <scheme/PolynomialReconstructionDescriptor.hpp>
 #include <utils/PugsTraits.hpp>
@@ -44,7 +42,8 @@ class CellByCellLimitation
   density_limiter(const MeshType& mesh,
                   const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
                   const DiscreteFunctionP0<const double>& rho,
-                  DiscreteFunctionDPk<Dimension, double>& rho_L) const
+                  DiscreteFunctionDPk<Dimension, double>& rho_L,
+                  const bool enableWeakBoundPositivityOnly = false) const
   {
     const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
     const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
@@ -95,9 +94,13 @@ class CellByCellLimitation
           for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
             const FaceId face_id = face_list[i_face];
 
-            const Rd& x0 = xr[face_to_node_matrix[face_id][0]];
-            const Rd& x2 = xr[face_to_node_matrix[face_id][1]];
-            const Rd& x1 = .5 * (x0 + x2);
+            const auto& face_to_node = face_to_node_matrix[face_id];
+            const NodeId& nodei      = face_to_node[0];
+            const NodeId& nodef      = face_to_node[face_to_node.size() - 1];
+
+            const Rd x0 = xr[nodei];
+            const Rd x1 = xl[face_id];
+            const Rd x2 = xr[nodef];
 
             const double rho_x0 = rho_L[cell_id](x0);
             const double rho_x1 = rho_L[cell_id](x1);
@@ -106,16 +109,120 @@ class CellByCellLimitation
             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, 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 {
+          throw NotImplementedError("not implement in 3D");
+        }
+        /*
+        const double eps = 1E-14;
+        double coef1     = 1;
+        if (std::abs(rho_bar_max - rhoj) > eps) {
+          coef1 = (rho_max - rhoj) / ((rho_bar_max - rhoj));
+        }
+
+        double coef2 = 1.;
+        if (std::abs(rho_bar_min - rhoj) > eps) {
+          coef2 = (rho_min - rhoj) / ((rho_bar_min - rhoj));
+        }
+*/
+        const double lambda =   // std::max(0., std::min(1., std::min(coef1, coef2)));
+          toolsLimitation::computeCoefLimitationBarthJespersen(rhoj, rho_bar_min, rho_bar_max, rho_min, rho_max,
+                                                               enableWeakBoundPositivityOnly);
+
+        auto coefficients = rho_L.coefficients(cell_id);
+
+        coefficients[0] = (1 - lambda) * rho[cell_id] + lambda * coefficients[0];
+
+        for (size_t i = 1; i < coefficients.size(); ++i) {
+          coefficients[i] *= lambda;
+        }
+      });
+  }
+
+  void
+  density_limiter(const MeshType& mesh,
+                  const FaceArray<TinyVector<2>>& QuadratureFace,
+                  const EdgeArray<TinyVector<2>>& QuadratureEdge,
+                  const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
+                  const DiscreteFunctionP0<const double>& rho,
+                  DiscreteFunctionDPk<Dimension, double>& rho_L,
+                  const bool enableWeakBoundPositivityOnly = false) const
+  {
+    const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+
+    const auto& xr = mesh.xr();
+    // const auto& xl = mesh.xl();
+
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    //    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();
+    const auto& xl = mesh_data.xl();
+
+    // const QuadratureFormula<1> qf =
+    //   QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        const double rhoj = rho[cell_id];
+
+        double rho_min = rhoj;
+        double rho_max = rhoj;
+
+        const auto cell_stencil = stencil[cell_id];
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          rho_min = std::min(rho_min, rho[cell_stencil[i_cell]]);
+          rho_max = std::max(rho_max, rho[cell_stencil[i_cell]]);
         }
 
+        double rho_bar_min = rhoj;
+        double rho_bar_max = rhoj;
+
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          const CellId cell_k_id = cell_stencil[i_cell];
+          const double rho_xk    = rho_L[cell_id](xj[cell_k_id]);
+
+          rho_bar_min = std::min(rho_bar_min, 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];
+
+          for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+            const FaceId face_id = face_list[i_face];
+
+            const auto& face_to_node = face_to_node_matrix[face_id];
+            const NodeId& nodei      = face_to_node[0];
+            const NodeId& nodef      = face_to_node[face_to_node.size() - 1];
+
+            const Rd& x0 = xr[nodei];
+            const Rd& x1 = xr[nodef];
+
+            const double rho_x0 = rho_L[cell_id](x0);
+            const double rho_x1 = rho_L[cell_id](x1);
+
+            rho_bar_min = std::min(rho_bar_min, std::min(rho_x0, rho_x1));
+            rho_bar_max = std::min(rho_bar_max, std::max(rho_x0, rho_x1));
+
+            for (size_t i = 0; i < QuadratureFace[face_id].size(); ++i) {
+              // const Rd& x_quad = xd + QuadratureFace[face][i][1] * (xf - xd);
+              const Rd& x_quad = xr[nodei] + QuadratureFace[face_id][i][1] * (xr[nodef] - xr[nodei]);
+
+              const double rho_xk = rho_L[cell_id](x_quad);
+
+              rho_bar_min = std::min(rho_bar_min, rho_xk);
+              rho_bar_max = std::max(rho_bar_max, rho_xk);
+            }
+          }
+        } else {
+          throw NotImplementedError("not implement in 3D");
+        }
+        /*
         const double eps = 1E-14;
         double coef1     = 1;
         if (std::abs(rho_bar_max - rhoj) > eps) {
@@ -126,8 +233,11 @@ class CellByCellLimitation
         if (std::abs(rho_bar_min - rhoj) > eps) {
           coef2 = (rho_min - rhoj) / ((rho_bar_min - rhoj));
         }
+*/
+        const double lambda =   // std::max(0., std::min(1., std::min(coef1, coef2)));
 
-        const double lambda = std::max(0., std::min(1., std::min(coef1, coef2)));
+          toolsLimitation::computeCoefLimitationBarthJespersen(rhoj, rho_bar_min, rho_bar_max, rho_min, rho_max,
+                                                               enableWeakBoundPositivityOnly);
 
         auto coefficients = rho_L.coefficients(cell_id);
 
@@ -139,26 +249,27 @@ class CellByCellLimitation
       });
   }
 
+  //
+  // Idem for epsilon
+  //
+
   void
   specific_internal_nrj_limiter(
     const MeshType& mesh,
     const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
-
-    // 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<double>& lambda_epsilon)   // const
-  //                                DiscreteFunctionDPk<Dimension, double>& epsilon_R) const
+    auto epsilon_R,   //(const CellId cell_id, const Rd& x),
+    CellValue<double>& lambda_epsilon,
+    const bool enableWeakBoundPositivityOnly = false) const
   {
     const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
     const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
 
     const auto& xr = mesh.xr();
-    const auto& xl = mesh.xl();
+    // const auto& xl = mesh.xl();
 
     MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    const auto& xl                = mesh_data.xl();
 
     //    auto stencil = StencilManager::instance()._getStencilArray(mesh.connectivity(), 1);
     auto stencil = StencilManager::instance()
@@ -216,7 +327,7 @@ class CellByCellLimitation
           // epsilon_R_max = std::max(epsilon_R_max, epsilon_xk);
           //}
         }
-
+        /*
         const double eps = 1E-14;
         double coef1     = 1;
         if (std::abs(epsilon_R_max - epsilonj) > eps) {
@@ -227,11 +338,114 @@ class CellByCellLimitation
         if (std::abs(epsilon_R_min - epsilonj) > eps) {
           coef2 = (epsilon_min - epsilonj) / ((epsilon_R_min - epsilonj));
         }
+*/
+        lambda_epsilon[cell_id] =   // std::max(0., std::min(1., std::min(coef1, coef2)));
+          toolsLimitation::computeCoefLimitationBarthJespersen(epsilonj, epsilon_R_min, epsilon_R_min, epsilon_min,
+                                                               epsilon_max, enableWeakBoundPositivityOnly);
+      });
+  }
+
+  //
+  //
+  //
+  void
+  specific_internal_nrj_limiter(
+    const MeshType& mesh,
+    const FaceArray<TinyVector<2>>& QuadratureFace,
+    const EdgeArray<TinyVector<2>>& QuadratureEdge,
+    const std::vector<std::shared_ptr<const IBoundaryDescriptor>>& symmetry_boundary_descriptor_list,
+    const DiscreteFunctionP0<const double>& epsilon,
+    auto epsilon_R,   //(const CellId cell_id, const Rd& x),
+    CellValue<double>& lambda_epsilon,
+    const bool enableWeakBoundPositivityOnly = false) const
+  {
+    const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+
+    const auto& xr = mesh.xr();
+    // const auto& xl = mesh.xl();
 
-        lambda_epsilon[cell_id] = std::max(0., std::min(1., std::min(coef1, coef2)));
+    MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    //    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();
+
+    //    const QuadratureFormula<1> qf =
+    // QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
+        const double epsilonj = epsilon[cell_id];
+
+        double epsilon_min = epsilonj;
+        double epsilon_max = epsilonj;
+
+        const auto cell_stencil = stencil[cell_id];
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          epsilon_min = std::min(epsilon_min, epsilon[cell_stencil[i_cell]]);
+          epsilon_max = std::max(epsilon_max, epsilon[cell_stencil[i_cell]]);
+        }
+
+        double epsilon_R_min = epsilonj;
+        double epsilon_R_max = epsilonj;
+
+        for (size_t i_cell = 0; i_cell < cell_stencil.size(); ++i_cell) {
+          const CellId cell_k_id  = cell_stencil[i_cell];
+          const double epsilon_xk = epsilon_R(cell_id, xj[cell_k_id]);
+
+          epsilon_R_min = std::min(epsilon_R_min, epsilon_xk);
+          epsilon_R_max = std::max(epsilon_R_max, epsilon_xk);
+        }
+
+        auto face_list = cell_to_face_matrix[cell_id];
+        for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+          const FaceId face_id = face_list[i_face];
+
+          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 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);
+          epsilon_R_min           = std::min(epsilon_R_min, std::min(epsilon_x0, epsilon_x2));
+          epsilon_R_max           = std::max(epsilon_R_max, std::max(epsilon_x0, epsilon_x2));
+
+          for (size_t i = 0; i < QuadratureFace[face_id].size(); ++i) {
+            // const Rd& x_quad = xd + QuadratureFace[face][i][1] * (xf - xd);
+            const Rd& x_quad = x0 + QuadratureFace[face_id][i][1] * (x2 - x0);
+
+            const double epsilon_xk = epsilon_R(cell_id, x_quad);
+
+            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;
+        double coef1     = 1;
+        if (std::abs(epsilon_R_max - epsilonj) > eps) {
+          coef1 = (epsilon_max - epsilonj) / ((epsilon_R_max - epsilonj));
+        }
+
+        double coef2 = 1.;
+        if (std::abs(epsilon_R_min - epsilonj) > eps) {
+          coef2 = (epsilon_min - epsilonj) / ((epsilon_R_min - epsilonj));
+        }
+*/
+        lambda_epsilon[cell_id] =   // std::max(0., std::min(1., std::min(coef1, coef2)));
+          toolsLimitation::computeCoefLimitationBarthJespersen(epsilonj, epsilon_R_min, epsilon_R_max, epsilon_min,
+                                                               epsilon_max, enableWeakBoundPositivityOnly);
       });
   }
 
+  //
+  // Other version (CellValue Coef for limitation )
+  //
+
   void
   computeLimitorVolumicScalarQuantityMinModDukowiczGradient(const MeshType& mesh,
                                                             const CellValue<double>& q,

src/scheme/RoeFluxFormEulerianCompositeSolver_v2.hpp

+#ifndef ROE_FLUX_FORM_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+#define ROE_FLUX_FORM_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+
+#include <mesh/MeshTraits.hpp>
+
+#include "JacobianAndStructuralInfoForSystemsofEquations.hpp"
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+
+#include <memory>
+#include <tuple>
+#include <vector>
+
+inline double
+signe(const double a)
+{
+  if (a < 0)
+    return -1.;
+  if (a > 0)
+    return 1;
+
+  return 0;
+}
+class Rp;
+class Rpxp;
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>,
+           std::shared_ptr<const DiscreteFunctionVariant>>
+roeFluxFormEulerianCompositeSolver_v2(
+  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,
+  const bool check = false);
+/*
+class DiscreteFunctionVariant;
+class IBoundaryConditionDescriptor;
+class MeshVariant;
+template <size_t Dimension>
+class Mesh;
+
+template <size_t Dimension, int number_of_unknowns>
+class RoeFluxFormEulerianCompositeSolver_v2
+{
+ private:
+  static constexpr int p = number_of_unknowns + (Dimension - 1);
+  using Rp               = TinyVector<p>;   // number_of_unknowns + (Dimension - 1)>;
+
+  using Rpxp = TinyMatrix<p, p, double>;
+  // private:
+ public:
+  JacobianStructuralInfoSystemsOfEquations<Dimension, number_of_unknowns> JacobianInfos;
+
+  // KN<std::vector<std::vector<Rn> > > _sign_jacobian_matrix_dof_in_cell;
+  // KN<std::vector<std::vector<Rn> > > _vp_jacobian_matrix_dof_in_cell;
+  // //KN<std::vector<Rnm> > _jac_jacobian_matrix_dof_in_cell;
+  // KN<std::vector<std::vector<Rnm> > > _right_eigenv_matrix_dof_in_cell;
+  // KN<std::vector<std::vector<Rnm> > >  _left_eigenv_matrix_dof_in_cell;
+  // //KN<std::vector<Rnm> >  _rotation_matrix_dof_in_cell;
+  // KN<std::vector<R2nm> >  _MatriceU_dof_in_cell;
+
+  NodeValuePerCell<const Rp> _vp_jacobian_matrix_node_in_cell;
+  EdgeValuePerCell<const Rp> _vp_jacobian_matrix_edge_in_cell;
+  FaceValuePerCell<const Rp> _vp_jacobian_matrix_face_in_cell;
+
+  NodeValuePerCell<const Rpxp> _right_eigenv_matrix_node_in_cell;
+  EdgeValuePerCell<const Rpxp> _right_eigenv_matrix_edge_in_cell;
+  FaceValuePerCell<const Rpxp> _right_eigenv_matrix_face_in_cell;
+
+  NodeValuePerCell<const Rpxp> _left_eigenv_matrix_node_in_cell;
+  EdgeValuePerCell<const Rpxp> _left_eigenv_matrix_edge_in_cell;
+  FaceValuePerCell<const Rpxp> _left_eigenv_matrix_face_in_cell;
+
+  NodeValuePerCell<const Rpxp> _MatriceU_node_in_cell;
+  EdgeValuePerCell<const Rpxp> _MatriceU_edge_in_cell;
+  FaceValuePerCell<const Rpxp> _MatriceU_face_in_cell;
+
+ public:
+  RoeFluxFormEulerianCompositeSolver_v2() {}
+
+  ~RoeFluxFormEulerianCompositeSolver_v2() {}
+};
+*/
+#endif