diff --git a/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp b/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp
index c76bd5ce7fb1b4b298e358f7472ed3163a9927cd..0f3196af245c5fa398cd2208a2bcd6db4f576cf6 100644
--- a/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/HybridHLLcRusanovEulerianCompositeSolver_v2.cpp
@@ -1150,93 +1150,86 @@ class HybridHLLcRusanovEulerianCompositeSolver_v2
const Rd& Cjf_loc = Cjf(j, l);
- CellId K=face_to_cell[0];
- unsigned int R =face_local_number_in_its_cells[0];
+ CellId K = face_to_cell[0];
+ unsigned int R = face_local_number_in_its_cells[0];
- 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 (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)
- {
+ if (j == K1) {
K = K2;
R = face_local_number_in_its_cells[1];
}
}
- 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;
- }
-
- 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;
- 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 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 Sl = MinVpNormjk / l2Norm(Cjf_loc);
+ const double Sl = MinVpNormjk / l2Norm(Cjf_loc);
const double MaxVpNormjk = std::max(MinMaxVpNormj.second, MinMaxVpNormk.second);
const double Sr = MaxVpNormjk / l2Norm(Cjf_loc);
-
- 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)
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ if (Sl >= 0) {
+ const Rd& uj_Cjf = Flux_qtmvtAtCellFace[j][l] * Cjf_loc;
+ Gjf[j][l][0] = dot(Flux_rhoAtCellFace[j][l], 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 (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)
+ Gjf[j][l][1 + Dimension] = dot(Flux_totnrjAtCellFace[j][l], Cjf_loc);
+ } else {
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ if (Sr <= 0) {
+ const Rd& uk_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc;
+ Gjf[j][l][0] = dot(Flux_rhoAtCellFace[K][R], 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{
- 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);
+ Gjf[j][l][1 + Dimension] = dot(Flux_totnrjAtCellFace[K][R], Cjf_loc);
+ } else {
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ const Rp& stateJ = StateAtFace[j][l];
+ const Rp& stateK = StateAtFace[K][R];
+
+ Rd unit_normal = Cjf_loc;
+ const double norm_Cjf_loc = l2Norm(Cjf_loc);
+ unit_normal *= 1. / norm_Cjf_loc;
+ Rd unit_tangent;
+ unit_tangent[0] = -unit_normal[1];
+ unit_tangent[1] = unit_normal[0];
+ // unit_tangent[2] = 0;
+
+ const double rhoj = stateJ[0];
+ const double rhoK = stateK[0];
+ Rd Uj;
+ Rd UK;
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ Uj[dim] = stateJ[dim + 1] / rhoj;
+ UK[dim] = stateK[dim + 1] / rhoK;
+ }
+
+ const double uL = dot(Uj, unit_normal);
+ const double uR = dot(UK, unit_normal);
- 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);
+ const double rhoEj = stateJ[Dimension + 1];
+ const double rhoEK = stateK[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);
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)
@@ -1245,42 +1238,52 @@ class HybridHLLcRusanovEulerianCompositeSolver_v2
}
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;
+ const Rp& U_left = MatriceRotI * stateJ;
+ const Rp& U_right = MatriceRotI * stateK;
- D_SC = transpose(MatriceRotI) * D_SC;
- D_SC *= PC;
+ 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);
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ if (SC >= 0) {
+ const Rd& uj_Cjf = Flux_qtmvtAtCellFace[j][l] * Cjf_loc;
+ Gjf[j][l][0] = dot(Flux_rhoAtCellFace[j][l], 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);
- if (SC >= 0){ // SL ≤ 0 ≤ SC
Rp DL;
DL[0] = rhoj;
DL[1] = rhoj * SC;
for (size_t d = 1; d < Dimension; ++d) {
- DL[d + 1] = rhoj * Uj[d];
- }
+ DL[d + 1] = U_left[d + 1];}
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 Rp U_star_L = (diffVelocityj / (Sl - SC)) * DL;
+
+ const Rp Delta_Ustar_U_L = norm_Cjf_loc * Sl * transpose(MatriceRotI) * (U_star_L - U_left);
+ Gjf[j][l] += Delta_Ustar_U_L;
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ } else {
+ const Rd& uk_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc;
+ Gjf[j][l][0] = dot(Flux_rhoAtCellFace[K][R], 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);
+
Rp DR;
DR[0] = rhoK;
DR[1] = rhoK * SC;
for (size_t d = 1; d < Dimension; ++d) {
- DR[d + 1] = rhoK * UK[d];
+ DR[d + 1] = U_right[d + 1];
}
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;
+ const Rp U_star_R = (diffVelocityK / (Sr - SC)) * DR;
+ const Rp Delta_Ustar_U_R = norm_Cjf_loc * Sr * transpose(MatriceRotI) * (U_star_R - U_right);
+ Gjf[j][l] += Delta_Ustar_U_R;
}
}
}
@@ -1402,7 +1405,7 @@ class HybridHLLcRusanovEulerianCompositeSolver_v2
} // dim 3
// Pour les assemblages
- double theta = 1; //.5; 2. / 3.
+ double theta = 2. / 3.; //.5; 2. / 3.
double eta = 0; //.2; 1. / 6.
if constexpr (Dimension == 2) {
eta = 0;