diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index a37461b77b51b4165a43eb2b4e2afe71b7d0cf86..5269981b3c1b00aa90f89444f2fd1a324322711b 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -43,6 +43,7 @@
#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
#include <scheme/OutflowBoundaryConditionDescriptor.hpp>
+#include <scheme/RoeViscousFormEulerianCompositeSolver_v2.hpp>
#include <scheme/RusanovEulerianCompositeSolver.hpp>
#include <scheme/RusanovEulerianCompositeSolver_v2.hpp>
#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
@@ -572,6 +573,39 @@ SchemeModule::SchemeModule()
std::shared_ptr<const DiscreteFunctionVariant>> {
return rusanovEulerianCompositeSolver_v2(rho, u, E, gamma, c, p, bc_descriptor_list, dt,
true);
+ }));
+ this->_addBuiltinFunction("roe_viscosityform_eulerian_composite_solver_version2",
+ std::function(
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return roeViscousFormEulerianCompositeSolver_v2(rho, u, E, gamma, c, p,
+ bc_descriptor_list, dt);
+ }
+
+ ));
+
+ this->_addBuiltinFunction("roe_viscosityform_eulerian_composite_solver_version2_with_checks",
+ std::function(
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+ const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+ bc_descriptor_list,
+ const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>,
+ std::shared_ptr<const DiscreteFunctionVariant>> {
+ return roeViscousFormEulerianCompositeSolver_v2(rho, u, E, gamma, c, p,
+ bc_descriptor_list, dt, true);
}
));
diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
index 646088154db01b557a029724d73315239240f7be..7650806336ffd4394a597b46463d8341aa80bb0f 100644
--- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
@@ -276,6 +276,8 @@ class RoeViscousFormEulerianCompositeSolver_v2
Rpxp LeftEigenVectors;
Rpxp RightEigenVectors;
Rp EigenValues;
+ Rpxp AbsJacobian;
+ double MaxErreurProd;
};
struct RoeAverageStateStructData
@@ -287,6 +289,15 @@ class RoeViscousFormEulerianCompositeSolver_v2
double gamma;
double p;
double c;
+ RoeAverageStateStructData(const double rhod,
+ const Rd Ud,
+ const double Ed,
+ const double Hd,
+ const double gammad,
+ const double pd,
+ const double cd)
+ : 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 )
@@ -301,14 +312,14 @@ class RoeViscousFormEulerianCompositeSolver_v2
const Rd& UD,
const double& ED,
const double gammaD,
- const double pD)
+ const double pD) const
{
double gamma = .5 * (gammaG + gammaD); // ou ponderation racine roG et roD
double RacineRoG = sqrt(rhoG);
double RacineRoD = sqrt(rhoD);
double rho_mean = RacineRoG * RacineRoD;
- Rd U_mean = (RacineRoG * UG + RacineRoD * UD) / (RacineRoG + RacineRoD);
+ Rd U_mean = 1. / (RacineRoG + RacineRoD) * (RacineRoG * UG + RacineRoD * UD);
double unorm2 = dot(U_mean, U_mean);
double NrjCin = .5 * unorm2;
const double TotEnthalpyG = (EG + pG / rhoG);
@@ -328,7 +339,7 @@ class RoeViscousFormEulerianCompositeSolver_v2
}
JacobianInformations
- JacobianFluxAlongUnitNormal(const RoeAverageStateStructData& RoeState, const Rd& normal)
+ JacobianFluxAlongUnitNormal(const RoeAverageStateStructData& RoeState, const Rd& normal) const
{
/* double rho;
Rd U;
@@ -340,7 +351,7 @@ class RoeViscousFormEulerianCompositeSolver_v2
*/
assert((l2norm(normal) - 1) < 1e-12);
- const double& rho = RoeState.rho;
+ // const double& rho = RoeState.rho;
const Rd& u_mean = RoeState.U;
const double H_mean = RoeState.H;
const double& cspeed = RoeState.c;
@@ -353,14 +364,62 @@ class RoeViscousFormEulerianCompositeSolver_v2
const double gm1 = gamma - 1;
const double K = c2 + gm1 * (dot(u_mean, u_mean) - H_mean);
- Rpxp Jacobian;
- Rdxd UdScaln(identity);
- UdScaln *= uscaln;
- Rdxd CentralF = tensorProduct(u_mean, normal) - gm1 * tensorProduct(normal, u_mean) + UdScaln;
- Jacobian[0] = Rp{0, normal, 0};
+ // Rdxd UdScaln(identity);
+ std::vector<Rd> UdScaln(Dimension);
+
+ // 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);
+
+ Jacobian(0, 0) = 0.;
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ Jacobian(0, dim + 1) = normal[dim];
+ }
+ 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;
for (size_t dim = 0; dim < Dimension; ++dim)
- Jacobian[1 + dim] = Rp{K * normal[dim] - uscaln * u_mean[dim], CentralF[dim], gm1 * normal[dim]};
- Jacobian[1 + Dimension] = Rp{(K - H_mean) * uscaln, (H_mean * normal - gm1 * uscaln * u_mean), gamma * uscaln};
+ 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
/*
@@ -414,75 +473,169 @@ class RoeViscousFormEulerianCompositeSolver_v2
*/
Rp EigenValues;
- Rpxp Left;
- Rpxp Right;
+ // Rpxp Left;
+ // Rpxp Right;
EigenValues[0] = uscaln - cspeed;
- for (int dim = 0; dim < Dimension; ++dim) // vp multidplicite d
+ for (size_t dim = 0; dim < Dimension; ++dim) // vp multiplicite d
EigenValues[1 + dim] = uscaln;
-
EigenValues[1 + Dimension] = uscaln + cspeed;
+
// Vecteur propres a droite et gauche
+
+ // hyper plan ortho à la normale
std::vector<Rd> ortho(Dimension - 1);
if constexpr (Dimension == 2) {
- ortho[0] = {normal[1], -normal[0]}; // aussi de norme 1
+ ortho[0] = Rd{normal[1], -normal[0]}; // aussi de norme 1
} else {
const double a = normal[0];
const double b = normal[1];
const double c = normal[2];
- if (a == b == c) {
+ if ((a == b) and (b == c)) {
static constexpr double invsqrt2 = 1. / sqrt(2.);
static constexpr double invsqrt6 = 1. / sqrt(6.);
- ortho[0] = {invsqrt2, -invsqrt2, 0};
- ortho[1] = {invsqrt6, invsqrt6, -2 * invsqrt6};
+ ortho[0] = Rd{invsqrt2, -invsqrt2, 0};
+ ortho[1] = Rd{invsqrt6, invsqrt6, -2 * invsqrt6}; // au signe pres
} else {
- ortho[0] = {b - c, -(a - c), a - b};
+ ortho[0] = Rd{b - c, -(a - c), a - b};
ortho[0] *= 1. / l2Norm(ortho[0]);
- ortho[1] = {a * (b + c) - b * b - c * c, b * (a + c) - a * a - c * c, c * (a + b) - a * a - b * b};
- ortho[1] *= 1. / l2Norm(ortho[1]);
+ ortho[1] = Rd{a * (b + c) - b * b - c * c, b * (a + c) - a * a - c * c, c * (a + b) - a * a - b * b};
+ ortho[1] *= 1. / l2Norm(ortho[1]); // au signe pres
}
}
//
- Rpxp RightT;
- RightT[0] = {1, u_mean - cspeed * normal, H_mean - uscaln * cspeed};
- RightT[1] = {1, u_mean, H_mean - c2 / gm1};
- for (size_t dim = 1; dim < Dimension; ++dim)
- RightT[1 + dim] = Rp{0, ortho[dim - 1], dot(u_mean, ortho[dim - 1])};
- RightT[1 + Dimension] = {1, u_mean + cspeed * normal, H_mean + uscaln * cspeed};
+ // 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])};
+ 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(1, 0) = 1; // Rp{0, ortho[0], dot(u_mean, ortho[0])};
+ for (size_t dim = 1; dim < Dimension + 1; ++dim)
+ RightTligne(1, dim) = u_mean[dim - 1];
+ RightTligne(1, Dimension + 1) = H_mean - c2 / gm1;
+
+ for (size_t dim = 1; dim < Dimension; ++dim) {
+ 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]);
+ }
- Right = RightT.transpose();
+ // 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};
+ 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();
const double invc2 = 1. / c2;
- Left[0] = .5 * invc2 * Rp{K + uscaln * cspeed, (-gm1 * u_mean - cspeed * normal), gm1};
- Left[1] = gm1 * invc2 * Rp{H_mean - u2, u_mean, -1};
+
+ 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)
- Left[1 + 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};
+ 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}));
+
+ 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};
+ for (size_t dim = 1; dim < Dimension + 1; ++dim)
+ Left(1, dim) = gm1 * invc2 * u_mean[dim - 1];
+ Left(1, 1 + Dimension) = -gm1 * invc2;
+
+ for (size_t dim = 1; dim < Dimension; ++dim) {
+ 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 + Dimension, 0) = .5 * invc2 * (K - uscaln * cspeed);
+ for (size_t dim = 1; dim < Dimension + 1; ++dim)
+ 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;
+ // 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;
- */
+ 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;
+ */
/*
Left(0,0)=((gamma-1)*NrjCin+uscaln*cspeed/lengthl)/2/(cspeed*cspeed);
@@ -501,8 +654,19 @@ class RoeViscousFormEulerianCompositeSolver_v2
Left(3,2)=-((gamma-1)*v-ny*cspeed/lengthl)/2/(cspeed*cspeed); Left(3,3)=(gamma-1)/2/(cspeed*cspeed);
*/
- Right = zero;
- Left = zero;
+
+ 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)));
+ }
+ // std::clog << " maxErr " << maxErr << " maxErrUnit " << maxErrLeftRightId;
+ // throw NormalError("STOP");
+
+ return JacobianInformations(Jacobian, Left, Right, EigenValues, AbsJacobian, maxErr);
}
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
@@ -613,7 +777,7 @@ class RoeViscousFormEulerianCompositeSolver_v2
Flux_rhoAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim];
}
});
-*/
+ */
NodeValuePerCell<Rdxd> Flux_qtmvtAtCellNode{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
EdgeValuePerCell<Rdxd> Flux_qtmvtAtCellEdge{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
FaceValuePerCell<Rdxd> Flux_qtmvtAtCellFace{p_mesh->connectivity()}; // = rhoUtensUPlusPid; // rhoUtensU + Pid;
@@ -641,7 +805,7 @@ class RoeViscousFormEulerianCompositeSolver_v2
Flux_qtmvtAtCellEdge[j][l][dim] = StateAtEdge[j][l][0] * StateAtEdge[j][l][dim];
}
});
-*/
+ */
// parallel_for(
// p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
// Flux_qtmvtAtCellNode[j] = Flux_qtmvtAtCellEdge[j] = Flux_qtmvtAtCellFace[j] = Flux_qtmvt[j];
@@ -778,6 +942,12 @@ class RoeViscousFormEulerianCompositeSolver_v2
Gje.fill(zero);
FaceValuePerCell<Rp> Gjf{p_mesh->connectivity()};
Gjf.fill(zero);
+ NodeValue<double> MaxErrNode{p_mesh->connectivity()};
+ MaxErrNode.fill(0);
+ FaceValue<double> MaxErrFace{p_mesh->connectivity()};
+ MaxErrFace.fill(0);
+ EdgeValue<double> MaxErrEdge{p_mesh->connectivity()};
+ MaxErrEdge.fill(0);
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
@@ -795,8 +965,8 @@ class RoeViscousFormEulerianCompositeSolver_v2
const size_t R = node_local_number_in_its_cells[k];
const Rd& Ckr_loc = Cjr(K, R);
- // Une moyenne entre les etats jk
-
+ /*
+ // Moyenne de 2 etats
Rd uNode = .5 * (u_n[j] + u_n[K]);
double cNode = .5 * (c_n[j] + c_n[K]);
@@ -809,6 +979,23 @@ class RoeViscousFormEulerianCompositeSolver_v2
Ckr_loc));
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]));
+ // Viscosity j k
+ const double nrmCkr = l2Norm(Ckr_loc);
+ const Rd& normal = 1. / nrmCkr * Ckr_loc;
+ const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, normal);
+
+ const double nrmCjr = l2Norm(Cjr_loc);
+ const Rd& normalJ = 1. / nrmCjr * Cjr_loc;
+ 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);
+
const Rd& u_Cjr = Flux_qtmvtAtCellNode[K][R] * Cjr_loc; // Flux_qtmvt[K] * Cjr_loc;
const Rp& statediff = StateAtNode[j][l] - StateAtNode[K][R]; // State[j] - State[K];
@@ -848,7 +1035,8 @@ class RoeViscousFormEulerianCompositeSolver_v2
MaxErrorConservationNode[l] = l2Norm(SumGjr);
}
});
- std::cout << " Max Error Node " << max(MaxErrorConservationNode) << "\n";
+ std::cout << " Max Error Node " << max(MaxErrorConservationNode) << " Max Error RoeMatrice " << max(MaxErrNode)
+ << "\n";
}
//
parallel_for(
@@ -867,8 +1055,8 @@ class RoeViscousFormEulerianCompositeSolver_v2
const CellId K = face_to_cell[k];
const unsigned int R = face_local_number_in_its_cells[k];
const Rd& Ckf_loc = Cjf(K, R);
- // Une moyenne entre les etats jk
-
+ /*
+ // Moyenne de 2 etats
Rd uFace = .5 * (u_n[j] + u_n[K]);
double cFace = .5 * (c_n[j] + c_n[K]);
@@ -881,6 +1069,21 @@ class RoeViscousFormEulerianCompositeSolver_v2
Ckf_loc));
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]));
+ // Viscosity j k
+ const double nrmCkf = l2Norm(Ckf_loc);
+ const Rd& normal = 1. / nrmCkf * Ckf_loc;
+ const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, normal);
+ const double nrmCjf = l2Norm(Cjf_loc);
+ const Rd& normalJ = 1. / nrmCjf * Cjf_loc;
+ 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);
const Rd& u_Cjf = Flux_qtmvtAtCellFace[K][R] * Cjf_loc; // Flux_qtmvt[K] * Cjf_loc;
@@ -921,7 +1124,8 @@ class RoeViscousFormEulerianCompositeSolver_v2
// MaxErrorConservationFace = std::max(MaxErrorConservationFace, l2Norm(SumGjf));
}
});
- std::cout << " Max Error Face " << max(MaxErrorConservationFace) << "\n";
+ std::cout << " Max Error Face " << max(MaxErrorConservationFace) << " Max Error RoeMatrice " << max(MaxErrFace)
+ << "\n";
}
if constexpr (Dimension == 3) {
@@ -948,9 +1152,8 @@ class RoeViscousFormEulerianCompositeSolver_v2
const size_t R = edge_local_number_in_its_cells[k];
const Rd& Cke_loc = Cje(K, R);
-
- // Une moyenne entre les etats jk
-
+ /*
+ // Moyenne de 2 etats
Rd uEdge = .5 * (u_n[j] + u_n[K]);
double cEdge = .5 * (c_n[j] + c_n[K]);
@@ -963,6 +1166,23 @@ class RoeViscousFormEulerianCompositeSolver_v2
Cke_loc));
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]));
+ // Viscosity j k
+ const double nrmCke = l2Norm(Cke_loc);
+ const Rd& normal = 1. / nrmCke * Cke_loc;
+ const JacobianInformations& JacInfoK = JacobianFluxAlongUnitNormal(RoeS, normal);
+
+ const double nrmCje = l2Norm(Cje_loc);
+ const Rd& normalJ = 1. / nrmCje * Cje_loc;
+ 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);
+
+ MaxErrEdge[edge] = std::max(MaxErrEdge[edge], JacInfoK.MaxErreurProd);
const Rd& u_Cje = Flux_qtmvtAtCellEdge[K][R] * Cje_loc; // Flux_qtmvt[K] * Cje_loc;
@@ -1003,7 +1223,8 @@ class RoeViscousFormEulerianCompositeSolver_v2
MaxErrorConservationEdge[l] = l2Norm(SumGje);
}
});
- std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << "\n";
+ std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << " Max Error RoeMatrice " << max(MaxErrEdge)
+ << "\n";
}
} // dim 3