correct SC and add /* Cjf_loc

src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp

@@ -1170,67 +1170,91 @@ class HybridHLLcRusanovEulerianCompositeSolver_v2
             }
           }
 
+          const double rhoj = StateAtFace[j][l][0];
+          const double rhoK = StateAtFace[K][R][0];
+          Rd Uj;
+          Rd UK;
+          for (size_t dim = 0; dim < Dimension; ++dim){
+            Uj[dim] = StateAtFace[j][l][dim + 1] / rhoj;
+            UK[dim] = StateAtFace[K][R][dim + 1] / rhoK;
+          }
+          const double uL = Uj[0];
+          const double uR = UK[0];
+          const double rhoEj = StateAtFace[j][l][Dimension + 1];
+          const double rhoEK = StateAtFace[K][R][Dimension + 1];
+          const double Ej        = rhoEj / rhoj;
+          const double EK        = rhoEK / rhoK;
+          const double epsilonj  = Ej - .5 * dot(Uj, Uj);
+          const double epsilonK  = EK - .5 * dot(UK, UK);
+          const double Pressionj = pression(rhoj, epsilonj, gamma);
+          const double PressionK = pression(rhoK, epsilonK, gamma);
+          
           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);   // SL
-          const double MaxVpNormjk = std::max(MinMaxVpNormj.second, MinMaxVpNormk.second); //  SR
-          Rp SC = zero;
-
-          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[K][R] - StateAtFace[j][l];   // State[K] - State[j];
-          const Rp& diff      =  MinVpNormjk * statediff;
+          const double MinVpNormjk = std::min(MinMaxVpNormj.first, MinMaxVpNormk.first) / l2Norm(Cjf_loc);   // SL
+          const double MaxVpNormjk = std::max(MinMaxVpNormj.second, MinMaxVpNormk.second) / l2Norm(Cjf_loc); //  SR
 
-          SC[0] = (dot(Flux_rhoAtCellFace[K][R], Cjf_loc) - dot(Flux_rhoAtCellFace[j][l], Cjf_loc) + diff[1 + Dimension]) / diff[0];
-          for (size_t d = 0; d < Dimension; ++d)
-            SC[1 + d] = (uk_Cjf[d] - uj_Cjf[d] + diff[d]) / diff[d];
-          SC[1 + Dimension] = (dot(Flux_totnrjAtCellFace[K][R], Cjf_loc) - dot(Flux_totnrjAtCellFace[j][l], Cjf_loc) + diff[1 + Dimension])  / diff[1 + Dimension];
+          const double diffP = PressionK - Pressionj;
+          const double diffVelocityj = MinVpNormjk - uL;
+          const double diffVelocityK = MaxVpNormjk - uR;
+          const double SC = (diffP + rhoj * uL * diffVelocityj - rhoK * uR * diffVelocityK)/(rhoj * diffVelocityj - rhoK * diffVelocityK); // SC
 
+          if (MinVpNormjk >= 0 ){ // SL ≥ 0
+            const Rd& uj_Cjf = Flux_qtmvtAtCellFace[j][l] * Cjf_loc;   // Flux_qtmvt[j] * Cjf_loc;
 
-          if (MinVpNormjk >= 0 ){ // SL >=0
             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
+            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] = uk_Cjf[d];
               Gjf[j][l][1 + Dimension] = dot(Flux_totnrjAtCellFace[K][R], Cjf_loc);    // dot(Flux_totnrj[K] , Cjf_loc)
             }
             else{
-              Rp diffL = zero;
-              Rp diffR = zero;
+              if (SC >= 0){ // SL ≤ 0 ≤ SC
+                const Rd& uj_Cjf = Flux_qtmvtAtCellFace[j][l] * Cjf_loc;   // Flux_qtmvt[j] * Cjf_loc;
     
-              if (SC[0] >= 0 ){// SL <= 0 <= SC
-                diffL[0] = MinVpNormjk * SC[0] * statediff[0];
-                Gjf[j][l][0] = (SC[0] * dot(Flux_rhoAtCellFace[j][l], Cjf_loc) - MinVpNormjk * dot(Flux_rhoAtCellFace[K][R], Cjf_loc) + diffL[0]) / (SC[0] - MinVpNormjk);
-              }
-              else{ // SC <= 0 <= SR
-                diffR[0] = SC[0] * MaxVpNormjk * statediff[0];
-                Gjf[j][l][0] = (MaxVpNormjk * dot(Flux_rhoAtCellFace[j][l], Cjf_loc) - SC[0] * dot(Flux_rhoAtCellFace[K][R], Cjf_loc) + diffR[0]) / (MaxVpNormjk - SC[0]);
-              }
+                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)
     
-              for (size_t d = 0; d < Dimension; ++d){
-                if (SC[1 + d] >= 0 ){// SL <= 0 <= SC
-                  diffL[1 + d] = MinVpNormjk * SC[1 + d] * statediff[1 + d];
-                  Gjf[j][l][1 + d] = (SC[1 + d] * uj_Cjf[d] - MinVpNormjk * uk_Cjf[d] + diffL[d])  / (SC[1 + d] - MinVpNormjk);
-                }
-                else{ // SC <= 0 <= SR
-                  diffR[1 + d] = SC[1 + d] * MaxVpNormjk * statediff[1 + d];
-                  Gjf[j][l][1 + d] = (MaxVpNormjk * uj_Cjf[d] - SC[1 + d] * uk_Cjf[d] + diffR[d])  / (MaxVpNormjk - SC[1 + d]);
+                Rp DL;
+                DL[0] = 1;
+                DL[1] = SC;
+                for (size_t d = 2; d < (Dimension + 1); ++d){
+                  DL[d] = Uj[d-1];
                 }
+                DL[Dimension + 1] = Ej + (SC - uL) * (SC + (Pressionj / (rhoj * diffVelocityj)));
+                  
+                const Rp UCL = rhoj * (diffVelocityj / (MinVpNormjk - SC)) * DL;
+                const Rp diffStates = MinVpNormjk * l2Norm(Cjf_loc) * (UCL - StateAtFace[j][l]);
+                Gjf[j][l] += diffStates;
               }
+              else{
+                // SC ≤ 0 ≤ SR
+                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] = uk_Cjf[d];
+                Gjf[j][l][1 + Dimension] = dot(Flux_totnrjAtCellFace[K][R], Cjf_loc);    // dot(Flux_totnrj[K] , Cjf_loc)
     
-              if (SC[1 + Dimension] >= 0 ){// SL <= 0 <= SC
-                diffL[1 + Dimension] = MinVpNormjk * SC[1 + Dimension] * statediff[1 + Dimension];
-                Gjf[j][l][1 + Dimension] = (SC[1 + Dimension] * dot(Flux_totnrjAtCellFace[j][l], Cjf_loc) - MinVpNormjk * dot(Flux_totnrjAtCellFace[K][R], Cjf_loc) + diffL[1 + Dimension]) / (SC[1 + Dimension] - MinVpNormjk);
+                Rp DR;
+                DR[0] = 1;
+                DR[1] = SC;
+                for (size_t d = 2; d < (Dimension + 1); ++d){
+                  DR[d] = UK[d-1];            
                 }
-              else{ // SC <= 0 <= SR
-                diffR[1 + Dimension] = SC[1 + Dimension] * MaxVpNormjk * statediff[1 + Dimension];
-                Gjf[j][l][1 + Dimension] = (MaxVpNormjk * dot(Flux_totnrjAtCellFace[j][l], Cjf_loc) - SC[1 + Dimension] * dot(Flux_totnrjAtCellFace[K][R], Cjf_loc) + diffR[1 + Dimension]) / (MaxVpNormjk - SC[1 + Dimension]); 
+                DR[Dimension + 1] = EK + (SC - uR) * (SC + (PressionK / (rhoK * diffVelocityK)));
+                    
+                const Rp UCR = rhoK * (diffVelocityK / (MaxVpNormjk - SC)) * DR;
+                const Rp diffStates = MaxVpNormjk * l2Norm(Cjf_loc) * (UCR - StateAtFace[K][R]);
+                Gjf[j][l] += diffStates;
               }
             }
           }