Select Git revision
RoeViscousFormEulerianCompositeSolver_v2.cpp
RoeViscousFormEulerianCompositeSolver_v2.cpp 88.35 KiB
#include <scheme/RoeViscousFormEulerianCompositeSolver_v2.hpp>
#include <language/utils/InterpolateItemArray.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/MeshEdgeBoundary.hpp>
#include <mesh/MeshFaceBoundary.hpp>
#include <mesh/MeshFlatEdgeBoundary.hpp>
#include <mesh/MeshFlatFaceBoundary.hpp>
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/MeshTraits.hpp>
#include <mesh/MeshVariant.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
#include <scheme/InflowListBoundaryConditionDescriptor.hpp>
#include <variant>
template <MeshConcept MeshTypeT>
class RoeViscousFormEulerianCompositeSolver_v2
{
private:
using MeshType = MeshTypeT;
static constexpr size_t Dimension = MeshType::Dimension;
using Rdxd = TinyMatrix<Dimension>;
using Rd = TinyVector<Dimension>;
using Rpxp = TinyMatrix<Dimension + 2>;
using Rp = TinyVector<Dimension + 2>;
using Rpxd = TinyMatrix<Dimension + 2, Dimension>;
class SymmetryBoundaryCondition;
class InflowListBoundaryCondition;
class OutflowBoundaryCondition;
class NeumannBoundaryCondition;
class NeumannflatBoundaryCondition;
using BoundaryCondition = std::variant<SymmetryBoundaryCondition,
InflowListBoundaryCondition,
OutflowBoundaryCondition,
NeumannflatBoundaryCondition,
NeumannBoundaryCondition>;
using BoundaryConditionList = std::vector<BoundaryCondition>;
BoundaryConditionList
_getBCList(const MeshType& mesh,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
{
BoundaryConditionList bc_list;
for (const auto& bc_descriptor : bc_descriptor_list) {
bool is_valid_boundary_condition = true;
switch (bc_descriptor->type()) {
case IBoundaryConditionDescriptor::Type::neumann: {
if constexpr (Dimension == 2) {
bc_list.emplace_back(
NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
} else {
static_assert(Dimension == 3);
bc_list.emplace_back(
NeumannBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor()))); }
break;
}
case IBoundaryConditionDescriptor::Type::symmetry: {
if constexpr (Dimension == 2) {
bc_list.emplace_back(
SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
} else {
static_assert(Dimension == 3);
bc_list.emplace_back(
SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFlatEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFlatFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
}
break;
}
case IBoundaryConditionDescriptor::Type::inflow_list: {
const InflowListBoundaryConditionDescriptor& inflow_list_bc_descriptor =
dynamic_cast<const InflowListBoundaryConditionDescriptor&>(*bc_descriptor);
if (inflow_list_bc_descriptor.functionSymbolIdList().size() != 2 + Dimension) {
std::ostringstream error_msg;
error_msg << "invalid number of functions for inflow boundary "
<< inflow_list_bc_descriptor.boundaryDescriptor() << ", found "
<< inflow_list_bc_descriptor.functionSymbolIdList().size() << ", expecting " << 2 + Dimension;
throw NormalError(error_msg.str());
}
if constexpr (Dimension == 2) {
auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
Table<const double> node_values =
InterpolateItemArray<double(Rd)>::template interpolate<ItemType::node>(inflow_list_bc_descriptor
.functionSymbolIdList(),
mesh.xr(), node_boundary.nodeList());
auto xl = MeshDataManager::instance().getMeshData(mesh).xl();
auto face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
Table<const double> face_values =
InterpolateItemArray<double(Rd)>::template interpolate<ItemType::face>(inflow_list_bc_descriptor
.functionSymbolIdList(),
xl, face_boundary.faceList());
bc_list.emplace_back(InflowListBoundaryCondition(node_boundary, face_boundary, node_values, face_values));
} else {
static_assert(Dimension == 3);
auto node_boundary = getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
Table<const double> node_values =
InterpolateItemArray<double(Rd)>::template interpolate<ItemType::node>(inflow_list_bc_descriptor
.functionSymbolIdList(),
mesh.xr(), node_boundary.nodeList());
auto xe = MeshDataManager::instance().getMeshData(mesh).xe();
auto edge_boundary = getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor());
Table<const double> edge_values =
InterpolateItemArray<double(Rd)>::template interpolate<ItemType::edge>(inflow_list_bc_descriptor
.functionSymbolIdList(),
xe, edge_boundary.edgeList());
auto xl = MeshDataManager::instance().getMeshData(mesh).xl();
auto face_boundary = getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
Table<const double> face_values =
InterpolateItemArray<double(Rd)>::template interpolate<ItemType::face>(inflow_list_bc_descriptor
.functionSymbolIdList(),
xl, face_boundary.faceList());
bc_list.emplace_back(InflowListBoundaryCondition(node_boundary, edge_boundary, face_boundary, node_values,
edge_values, face_values)); }
break;
}
case IBoundaryConditionDescriptor::Type::outflow: {
if constexpr (Dimension == 2) {
bc_list.emplace_back(
OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
} else {
static_assert(Dimension == 3);
bc_list.emplace_back(
OutflowBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshEdgeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor())));
}
break;
// std::cout << "outflow not implemented yet\n";
// break;
}
default: {
is_valid_boundary_condition = false;
}
}
if (not is_valid_boundary_condition) {
std::ostringstream error_msg;
error_msg << *bc_descriptor << " is an invalid boundary condition for Rusanov v2 Eulerian Composite solver";
throw NormalError(error_msg.str());
}
}
return bc_list;
}
public:
void _applyNeumannBoundaryCondition(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
FaceValuePerCell<Rp>& stateFace) const;
void _applyNeumannflatBoundaryCondition(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
FaceValuePerCell<Rp>& stateFace) const;
void _applyInflowBoundaryCondition(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
FaceValuePerCell<Rp>& stateFace) const;
void _applyOutflowBoundaryCondition(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
FaceValuePerCell<Rp>& stateFace) const;
void _applySymmetricBoundaryCondition(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
FaceValuePerCell<Rp>& stateFace) const;
public:
double
EvaluateMaxEigenValueInGivenUnitDirection(const double& rho_mean,
const Rd& U_mean,
const double& E_mean,
const double& c_mean, const Rd& normal) const
{
const double norme_normal = l2norm(normal);
Rd unit_normal = normal;
unit_normal *= 1. / norme_normal;
const double uscaln = dot(U_mean, unit_normal);
return std::max(std::fabs(uscaln - c_mean), std::fabs(uscaln + c_mean));
}
inline double
pression(const double rho, const double epsilon, const double gam) const
{
return (gam - 1) * rho * epsilon;
}
inline Rpxd
Flux(const double& rho, const Rd& U, const double& E, const double gam) const
{
// const R2 flux_rho = rhoU;
// const R22 flux_rhoU = R22(rhoU.x1() * rhoU.x1() / rho + P, rhoU.x1() * rhoU.x2() / rho, rhoU.x2() * rhoU.x1()
// rho,rhoU.x2() * rhoU.x2() / rho + P);
// const R2 flux_rhoE = ((rhoE + P) / rho) * rhoU;
/* CellValue<Rdxd> rhoUtensU{p_mesh->connectivity()};
CellValue<Rdxd> Pid(p_mesh->connectivity());
Pid.fill(identity);
CellValue<Rdxd> rhoUtensUPlusPid{p_mesh->connectivity()};
rhoUtensUPlusPid.fill(zero);
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
rhoUtensU[j] = tensorProduct(rhoU[j], u[j]);
Pid[j] *= p_n[j];
rhoUtensUPlusPid[j] = rhoUtensU[j] + Pid[j];
});
auto Flux_rho = rhoU;
auto Flux_qtmvt = rhoUtensUPlusPid; // rhoUtensU + Pid;
auto Flux_totnrj = (rhoE + p_n) * u;
*/
const Rd& rhoU = rho * U;
const Rdxd rhoUTensU = tensorProduct(rhoU, U);
const double p = pression(rho, E - .5 * dot(U, U), gam);
Rdxd pid(identity);
pid *= p;
const Rdxd rhoUTensUPlusPid = rhoUTensU + pid;
const double rhoEPlusP = rho * E + p;
const Rd& rhoEPlusPtimesU = rhoEPlusP * U;
Rpxd Fluxx; // En ligne ci dessous
Fluxx[0] = rhoU;
for (size_t dim = 0; dim < Dimension; ++dim)
// for (size_t dim2 = 0; dim2 < Dimension; ++dim2)
Fluxx[1 + dim] = rhoUTensUPlusPid[dim];
// Fluxx[1 + dim][dim2] = rhoUTensUPlusPid[dim][dim2];
Fluxx[1 + Dimension] = rhoEPlusPtimesU;
return Fluxx;
}
struct JacobianInformations
{
Rpxp Jacobian;
Rpxp LeftEigenVectors;
Rpxp RightEigenVectors;
Rp EigenValues;
Rpxp AbsJacobian;
double MaxErreurProd; };
struct RoeAverageStateStructData
{
double rho;
Rd U;
double E;
double H;
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 )
// Rp
RoeAverageStateStructData
RoeAverageState(const double& rhoG,
const Rd& UG,
const double& EG,
const double& gammaG,
const double& pG,
const double& rhoD,
const Rd& UD,
const double& ED,
const double gammaD,
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 = 1. / (RacineRoG + RacineRoD) * (RacineRoG * UG + RacineRoD * UD);
double unorm2 = dot(U_mean, U_mean);
double NrjCin = .5 * unorm2;
const double TotEnthalpyG = (EG + pG / rhoG);
const double TotEnthalpyD = (ED + pD / rhoD);
double H_mean = (RacineRoG * TotEnthalpyG + RacineRoD * TotEnthalpyD) / (RacineRoG + RacineRoD);
double E_mean = H_mean / gamma + ((gamma - 1) / (gamma)) * (NrjCin);
double P_mean = rho_mean * (H_mean - E_mean);
double c2 = gamma * P_mean / rho_mean; // cspeed*cspeed;
// assert(fabs((gamma - 1) * rho_mean * (E_mean - .5 * (u_mean, u_mean)) - P_mean) < 1e-13); // equilibre GP
double c_mean = sqrt(c2); // cspeed_meandof/Area;
return RoeAverageStateStructData(rho_mean, U_mean, E_mean, H_mean, gamma, P_mean, c_mean);
}
JacobianInformations
JacobianFluxAlongUnitNormal(const RoeAverageStateStructData& RoeState, const Rd& normal) const
{
/* double rho;
Rd U;
double E;
double H;
double gamma;
double p;
double c; */
assert((l2norm(normal) - 1) < 1e-12);
// const double& rho = RoeState.rho;
const Rd& u_mean = RoeState.U;
const double H_mean = RoeState.H;
const double& cspeed = RoeState.c;
const double& gamma = RoeState.gamma;
const double& uscaln = dot(u_mean, normal);
const double& u2 = dot(u_mean, u_mean);
// const R NrjCin=.5*unorm2;
const double c2 = cspeed * cspeed;
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);
// 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(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;
*/
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;
// Vecteur propres a droite et gauche
// hyper plan ortho à la normale
std::vector<Rd> ortho(Dimension - 1);
if constexpr (Dimension == 2) {
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) and (b == c)) {
static constexpr double invsqrt2 = 1. / sqrt(2.);
static constexpr double invsqrt6 = 1. / sqrt(6.);
ortho[0] = Rd{invsqrt2, -invsqrt2, 0};
ortho[1] = Rd{invsqrt6, invsqrt6, -2 * invsqrt6}; // au signe pres
} else {
ortho[0] = Rd{b - c, -(a - c), a - b};
ortho[0] *= 1. / l2Norm(ortho[0]);
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;
// 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]);
}
// 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;
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}));
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;
*/
/*
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)));
}
// std::clog << " maxErr " << maxErr << " maxErrUnit " << maxErrLeftRightId;
// throw NormalError("STOP");
return JacobianInformations(Jacobian, Left, Right, EigenValues, AbsJacobian, maxErr);
}
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
solve(const std::shared_ptr<const MeshType>& p_mesh,
const DiscreteFunctionP0<const double>& rho_n,
const DiscreteFunctionP0<const Rd>& u_n,
const DiscreteFunctionP0<const double>& E_n,
const double& gamma,
const DiscreteFunctionP0<const double>& c_n,
const DiscreteFunctionP0<const double>& p_n,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt,
const bool checkLocalConservation) const
{
auto rho = copy(rho_n);
auto u = copy(u_n);
auto E = copy(E_n);
// auto c = copy(c_n);
// auto p = copy(p_n);
auto bc_list = this->_getBCList(*p_mesh, bc_descriptor_list);
auto rhoE = rho * E;
auto rhoU = rho * u;
// Construction du vecteur des variables conservatives
//
// Ci dessous juste ordre 1
// CellValue<Rp> State{p_mesh->connectivity()};
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
State[j][0] = rho[j];
for (size_t d = 0; d < Dimension; ++d)
State[j][1 + d] = rhoU[j][d];
State[j][1 + Dimension] = rhoE[j];
});
// CellValue<Rp> State{p_mesh->connectivity()};
NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()};
StateAtNode.fill(zero);
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtNode[j].fill(State[j]); });
EdgeValuePerCell<Rp> StateAtEdge{p_mesh->connectivity()};
StateAtEdge.fill(zero);
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtEdge[j].fill(State[j]); });
FaceValuePerCell<Rp> StateAtFace{p_mesh->connectivity()};
StateAtFace.fill(zero);
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtFace[j].fill(State[j]); });
// Conditions aux limites des etats conservatifs
_applyInflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
//_applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
_applySymmetricBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
_applyNeumannflatBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace);
//
// Construction du flux .. ok pour l'ordre 1
//
CellValue<Rdxd> rhoUtensU{p_mesh->connectivity()};
CellValue<Rdxd> Pid(p_mesh->connectivity());
Pid.fill(identity);
CellValue<Rdxd> rhoUtensUPlusPid{p_mesh->connectivity()};
rhoUtensUPlusPid.fill(zero);
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
rhoUtensU[j] = tensorProduct(rhoU[j], u[j]);
Pid[j] *= p_n[j];
rhoUtensUPlusPid[j] = rhoUtensU[j] + Pid[j];
});
auto Flux_rho = rhoU;
auto Flux_qtmvt = rhoUtensUPlusPid; // rhoUtensU + Pid;
auto Flux_totnrj = (rhoE + p_n) * u;
// Flux avec prise en compte des states at Node/Edge/Face
NodeValuePerCell<Rd> Flux_rhoAtCellNode{p_mesh->connectivity()};
EdgeValuePerCell<Rd> Flux_rhoAtCellEdge{p_mesh->connectivity()};
FaceValuePerCell<Rd> Flux_rhoAtCellFace{p_mesh->connectivity()};
/*
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_to_node = cell_to_node_matrix[j];
for (size_t l = 0; l < cell_to_node.size(); ++l) {
for (size_t dim = 0; dim < Dimension; ++dim)
Flux_rhoAtCellNode[j][l][dim] = StateAtNode[j][l][0] * StateAtNode[j][l][dim];
}
const auto& cell_to_face = cell_to_face_matrix[j];
for (size_t l = 0; l < cell_to_face.size(); ++l) {
for (size_t dim = 0; dim < Dimension; ++dim)
Flux_rhoAtCellFace[j][l][dim] = StateAtFace[j][l][0] * StateAtFace[j][l][dim]; }
const auto& cell_to_edge = cell_to_edge_matrix[j];
for (size_t l = 0; l < cell_to_edge.size(); ++l) {
for (size_t dim = 0; dim < Dimension; ++dim)
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;
/*
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_to_node = cell_to_node_matrix[j];
for (size_t l = 0; l < cell_to_node.size(); ++l) {
for (size_t dim = 0; dim < Dimension; ++dim)
Flux_qtmvtAtCellNode[j][l][dim] = StateAtNode[j][l][0] * StateAtNode[j][l][dim];
}
const auto& cell_to_face = cell_to_face_matrix[j];
for (size_t l = 0; l < cell_to_face.size(); ++l) {
for (size_t dim = 0; dim < Dimension; ++dim)
Flux_qtmvtAtCellFace[j][l][dim] = StateAtFace[j][l][0] * StateAtFace[j][l][dim];
}
const auto& cell_to_edge = cell_to_edge_matrix[j];
for (size_t l = 0; l < cell_to_edge.size(); ++l) {
for (size_t dim = 0; dim < Dimension; ++dim)
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];
// });
NodeValuePerCell<Rd> Flux_totnrjAtCellNode{p_mesh->connectivity()};
EdgeValuePerCell<Rd> Flux_totnrjAtCellEdge{p_mesh->connectivity()};
FaceValuePerCell<Rd> Flux_totnrjAtCellFace{p_mesh->connectivity()};
const auto& cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
const auto& cell_to_edge_matrix = p_mesh->connectivity().cellToEdgeMatrix();
const auto& cell_to_face_matrix = p_mesh->connectivity().cellToFaceMatrix();
Flux_rhoAtCellEdge.fill(zero);
Flux_totnrjAtCellEdge.fill(zero);
Flux_qtmvtAtCellEdge.fill(zero);
// Les flux aux nodes/edge/faces
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_to_node = cell_to_node_matrix[j];
for (size_t l = 0; l < cell_to_node.size(); ++l) {
// Etats conservatifs aux noeuds
const double rhonode = StateAtNode[j][l][0];
Rd Unode;
for (size_t dim = 0; dim < Dimension; ++dim)
Unode[dim] = StateAtNode[j][l][dim + 1] / rhonode;
const double rhoEnode = StateAtNode[j][l][Dimension + 1];
//
const double Enode = rhoEnode / rhonode;
const double epsilonnode = Enode - .5 * dot(Unode, Unode); const Rd rhoUnode = rhonode * Unode;
const Rdxd rhoUtensUnode = tensorProduct(rhoUnode, Unode);
const double Pressionnode = pression(rhonode, epsilonnode, gamma);
const double rhoEnodePlusP = rhoEnode + Pressionnode;
Rdxd rhoUtensUPlusPidnode(identity);
rhoUtensUPlusPidnode *= Pressionnode;
rhoUtensUPlusPidnode += rhoUtensUnode;
Flux_rhoAtCellNode[j][l] = rhoUnode;
Flux_qtmvtAtCellNode[j][l] = rhoUtensUPlusPidnode;
Flux_totnrjAtCellNode[j][l] = rhoEnodePlusP * Unode;
}
const auto& cell_to_face = cell_to_face_matrix[j];
for (size_t l = 0; l < cell_to_face.size(); ++l) {
const double rhoface = StateAtFace[j][l][0];
Rd Uface;
for (size_t dim = 0; dim < Dimension; ++dim)
Uface[dim] = StateAtFace[j][l][dim + 1] / rhoface;
const double rhoEface = StateAtFace[j][l][Dimension + 1];
//
const double Eface = rhoEface / rhoface;
const double epsilonface = Eface - .5 * dot(Uface, Uface);
const Rd rhoUface = rhoface * Uface;
const Rdxd rhoUtensUface = tensorProduct(rhoUface, Uface);
const double Pressionface = pression(rhoface, epsilonface, gamma);
const double rhoEfacePlusP = rhoEface + Pressionface;
Rdxd rhoUtensUPlusPidface(identity);
rhoUtensUPlusPidface *= Pressionface;
rhoUtensUPlusPidface += rhoUtensUface;
Flux_rhoAtCellFace[j][l] = rhoUface;
Flux_qtmvtAtCellFace[j][l] = rhoUtensUPlusPidface;
Flux_totnrjAtCellFace[j][l] = rhoEfacePlusP * Uface;
}
if constexpr (Dimension == 3) {
const auto& cell_to_edge = cell_to_edge_matrix[j];
for (size_t l = 0; l < cell_to_edge.size(); ++l) {
const double rhoedge = StateAtEdge[j][l][0];
Rd Uedge;
for (size_t dim = 0; dim < Dimension; ++dim)
Uedge[dim] = StateAtEdge[j][l][dim + 1] / rhoedge;
const double rhoEedge = StateAtEdge[j][l][Dimension + 1];
//
const double Eedge = rhoEedge / rhoedge;
const double epsilonedge = Eedge - .5 * dot(Uedge, Uedge);
const Rd rhoUedge = rhoedge * Uedge;
const Rdxd rhoUtensUedge = tensorProduct(rhoUedge, Uedge);
const double Pressionedge = pression(rhoedge, epsilonedge, gamma);
const double rhoEedgePlusP = rhoEedge + Pressionedge;
Rdxd rhoUtensUPlusPidedge(identity);
rhoUtensUPlusPidedge *= Pressionedge;
rhoUtensUPlusPidedge += rhoUtensUedge;
Flux_rhoAtCellEdge[j][l] = rhoUedge;
Flux_qtmvtAtCellEdge[j][l] = rhoUtensUPlusPidedge;
Flux_totnrjAtCellEdge[j][l] = rhoEedgePlusP * Uedge;
} }
});
// parallel_for(
// p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
// Flux_totnrjAtCellNode[j] = Flux_totnrjAtCellEdge[j] = Flux_totnrjAtCellFace[j] = Flux_totnrj[j];
// });
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(*p_mesh);
auto volumes_cell = mesh_data.Vj();
// Calcul des matrices de viscosité aux node/edge/face
const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
const NodeValuePerCell<const Rd> njr = mesh_data.njr();
const FaceValuePerCell<const Rd> Cjf = mesh_data.Cjf();
const FaceValuePerCell<const Rd> njf = mesh_data.njf();
const auto& node_to_cell_matrix = p_mesh->connectivity().nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells = p_mesh->connectivity().nodeLocalNumbersInTheirCells();
const auto& face_to_cell_matrix = p_mesh->connectivity().faceToCellMatrix();
const auto& face_local_numbers_in_their_cells = p_mesh->connectivity().faceLocalNumbersInTheirCells();
// We compute Gjr, Gjf
NodeValuePerCell<Rp> Gjr{p_mesh->connectivity()};
Gjr.fill(zero);
EdgeValuePerCell<Rp> Gje{p_mesh->connectivity()};
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) {
const auto& cell_to_node = cell_to_node_matrix[j];
for (size_t l = 0; l < cell_to_node.size(); ++l) {
const NodeId& node = cell_to_node[l];
const auto& node_to_cell = node_to_cell_matrix[node];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node);
const Rd& Cjr_loc = Cjr(j, l);
for (size_t k = 0; k < node_to_cell.size(); ++k) {
const CellId K = node_to_cell[k];
const size_t R = node_local_number_in_its_cells[k];
const Rd& Ckr_loc = Cjr(K, R);
/*
// Moyenne de 2 etats
Rd uNode = .5 * (u_n[j] + u_n[K]);
double cNode = .5 * (c_n[j] + c_n[K]);
// Viscosity j k
Rpxp ViscosityMatrixJK(identity);
const double MaxmaxabsVpNormjk =
std::max(toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uNode, cNode,
Cjr_loc),
toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uNode, cNode,
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];
const Rp& diff = ViscosityMatrixJK * statediff;
Gjr[j][l][0] += dot(Flux_rhoAtCellNode[K][R], Cjr_loc); // dot(Flux_rho[K], Cjr_loc);
for (size_t d = 0; d < Dimension; ++d)
Gjr[j][l][1 + d] += u_Cjr[d];
Gjr[j][l][1 + Dimension] += dot(Flux_totnrjAtCellNode[K][R], Cjr_loc); // dot(Flux_totnrj[K], Cjr_loc);
Gjr[j][l] += diff;
}
Gjr[j][l] *= 1. / node_to_cell.size();
}
});
if (checkLocalConservation) {
auto is_boundary_node = p_mesh->connectivity().isBoundaryNode();
NodeValue<double> MaxErrorConservationNode(p_mesh->connectivity());
MaxErrorConservationNode.fill(0.);
// double MaxErrorConservationNode = 0;
parallel_for(
p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId l) {
const auto& node_to_cell = node_to_cell_matrix[l];
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(l);
if (not is_boundary_node[l]) {
Rp SumGjr(zero);
for (size_t k = 0; k < node_to_cell.size(); ++k) {
const CellId K = node_to_cell[k];
const unsigned int R = node_local_number_in_its_cells[k];
SumGjr += Gjr[K][R];
}
// MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr));
MaxErrorConservationNode[l] = l2Norm(SumGjr);
}
});
std::cout << " Max Error Node " << max(MaxErrorConservationNode) << " Max Error RoeMatrice " << max(MaxErrNode)
<< "\n";
}
//
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
// Edge
const auto& cell_to_face = cell_to_face_matrix[j];
for (size_t l = 0; l < cell_to_face.size(); ++l) {
const FaceId& face = cell_to_face[l];
const auto& face_to_cell = face_to_cell_matrix[face];
const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(face);
const Rd& Cjf_loc = Cjf(j, l);
for (size_t k = 0; k < face_to_cell.size(); ++k) {
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);
/*
// Moyenne de 2 etats
Rd uFace = .5 * (u_n[j] + u_n[K]);
double cFace = .5 * (c_n[j] + c_n[K]);
// Viscosity j k
Rpxp ViscosityMatrixJK(identity);
const double MaxmaxabsVpNormjk =
std::max(toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uFace, cFace,
Cjf_loc),
toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uFace, cFace,
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;
const Rp& statediff = StateAtFace[j][l] - StateAtFace[K][R]; // State[j] - State[K];
const Rp& diff = ViscosityMatrixJK * statediff;
Gjf[j][l][0] += dot(Flux_rhoAtCellFace[K][R], Cjf_loc); // dot(Flux_rho[K], Cjf_loc);
for (size_t d = 0; d < Dimension; ++d)
Gjf[j][l][1 + d] += u_Cjf[d];
Gjf[j][l][1 + Dimension] += dot(Flux_totnrjAtCellFace[K][R], Cjf_loc); // dot(Flux_totnrj[K], Cjf_loc);
Gjf[j][l] += diff;
}
Gjf[j][l] *= 1. / face_to_cell.size();
}
});
if (checkLocalConservation) {
auto is_boundary_face = p_mesh->connectivity().isBoundaryFace();
FaceValue<double> MaxErrorConservationFace(p_mesh->connectivity());
MaxErrorConservationFace.fill(0.);
parallel_for(
p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) {
const auto& face_to_cell = face_to_cell_matrix[l];
const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(l);
if (not is_boundary_face[l]) {
Rp SumGjf(zero);
for (size_t k = 0; k < face_to_cell.size(); ++k) {
const CellId K = face_to_cell[k];
const unsigned int R = face_local_number_in_its_cells[k]; SumGjf += Gjf[K][R];
}
MaxErrorConservationFace[l] = l2Norm(SumGjf);
// MaxErrorConservationFace = std::max(MaxErrorConservationFace, l2Norm(SumGjf));
}
});
std::cout << " Max Error Face " << max(MaxErrorConservationFace) << " Max Error RoeMatrice " << max(MaxErrFace)
<< "\n";
}
if constexpr (Dimension == 3) {
const auto& edge_to_cell_matrix = p_mesh->connectivity().edgeToCellMatrix();
const auto& edge_local_numbers_in_their_cells = p_mesh->connectivity().edgeLocalNumbersInTheirCells();
const EdgeValuePerCell<const Rd> Cje = mesh_data.Cje();
const EdgeValuePerCell<const Rd> nje = mesh_data.nje();
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
// Edge
const auto& cell_to_edge = cell_to_edge_matrix[j];
for (size_t l = 0; l < cell_to_edge.size(); ++l) {
const EdgeId& edge = cell_to_edge[l];
const auto& edge_to_cell = edge_to_cell_matrix[edge];
const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge);
const Rd& Cje_loc = Cje(j, l);
for (size_t k = 0; k < edge_to_cell.size(); ++k) {
const CellId K = edge_to_cell[k];
const size_t R = edge_local_number_in_its_cells[k];
const Rd& Cke_loc = Cje(K, R);
/*
// Moyenne de 2 etats
Rd uEdge = .5 * (u_n[j] + u_n[K]);
double cEdge = .5 * (c_n[j] + c_n[K]);
// Viscosity j k
Rpxp ViscosityMatrixJK(identity);
const double MaxmaxabsVpNormjk =
std::max(toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uEdge, cEdge,
Cje_loc),
toolsCompositeSolver::EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(uEdge, cEdge,
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;
const Rp& statediff = StateAtEdge[j][l] - StateAtEdge[K][R]; // State[j] - State[K];
const Rp& diff = ViscosityMatrixJK * statediff;
Gje[j][l][0] += dot(Flux_rhoAtCellEdge[K][R], Cje_loc); // dot(Flux_rho[K], Cje_loc);
for (size_t d = 0; d < Dimension; ++d)
Gje[j][l][1 + d] += u_Cje[d];
Gje[j][l][1 + Dimension] += dot(Flux_totnrjAtCellEdge[K][R], Cje_loc); // dot(Flux_totnrj[K], Cje_loc);
Gje[j][l] += diff;
}
Gje[j][l] *= 1. / edge_to_cell.size();
}
});
if (checkLocalConservation) {
auto is_boundary_edge = p_mesh->connectivity().isBoundaryEdge();
EdgeValue<double> MaxErrorConservationEdge(p_mesh->connectivity());
MaxErrorConservationEdge.fill(0.);
// double MaxErrorConservationEdge = 0;
parallel_for(
p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId l) {
const auto& edge_to_cell = edge_to_cell_matrix[l];
const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(l);
if (not is_boundary_edge[l]) {
Rp SumGje(zero);
for (size_t k = 0; k < edge_to_cell.size(); ++k) {
const CellId K = edge_to_cell[k];
const unsigned int R = edge_local_number_in_its_cells[k];
SumGje += Gje[K][R];
}
// MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2norm(SumGje));
MaxErrorConservationEdge[l] = l2Norm(SumGje);
}
});
std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << " Max Error RoeMatrice " << max(MaxErrEdge)
<< "\n";
}
} // dim 3
// Pour les assemblages
double theta = .5;
double eta = .2;
if constexpr (Dimension == 2) {
eta = 0;
} else {
theta = 1. / 3.;
eta = 1. / 3.;
// theta = .5;
// eta = 0;
}
//
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
const auto& cell_to_node = cell_to_node_matrix[j];
Rp SumFluxesNode(zero);
for (size_t l = 0; l < cell_to_node.size(); ++l) {
SumFluxesNode += Gjr[j][l];
}
// SumFluxesNode *= (1 - theta);
const auto& cell_to_edge = cell_to_edge_matrix[j];
Rp SumFluxesEdge(zero);
for (size_t l = 0; l < cell_to_edge.size(); ++l) {
SumFluxesEdge += Gje[j][l];
}
const auto& cell_to_face = cell_to_face_matrix[j];
Rp SumFluxesFace(zero);
for (size_t l = 0; l < cell_to_face.size(); ++l) {
SumFluxesFace += Gjf[j][l];
}
// SumFluxesEdge *= (theta);
const Rp SumFluxes = (1 - theta - eta) * SumFluxesNode + eta * SumFluxesEdge + theta * SumFluxesFace;
State[j] -= dt / volumes_cell[j] * SumFluxes;
rho[j] = State[j][0];
for (size_t d = 0; d < Dimension; ++d)
u[j][d] = State[j][1 + d] / rho[j];
E[j] = State[j][1 + Dimension] / rho[j];
});
return std::make_tuple(std::make_shared<const DiscreteFunctionVariant>(rho),
std::make_shared<const DiscreteFunctionVariant>(u),
std::make_shared<const DiscreteFunctionVariant>(E));
}
RoeViscousFormEulerianCompositeSolver_v2() = default;
~RoeViscousFormEulerianCompositeSolver_v2() = default;
};
template <MeshConcept MeshType>
class RoeViscousFormEulerianCompositeSolver_v2<MeshType>::NeumannBoundaryCondition
{
};
template <>
class RoeViscousFormEulerianCompositeSolver_v2<Mesh<2>>::NeumannBoundaryCondition
{
using Rd = TinyVector<Dimension, double>;
private:
const MeshNodeBoundary m_mesh_node_boundary;
const MeshFaceBoundary m_mesh_face_boundary;
// const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
// const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;
public:
size_t
numberOfNodes() const
{
return m_mesh_node_boundary.nodeList().size();
}
size_t
numberOfFaces() const
{
return m_mesh_face_boundary.faceList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
const Array<const FaceId>&
faceList() const
{
return m_mesh_face_boundary.faceList();
}
NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary) : m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary)
{
;
}
};
template <>
class RoeViscousFormEulerianCompositeSolver_v2<Mesh<3>>::NeumannBoundaryCondition
{
using Rd = TinyVector<Dimension, double>;
private:
const MeshNodeBoundary m_mesh_node_boundary;
const MeshEdgeBoundary m_mesh_edge_boundary;
const MeshFaceBoundary m_mesh_face_boundary;
public:
size_t
numberOfNodes() const
{
return m_mesh_node_boundary.nodeList().size();
}
size_t
numberOfEdges() const
{
return m_mesh_edge_boundary.edgeList().size();
}
size_t
numberOfFaces() const
{
return m_mesh_face_boundary.faceList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
const Array<const EdgeId>&
edgeList() const
{
return m_mesh_edge_boundary.edgeList();
}
const Array<const FaceId>&
faceList() const
{
return m_mesh_face_boundary.faceList();
}
NeumannBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
const MeshEdgeBoundary& mesh_edge_boundary,
const MeshFaceBoundary& mesh_face_boundary)
: m_mesh_node_boundary(mesh_node_boundary),
m_mesh_edge_boundary(mesh_edge_boundary),
m_mesh_face_boundary(mesh_face_boundary)
{
;
}
};
template <MeshConcept MeshType>
class RoeViscousFormEulerianCompositeSolver_v2<MeshType>::NeumannflatBoundaryCondition
{
};
template <>class RoeViscousFormEulerianCompositeSolver_v2<Mesh<2>>::NeumannflatBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
private:
const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;
public:
const Rd&
outgoingNormal() const
{
return m_mesh_flat_node_boundary.outgoingNormal();
}
size_t
numberOfNodes() const
{
return m_mesh_flat_node_boundary.nodeList().size();
}
size_t
numberOfFaces() const
{
return m_mesh_flat_face_boundary.faceList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_flat_node_boundary.nodeList();
}
const Array<const FaceId>&
faceList() const
{
return m_mesh_flat_face_boundary.faceList();
}
NeumannflatBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
: m_mesh_flat_node_boundary(mesh_flat_node_boundary), m_mesh_flat_face_boundary(mesh_flat_face_boundary)
{
;
}
~NeumannflatBoundaryCondition() = default;
};
template <>
class RoeViscousFormEulerianCompositeSolver_v2<Mesh<3>>::NeumannflatBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
private:
const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
const MeshFlatEdgeBoundary<MeshType> m_mesh_flat_edge_boundary;
const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;
public:
const Rd&
outgoingNormal() const
{
return m_mesh_flat_node_boundary.outgoingNormal();
}
size_t
numberOfNodes() const {
return m_mesh_flat_node_boundary.nodeList().size();
}
size_t
numberOfEdges() const
{
return m_mesh_flat_edge_boundary.edgeList().size();
}
size_t
numberOfFaces() const
{
return m_mesh_flat_face_boundary.faceList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_flat_node_boundary.nodeList();
}
const Array<const EdgeId>&
edgeList() const
{
return m_mesh_flat_edge_boundary.edgeList();
}
const Array<const FaceId>&
faceList() const
{
return m_mesh_flat_face_boundary.faceList();
}
NeumannflatBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
const MeshFlatEdgeBoundary<MeshType>& mesh_flat_edge_boundary,
const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
: m_mesh_flat_node_boundary(mesh_flat_node_boundary),
m_mesh_flat_edge_boundary(mesh_flat_edge_boundary),
m_mesh_flat_face_boundary(mesh_flat_face_boundary)
{
;
}
~NeumannflatBoundaryCondition() = default;
};
template <MeshConcept MeshType>
class RoeViscousFormEulerianCompositeSolver_v2<MeshType>::SymmetryBoundaryCondition
{
};
template <>
class RoeViscousFormEulerianCompositeSolver_v2<Mesh<2>>::SymmetryBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
private:
const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;
public:
const Rd&
outgoingNormal() const
{
return m_mesh_flat_node_boundary.outgoingNormal();
}
size_t numberOfNodes() const
{
return m_mesh_flat_node_boundary.nodeList().size();
}
size_t
numberOfFaces() const
{
return m_mesh_flat_face_boundary.faceList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_flat_node_boundary.nodeList();
}
const Array<const FaceId>&
faceList() const
{
return m_mesh_flat_face_boundary.faceList();
}
SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
: m_mesh_flat_node_boundary(mesh_flat_node_boundary), m_mesh_flat_face_boundary(mesh_flat_face_boundary)
{
;
}
~SymmetryBoundaryCondition() = default;
};
template <>
class RoeViscousFormEulerianCompositeSolver_v2<Mesh<3>>::SymmetryBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
private:
const MeshFlatNodeBoundary<MeshType> m_mesh_flat_node_boundary;
const MeshFlatEdgeBoundary<MeshType> m_mesh_flat_edge_boundary;
const MeshFlatFaceBoundary<MeshType> m_mesh_flat_face_boundary;
public:
const Rd&
outgoingNormal() const
{
return m_mesh_flat_node_boundary.outgoingNormal();
}
size_t
numberOfNodes() const
{
return m_mesh_flat_node_boundary.nodeList().size();
}
size_t
numberOfEdges() const
{
return m_mesh_flat_edge_boundary.edgeList().size();
}
size_t
numberOfFaces() const
{
return m_mesh_flat_face_boundary.faceList().size();
}
const Array<const NodeId>& nodeList() const
{
return m_mesh_flat_node_boundary.nodeList();
}
const Array<const EdgeId>&
edgeList() const
{
return m_mesh_flat_edge_boundary.edgeList();
}
const Array<const FaceId>&
faceList() const
{
return m_mesh_flat_face_boundary.faceList();
}
SymmetryBoundaryCondition(const MeshFlatNodeBoundary<MeshType>& mesh_flat_node_boundary,
const MeshFlatEdgeBoundary<MeshType>& mesh_flat_edge_boundary,
const MeshFlatFaceBoundary<MeshType>& mesh_flat_face_boundary)
: m_mesh_flat_node_boundary(mesh_flat_node_boundary),
m_mesh_flat_edge_boundary(mesh_flat_edge_boundary),
m_mesh_flat_face_boundary(mesh_flat_face_boundary)
{
;
}
~SymmetryBoundaryCondition() = default;
};
template <MeshConcept MeshType>
class RoeViscousFormEulerianCompositeSolver_v2<MeshType>::InflowListBoundaryCondition
{
};
template <>
class RoeViscousFormEulerianCompositeSolver_v2<Mesh<2>>::InflowListBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
private:
const MeshNodeBoundary m_mesh_node_boundary;
const MeshFaceBoundary m_mesh_face_boundary;
const Table<const double> m_node_array_list;
const Table<const double> m_face_array_list;
public:
size_t
numberOfNodes() const
{
return m_mesh_node_boundary.nodeList().size();
}
size_t
numberOfFaces() const
{
return m_mesh_face_boundary.faceList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
const Array<const FaceId>&
faceList() const
{
return m_mesh_face_boundary.faceList(); }
const Table<const double>&
nodeArrayList() const
{
return m_node_array_list;
}
const Table<const double>&
faceArrayList() const
{
return m_face_array_list;
}
InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
const MeshFaceBoundary& mesh_face_boundary,
const Table<const double>& node_array_list,
const Table<const double>& face_array_list)
: m_mesh_node_boundary(mesh_node_boundary),
m_mesh_face_boundary(mesh_face_boundary),
m_node_array_list(node_array_list),
m_face_array_list(face_array_list)
{
;
}
~InflowListBoundaryCondition() = default;
};
template <>
class RoeViscousFormEulerianCompositeSolver_v2<Mesh<3>>::InflowListBoundaryCondition
{
public:
using Rd = TinyVector<Dimension, double>;
private:
const MeshNodeBoundary m_mesh_node_boundary;
const MeshEdgeBoundary m_mesh_edge_boundary;
const MeshFaceBoundary m_mesh_face_boundary;
const Table<const double> m_node_array_list;
const Table<const double> m_edge_array_list;
const Table<const double> m_face_array_list;
public:
size_t
numberOfNodes() const
{
return m_mesh_node_boundary.nodeList().size();
}
size_t
numberOfEdges() const
{
return m_mesh_edge_boundary.edgeList().size();
}
size_t
numberOfFaces() const
{
return m_mesh_face_boundary.faceList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
const Array<const EdgeId>&
edgeList() const {
return m_mesh_edge_boundary.edgeList();
}
const Array<const FaceId>&
faceList() const
{
return m_mesh_face_boundary.faceList();
}
const Table<const double>&
nodeArrayList() const
{
return m_node_array_list;
}
const Table<const double>&
edgeArrayList() const
{
return m_edge_array_list;
}
const Table<const double>&
faceArrayList() const
{
return m_face_array_list;
}
InflowListBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary,
const MeshEdgeBoundary& mesh_edge_boundary,
const MeshFaceBoundary& mesh_face_boundary,
const Table<const double>& node_array_list,
const Table<const double>& edge_array_list,
const Table<const double>& face_array_list)
: m_mesh_node_boundary(mesh_node_boundary),
m_mesh_edge_boundary(mesh_edge_boundary),
m_mesh_face_boundary(mesh_face_boundary),
m_node_array_list(node_array_list),
m_edge_array_list(edge_array_list),
m_face_array_list(face_array_list)
{
;
}
~InflowListBoundaryCondition() = default;
};
template <MeshConcept MeshType>
class RoeViscousFormEulerianCompositeSolver_v2<MeshType>::OutflowBoundaryCondition
{
};
template <>
class RoeViscousFormEulerianCompositeSolver_v2<Mesh<2>>::OutflowBoundaryCondition
{
using Rd = TinyVector<Dimension, double>;
private:
const MeshNodeBoundary m_mesh_node_boundary;
const MeshFaceBoundary m_mesh_face_boundary;
public:
size_t
numberOfNodes() const
{
return m_mesh_node_boundary.nodeList().size();
}
size_t
numberOfFaces() const {
return m_mesh_face_boundary.faceList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
const Array<const FaceId>&
faceList() const
{
return m_mesh_face_boundary.faceList();
}
OutflowBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshFaceBoundary& mesh_face_boundary)
: m_mesh_node_boundary(mesh_node_boundary), m_mesh_face_boundary(mesh_face_boundary)
{
;
}
};
template <>
class RoeViscousFormEulerianCompositeSolver_v2<Mesh<3>>::OutflowBoundaryCondition
{
using Rd = TinyVector<Dimension, double>;
private:
const MeshNodeBoundary m_mesh_node_boundary;
const MeshEdgeBoundary m_mesh_edge_boundary;
const MeshFaceBoundary m_mesh_face_boundary;
public:
size_t
numberOfNodes() const
{
return m_mesh_node_boundary.nodeList().size();
}
size_t
numberOfEdges() const
{
return m_mesh_edge_boundary.edgeList().size();
}
size_t
numberOfFaces() const
{
return m_mesh_face_boundary.faceList().size();
}
const Array<const NodeId>&
nodeList() const
{
return m_mesh_node_boundary.nodeList();
}
const Array<const EdgeId>&
edgeList() const
{
return m_mesh_edge_boundary.edgeList();
}
const Array<const FaceId>&
faceList() const
{
return m_mesh_face_boundary.faceList();
}
OutflowBoundaryCondition(const MeshNodeBoundary& mesh_node_boundary, const MeshEdgeBoundary& mesh_edge_boundary,
const MeshFaceBoundary& mesh_face_boundary)
: m_mesh_node_boundary(mesh_node_boundary),
m_mesh_edge_boundary(mesh_edge_boundary),
m_mesh_face_boundary(mesh_face_boundary)
{
;
}
};
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
roeViscousFormEulerianCompositeSolver_v2(
const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
const std::shared_ptr<const DiscreteFunctionVariant>& E_v,
const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt,
const bool check)
{
std::shared_ptr mesh_v = getCommonMesh({rho_v, u_v, E_v, c_v, p_v});
if (not mesh_v) {
throw NormalError("discrete functions are not defined on the same mesh");
}
if (not checkDiscretizationType({rho_v, u_v, E_v}, DiscreteFunctionType::P0)) {
throw NormalError("acoustic solver expects P0 functions");
}
return std::visit(
PUGS_LAMBDA(auto&& p_mesh)
->std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
using MeshType = mesh_type_t<decltype(p_mesh)>;
static constexpr size_t Dimension = MeshType::Dimension;
using Rd = TinyVector<Dimension>;
if constexpr (Dimension == 1) {
throw NormalError("RoeViscousFormEulerianCompositeSolver v2 is not available in 1D");
} else {
if constexpr (is_polygonal_mesh_v<MeshType>) {
return RoeViscousFormEulerianCompositeSolver_v2<MeshType>{}
.solve(p_mesh, rho_v->get<DiscreteFunctionP0<const double>>(), u_v->get<DiscreteFunctionP0<const Rd>>(),
E_v->get<DiscreteFunctionP0<const double>>(), gamma, c_v->get<DiscreteFunctionP0<const double>>(),
p_v->get<DiscreteFunctionP0<const double>>(), bc_descriptor_list, dt, check);
} else {
throw NormalError("RoeViscousFormEulerianCompositeSolver v2 is only defined on polygonal meshes");
}
}
},
mesh_v->variant());
}
template <MeshConcept MeshType>
void
RoeViscousFormEulerianCompositeSolver_v2<MeshType>::_applyOutflowBoundaryCondition(
const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
FaceValuePerCell<Rp>& stateFace) const
{
for (const auto& boundary_condition : bc_list) {
std::visit( [&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<OutflowBoundaryCondition, T>) {
std::cout << " Traitement Outflow \n";
// const Rd& normal = bc.outgoingNormal();
/*
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
// const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
const auto& face_list = bc.faceList();
const auto& node_list = bc.nodeList();
const auto xj = mesh.xj();
const auto xr = mesh.xr();
const auto xf = mesh.xl();
const auto xe = mesh.xe();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_cell_list = node_to_cell_matrix[node_id];
// Assert(face_cell_list.size() == 1);
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
CellId node_cell_id = node_cell_list[i_cell];
size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];
for (size_t dim = 0; dim < Dimension + 2; ++dim)
stateNode[node_cell_id][node_local_number_in_cell][dim] += vectorSym[dim];
Rd vectorSym(zero);
for (size_t dim = 0; dim < Dimension; ++dim)
vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
vectorSym = MatriceProj * vectorSym;
for (size_t dim = 0; dim < Dimension; ++dim)
stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim];
// stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim];
}
}
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_cell_list = face_to_cell_matrix[face_id];
Assert(face_cell_list.size() == 1);
CellId face_cell_id = face_cell_list[0];
size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
Rd vectorSym(zero);
for (size_t dim = 0; dim < Dimension; ++dim)
vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
vectorSym = MatriceProj * vectorSym;
for (size_t dim = 0; dim < Dimension; ++dim)
stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim];
}
if constexpr (Dimension == 3) {
const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();
const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();
const auto& edge_list = bc.edgeList();
for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
const EdgeId edge_id = edge_list[i_edge];
const auto& edge_cell_list = edge_to_cell_matrix[edge_id];
// Assert(face_cell_list.size() == 1);
const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
CellId edge_cell_id = edge_cell_list[i_cell];
size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];
Rd vectorSym(zero);
for (size_t dim = 0; dim < Dimension; ++dim)
vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
vectorSym = MatriceProj * vectorSym;
for (size_t dim = 0; dim < Dimension; ++dim)
stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
}
}
// throw NormalError("Not implemented");
}
*/
}
},
boundary_condition);
}
}
template <MeshConcept MeshType>
void
RoeViscousFormEulerianCompositeSolver_v2<MeshType>::_applySymmetricBoundaryCondition(
const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
FaceValuePerCell<Rp>& stateFace) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
// MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
std::cout << " Traitement SYMMETRY \n";
const Rd& normal = bc.outgoingNormal();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
// const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
const auto& face_list = bc.faceList();
const auto& node_list = bc.nodeList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_cell_list = node_to_cell_matrix[node_id];
// Assert(face_cell_list.size() == 1);
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
CellId node_cell_id = node_cell_list[i_cell];
size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];
Rd vectorSym(zero);
for (size_t dim = 0; dim < Dimension; ++dim)
vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
vectorSym = MatriceProj * vectorSym;
for (size_t dim = 0; dim < Dimension; ++dim)
stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim];
// stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim];
}
}
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_cell_list = face_to_cell_matrix[face_id];
Assert(face_cell_list.size() == 1);
CellId face_cell_id = face_cell_list[0];
size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
Rd vectorSym(zero);
for (size_t dim = 0; dim < Dimension; ++dim)
vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
vectorSym = MatriceProj * vectorSym;
for (size_t dim = 0; dim < Dimension; ++dim)
stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim];
}
if constexpr (Dimension == 3) {
const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();
const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();
const auto& edge_list = bc.edgeList();
for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
const EdgeId edge_id = edge_list[i_edge];
const auto& edge_cell_list = edge_to_cell_matrix[edge_id];
// Assert(face_cell_list.size() == 1);
const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
CellId edge_cell_id = edge_cell_list[i_cell];
size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];
Rd vectorSym(zero);
for (size_t dim = 0; dim < Dimension; ++dim)
vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim];
Rdxd MatriceProj(identity);
MatriceProj -= tensorProduct(normal, normal);
vectorSym = MatriceProj * vectorSym;
for (size_t dim = 0; dim < Dimension; ++dim)
stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
}
}
}
}
},
boundary_condition);
}
}
template <MeshConcept MeshType>
void
RoeViscousFormEulerianCompositeSolver_v2<MeshType>::_applyNeumannflatBoundaryCondition(
const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
FaceValuePerCell<Rp>& stateFace) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<NeumannflatBoundaryCondition, T>) {
// MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
std::cout << " Traitement WALL \n";
const Rd& normal = bc.outgoingNormal();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
// const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
const auto& face_list = bc.faceList();
const auto& node_list = bc.nodeList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_cell_list = node_to_cell_matrix[node_id];
// Assert(face_cell_list.size() == 1);
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
CellId node_cell_id = node_cell_list[i_cell];
size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];
Rd vectorSym(zero);
for (size_t dim = 0; dim < Dimension; ++dim)
vectorSym[dim] = stateNode[node_cell_id][node_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
for (size_t dim = 0; dim < Dimension; ++dim)
stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = vectorSym[dim];
// stateNode[node_cell_id][node_local_number_in_cell][dim] = 0; // node_array_list[i_node][dim];
}
}
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_cell_list = face_to_cell_matrix[face_id];
Assert(face_cell_list.size() == 1);
CellId face_cell_id = face_cell_list[0];
size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
Rd vectorSym(zero);
for (size_t dim = 0; dim < Dimension; ++dim)
vectorSym[dim] = stateEdge[face_cell_id][face_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
for (size_t dim = 0; dim < Dimension; ++dim)
stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim];
}
if constexpr (Dimension == 3) {
const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();
const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();
const auto& edge_list = bc.edgeList();
for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
const EdgeId edge_id = edge_list[i_edge];
const auto& edge_cell_list = edge_to_cell_matrix[edge_id];
// Assert(face_cell_list.size() == 1);
const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
CellId edge_cell_id = edge_cell_list[i_cell];
size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];
Rd vectorSym(zero);
for (size_t dim = 0; dim < Dimension; ++dim)
vectorSym[dim] = stateEdge[edge_cell_id][edge_local_number_in_cell][1 + dim];
vectorSym -= dot(vectorSym, normal) * normal;
for (size_t dim = 0; dim < Dimension; ++dim)
stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim];
}
}
}
}
},
boundary_condition);
}
}
template <MeshConcept MeshType>
void
RoeViscousFormEulerianCompositeSolver_v2<MeshType>::_applyInflowBoundaryCondition(const BoundaryConditionList& bc_list,
const MeshType& mesh,
NodeValuePerCell<Rp>& stateNode,
EdgeValuePerCell<Rp>& stateEdge,
FaceValuePerCell<Rp>& stateFace) const
{
for (const auto& boundary_condition : bc_list) {
std::visit(
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<InflowListBoundaryCondition, T>) {
// MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
std::cout << " Traitement INFLOW \n";
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
// const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
const auto& face_list = bc.faceList(); const auto& node_list = bc.nodeList();
const auto& face_array_list = bc.faceArrayList();
const auto& node_array_list = bc.nodeArrayList();
for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
const NodeId node_id = node_list[i_node];
const auto& node_cell_list = node_to_cell_matrix[node_id];
// Assert(face_cell_list.size() == 1);
const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
CellId node_cell_id = node_cell_list[i_cell];
size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell];
for (size_t dim = 0; dim < Dimension + 2; ++dim)
stateNode[node_cell_id][node_local_number_in_cell][dim] = node_array_list[i_node][dim];
}
}
for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
const FaceId face_id = face_list[i_face];
const auto& face_cell_list = face_to_cell_matrix[face_id];
Assert(face_cell_list.size() == 1);
CellId face_cell_id = face_cell_list[0];
size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
for (size_t dim = 0; dim < Dimension + 2; ++dim)
stateFace[face_cell_id][face_local_number_in_cell][dim] = face_array_list[i_face][dim];
}
if constexpr (Dimension == 3) {
const auto& edge_to_cell_matrix = mesh.connectivity().edgeToCellMatrix();
const auto& edge_local_numbers_in_their_cells = mesh.connectivity().edgeLocalNumbersInTheirCells();
// const auto& face_cell_is_reversed = mesh.connectivity().cellFaceIsReversed();
const auto& edge_list = bc.edgeList();
const auto& edge_array_list = bc.edgeArrayList();
for (size_t i_edge = 0; i_edge < edge_list.size(); ++i_edge) {
const EdgeId edge_id = edge_list[i_edge];
const auto& edge_cell_list = edge_to_cell_matrix[edge_id];
const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(edge_id);
// Assert(face_cell_list.size() == 1);
for (size_t i_cell = 0; i_cell < edge_cell_list.size(); ++i_cell) {
CellId edge_cell_id = edge_cell_list[i_cell];
size_t edge_local_number_in_cell = edge_local_number_in_its_cells[i_cell];
for (size_t dim = 0; dim < Dimension + 2; ++dim)
stateEdge[edge_cell_id][edge_local_number_in_cell][dim] = edge_array_list[i_edge][dim];
}
}
}
}
},
boundary_condition);
}
}