Intro of any order induced limitation

src/scheme/CellbyCellLimitation.hpp

@@ -37,6 +37,168 @@ class CellByCellLimitation
   using Rd   = TinyVector<Dimension>;
 
  public:
+  void
+  density_limiter(const MeshType& mesh,
+                  const DiscreteFunctionP0<const double>& rho,
+                  DiscreteFunctionDPk<Dimension, double>& rho_L) 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().getStencil(mesh.connectivity(), 1);
+
+    auto xj = mesh_data.xj();
+
+    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);
+        }
+
+        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 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);
+          }
+        }
+
+        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)));
+
+        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
+  specific_internal_nrj_limiter(const MeshType& mesh,
+                                const DiscreteFunctionP0<const double>& rho,
+                                const DiscreteFunctionDPk<Dimension, double>& rho_L
+                                  DiscreteFunctionP0<const double>& epsilon,
+                                DiscreteFunctionDPk<Dimension, double>& epsilon_R) 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().getStencil(mesh.connectivity(), 1);
+
+    auto xj = mesh_data.xj();
+
+    const QuadratureFormula<1> qf =
+      QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_quadrature_degree));
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(const 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 LineParabolicTransformation<Dimension> t(x0, x1, x2);
+          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;
+        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)));
+      });
+  }
+
   void
   computeLimitorVolumicScalarQuantityMinModDukowiczGradient(const MeshType& mesh,
                                                             const CellValue<double>& q,