From 495a411d73264a6018f8ce158f0c473fb65f5847 Mon Sep 17 00:00:00 2001
From: Dylan Cormet <cormet.dylan@gmail.com>
Date: Wed, 16 Apr 2025 14:20:21 +0200
Subject: [PATCH] correct MinVpNormjk and MaxVpNormjk implementation, add
 missing templates for EvaluateMinMax

---
 .../CompositeSchemeOtherFluxesModule.cpp      | 21 ++++++++
 src/scheme/CMakeLists.txt                     |  1 +
 ...idHLLRusanovEulerianCompositeSolver_v2.cpp | 53 +++++++++----------
 ...idHLLRusanovEulerianCompositeSolver_v2.hpp |  8 +--
 .../RusanovEulerianCompositeSolverTools.cpp   | 10 ++++
 5 files changed, 62 insertions(+), 31 deletions(-)

diff --git a/src/language/modules/CompositeSchemeOtherFluxesModule.cpp b/src/language/modules/CompositeSchemeOtherFluxesModule.cpp
index dfc28723a..b6e84023d 100644
--- a/src/language/modules/CompositeSchemeOtherFluxesModule.cpp
+++ b/src/language/modules/CompositeSchemeOtherFluxesModule.cpp
@@ -19,6 +19,8 @@
 #include <scheme/RusanovEulerianCompositeSolver_v2_order_n.hpp>
 #include <scheme/HybridRoeRusanovViscousFormEulerianCompositeSolver_v2.hpp>
 #include <scheme/HybridRusanovRoeViscousFormEulerianCompositeSolver_v2.hpp>
+#include <scheme/HybridHLLRusanovEulerianCompositeSolver_v2.hpp>
+
 
 CompositeSchemeOtherFluxesModule::CompositeSchemeOtherFluxesModule() {
 
@@ -58,6 +60,25 @@ CompositeSchemeOtherFluxesModule::CompositeSchemeOtherFluxesModule() {
                                                   bc_descriptor_list, dt);
           }
     
+          )); 
+
+    this->_addBuiltinFunction("hybrid_hll_rusanov_eulerian_composite_solver_version2",
+        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 hybridHLLRusanovEulerianCompositeSolver_v2(rho, u, E, gamma, c, p,   // degree,
+                                                  bc_descriptor_list, dt);
+          }
+    
           )); 
 }
 
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index c3d8d6050..74f878fcc 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -34,6 +34,7 @@ add_library(
   RoeViscousFormEulerianCompositeSolver_v2_o2.cpp
   HybridRoeRusanovViscousFormEulerianCompositeSolver_v2.cpp
   HybridRusanovRoeViscousFormEulerianCompositeSolver_v2.cpp
+  HybridHLLRusanovEulerianCompositeSolver_v2.cpp
 )
 
 target_link_libraries(
diff --git a/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp b/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp
index 9caaf85c6..da6bb44b7 100644
--- a/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp
@@ -1,3 +1,4 @@
+#include <scheme/HybridHLLRusanovEulerianCompositeSolver_v2.hpp>
 #include <scheme/RusanovEulerianCompositeSolver_v2.hpp>
 
 #include <language/utils/InterpolateItemArray.hpp>
@@ -18,7 +19,7 @@
 #include <variant>
 
 template <MeshConcept MeshTypeT>
-class RusanovEulerianCompositeSolver_v2
+class HybridHLLRusanovEulerianCompositeSolver_v2
 {
  private:
   using MeshType = MeshTypeT;
@@ -1158,13 +1159,11 @@ class RusanovEulerianCompositeSolver_v2
 
             Rd uFace     = .5 * (u_n[j] + u_n[K]);
             double cFace = .5 * (c_n[j] + c_n[K]);
-
-            const std::pair<double, double> (MinVpNormjk1, MaxmVpNormjk1) =
-              toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(uFace, cFace, Cjf_loc)
-            const std::pair<double, double> (MinVpNormjk2, MaxmVpNormjk2) =
-              toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(uFace, cFace, Ckf_loc);
-            const double MinVpNormjk = std::min(MinVpNormjk1, MinVpNormjk2);
-            const double MaxVpNormjk = std::max(MaxmVpNormjk1, MaxmVpNormjk2);
+            
+            const std::pair<double, double> MinMaxVpNormj = toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(uFace, cFace, Cjf_loc);
+            const std::pair<double, double> MinMaxVpNormk = toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(uFace, cFace, Ckf_loc);
+            const double MinVpNormjk = std::min(MinMaxVpNormj.first, MinMaxVpNormk.first);
+            const double MaxVpNormjk = std::max(MinMaxVpNormj.second, MinMaxVpNormk.second);
 
             const Rd& u_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc;   // Flux_qtmvt[K] * Cjf_loc;
 
@@ -1362,17 +1361,17 @@ class RusanovEulerianCompositeSolver_v2
                            std::make_shared<const DiscreteFunctionVariant>(E));
   }
 
-  RusanovEulerianCompositeSolver_v2()  = default;
-  ~RusanovEulerianCompositeSolver_v2() = default;
+  HybridHLLRusanovEulerianCompositeSolver_v2()  = default;
+  ~HybridHLLRusanovEulerianCompositeSolver_v2() = default;
 };
 
 template <MeshConcept MeshType>
-class RusanovEulerianCompositeSolver_v2<MeshType>::WallBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<MeshType>::WallBoundaryCondition
 {
 };
 
 template <>
-class RusanovEulerianCompositeSolver_v2<Mesh<2>>::WallBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<Mesh<2>>::WallBoundaryCondition
 {
   using Rd = TinyVector<Dimension, double>;
 
@@ -1415,7 +1414,7 @@ class RusanovEulerianCompositeSolver_v2<Mesh<2>>::WallBoundaryCondition
 };
 
 template <>
-class RusanovEulerianCompositeSolver_v2<Mesh<3>>::WallBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<Mesh<3>>::WallBoundaryCondition
 {
   using Rd = TinyVector<Dimension, double>;
 
@@ -1474,11 +1473,11 @@ class RusanovEulerianCompositeSolver_v2<Mesh<3>>::WallBoundaryCondition
 };
 
 template <MeshConcept MeshType>
-class RusanovEulerianCompositeSolver_v2<MeshType>::NeumannflatBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<MeshType>::NeumannflatBoundaryCondition
 {
 };
 template <>
-class RusanovEulerianCompositeSolver_v2<Mesh<2>>::NeumannflatBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<Mesh<2>>::NeumannflatBoundaryCondition
 {
  public:
   using Rd = TinyVector<Dimension, double>;
@@ -1529,7 +1528,7 @@ class RusanovEulerianCompositeSolver_v2<Mesh<2>>::NeumannflatBoundaryCondition
 };
 
 template <>
-class RusanovEulerianCompositeSolver_v2<Mesh<3>>::NeumannflatBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<Mesh<3>>::NeumannflatBoundaryCondition
 {
  public:
   using Rd = TinyVector<Dimension, double>;
@@ -1596,12 +1595,12 @@ class RusanovEulerianCompositeSolver_v2<Mesh<3>>::NeumannflatBoundaryCondition
 };
 
 template <MeshConcept MeshType>
-class RusanovEulerianCompositeSolver_v2<MeshType>::SymmetryBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<MeshType>::SymmetryBoundaryCondition
 {
 };
 
 template <>
-class RusanovEulerianCompositeSolver_v2<Mesh<2>>::SymmetryBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<Mesh<2>>::SymmetryBoundaryCondition
 {
  public:
   using Rd = TinyVector<Dimension, double>;
@@ -1652,7 +1651,7 @@ class RusanovEulerianCompositeSolver_v2<Mesh<2>>::SymmetryBoundaryCondition
 };
 
 template <>
-class RusanovEulerianCompositeSolver_v2<Mesh<3>>::SymmetryBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<Mesh<3>>::SymmetryBoundaryCondition
 {
  public:
   using Rd = TinyVector<Dimension, double>;
@@ -1719,12 +1718,12 @@ class RusanovEulerianCompositeSolver_v2<Mesh<3>>::SymmetryBoundaryCondition
 };
 
 template <MeshConcept MeshType>
-class RusanovEulerianCompositeSolver_v2<MeshType>::InflowListBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<MeshType>::InflowListBoundaryCondition
 {
 };
 
 template <>
-class RusanovEulerianCompositeSolver_v2<Mesh<2>>::InflowListBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<Mesh<2>>::InflowListBoundaryCondition
 {
  public:
   using Rd = TinyVector<Dimension, double>;
@@ -1788,7 +1787,7 @@ class RusanovEulerianCompositeSolver_v2<Mesh<2>>::InflowListBoundaryCondition
 };
 
 template <>
-class RusanovEulerianCompositeSolver_v2<Mesh<3>>::InflowListBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<Mesh<3>>::InflowListBoundaryCondition
 {
  public:
   using Rd = TinyVector<Dimension, double>;
@@ -1876,12 +1875,12 @@ class RusanovEulerianCompositeSolver_v2<Mesh<3>>::InflowListBoundaryCondition
 };
 
 template <MeshConcept MeshType>
-class RusanovEulerianCompositeSolver_v2<MeshType>::OutflowBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<MeshType>::OutflowBoundaryCondition
 {
 };
 
 template <>
-class RusanovEulerianCompositeSolver_v2<Mesh<2>>::OutflowBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<Mesh<2>>::OutflowBoundaryCondition
 {
   using Rd = TinyVector<Dimension, double>;
 
@@ -1922,7 +1921,7 @@ class RusanovEulerianCompositeSolver_v2<Mesh<2>>::OutflowBoundaryCondition
 };
 
 template <>
-class RusanovEulerianCompositeSolver_v2<Mesh<3>>::OutflowBoundaryCondition
+class HybridHLLRusanovEulerianCompositeSolver_v2<Mesh<3>>::OutflowBoundaryCondition
 {
   using Rd = TinyVector<Dimension, double>;
 
@@ -1983,7 +1982,7 @@ class RusanovEulerianCompositeSolver_v2<Mesh<3>>::OutflowBoundaryCondition
 std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
            std::shared_ptr<const DiscreteFunctionVariant>,
            std::shared_ptr<const DiscreteFunctionVariant>>
-rusanovEulerianCompositeSolver_v2(
+hybridHLLRusanovEulerianCompositeSolver_v2(
   const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
   const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
   const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
@@ -2016,7 +2015,7 @@ rusanovEulerianCompositeSolver_v2(
           throw NormalError("RusanovEulerianCompositeSolver v2 is not available in 1D");
         } else {
           if constexpr (is_polygonal_mesh_v<MeshType>) {
-            return RusanovEulerianCompositeSolver_v2<MeshType>{}
+            return HybridHLLRusanovEulerianCompositeSolver_v2<MeshType>{}
               .solve(p_mesh, rho_v->get<DiscreteFunctionP0<const double>>(), u_v->get<DiscreteFunctionP0<const Rd>>(),
                      E_v->get<DiscreteFunctionP0<const double>>(), gamma, c_v->get<DiscreteFunctionP0<const double>>(),
                      p_v->get<DiscreteFunctionP0<const double>>(), bc_descriptor_list, dt, check);
diff --git a/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.hpp b/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.hpp
index 834fe05e4..f1541b87b 100644
--- a/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.hpp
+++ b/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.hpp
@@ -1,5 +1,5 @@
-#ifndef RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
-#define RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+#ifndef HYBRID_HLL_RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+#define HYBRID_HLL_RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
 
 #include <mesh/MeshVariant.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
@@ -15,7 +15,7 @@
 std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
            std::shared_ptr<const DiscreteFunctionVariant>,
            std::shared_ptr<const DiscreteFunctionVariant>>
-rusanovEulerianCompositeSolver_v2(
+hybridHLLRusanovEulerianCompositeSolver_v2(
   const std::shared_ptr<const DiscreteFunctionVariant>& rho,
   const std::shared_ptr<const DiscreteFunctionVariant>& u,
   const std::shared_ptr<const DiscreteFunctionVariant>& E,
@@ -27,4 +27,4 @@ rusanovEulerianCompositeSolver_v2(
   const double& dt,
   const bool check = false);
 
-#endif   // RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+#endif   // HYBRID_HLL_RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
diff --git a/src/scheme/RusanovEulerianCompositeSolverTools.cpp b/src/scheme/RusanovEulerianCompositeSolverTools.cpp
index 37ec39c59..85ccb357b 100644
--- a/src/scheme/RusanovEulerianCompositeSolverTools.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolverTools.cpp
@@ -109,3 +109,13 @@ template double toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGi
   const TinyVector<3>& U_mean,
   const double& c_mean,
   const TinyVector<3>& normal);
+
+template std::pair<double, double> toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(
+  const TinyVector<2>& U_mean,
+  const double& c_mean,
+  const TinyVector<2>& normal);
+
+template std::pair<double, double> toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(
+  const TinyVector<3>& U_mean,
+  const double& c_mean,
+  const TinyVector<3>& normal);
-- 
GitLab