diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp index 7650806336ffd4394a597b46463d8341aa80bb0f..b2c4d1deb242be0363af2662967b2cffe97a7a32 100644 --- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp +++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp @@ -277,9 +277,31 @@ class RoeViscousFormEulerianCompositeSolver_v2 Rpxp RightEigenVectors; Rp EigenValues; Rpxp AbsJacobian; - double MaxErreurProd; + double maxabsVpNormal; + // bool changingSign; + double MaxErrorProd; }; + Rp + EvaluateEigenValuesNormal( // const double rhod, // rhoJ, uJ, EJ, gammaJ, cJ, pJ, normal + const Rd& Ud, + // const double Ed, + // const double gammad, + const double cd, + // const double pd, + const Rd& normal) const + { + Rp Eigen; + const double uscaln = dot(Ud, normal); + + Eigen[0] = uscaln - cd; + for (size_t dim = 0; dim < Dimension; ++dim) + Eigen[dim] = uscaln; + Eigen[1 + Dimension] = uscaln + cd; + + return Eigen; + } + struct RoeAverageStateStructData { double rho; @@ -299,9 +321,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 : rho(rhod), U(Ud), E(Ed), H(Hd), gamma(gammad), p(pd), c(cd) {} }; - // Rp RoeAverageState(const Rp& val1, const Rp& val2 ) - // Rp RoeAverageState(const Rp& val1, const Rp& val2 ) - // Rp + RoeAverageStateStructData RoeAverageState(const double& rhoG, const Rd& UG, @@ -338,17 +358,63 @@ class RoeViscousFormEulerianCompositeSolver_v2 return RoeAverageStateStructData(rho_mean, U_mean, E_mean, H_mean, gamma, P_mean, c_mean); } + bool + EvaluateChangingSignVpAlongNormal( // const double rhoJ, + const Rd& uJ, + // const double EJ, + // const double gammaJ, + const double cJ, + // const double pJ, + + // const double rhoK, + const Rd& uK, + // const double EK, + // const double gammaK, + const double cK, + // const double pK, + const Rd& normal) const + { + const Rp& EigenJ = EvaluateEigenValuesNormal( // rhoJ, + uJ, // EJ, gammaJ, + cJ, // pJ, + normal); + + const Rp& EigenK = EvaluateEigenValuesNormal( // rhoK, + uK, // EK, gammaK, + cK, // pK, + normal); + + if (EigenJ[0] * EigenK[0] < -1e-12) // <=0) //< -1e-12) + return true; + for (size_t dim = 1; dim < Dimension + 1; ++dim) // vp multiplicite d + if (EigenJ[dim] * EigenK[dim] < -1e-12) // <=0) //< -1e-12) + return true; + if (EigenJ[1 + Dimension] * EigenK[1 + Dimension] < -1e-12) // <=0) //< -1e-12) + return true; + + return false; + } + JacobianInformations - JacobianFluxAlongUnitNormal(const RoeAverageStateStructData& RoeState, const Rd& normal) const - { - /* double rho; - Rd U; - double E; - double H; - double gamma; - double p; - double c; - */ + JacobianFluxAlongUnitNormal(const RoeAverageStateStructData& RoeState, + /* + const double rhoJ, + const Rd& uJ, + const double EJ, + const double gammaJ, + const double cJ, + const double pJ, + + const double rhoK, + const Rd& uK, + const double EK, + const double gammaK, + const double cK, + const double pK, +*/ + const Rd& normal, + const bool check = false) const + { assert((l2norm(normal) - 1) < 1e-12); // const double& rho = RoeState.rho; @@ -364,24 +430,13 @@ class RoeViscousFormEulerianCompositeSolver_v2 const double gm1 = gamma - 1; const double K = c2 + gm1 * (dot(u_mean, u_mean) - H_mean); - // Rdxd UdScaln(identity); - std::vector<Rd> UdScaln(Dimension); + // Le jacobien est lineaire par rapport a la normale + // Ref : PAr ex. le papier de D. Chauvheid sur la tension de surface + Rpxp Jacobian; - // UdScaln *= uscaln; - for (size_t dim = 0; dim < Dimension; ++dim) { - Rd e(zero); - e[dim] = 1; - UdScaln[dim] = uscaln * e; - } - // Rdxd CentralT = g; Rdxd CentralT = identity; CentralT *= uscaln; - CentralT += tensorProduct(u_mean, normal) - gm1 * tensorProduct(normal, u_mean); // + UdScaln; - // std::vector<Rd> CentralTerm(Dimension); - // for (size_t dim = 0; dim < Dimension; ++dim) - // CentralTerm[dim] = Rd{u_mean[dim] * normal - gm1 * normal[dim] * u_mean + UdScaln[dim]}; - - Rpxp Jacobian(zero); + CentralT += tensorProduct(u_mean, normal) - gm1 * tensorProduct(normal, u_mean); Jacobian(0, 0) = 0.; for (size_t dim = 0; dim < Dimension; ++dim) { @@ -390,16 +445,10 @@ class RoeViscousFormEulerianCompositeSolver_v2 Jacobian(0, Dimension + 1) = 0; for (size_t dim = 0; dim < Dimension; ++dim) { - // const double& cp1 = K * normal[dim] - uscaln * u_mean[dim]; - // const Rd& cpd = CentralTerm[dim]; - // const double& cp2 = gm1 * normal[dim]; - // Rp c; Jacobian(dim + 1, 0) = K * normal[dim] - uscaln * u_mean[dim]; for (size_t dim2 = 0; dim2 < Dimension; ++dim2) Jacobian(dim + 1, dim2 + 1) = CentralT(dim, dim2); Jacobian(dim + 1, Dimension + 1) = gm1 * normal[dim]; - // Jacobian(dim + 1, dim) = c; - // Rp{K * normal[dim] - uscaln * u_mean[dim], Rd(CentralTerm[dim]), gm1 * normal[dim]}; } Jacobian(Dimension + 1, 0) = (K - H_mean) * uscaln; @@ -407,80 +456,16 @@ class RoeViscousFormEulerianCompositeSolver_v2 Jacobian(Dimension + 1, dim + 1) = (H_mean * normal[dim] - gm1 * uscaln * u_mean[dim]); Jacobian(Dimension + 1, Dimension + 1) = gamma * uscaln; - // Rp{(K - H_mean) * uscaln, (H_mean * normal - gm1 * uscaln * u_mean), gamma * uscaln}; - - // l[Dimension + 1] = Rp((K - H_mean) * uscaln, (H_mean * normal - gm1 * uscaln * u_mean), gamma * uscaln); - // std::vector<Rp> ll(Dimension + 2); - // TinyVector<Rp> ll(Dimension); - // Rpxp Jacobian(Rp(0, normal, 0), l, - // Rp((K - H_mean) * uscaln, (H_mean * normal - gm1 * uscaln * u_mean), gamma * uscaln)); - // for (size_t i = 0; i < Dimension + 2;) - // JacobianL[0] = Rp{0, normal, 0}; - //, for (size_t dim = 0; dim < Dimension; ++dim) JacobianL[1 + dim] = - // Rp{K * normal[dim] - uscaln * u_mean[dim], CentralF[dim], gm1 * normal[dim]}; - // JacobianL[1 + Dimension] = Rp{(K - H_mean) * uscaln, (H_mean * normal - gm1 * uscaln * u_mean), gamma * - // uscaln}; - - // Le jacobien est lineaire par rapport a la normale - /* - // R2 - Jac(0, 0) = 0; - Jac(0, 1) = nx; - Jac(0, 2) = ny; - Jac(0, 3) = 0; - Jac(1, 0) = (gamma - 1) * NrjCin * nx - uscaln * u; - Jac(1, 1) = uscaln + (2 - gamma) * u * nx; - Jac(1, 2) = u * ny + (1 - gamma) * v * nx; - Jac(1, 3) = (gamma - 1) * nx; - Jac(2, 0) = (gamma - 1) * NrjCin * ny - uscaln * v; - Jac(2, 1) = v * nx + (1 - gamma) * u * ny; - Jac(2, 2) = uscaln + (2 - gamma) * v * ny; - Jac(2, 3) = (gamma - 1) * ny; - Jac(3, 0) = uscaln * ((gamma - 1) * NrjCin - H_meandof); - Jac(3, 1) = H_meandof * nx + (1 - gamma) * u * uscaln; - Jac(3, 2) = H_meandof * ny + (1 - gamma) * v * uscaln; - Jac(3, 3) = gamma * uscaln; - - // R3 - Jac(0, 0) = 0; - Jac(0, 1) = nx; - Jac(0, 2) = ny; - Jac(0, 3) = nz - Jac(0, 4) = 0; - - Jac(1, 0) = (gamma - 1) * NrjCin * nx - uscaln * u; - Jac(1, 1) = uscaln + (2 - gamma) * u * nx; - - Jac(1, 4) = (gamma - 1) * nx; - - Jac(2, 0) = (gamma - 1) * NrjCin * ny - uscaln * v; - - Jac(2, 2) = uscaln + (2 - gamma) * v * ny; - - Jac(2, 4) = (gamma - 1) * ny; - - Jac(3, 0) = (gamma - 1) * NrjCin * nz - uscaln * w; - - Jac(3, 3) = uscaln + (2 - gamma) * w * nz; - Jac(3, 4) = (gamma - 1) * nz; - - Jac(4, 0) = uscaln * ((gamma - 1) * NrjCin - H_meandof); - Jac(4, 1) = H_meandof * nx + (1 - gamma) * u * uscaln; - Jac(4, 2) = H_meandof * ny + (1 - gamma) * v * uscaln; - Jac(4, 3) = H_meandof * nz + (1 - gamma) * w * uscaln; - Jac(4, 4) = gamma * uscaln; - - */ - + // Valeur propres.. Rp EigenValues; - // Rpxp Left; - // Rpxp Right; EigenValues[0] = uscaln - cspeed; for (size_t dim = 0; dim < Dimension; ++dim) // vp multiplicite d EigenValues[1 + dim] = uscaln; EigenValues[1 + Dimension] = uscaln + cspeed; + const double maxabsVpNormal = std::max(fabs(EigenValues[0]), fabs(EigenValues[1 + Dimension])); + // Vecteur propres a droite et gauche // hyper plan ortho à la normale @@ -491,6 +476,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 const double a = normal[0]; const double b = normal[1]; const double c = normal[2]; + if ((a == b) and (b == c)) { static constexpr double invsqrt2 = 1. / sqrt(2.); static constexpr double invsqrt6 = 1. / sqrt(6.); @@ -506,24 +492,13 @@ class RoeViscousFormEulerianCompositeSolver_v2 } } - // - // Rpxp RightT; - // std::vector<Rp> RightTligne(Dimension); - // std::vector<Rp> l(Dimension); - - // RightTligne[0] = Rp{1, u_mean - cspeed * normal, H_mean - uscaln * cspeed}; - // RightTligne[0] = Rp(1, u_mean, H_mean - c2 / gm1); - - // Rp(1, u_mean, H_mean - c2 / gm1); - // Rp s = Rp{1, u_mean, H_mean - c2 / gm1}; - Rpxp RightTligne; - RightTligne(0, 0) = 1; // Rp{0, ortho[0], dot(u_mean, ortho[0])}; + RightTligne(0, 0) = 1; for (size_t dim = 1; dim < Dimension + 1; ++dim) - RightTligne(0, dim) = u_mean[dim - 1] - cspeed * normal[dim - 1]; // ortho[0][dim - 1]; - RightTligne(0, Dimension + 1) = H_mean - uscaln * cspeed; // dot(u_mean, ortho[0]); + RightTligne(0, dim) = u_mean[dim - 1] - cspeed * normal[dim - 1]; + RightTligne(0, Dimension + 1) = H_mean - uscaln * cspeed; - RightTligne(1, 0) = 1; // Rp{0, ortho[0], dot(u_mean, ortho[0])}; + RightTligne(1, 0) = 1; for (size_t dim = 1; dim < Dimension + 1; ++dim) RightTligne(1, dim) = u_mean[dim - 1]; RightTligne(1, Dimension + 1) = H_mean - c2 / gm1; @@ -532,48 +507,27 @@ class RoeViscousFormEulerianCompositeSolver_v2 RightTligne(dim + 1, 0) = 0.; for (size_t dim2 = 1; dim2 < Dimension + 1; ++dim2) { RightTligne(dim + 1, dim2) = ortho[dim - 1][dim2 - 1]; - // Rp{0, ortho[dim - 1], dot(u_mean, ortho[dim - 1])}; } RightTligne(dim + 1, Dimension + 1) = dot(u_mean, ortho[dim - 1]); } - // for (size_t dim = 1; dim < Dimension; ++dim) - // RightTligne[dim] = Rp{0, ortho[dim - 1], dot(u_mean, ortho[dim - 1])}; - RightTligne(Dimension + 1, 0) = 1; // = Rp{1, u_mean + cspeed * normal, H_mean + uscaln * cspeed}; + RightTligne(Dimension + 1, 0) = 1; for (size_t dim = 1; dim < Dimension + 1; ++dim) { RightTligne(Dimension + 1, dim) = u_mean[dim - 1] + cspeed * normal[dim - 1]; } RightTligne(Dimension + 1, Dimension + 1) = H_mean + uscaln * cspeed; - // Rpxp RightT(Rp(1, u_mean - cspeed * normal, H_mean - uscaln * cspeed), RightTligne, - // Rp{1, u_mean + cspeed * normal, H_mean + uscaln * cspeed}); - - Rpxp Right = transpose(RightTligne); //.transpose(); + Rpxp Right = transpose(RightTligne); const double invc2 = 1. / c2; - Rpxp Left; //(zero); // std::vector<Rp> LeftLigne(Dimension); - /* - // LeftL[0] = .5 * invc2 * Rp{K + uscaln * cspeed, (-gm1 * u_mean - cspeed * normal), gm1}; - Left(0, 0) = gm1 * invc2 * H_mean - u2; // gm1 * invc2 * Rp{H_mean - u2, u_mean, -1}; - for (size_t dim = 1; dim < Dimension + 1; ++dim) - Left(0, dim) = gm1 * invc2 * u_mean[dim - 1]; - Left(0, Dimension + 1) = -gm1 * invc2; - - for (size_t dim = 1; dim < Dimension; ++dim) - LeftLigne[dim] = Rp{-dot(u_mean, ortho[dim - 1]), ortho[dim - 1], 0}; - // Left[1 + Dimension] = .5 * invc2 * Rp{K - uscaln * cspeed, (-gm1 * u_mean + cspeed * normal), gm1}; - - */ - // Rpxp Left(Rp(.5 * invc2 * Rp{K + uscaln * cspeed, (-gm1 * u_mean - cspeed * normal), gm1}), LeftLigne, - // Rp(.5 * invc2 * Rp{K - uscaln * cspeed, (-gm1 * u_mean + cspeed * normal), gm1})); - + Rpxp Left; //(zero); Left(0, 0) = .5 * invc2 * (K + uscaln * cspeed); for (size_t dim = 1; dim < Dimension + 1; ++dim) Left(0, dim) = .5 * invc2 * (-gm1 * u_mean[dim - 1] - cspeed * normal[dim - 1]); Left(0, Dimension + 1) = .5 * invc2 * gm1; - Left(1, 0) = gm1 * invc2 * (H_mean - u2); // gm1 * invc2 * Rp{H_mean - u2, u_mean, -1}; + Left(1, 0) = gm1 * invc2 * (H_mean - u2); for (size_t dim = 1; dim < Dimension + 1; ++dim) Left(1, dim) = gm1 * invc2 * u_mean[dim - 1]; Left(1, 1 + Dimension) = -gm1 * invc2; @@ -582,7 +536,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 Left(1 + dim, 0) = -dot(u_mean, ortho[dim - 1]); for (size_t dim2 = 0; dim2 < Dimension; ++dim2) Left(1 + dim, dim2 + 1) = ortho[dim - 1][dim2]; - Left(1 + dim, 1 + Dimension) = 0; // Rp{-dot(u_mean, ortho[dim - 1]), ortho[dim - 1], 0}; + Left(1 + dim, 1 + Dimension) = 0; } Left(1 + Dimension, 0) = .5 * invc2 * (K - uscaln * cspeed); @@ -590,83 +544,39 @@ class RoeViscousFormEulerianCompositeSolver_v2 Left(1 + Dimension, dim) = .5 * invc2 * (-gm1 * u_mean[dim - 1] + cspeed * normal[dim - 1]); Left(1 + Dimension, Dimension + 1) = .5 * invc2 * gm1; - // Left(1 + Dimension, 0) = -dot(u_mean, ortho[dim - 1]); - // for (size_t dim2 = 0; dim2 < Dimension; ++dim2) - // Left(1 + Dimension, dim2 + 1) = ortho[dim - 1][dim2]; - // Left(1 + Dimension, 1 + Dimension) = 0; // Rp{-dot(u_mean, ortho[dim - 1]), ortho[dim - 1], 0}; - // check - - Rpxp prod = Right * Left; + // Rpxp prod = Right * Left; // std::cout << prod << "\n"; Rpxp EigenMatAbs(identity); - Rpxp EigenAbs(identity); for (size_t dim = 0; dim < Dimension + 2; ++dim) { EigenMatAbs(dim, dim) *= fabs(EigenValues[dim]); - EigenAbs(dim, dim) *= EigenValues[dim]; } Rpxp ProdLefAbs = Right * EigenMatAbs; Rpxp AbsJacobian = ProdLefAbs * Left; - Rpxp ProdLeft = Right * EigenAbs; - Rpxp JacobianR = ProdLeft * Left; - - // Rpxp prodI = Left * Right; - - // std::cout << prodI << "\n"; - - /* - Right(0, 0) = 1; - Right(0, 1) = 1; - Right(0, 2) = 0; - Right(0, 3) = 1; - - Right(1, 0) = u - cspeed * nx / length; - Right(1, 1) = u; - Right(1, 2) = -ny / length; - Right(1, 3) = u + cspeed * nx / length; - - Right(2, 0) = v - cspeed * ny / length; - Right(2, 1) = v; - Right(2, 2) = nx / length; - Right(2, 3) = v + cspeed * ny / length; - - Right(3, 0) = H_meandof - uscaln * cspeed / length; - Right(3, 1) = NrjCin; - Right(3, 2) = (u_meandof, R2(-ny, nx)) / length; - Right(3, 3) = H_meandof + uscaln * cspeed / length; - */ + double maxErr = 0; - /* - Left(0,0)=((gamma-1)*NrjCin+uscaln*cspeed/lengthl)/2/(cspeed*cspeed); - Left(0,1)=-((gamma-1)*u+nx*cspeed/lengthl)/2/(cspeed*cspeed); - Left(0,2)=-((gamma-1)*v+ny*cspeed/lengthl)/2/(cspeed*cspeed); Left(0,3)=(gamma-1)/2/(cspeed*cspeed); - - Left(1,0)=1-(gamma-1)*NrjCin/(cspeed*cspeed); - Left(1,1)=(gamma-1)*u/(cspeed*cspeed);Left(1,2)=(gamma-1)*v/(cspeed*cspeed);Left(1,3)=(1-gamma)/(cspeed*cspeed); - - //Left(2,0)=-(u_meandof,R2(-ny,nx))/lengthl/cspeed; Left(2,1)=-ny/lengthl/cspeed; - Left(2,2)=nx/lengthl/cspeed; Left(2,3)=0; Left(2,0)=-(u_meandof,R2(-ny,nx))/lengthl; Left(2,1)=-ny/lengthl; - Left(2,2)=nx/lengthl; Left(2,3)=0; - //Left(2,0)=-(u_meandof,R2(-ny,nx)); Left(2,1)=-ny; Left(2,2)=nx; Left(2,3)=0; - Left(3,0)=((gamma-1)*NrjCin-uscaln*cspeed/lengthl)/2/(cspeed*cspeed); - Left(3,1)=-((gamma-1)*u-nx*cspeed/lengthl)/2/(cspeed*cspeed); - Left(3,2)=-((gamma-1)*v-ny*cspeed/lengthl)/2/(cspeed*cspeed); Left(3,3)=(gamma-1)/2/(cspeed*cspeed); - - */ - - Rpxp id = identity; - double maxErr = 0; - double maxErrLeftRightId = 0; - for (size_t i = 0; i < Dimension + 2; ++i) - for (size_t j = 0; j < Dimension + 2; ++j) { - maxErr = std::max(maxErr, fabs(Jacobian(i, j) - JacobianR(i, j))); - // maxErrLeftRightId = std::max(maxErrLeftRightId, fabs(id(i, j) - prod(i, j))); + if (check) { + Rpxp EigenAbs(identity); + for (size_t dim = 0; dim < Dimension + 2; ++dim) { + EigenAbs(dim, dim) *= EigenValues[dim]; } - // std::clog << " maxErr " << maxErr << " maxErrUnit " << maxErrLeftRightId; - // throw NormalError("STOP"); + Rpxp ProdLeft = Right * EigenAbs; + Rpxp JacobianR = ProdLeft * Left; + + Rpxp id = identity; + maxErr = 0; + // double maxErrLeftRightId = 0; + for (size_t i = 0; i < Dimension + 2; ++i) + for (size_t j = 0; j < Dimension + 2; ++j) { + maxErr = std::max(maxErr, fabs(Jacobian(i, j) - JacobianR(i, j))); + // maxErrLeftRightId = std::max(maxErrLeftRightId, fabs(id(i, j) - prod(i, j))); + } + // std::clog << " maxErr " << maxErr << " maxErrUnit " << maxErrLeftRightId; + // throw NormalError("STOP"); + } - return JacobianInformations(Jacobian, Left, Right, EigenValues, AbsJacobian, maxErr); + return JacobianInformations(Jacobian, Left, Right, EigenValues, AbsJacobian, maxabsVpNormal, maxErr); } std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, @@ -949,6 +859,8 @@ class RoeViscousFormEulerianCompositeSolver_v2 EdgeValue<double> MaxErrEdge{p_mesh->connectivity()}; MaxErrEdge.fill(0); + NodeValue<int> nbChangingSignNode{p_mesh->connectivity()}; + nbChangingSignNode.fill(0); parallel_for( p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { const auto& cell_to_node = cell_to_node_matrix[j]; @@ -980,6 +892,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 ViscosityMatrixJK *= MaxmaxabsVpNormjk; */ + // La moyenne de Roe entre les etats jk RoeAverageStateStructData RoeS( RoeAverageState(rho_n[j], u_n[j], E_n[j], gamma, p_n[j], rho_n[K], u_n[K], E_n[K], gamma, p_n[K])); @@ -993,8 +906,36 @@ class RoeViscousFormEulerianCompositeSolver_v2 const JacobianInformations& JacInfoJ = JacobianFluxAlongUnitNormal(RoeS, normalJ); // Matrice de Viscosite : Valeur Abs de la jacobienne en l'etat de Roe - const Rpxp& ViscosityMatrixJK = .5 * (nrmCjr * JacInfoJ.AbsJacobian + nrmCkr * JacInfoK.AbsJacobian); - MaxErrNode[node] = std::max(MaxErrNode[node], JacInfoK.MaxErreurProd); + Rpxp ViscosityMatrixJK = .5 * (nrmCjr * JacInfoJ.AbsJacobian + nrmCkr * JacInfoK.AbsJacobian); + + // Test Rajout Viscosité type Rusanov v2 si chgt signe vp .. + bool anySignChgJ = EvaluateChangingSignVpAlongNormal( // rho_n[j], + u_n[j], // E_n[j], gamma, + c_n[j], // p_n[j], rho_n[K], + u_n[K], + // E_n[K], gamma, + c_n[K], // p_n[K], + normalJ); + bool anySignChgK = EvaluateChangingSignVpAlongNormal( // rho_n[j], + u_n[j], // E_n[j], gamma, + c_n[j], // p_n[j], rho_n[K], + u_n[K], + // E_n[K], gamma, + c_n[K], // p_n[K], + normal); + + bool anySignDif = false; + if (anySignChgJ or anySignChgK) + anySignDif = true; + + if (anySignDif) { + Rpxp AddedViscosity(identity); + AddedViscosity *= std::max(JacInfoJ.maxabsVpNormal * nrmCjr, JacInfoK.maxabsVpNormal * nrmCkr); + ViscosityMatrixJK += AddedViscosity; + nbChangingSignNode[node] = 1; + } + + MaxErrNode[node] = std::max(MaxErrNode[node], JacInfoK.MaxErrorProd); const Rd& u_Cjr = Flux_qtmvtAtCellNode[K][R] * Cjr_loc; // Flux_qtmvt[K] * Cjr_loc; @@ -1036,9 +977,11 @@ class RoeViscousFormEulerianCompositeSolver_v2 } }); std::cout << " Max Error Node " << max(MaxErrorConservationNode) << " Max Error RoeMatrice " << max(MaxErrNode) - << "\n"; + << " nb Chg Sign " << sum(nbChangingSignNode) << "\n"; } // + FaceValue<int> nbChangingSignFace{p_mesh->connectivity()}; + nbChangingSignFace.fill(0); parallel_for( p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { // Edge @@ -1082,8 +1025,37 @@ class RoeViscousFormEulerianCompositeSolver_v2 const JacobianInformations& JacInfoJ = JacobianFluxAlongUnitNormal(RoeS, normalJ); // Matrice de Viscosite : Valeur Abs de la jacobienne en l'etat de Roe - const Rpxp& ViscosityMatrixJK = .5 * (nrmCjf * JacInfoJ.AbsJacobian + nrmCkf * JacInfoK.AbsJacobian); - MaxErrFace[face] = std::max(MaxErrFace[face], JacInfoK.MaxErreurProd); + Rpxp ViscosityMatrixJK = .5 * (nrmCjf * JacInfoJ.AbsJacobian + nrmCkf * JacInfoK.AbsJacobian); + + // Test Rajout Viscosité type Rusanov v2 si chgt signe vp .. + + bool anySignChgJ = EvaluateChangingSignVpAlongNormal( // rho_n[j], + u_n[j], // E_n[j], gamma, + c_n[j], // p_n[j], rho_n[K], + u_n[K], + // E_n[K], gamma, + c_n[K], // p_n[K], + normalJ); + bool anySignChgK = EvaluateChangingSignVpAlongNormal( // rho_n[j], + u_n[j], // E_n[j], gamma, + c_n[j], // p_n[j], rho_n[K], + u_n[K], + // E_n[K], gamma, + c_n[K], // p_n[K], + normal); + + bool anySignDif = false; + if (anySignChgJ or anySignChgK) + anySignDif = true; + + if (anySignDif) { + Rpxp AddedViscosity(identity); + AddedViscosity *= std::max(JacInfoJ.maxabsVpNormal * nrmCjf, JacInfoK.maxabsVpNormal * nrmCkf); + ViscosityMatrixJK += AddedViscosity; + nbChangingSignFace[face] = 1; + } + + MaxErrFace[face] = std::max(MaxErrFace[face], JacInfoK.MaxErrorProd); const Rd& u_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc; // Flux_qtmvt[K] * Cjf_loc; @@ -1125,7 +1097,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 } }); std::cout << " Max Error Face " << max(MaxErrorConservationFace) << " Max Error RoeMatrice " << max(MaxErrFace) - << "\n"; + << " nb Chg Sign " << sum(nbChangingSignFace) << "\n"; } if constexpr (Dimension == 3) { @@ -1135,6 +1107,9 @@ class RoeViscousFormEulerianCompositeSolver_v2 const EdgeValuePerCell<const Rd> Cje = mesh_data.Cje(); const EdgeValuePerCell<const Rd> nje = mesh_data.nje(); + EdgeValue<int> nbChangingSignEdge(p_mesh->connectivity()); + nbChangingSignEdge.fill(0); + parallel_for( p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { // Edge @@ -1180,9 +1155,37 @@ class RoeViscousFormEulerianCompositeSolver_v2 const JacobianInformations& JacInfoJ = JacobianFluxAlongUnitNormal(RoeS, normalJ); // Matrice de Viscosite : Valeur Abs de la jacobienne en l'etat de Roe - const Rpxp& ViscosityMatrixJK = .5 * (nrmCje * JacInfoJ.AbsJacobian + nrmCke * JacInfoK.AbsJacobian); + Rpxp ViscosityMatrixJK = .5 * (nrmCje * JacInfoJ.AbsJacobian + nrmCke * JacInfoK.AbsJacobian); + + // Test Rajout Viscosité type Rusanov v2 si chgt signe vp .. + + bool anySignChgJ = EvaluateChangingSignVpAlongNormal( // rho_n[j], + u_n[j], // E_n[j], gamma, + c_n[j], // p_n[j], rho_n[K], + u_n[K], + // E_n[K], gamma, + c_n[K], // p_n[K], + normalJ); + bool anySignChgK = EvaluateChangingSignVpAlongNormal( // rho_n[j], + u_n[j], // E_n[j], gamma, + c_n[j], // p_n[j], rho_n[K], + u_n[K], + // E_n[K], gamma, + c_n[K], // p_n[K], + normal); + + bool anySignDif = false; + if (anySignChgJ or anySignChgK) + anySignDif = true; + + if (anySignDif) { + Rpxp AddedViscosity(identity); + AddedViscosity *= std::max(JacInfoJ.maxabsVpNormal * nrmCje, JacInfoK.maxabsVpNormal * nrmCke); + ViscosityMatrixJK += AddedViscosity; + nbChangingSignEdge[edge] = 1; + } - MaxErrEdge[edge] = std::max(MaxErrEdge[edge], JacInfoK.MaxErreurProd); + MaxErrEdge[edge] = std::max(MaxErrEdge[edge], JacInfoK.MaxErrorProd); const Rd& u_Cje = Flux_qtmvtAtCellEdge[K][R] * Cje_loc; // Flux_qtmvt[K] * Cje_loc; @@ -1224,7 +1227,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 } }); std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << " Max Error RoeMatrice " << max(MaxErrEdge) - << "\n"; + << " nb Chg Sign " << sum(nbChangingSignEdge) << "\n"; } } // dim 3