From a66fe835213e6932221601e53803fd794faa5d15 Mon Sep 17 00:00:00 2001
From: Dylan Cormet <cormet.dylan@gmail.com>
Date: Thu, 17 Apr 2025 14:56:19 +0200
Subject: [PATCH] correct HLL scheme

---
 ...idHLLRusanovEulerianCompositeSolver_v2.cpp | 66 ++++++++++++-------
 1 file changed, 44 insertions(+), 22 deletions(-)

diff --git a/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp b/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp
index da6bb44b..f5025f22 100644
--- a/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp
@@ -1,5 +1,4 @@
 #include <scheme/HybridHLLRusanovEulerianCompositeSolver_v2.hpp>
-#include <scheme/RusanovEulerianCompositeSolver_v2.hpp>
 
 #include <language/utils/InterpolateItemArray.hpp>
 #include <mesh/Mesh.hpp>
@@ -1151,43 +1150,66 @@ class HybridHLLRusanovEulerianCompositeSolver_v2
 
           const Rd& Cjf_loc = Cjf(j, l);
 
-          for (size_t k = 0; k < face_to_cell.size(); ++k) {
-            const CellId K       = face_to_cell[k];
-            const unsigned int R = face_local_number_in_its_cells[k];
-            const Rd& Ckf_loc    = Cjf(K, R);
-            // Une moyenne entre les etats jk
+          CellId K=face_to_cell[0];
+	        unsigned int R =face_local_number_in_its_cells[0];
 
-            Rd uFace     = .5 * (u_n[j] + u_n[K]);
-            double cFace = .5 * (c_n[j] + c_n[K]);
-            
-            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);
+          if (face_to_cell.size()==1)
+          {
+            K=j;
+            R=l;
+          }
+          else
+          {
+            const CellId K1       = face_to_cell[0];
+  	        const CellId K2       = face_to_cell[1];
+
+            if (j==K1)
+            {
+              K = K2;
+              R = face_local_number_in_its_cells[1];
+            }
+          }
 
-            const Rd& u_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc;   // Flux_qtmvt[K] * Cjf_loc;
+          //const Rd& Ckf_loc    = Cjf(K, R);
+          // Une moyenne entre les etats jk
 
+          //Rd uFace     = .5 * (u_n[j] + u_n[K]);
+          //double cFace = .5 * (c_n[j] + c_n[K]);
 
-            if (MinVpNormjk >= 0 ){ // SL >=0
-              // Gjf does not change
-            }
+          const std::pair<double, double> MinMaxVpNormj = toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(u_n[j], c_n[j], Cjf_loc);
+          const std::pair<double, double> MinMaxVpNormk = toolsCompositeSolver::EvaluateMinMaxEigenValueTimesNormalLengthInGivenDirection(u_n[K], c_n[K], Cjf_loc); 
+          const double MinVpNormjk = std::min(MinMaxVpNormj.first, MinMaxVpNormk.first);
+          const double MaxVpNormjk = std::max(MinMaxVpNormj.second, MinMaxVpNormk.second);
+
+          if (MinVpNormjk >= 0 ){ // SL >=0
+            const Rd& uj_Cjf = Flux_qtmvtAtCellFace[j][l] * Cjf_loc;   // Flux_qtmvt[j] * Cjf_loc;
+
+            Gjf[j][l][0] = dot(Flux_rhoAtCellFace[j][l], Cjf_loc);   // dot(Flux_roh[j] , Cjf_loc)
+            for (size_t d = 0; d < Dimension; ++d)
+              Gjf[j][l][1 + d] = uj_Cjf[d];
+            Gjf[j][l][1 + Dimension] = dot(Flux_totnrjAtCellFace[j][l], Cjf_loc);    // dot(Flux_totnrj[K] , Cjf_loc)
+          }
+          else{
             if (MaxVpNormjk <= 0 ){ // SR <= 0
+              const Rd& uk_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc;   // Flux_qtmvt[K] * Cjf_loc;
               Gjf[j][l][0] = dot(Flux_rhoAtCellFace[K][R], Cjf_loc);   // dot(Flux_roh[K] , Cjf_loc)
               for (size_t d = 0; d < Dimension; ++d)
-                Gjf[j][l][1 + d] = u_Cjf[d];
+                Gjf[j][l][1 + d] = uk_Cjf[d];
               Gjf[j][l][1 + Dimension] = dot(Flux_totnrjAtCellFace[K][R], Cjf_loc);    // dot(Flux_totnrj[K] , Cjf_loc)
             }
             else{ // SL <= 0 <= < SR
+              const Rd& uk_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc;   // Flux_qtmvt[K] * Cjf_loc;
+              const Rd& uj_Cjf = Flux_qtmvtAtCellFace[j][l] * Cjf_loc;   // Flux_qtmvt[j] * Cjf_loc;
               const Rp& statediff = StateAtFace[j][l] - StateAtFace[K][R];   // State[j] - State[K];
               const Rp& diff      =  (MinVpNormjk * MaxVpNormjk) * statediff;
 
-              Gjf[j][l][0] = MaxVpNormjk * Gjf[j][l][0] - MinVpNormjk * dot(Flux_rhoAtCellFace[K][R], Cjf_loc);   // SR * Flux_roh[j] - SL * Flux_roh[K]
+              Gjf[j][l][0] = MaxVpNormjk * dot(Flux_rhoAtCellFace[j][l], Cjf_loc) - MinVpNormjk * dot(Flux_rhoAtCellFace[K][R], Cjf_loc);   // SR * Flux_roh[j] - SL * Flux_roh[K]
               for (size_t d = 0; d < Dimension; ++d)
-                Gjf[j][l][1 + d] = MaxVpNormjk * Gjf[j][l][1 + Dimension] - MinVpNormjk * u_Cjf[d];
-              Gjf[j][l][1 + Dimension] = MaxVpNormjk * Gjf[j][l][1 + Dimension] - MinVpNormjk * dot(Flux_totnrjAtCellFace[K][R], Cjf_loc); // SR * Flux_totnrj[j] - SL * Flux_totnrj[K]
+                Gjf[j][l][1 + d] = MaxVpNormjk * uj_Cjf[d] - MinVpNormjk * uk_Cjf[d];
+              Gjf[j][l][1 + Dimension] = MaxVpNormjk * dot(Flux_totnrjAtCellFace[j][l], Cjf_loc) - MinVpNormjk * dot(Flux_totnrjAtCellFace[K][R], Cjf_loc); // SR * Flux_totnrj[j] - SL * Flux_totnrj[K]
 
               Gjf[j][l] -= diff;
-              Gjf[j][l] *= 1. / (MaxVpNormjk - MinVpNormjk); 
+              Gjf[j][l] *= 1. / (MaxVpNormjk - MinVpNormjk);
             }
           }
         }
-- 
GitLab