diff --git a/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp b/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp index f5025f2284c8125105967d9bd91ec6ae11147cba..e0e57a6856ccd912b134b1ecd7a571b4877324c3 100644 --- a/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp +++ b/src/scheme/HybridHLLRusanovEulerianCompositeSolver_v2.cpp @@ -1330,8 +1330,8 @@ class HybridHLLRusanovEulerianCompositeSolver_v2 } // dim 3 // Pour les assemblages - double theta = 2. / 3.; //.5; 2. / 3. - double eta = 1. / 6.; //.2; 1. / 6. + double theta = 1.; //.5; 2. / 3. + double eta = 0; //.2; 1. / 6. if constexpr (Dimension == 2) { eta = 0; } diff --git a/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp b/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp index 1a30528b9e6a38b63fe39841ed244cd053977f86..36dae3df541f5d10f36839a780745f610bc7dea5 100644 --- a/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp +++ b/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp @@ -1178,8 +1178,17 @@ class HybridHLLcRusanovEulerianCompositeSolver_v2 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]; + + Rd unit_normal = Cjf_loc; + unit_normal *= 1. / l2Norm(Cjf_loc); + Rd unit_tangent; + unit_tangent[0] = -unit_normal[1]; + unit_tangent[1] = unit_normal[0]; + // unit_tangent[2] = 0; + + const double uL = dot(Uj, unit_normal); + const double uR = dot(UK, unit_normal); + const double rhoEj = StateAtFace[j][l][Dimension + 1]; const double rhoEK = StateAtFace[K][R][Dimension + 1]; const double Ej = rhoEj / rhoj; @@ -1191,15 +1200,12 @@ class HybridHLLcRusanovEulerianCompositeSolver_v2 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) / l2Norm(Cjf_loc); // SL - const double MaxVpNormjk = std::max(MinMaxVpNormj.second, MinMaxVpNormk.second) / l2Norm(Cjf_loc); // SR + const double MinVpNormjk = std::min(MinMaxVpNormj.first, MinMaxVpNormk.first); + const double Sl = MinVpNormjk / l2Norm(Cjf_loc); + const double MaxVpNormjk = std::max(MinMaxVpNormj.second, MinMaxVpNormk.second); + const double Sr = MaxVpNormjk / l2Norm(Cjf_loc); - 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 + if (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) @@ -1208,7 +1214,7 @@ class HybridHLLcRusanovEulerianCompositeSolver_v2 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 (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) @@ -1216,45 +1222,65 @@ class HybridHLLcRusanovEulerianCompositeSolver_v2 Gjf[j][l][1 + Dimension] = dot(Flux_totnrjAtCellFace[K][R], Cjf_loc); // dot(Flux_totnrj[K] , Cjf_loc) } else{ + const double diffP = PressionK - Pressionj; + const double diffVelocityj = Sl - uL; + const double diffVelocityK = Sr - uR; + const double SC = (diffP + rhoj * uL * diffVelocityj - rhoK * uR * diffVelocityK)/(rhoj * diffVelocityj - rhoK * diffVelocityK); + + const double PC_R = PressionK + rhoK * diffVelocityK * (SC - uR); + const double PC_L = Pressionj + rhoj * diffVelocityj * (SC - uL); + const double PC = .5 * (PC_R + PC_L); + + Rdxd MatriceRot(identity); + MatriceRot(0, 0) = unit_normal[0]; + MatriceRot(0, 1) = unit_normal[1]; + MatriceRot(1, 0) = unit_tangent[0]; + MatriceRot(1, 1) = unit_tangent[1]; + + Rpxp MatriceRotI(identity); + MatriceRotI(0, 0) = 1; + for (size_t d1 = 0; d1 < Dimension; ++d1) + for (size_t d2 = 0; d2 < Dimension; ++d2) { + MatriceRotI(d1 + 1, d2 + 1) = MatriceRot(d1, d2); + } + MatriceRotI(Dimension + 1, Dimension + 1) = 1; + + Rp D_SC; + D_SC[0] = 0; + D_SC[1] = 1; + for (size_t d = 1; d < Dimension; ++d) + D_SC[d + 1] = 0; + D_SC[Dimension + 1] = SC; + + D_SC = transpose(MatriceRotI) * D_SC; + D_SC *= PC; + if (SC >= 0){ // SL ≤ 0 ≤ SC - 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) - Rp DL; - DL[0] = 1; - DL[1] = SC; - for (size_t d = 2; d < (Dimension + 1); ++d){ - DL[d] = Uj[d-1]; + DL[0] = rhoj; + DL[1] = rhoj * SC; + for (size_t d = 1; d < Dimension; ++d) { + DL[d + 1] = rhoj * Uj[d]; } - DL[Dimension + 1] = Ej + rhoj * (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; + DL[Dimension + 1] = rhoj * Ej + (SC - uL) * (rhoj * SC + Pressionj / diffVelocityj); + + const Rp UCL = (diffVelocityj / (Sl - SC)) * transpose(MatriceRotI) * DL; + + Gjf[j][l] = SC * UCL + D_SC; } 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) - + // SC ≤ 0 ≤ SR Rp DR; - DR[0] = 1; - DR[1] = SC; - for (size_t d = 2; d < (Dimension + 1); ++d){ - DR[d] = UK[d-1]; + DR[0] = rhoK; + DR[1] = rhoK * SC; + for (size_t d = 1; d < Dimension; ++d) { + DR[d + 1] = rhoK * UK[d]; } - DR[Dimension + 1] = EK + rhoK * (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; + DR[Dimension + 1] = rhoK * EK + (SC - uR) * (rhoK * SC + PressionK / diffVelocityK); + + const Rp UCR = (diffVelocityK / (Sr - SC)) * transpose(MatriceRotI) * DR; + + Gjf[j][l] = SC * UCR + D_SC; } } } @@ -1376,8 +1402,8 @@ class HybridHLLcRusanovEulerianCompositeSolver_v2 } // dim 3 // Pour les assemblages - double theta = 2. / 3.; //.5; 2. / 3. - double eta = 1. / 6.; //.2; 1. / 6. + double theta = 1; //.5; 2. / 3. + double eta = 0; //.2; 1. / 6. if constexpr (Dimension == 2) { eta = 0; }