From e3246768de8b49d3a4bf771dd2526fa0f94044e5 Mon Sep 17 00:00:00 2001
From: HOCH PHILIPPE <philippe.hoch@gmail.com>
Date: Wed, 7 Aug 2024 23:58:52 +0200
Subject: [PATCH] Codage Rusanov 2D 3D composite version 1
---
src/language/modules/SchemeModule.cpp | 25 +
src/scheme/CMakeLists.txt | 10 +
src/scheme/RusanovEulerianCompositeSolver.cpp | 545 +++++++++++++++++-
src/scheme/RusanovEulerianCompositeSolver.hpp | 6 +-
4 files changed, 580 insertions(+), 6 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 43b6e5c3d..3960f2667 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -520,6 +520,24 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("rusanov_eulerian_composite_solver_version1_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 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 rusanovEulerianCompositeSolver(rho, u, E, c, p, bc_descriptor_list, dt, true);
+ }
+
+ ));
+
this->_addBuiltinFunction("eucclhyd_fluxes",
std::function(
@@ -775,6 +793,13 @@ SchemeModule::SchemeModule()
));
+ this->_addBuiltinFunction("compute_dt", std::function(
+
+ [](const std::shared_ptr<const DiscreteFunctionVariant>& u,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
+ return compute_dt(u, c);
+ }));
+
this->_addBuiltinFunction("cell_volume",
std::function(
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index e285dbde1..260582273 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -4,6 +4,16 @@ add_library(
PugsScheme
AcousticSolver.cpp
AcousticCompositeSolver.cpp
+# JacobianAndStructuralInfoForSystemsofEquations.cpp
+# ApproximateRiemannCompositeSolver.cpp
+# ApproximateRiemannCompositeSolver_FluxForm.cpp
+# ApproximateRiemannCompositeSolver_ViscousForm.cpp
+# Roe_FluxForm_v2_CompositeSolver.cpp
+# Roe_ViscousForm_v2_CompositeSolver.cpp
+# VFFC_FluxForm_v1_CompositeSolver.cpp
+# VFFC_ViscousForm_v1_CompositeSolver.cpp
+# VFFC_FluxForm_v2_CompositeSolver.cpp
+# VFFC_ViscousForm_v2_CompositeSolver.cpp
HyperelasticSolver.cpp
DiscreteFunctionIntegrator.cpp
DiscreteFunctionInterpoler.cpp
diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp
index 633c324cd..e2c9f5af5 100644
--- a/src/scheme/RusanovEulerianCompositeSolver.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver.cpp
@@ -17,6 +17,94 @@
#include <variant>
+template <class Rd>
+double
+EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // 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) * norme_normal, std::fabs(uscaln + c_mean) * norme_normal);
+}
+
+double
+compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
+{
+ const auto& c = c_v->get<DiscreteFunctionP0<const double>>();
+
+ return std::visit(
+ [&](auto&& p_mesh) -> double {
+ const auto& mesh = *p_mesh;
+
+ using MeshType = mesh_type_t<decltype(p_mesh)>;
+ static constexpr size_t Dimension = MeshType::Dimension;
+ using Rd = TinyVector<Dimension>;
+
+ const auto& u = u_v->get<DiscreteFunctionP0<const Rd>>();
+
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+ // const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
+
+ const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
+ const NodeValuePerCell<const Rd> njr = MeshDataManager::instance().getMeshData(mesh).njr();
+
+ const FaceValuePerCell<const Rd> Cjf = MeshDataManager::instance().getMeshData(mesh).Cjf();
+ const FaceValuePerCell<const Rd> njf = MeshDataManager::instance().getMeshData(mesh).njf();
+
+ const EdgeValuePerCell<const Rd> Cje = MeshDataManager::instance().getMeshData(mesh).Cje();
+ const EdgeValuePerCell<const Rd> nje = MeshDataManager::instance().getMeshData(mesh).nje();
+
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ const auto& cell_to_edge_matrix = mesh.connectivity().cellToEdgeMatrix();
+ const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+
+ CellValue<double> local_dt{mesh.connectivity()};
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ double maxm(0);
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ const Rd normalCjr = Cjr(j, l);
+ // maxm = std::max(maxm, EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u, c, normalCjr));
+ maxm += EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u[j], c[j], normalCjr);
+ }
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const Rd normalCjf = Cjf(j, l);
+ // maxm = std::max(maxm, EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u, c, normalCjr));
+ maxm += EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u[j], c[j], normalCjf);
+ }
+
+ if constexpr (MeshType::Dimension == 3) {
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ const Rd normalCje = Cje(j, l);
+ // maxm = std::max(maxm, EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u, c, normalCjr));
+ maxm += EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(u[j], c[j], normalCje);
+ }
+ }
+
+ local_dt[j] = Vj[j] / maxm; //(Sj[j] * c[j]);
+ });
+
+ return min(local_dt);
+ } else {
+ throw NormalError("unexpected mesh type");
+ }
+ },
+ c.meshVariant()->variant());
+}
+
template <MeshConcept MeshTypeT>
class RusanovEulerianCompositeSolver
{
@@ -28,6 +116,9 @@ class RusanovEulerianCompositeSolver
using Rdxd = TinyMatrix<Dimension>;
using Rd = TinyVector<Dimension>;
+ using Rpxp = TinyMatrix<Dimension + 2>;
+ using Rp = TinyVector<Dimension + 2>;
+
class SymmetryBoundaryCondition;
class InflowListBoundaryCondition;
class OutflowBoundaryCondition;
@@ -136,6 +227,21 @@ class RusanovEulerianCompositeSolver
}
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));
+ }
+
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
@@ -146,16 +252,444 @@ class RusanovEulerianCompositeSolver
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
+ 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 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()};
+
+ //
+ // 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;
+
+ 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();
+
+ NodeValue<Rpxp> ViscosityNodeMatrix{p_mesh->connectivity()};
+ ViscosityNodeMatrix.fill(identity);
+
+ parallel_for(
+ p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId l) {
+ double maxabsVpNormCjr = 0;
+ 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);
+ // double rhoNode = 0;
+ Rd uNode = zero;
+ // double ENode = 0;
+ double cNode = 0;
+
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ // rhoNode += rho_n[J];
+ uNode += u_n[J];
+ // ENode += E_n[J];
+ cNode += c_n[J];
+ }
+ // rhoNode /= node_to_cell.size();
+ uNode *= 1. / node_to_cell.size();
+ // ENode /= node_to_cell.size();
+ cNode /= node_to_cell.size();
+
+ // Boundary conditions ..
+
+ for (size_t j = 0; j < node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ const unsigned int R = node_local_number_in_its_cells[j];
+
+ const Rd& Cjr_loc = Cjr(J, R);
+ maxabsVpNormCjr = std::max(maxabsVpNormCjr,
+ EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // rhoNode,
+ uNode, // ENode,
+ cNode, Cjr_loc));
+ }
+ ViscosityNodeMatrix[l] *= maxabsVpNormCjr;
+ });
+
+ FaceValue<Rpxp> ViscosityFaceMatrix{p_mesh->connectivity()};
+ ViscosityFaceMatrix.fill(identity);
+
+ parallel_for(
+ p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId l) {
+ double maxabsVpNormCjf = 0;
+ 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);
+ // double rhoFace = 0;
+ Rd uFace = zero;
+ // double EFace = 0;
+ double cFace = 0;
+
+ for (size_t j = 0; j < face_to_cell.size(); ++j) {
+ const CellId J = face_to_cell[j];
+ // rhoFace += rho_n[J];
+ uFace += u_n[J];
+ // EFace += E_n[J];
+ cFace += c_n[J];
+ }
+ // rhoFace /= face_to_cell.size();
+ uFace *= 1. / face_to_cell.size();
+ // EFace /= face_to_cell.size();
+ cFace /= face_to_cell.size();
+
+ // Boundary conditions ..
+
+ for (size_t j = 0; j < face_to_cell.size(); ++j) {
+ const CellId J = face_to_cell[j];
+ const unsigned int R = face_local_number_in_its_cells[j];
+
+ const Rd& Cjf_loc = Cjf(J, R);
+
+ maxabsVpNormCjf = std::max(maxabsVpNormCjf,
+ EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // rhoEdge,
+ uFace, // EEdge,
+ cFace, Cjf_loc));
+ }
+ ViscosityFaceMatrix[l] *= maxabsVpNormCjf;
+ });
+
+ // Computation of Viscosity Matrices at nodes and edges are done
+
+ // We may compute Gjr, Gje
+
+ 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);
+
+ 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();
+
+ 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 Rd& Cjr_loc = Cjr(j, l);
+ Rp statediff(zero);
+
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const CellId K = node_to_cell[k];
+
+ const Rd& u_Cjr = Flux_qtmvt[K] * Cjr_loc;
+
+ // const Rp&
+ statediff += State[j] - State[K];
+ // const Rp& diff = ViscosityNodeMatrix[node] * statediff;
+
+ Gjr[j][l][0] += 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_totnrj[K], Cjr_loc);
+
+ // Gjr[j][l] += diff;
+ }
+ Gjr[j][l] += ViscosityNodeMatrix[node] * statediff;
+
+ 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) << "\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 Rd& Cjf_loc = Cjf(j, l);
+ Rp statediff(zero);
+
+ for (size_t k = 0; k < face_to_cell.size(); ++k) {
+ const CellId K = face_to_cell[k];
+
+ const Rd& u_Cjf = Flux_qtmvt[K] * Cjf_loc;
+
+ // const Rp&
+ statediff += State[j] - State[K];
+ // const Rp& diff = ViscosityFaceMatrix[face] * statediff;
+
+ Gjf[j][l][0] += 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_totnrj[K], Cjf_loc);
+
+ // Gjf[j][l] += diff;
+ }
+
+ Gjf[j][l] += ViscosityFaceMatrix[face] * statediff;
+
+ 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) << "\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();
+
+ EdgeValue<Rpxp> ViscosityEdgeMatrix{p_mesh->connectivity()};
+ ViscosityEdgeMatrix.fill(identity);
+
+ parallel_for(
+ p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId l) {
+ double maxabsVpNormCje = 0;
+ 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);
+ // double rhoEdge = 0;
+ Rd uEdge = zero;
+ // double EEdge = 0;
+ double cEdge = 0;
+
+ for (size_t j = 0; j < edge_to_cell.size(); ++j) {
+ const CellId J = edge_to_cell[j];
+ // rhoEdge += rho_n[J];
+ uEdge += u_n[J];
+ // EEdge += E_n[J];
+ cEdge += c_n[J];
+ }
+ // rhoEdge /= edge_to_cell.size();
+ uEdge *= 1. / edge_to_cell.size();
+ // EEdge /= edge_to_cell.size();
+ cEdge /= edge_to_cell.size();
+
+ // Boundary conditions ..
+
+ for (size_t j = 0; j < edge_to_cell.size(); ++j) {
+ const CellId J = edge_to_cell[j];
+ const unsigned int R = edge_local_number_in_its_cells[j];
+
+ const Rd& Cje_loc = Cje(J, R);
+
+ maxabsVpNormCje = std::max(maxabsVpNormCje,
+ EvaluateMaxEigenValueTimesNormalLengthInGivenDirection( // rhoEdge,
+ uEdge, // EEdge,
+ cEdge, Cje_loc));
+ }
+ ViscosityEdgeMatrix[l] *= maxabsVpNormCje;
+ });
+
+ 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 Rd& Cje_loc = Cje(j, l);
+ Rp statediff(zero);
+
+ for (size_t k = 0; k < edge_to_cell.size(); ++k) {
+ const CellId K = edge_to_cell[k];
+
+ const Rd& u_Cje = Flux_qtmvt[K] * Cje_loc;
+
+ // const Rp&
+ statediff += State[j] - State[K];
+ // const Rp& diff = ViscosityEdgeMatrix[edge] * statediff;
+
+ Gje[j][l][0] += 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_totnrj[K], Cje_loc);
+
+ // Gje[j][l] += diff;
+ }
+ Gje[j][l] += ViscosityEdgeMatrix[edge] * statediff;
+
+ 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) << "\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));
@@ -460,7 +994,8 @@ rusanovEulerianCompositeSolver(
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 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) {
@@ -489,7 +1024,7 @@ rusanovEulerianCompositeSolver(
E_v->get<DiscreteFunctionP0<const double>>(),
c_v->get<DiscreteFunctionP0<const double>>(),
p_v->get<DiscreteFunctionP0<const double>>(),
- bc_descriptor_list, dt);
+ bc_descriptor_list, dt, check);
} else {
throw NormalError("RusanovEulerianCompositeSolver is only defined on polygonal meshes");
}
diff --git a/src/scheme/RusanovEulerianCompositeSolver.hpp b/src/scheme/RusanovEulerianCompositeSolver.hpp
index 69dfb6be8..357adf245 100644
--- a/src/scheme/RusanovEulerianCompositeSolver.hpp
+++ b/src/scheme/RusanovEulerianCompositeSolver.hpp
@@ -8,6 +8,9 @@
#include <memory>
#include <vector>
+double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+ const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>>
@@ -18,6 +21,7 @@ rusanovEulerianCompositeSolver(
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);
+ const double& dt,
+ const bool check = false);
#endif // RUSANOV_EULERIAN_COMPOSITE_SOLVER_HPP
--
GitLab