From f81ac89d7dea1d743164e1a05f54a7c54fd32f1b Mon Sep 17 00:00:00 2001
From: HOCH PHILIPPE <philippe.hoch@gmail.com>
Date: Wed, 21 Aug 2024 17:27:21 +0200
Subject: [PATCH] Rajout degree reconstruction polynome dans les solve et
codage effectif pour Rusanov_v2
---
src/language/modules/SchemeModule.cpp | 31 ++--
...eViscousFormEulerianCompositeSolver_v2.cpp | 4 +-
...eViscousFormEulerianCompositeSolver_v2.hpp | 1 +
src/scheme/RusanovEulerianCompositeSolver.cpp | 4 +-
src/scheme/RusanovEulerianCompositeSolver.hpp | 1 +
.../RusanovEulerianCompositeSolver_v2.cpp | 155 +++++++++++++-----
.../RusanovEulerianCompositeSolver_v2.hpp | 1 +
7 files changed, 143 insertions(+), 54 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 3c8fe1fcd..32f689003 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -512,13 +512,14 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
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, gamma, c, p, bc_descriptor_list, dt);
+ return rusanovEulerianCompositeSolver(rho, u, E, gamma, c, p, degree,
+ bc_descriptor_list, dt);
}
));
@@ -530,14 +531,14 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
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, gamma, c, p, bc_descriptor_list, dt,
- true);
+ return rusanovEulerianCompositeSolver(rho, u, E, gamma, c, p, degree,
+ bc_descriptor_list, dt, true);
}
));
@@ -548,14 +549,14 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
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_v2(rho, u, E, gamma, c, p, bc_descriptor_list,
- dt);
+ return rusanovEulerianCompositeSolver_v2(rho, u, E, gamma, c, p, degree,
+ bc_descriptor_list, dt);
}
));
@@ -566,14 +567,14 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
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_v2(rho, u, E, gamma, c, p, bc_descriptor_list, dt,
- true);
+ return rusanovEulerianCompositeSolver_v2(rho, u, E, gamma, c, p, degree,
+ bc_descriptor_list, dt, true);
}));
this->_addBuiltinFunction("roe_viscosityform_eulerian_composite_solver_version2",
std::function(
@@ -581,13 +582,13 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
- return roeViscousFormEulerianCompositeSolver_v2(rho, u, E, gamma, c, p,
+ return roeViscousFormEulerianCompositeSolver_v2(rho, u, E, gamma, c, p, degree,
bc_descriptor_list, dt);
}
@@ -599,13 +600,13 @@ SchemeModule::SchemeModule()
const std::shared_ptr<const DiscreteFunctionVariant>& u,
const std::shared_ptr<const DiscreteFunctionVariant>& E, const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
- const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const std::shared_ptr<const DiscreteFunctionVariant>& p, const size_t& degree,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
bc_descriptor_list,
const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
- return roeViscousFormEulerianCompositeSolver_v2(rho, u, E, gamma, c, p,
+ return roeViscousFormEulerianCompositeSolver_v2(rho, u, E, gamma, c, p, degree,
bc_descriptor_list, dt, true);
}
diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
index cfe2edb3f..3490fa4f1 100644
--- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp
@@ -983,6 +983,7 @@ class RoeViscousFormEulerianCompositeSolver_v2
const double& gamma,
const DiscreteFunctionP0<const double>& c_n,
const DiscreteFunctionP0<const double>& p_n,
+ const size_t& degree,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt,
const bool checkLocalConservation) const
@@ -2307,6 +2308,7 @@ roeViscousFormEulerianCompositeSolver_v2(
const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
+ const size_t& degree,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt,
const bool check)
@@ -2335,7 +2337,7 @@ roeViscousFormEulerianCompositeSolver_v2(
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);
+ p_v->get<DiscreteFunctionP0<const double>>(), degree, bc_descriptor_list, dt, check);
} else {
throw NormalError("RoeViscousFormEulerianCompositeSolver v2 is only defined on polygonal meshes");
}
diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.hpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.hpp
index a277ae26d..359b8cc9d 100644
--- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.hpp
+++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.hpp
@@ -22,6 +22,7 @@ roeViscousFormEulerianCompositeSolver_v2(
const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const size_t& degree,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt,
const bool check = false);
diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp
index dd0f65402..05ce4a9aa 100644
--- a/src/scheme/RusanovEulerianCompositeSolver.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver.cpp
@@ -608,6 +608,7 @@ class RusanovEulerianCompositeSolver
const double& gamma,
const DiscreteFunctionP0<const double>& c_n,
const DiscreteFunctionP0<const double>& p_n,
+ const size_t& degree,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt,
const bool checkLocalConservation) const
@@ -1876,6 +1877,7 @@ rusanovEulerianCompositeSolver(
const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
+ const size_t& degree,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt,
const bool check)
@@ -1907,7 +1909,7 @@ rusanovEulerianCompositeSolver(
E_v->get<DiscreteFunctionP0<const double>>(), gamma,
c_v->get<DiscreteFunctionP0<const double>>(),
p_v->get<DiscreteFunctionP0<const double>>(),
- bc_descriptor_list, dt, check);
+ degree, 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 593b6f763..9f1c131b3 100644
--- a/src/scheme/RusanovEulerianCompositeSolver.hpp
+++ b/src/scheme/RusanovEulerianCompositeSolver.hpp
@@ -22,6 +22,7 @@ rusanovEulerianCompositeSolver(
const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const size_t& degree,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt,
const bool check = false);
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
index 5111335a2..5e8bc7ffe 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp
@@ -1123,6 +1123,7 @@ class RusanovEulerianCompositeSolver_v2
const double& gamma,
const DiscreteFunctionP0<const double>& c_n,
const DiscreteFunctionP0<const double>& p_n,
+ const size_t degree,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt,
const bool checkLocalConservation) const
@@ -1143,7 +1144,7 @@ class RusanovEulerianCompositeSolver_v2
symmetry_boundary_descriptor_list.push_back(bc_descriptor->boundaryDescriptor_shared());
}
}
- static const int degree = 1;
+ // static const int degree = 1;
PolynomialReconstructionDescriptor
reconstruction_descriptor(IntegrationMethodType::cell_center, // IntegrationMethodType::boundary,
degree, symmetry_boundary_descriptor_list);
@@ -1176,9 +1177,6 @@ class RusanovEulerianCompositeSolver_v2
// DiscreteFunctionDPk rho_u_bar = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const Rd>>();
// DiscreteFunctionDPk rho_E_bar = reconstructions[2]->template get<DiscreteFunctionDPk<Dimension, const double>>();
- auto rho_bar_lim = rho_bar;
- auto rho_u_bar_lim = u_bar;
- auto rho_E_bar_lim = eps_bar;
CellValue<double> Limitor_rho(p_mesh->connectivity());
CellValue<double> Limitor_u(p_mesh->connectivity());
@@ -1204,14 +1202,55 @@ class RusanovEulerianCompositeSolver_v2
});
// Appels des limiteurs cell by cell
+ // Pour rho (plus ou moins classique)
computeLimitorVolumicScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, Limitor_rho);
+ // Pour les variables d'interet (physique)
+ // Pour eps ici G.P. (avec Leibniz et toute l artillerie)
computeLimitorInternalNrjScalarQuantityMinModDukowiczGradient(*p_mesh, valrho, grad_rho_cell, valeps, grad_eps_cell,
// valu,
grad_u_cell, Limitor_eps);
+ // Pour u : relation deduite du principe d'induction
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Limitor_u[j] = sqrt(Limitor_eps[j]); });
+ //
+ // Creation des variables conservatives limitées
+ //
+ auto rho_bar_lim = copy(rho_bar);
+ auto rho_u_bar_lim = copy(u_bar); // a modifier .. du a la densite
+ auto rho_E_bar_lim = copy(eps_bar); // a modifier et a completer : densité et energie cinetique
+
+ // Assemblage des ctes et ordres élevés
+ parallel_for(
+ p_mesh->numberOfCells(),
+ PUGS_LAMBDA(CellId j) { // rho_bar_lim[j] = rho[j] + Limitor_rho[j] * (rho_bar[j](x) - rho[j]);
+ // rho_u_bar_lim.coefficients(j)[0] = valrho[j] * valu[j]; // * valu[j];
+ // rho_E_bar_lim.coefficients(j)[0] = valrho[j] * (valeps[j] + .5 * dot(valu[j], valu[j]));
+
+ // rho_bar_lim.coefficients(j)[0] // cte inchangée
+ rho_u_bar_lim.coefficients(j)[0] *= valrho[j]; // du coup (rho u)
+ rho_E_bar_lim.coefficients(j)[0] *= valrho[j]; // du coup (rho*epsilon)
+ rho_E_bar_lim.coefficients(j)[0] += valrho[j] * .5 * dot(valu[j], valu[j]);
+ // Rdxd grad_rhou_lim=
+ // valrho[j]*Limitor_u[j]*grad_u_cell[j]+Limitor_rho[j]*tensorProduct(valu[j],grad_rho_cell[j]);
+ //
+ Rdxd gradU_lim_transpose = Limitor_u[j] * transpose(grad_u_cell[j]);
+ const Rd gradUtu_lim = gradU_lim_transpose * valu[j];
+
+ for (size_t dim = 0; dim < Dimension; ++dim) {
+ rho_bar_lim.coefficients(j)[1 + dim] *= Limitor_rho[j];
+ rho_u_bar_lim.coefficients(j)[1 + dim] *= Limitor_u[j] * valrho[j]; // c'est la partie gradu limite
+ rho_u_bar_lim.coefficients(j)[1 + dim] += Limitor_rho[j] * grad_rho_cell[j][dim] * valu[j];
+
+ rho_E_bar_lim.coefficients(j)[1 + dim] *= Limitor_eps[j] * valrho[j];
+ rho_E_bar_lim.coefficients(j)[1 + dim] += gradUtu_lim[dim] * valrho[j];
+
+ rho_E_bar_lim.coefficients(j)[1 + dim] +=
+ (valeps[j] + .5 * dot(valu[j], valu[j])) * Limitor_rho[j] * grad_rho_cell[j][dim];
+ }
+ });
+
// auto c = copy(c_n);
// auto p = copy(p_n);
@@ -1237,18 +1276,64 @@ class RusanovEulerianCompositeSolver_v2
//
// Remplir ici les reconstructions d'ordre élevé
+ //
+ 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();
+
+ const auto xr = p_mesh->xr();
+ auto xl = MeshDataManager::instance().getMeshData(*p_mesh).xl();
+ auto xe = MeshDataManager::instance().getMeshData(*p_mesh).xe();
+
NodeValuePerCell<Rp> StateAtNode{p_mesh->connectivity()};
StateAtNode.fill(zero);
parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtNode[j].fill(State[j]); });
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // StateAtNode[j].fill(State[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];
+
+ StateAtNode[j][l][0] = rho_bar_lim[j](xr[node]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtNode[j][l][1 + dim] = rho_u_bar_lim[j](xr[node])[dim];
+ StateAtNode[j][l][1 + Dimension] = rho_E_bar_lim[j](xr[node]);
+ }
+ });
+
EdgeValuePerCell<Rp> StateAtEdge{p_mesh->connectivity()};
StateAtEdge.fill(zero);
parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtEdge[j].fill(State[j]); });
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // eStateAtEdge[j].fill(State[j]);
+ 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];
+
+ StateAtEdge[j][l][0] = rho_bar_lim[j](xe[edge]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtEdge[j][l][1 + dim] = rho_u_bar_lim[j](xe[edge])[dim];
+ StateAtEdge[j][l][1 + Dimension] = rho_E_bar_lim[j](xe[edge]);
+ }
+ });
FaceValuePerCell<Rp> StateAtFace{p_mesh->connectivity()};
StateAtFace.fill(zero);
parallel_for(
- p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { StateAtFace[j].fill(State[j]); });
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ // StateAtFace[j].fill(State[j]);
+ 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];
+
+ StateAtFace[j][l][0] = rho_bar_lim[j](xl[face]);
+ for (size_t dim = 0; dim < Dimension; ++dim)
+ StateAtFace[j][l][1 + dim] = rho_u_bar_lim[j](xl[face])[dim];
+ StateAtFace[j][l][1 + Dimension] = rho_E_bar_lim[j](xl[face]);
+ }
+ });
// Conditions aux limites des etats conservatifs
@@ -1310,31 +1395,30 @@ class RusanovEulerianCompositeSolver_v2
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) {
+ 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];
@@ -1344,10 +1428,6 @@ class RusanovEulerianCompositeSolver_v2
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);
@@ -2382,6 +2462,7 @@ rusanovEulerianCompositeSolver_v2(
const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c_v,
const std::shared_ptr<const DiscreteFunctionVariant>& p_v,
+ const size_t& degree,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt,
const bool check)
@@ -2410,7 +2491,7 @@ rusanovEulerianCompositeSolver_v2(
return RusanovEulerianCompositeSolver_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);
+ p_v->get<DiscreteFunctionP0<const double>>(), degree, bc_descriptor_list, dt, check);
} else {
throw NormalError("RusanovEulerianCompositeSolver v2 is only defined on polygonal meshes");
}
diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.hpp b/src/scheme/RusanovEulerianCompositeSolver_v2.hpp
index 749c09a18..4b5864abe 100644
--- a/src/scheme/RusanovEulerianCompositeSolver_v2.hpp
+++ b/src/scheme/RusanovEulerianCompositeSolver_v2.hpp
@@ -22,6 +22,7 @@ rusanovEulerianCompositeSolver_v2(
const double& gamma,
const std::shared_ptr<const DiscreteFunctionVariant>& c,
const std::shared_ptr<const DiscreteFunctionVariant>& p,
+ const size_t& degree,
const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
const double& dt,
const bool check = false);
--
GitLab