diff --git a/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.cpp
index e5eab006bc109462784b1c592131fbe15ae1d7e0..8a3a96653a3b2648a08f32230c092b4c2b7127a3 100644
--- a/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.cpp
+++ b/src/scheme/RoeFluxFormEulerianCompositeSolver_v2.cpp
@@ -202,7 +202,7 @@ class RoeFluxFormEulerianCompositeSolver_v2
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<OutflowBoundaryCondition, T>) {
- std::cout << " Traitement Outflow \n";
+ // std::cout << " Traitement Outflow \n";
// const Rd& normal = bc.outgoingNormal();
/*
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
@@ -322,7 +322,7 @@ class RoeFluxFormEulerianCompositeSolver_v2
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";
+ // std::cout << " Traitement SYMMETRY \n";
const Rd& normal = bc.outgoingNormal();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
@@ -431,7 +431,7 @@ class RoeFluxFormEulerianCompositeSolver_v2
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";
+ // std::cout << " Traitement WALL \n";
const Rd& normal = bc.outgoingNormal();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
@@ -534,7 +534,7 @@ class RoeFluxFormEulerianCompositeSolver_v2
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<WallBoundaryCondition, T>) {
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- std::cout << " Traitement WALL local (non flat) \n";
+ // std::cout << " Traitement WALL local (non flat) \n";
// const Rd& normal = bc.outgoingNormal();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
@@ -719,7 +719,7 @@ class RoeFluxFormEulerianCompositeSolver_v2
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";
+ // std::cout << " Traitement INFLOW \n";
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
@@ -3081,7 +3081,7 @@ class RoeFluxFormEulerianCompositeSolver_v2_with_source
[&](auto&& bc) {
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<OutflowBoundaryCondition, T>) {
- std::cout << " Traitement Outflow \n";
+ // std::cout << " Traitement Outflow \n";
// const Rd& normal = bc.outgoingNormal();
/*
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
@@ -3201,7 +3201,7 @@ class RoeFluxFormEulerianCompositeSolver_v2_with_source
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";
+ // std::cout << " Traitement SYMMETRY \n";
const Rd& normal = bc.outgoingNormal();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
@@ -3310,7 +3310,7 @@ class RoeFluxFormEulerianCompositeSolver_v2_with_source
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";
+ // std::cout << " Traitement WALL \n";
const Rd& normal = bc.outgoingNormal();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
@@ -3413,7 +3413,7 @@ class RoeFluxFormEulerianCompositeSolver_v2_with_source
using T = std::decay_t<decltype(bc)>;
if constexpr (std::is_same_v<WallBoundaryCondition, T>) {
MeshData<MeshType>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
- std::cout << " Traitement WALL local (non flat) \n";
+ // std::cout << " Traitement WALL local (non flat) \n";
// const Rd& normal = bc.outgoingNormal();
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
@@ -3598,7 +3598,7 @@ class RoeFluxFormEulerianCompositeSolver_v2_with_source
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";
+ // std::cout << " Traitement INFLOW \n";
const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
@@ -5051,17 +5051,250 @@ class RoeFluxFormEulerianCompositeSolver_v2_with_source
});
+/*
+ const auto xr = p_mesh->xr(); // coordonnée noeud
+ const auto xj = mesh_data.xj(); // coordonnée centre maille
+ const auto xf = mesh_data.xl(); // coordonnée centre face
+ const auto xe = mesh_data.xe(); // coordonnée centre ?
+
+ CellValue<Rp> Sources{p_mesh->connectivity()};
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+
+ Sources[j][0] -= rho_source[j];
+ for (size_t d = 0; d < Dimension; ++d)
+ Sources[j][d + 1]-= u_source[j][d];
+ Sources[j][1 + Dimension] -= E_source[j];
+ });
+
+
+
+ const auto& edge_to_cell_matrix = p_mesh->connectivity().edgeToCellMatrix();
+ const auto& edge_local_numbers_in_their_cells = p_mesh->connectivity().edgeLocalNumbersInTheirCells();
+
+ CellValue<Rp> Sj{p_mesh->connectivity()};
+ CellValue<Rp> Sj_gauche{p_mesh->connectivity()};
+ CellValue<Rp> Sj_droite{p_mesh->connectivity()};
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ std::cout << "Traitement cellule " << j << std::endl;
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ // double Somme_1_Nr = 0.;
+
+
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ std::cout << "Traitement Noeud" << std::endl;
+ CellValue<Rp> temp_gauche{p_mesh->connectivity()};
+
+ const NodeId& node = cell_to_node[l];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node);
+ const auto& node_to_cell = node_to_cell_matrix[node];
+
+ for (size_t k = 0; k < node_to_cell.size(); ++k) {
+ const size_t R = node_local_number_in_its_cells[k];
+ const CellId& i_cell = node_to_cell[k];
+
+ const Rpxp Lambda_rl = Lambda_rj[node][i_cell];
+ Rpxp MatriceIDPlusLambda_ri = identity;
+ MatriceIDPlusLambda_ri += Lambda_rl;
+
+ const auto Sub_volume_ir = (1 - theta)/2 * dot((xr[node] - xj[i_cell]),Cjr(i_cell,R));
+
+ temp_gauche[j] += Sub_volume_ir * MatriceIDPlusLambda_ri * Sources[i_cell];
+ }
+
+
+ // Somme_1_Nr += 1/node_to_cell_matrix[node].size();
+ Sj_gauche[j] += 1/node_to_cell_matrix[node].size() * temp_gauche[j];
+ }
+
+
+
+ const auto& cell_to_edge = cell_to_edge_matrix[j];
+ // double Somme_1_2 = 0.;
+
+ for (size_t l = 0; l < cell_to_edge.size(); ++l) {
+ CellValue<Rp> temp_droite{p_mesh->connectivity()};
+
+ 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);
+
+ for (size_t k = 0; k < edge_to_cell.size(); ++k) {
+ const size_t Ed = edge_local_number_in_its_cells[k];
+ const CellId& i_cell = edge_to_cell[k];
+
+ const Rpxp Lambda_el = Lambda_ej[edge][i_cell];
+ Rpxp MatriceIDPlusLambda_ei = identity;
+ MatriceIDPlusLambda_ei += Lambda_el;
+
+ const auto Sub_volume_ie = theta/2 * dot((xe[edge] - xj[i_cell]),Cje(i_cell,Ed));
+
+ temp_droite[j] += Sub_volume_ie * MatriceIDPlusLambda_ei * Sources[i_cell];
+
+ }
+
+ // Somme_1_2 += 1/2;
+
+ Sj_droite[j] += 1/2 * temp_droite[j];
+ }
- std::cout << rho_source[j] << " / " << u_source[j] << " / " << E_source[j] << std::endl;
+
- SumFluxes[j][0] -= rho_source[j] * volumes_cell[j];//* volumes_cell[j];
+ Sj[j] = Sj_gauche[j] + Sj_droite[j];
+ }
+ );
+
+
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+
+ SumFluxes[j][0] -= Sj[j][0] * volumes_cell[j];
for (size_t d = 0; d < Dimension; ++d)
- SumFluxes[j][d + 1]-= u_source[j][d] * volumes_cell[j];//* volumes_cell[j];
- SumFluxes[j][1 + Dimension] -= E_source[j] * volumes_cell[j];//* volumes_cell[j];
+ SumFluxes[j][d + 1]-= Sj[j][d+1] * volumes_cell[j];
+ SumFluxes[j][1 + Dimension] -= Sj[j][1 + Dimension] * volumes_cell[j];
});
-
+
+
+*/
+
+/*
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+
+ SumFluxes[j][0] -= rho_source[j] * volumes_cell[j];
+ for (size_t d = 0; d < Dimension; ++d)
+ SumFluxes[j][d + 1]-= u_source[j][d] * volumes_cell[j];
+ SumFluxes[j][1 + Dimension] -= E_source[j] * volumes_cell[j];
+ });
+*/
+
+CellValue<Rp> Sj{p_mesh->connectivity()};
+
+// Classical discretisation
+parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ Sj[j][0] = rho_source[j];
+ for (size_t d = 0; d < Dimension; ++d)
+ Sj[j][d + 1] = u_source[j][d];
+ Sj[j][1 + Dimension] = E_source[j];
+ });
+
+bool SubVolumeUpwinding = false; // classical
+// SubVolumeUpwinding = true;
+
+if (SubVolumeUpwinding) {
+ // Discretisation of subvolumes and decentering
+ NodeValue<Rp> S_node{p_mesh->connectivity()};
+ FaceValue<Rp> S_face{p_mesh->connectivity()};
+ EdgeValue<Rp> S_edge{p_mesh->connectivity()};
+
+ const auto xr = p_mesh->xr();
+
+ const auto xj = mesh_data.xj();
+ const auto xf = mesh_data.xl();
+ const auto xe = mesh_data.xe();
+
+ parallel_for(
+ p_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
+ const auto& node_to_cell = node_to_cell_matrix[r];
+ S_node[r] = zero;
+
+ for (size_t l = 0; l < node_to_cell.size(); ++l) {
+ const CellId& cell = node_to_cell[l];
+ const size_t R = node_local_number_in_its_cells[l];
+
+ const Rpxp Lambda_rl = Lambda_rj[r][l];
+ Rpxp MatriceIDPlusLambda_ri = identity;
+ MatriceIDPlusLambda_ri += Lambda_rl;
+
+ const auto Sub_volume_ir = (1 - theta) / Dimension * dot((xr[r] - xj[cell]), Cjr(cell, R));
+
+ S_node[r] += Sub_volume_ir * MatriceIDPlusLambda_ri * Sj[cell];
+ }
+ S_node[r] *= 1. / node_to_cell.size();
+ });
+
+ parallel_for(
+ p_mesh->numberOfFaces(), PUGS_LAMBDA(FaceId f) {
+ const auto& face_local_number_in_its_cells = face_local_numbers_in_their_cells.itemArray(f);
+ const auto& face_to_cell = face_to_cell_matrix[f];
+ S_face[f] = zero;
+
+ for (size_t l = 0; l < face_to_cell.size(); ++l) {
+ const CellId& cell = face_to_cell[l];
+ const size_t R = face_local_number_in_its_cells[l];
+
+ const Rpxp Lambda_fl = Lambda_jf[cell][R];
+ Rpxp MatriceIDPlusLambda_fi = identity;
+ MatriceIDPlusLambda_fi += Lambda_fl;
+
+ const auto Sub_volume_if = (theta) / Dimension * dot((xf[f] - xj[cell]), Cjf(cell, R));
+
+ S_face[f] += Sub_volume_if * MatriceIDPlusLambda_fi * Sj[cell];
+ }
+ S_face[f] *= 1. / face_to_cell.size();
+ });
+
+ 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();
+ parallel_for(
+ p_mesh->numberOfEdges(), PUGS_LAMBDA(EdgeId e) {
+ const auto& edge_local_number_in_its_cells = edge_local_numbers_in_their_cells.itemArray(e);
+ const auto& edge_to_cell = edge_to_cell_matrix[e];
+ S_edge[e] = zero;
+
+ for (size_t l = 0; l < edge_to_cell.size(); ++l) {
+ const CellId& cell = edge_to_cell[l];
+ const size_t R = edge_local_number_in_its_cells[l];
+
+ const Rpxp Lambda_el = Lambda_ej[e][l];
+ Rpxp MatriceIDPlusLambda_ei = identity;
+ MatriceIDPlusLambda_ei += Lambda_el;
+
+ const auto Sub_volume_ie = (eta) / Dimension * dot((xe[e] - xj[cell]), Cje(cell, R));
+
+ S_edge[e] += Sub_volume_ie * MatriceIDPlusLambda_ei * Sj[cell];
+ }
+ S_edge[e] *= 1. / edge_to_cell.size();
+ });
+ }
+
+ // Discretization of source terms as a balance of (continuous) flux sources.
+
+ parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const auto& cell_to_node = cell_to_node_matrix[j];
+ const auto& cell_to_face = cell_to_face_matrix[j];
+
+ Sj[j] = zero;
+ for (size_t l = 0; l < cell_to_node.size(); ++l) {
+ const NodeId& node = cell_to_node[l];
+ Sj[j] += S_node[node];
+ }
+
+ for (size_t l = 0; l < cell_to_face.size(); ++l) {
+ const FaceId& face = cell_to_face[l];
+ Sj[j] += S_face[face];
+ }
+
+ 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 EdgeId& edge = cell_to_edge[l];
+ Sj[j] += S_edge[edge];
+ }
+ }
+ Sj[j] *= 1. / volumes_cell[j];
+ });
+}
+
+parallel_for(
+ p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { SumFluxes[j] -= volumes_cell[j] * Sj[j]; });
parallel_for(
p_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {