From 5893b48876f113537cad44d0e8c1f8ee9172cbda Mon Sep 17 00:00:00 2001 From: Dylan Cormet <cormet.dylan@gmail.com> Date: Mon, 28 Apr 2025 09:43:00 +0200 Subject: [PATCH] correct SC and add /* Cjf_loc --- ...dHLLcRusanovEulerianCompositeSolver_v2.cpp | 108 +++++++++++------- 1 file changed, 66 insertions(+), 42 deletions(-) diff --git a/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp b/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp index 0a7c00c6..373dc1c2 100644 --- a/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp +++ b/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] >= 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]); - } - - 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]); + 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[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; } - - 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); - } - 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]); + 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) + + Rp DR; + DR[0] = 1; + DR[1] = SC; + for (size_t d = 2; d < (Dimension + 1); ++d){ + DR[d] = UK[d-1]; + } + 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; } } } -- GitLab