diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp index b6d8be8b0a5ddd73e569d63277211108dade575a..ce30fc3520ff3d9239a8cbbec3f539fcdfaf6466 100644 --- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp +++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2.cpp @@ -1162,6 +1162,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 _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); _applyWallBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); // @@ -1435,6 +1436,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 Rpxp ViscosityMatrixJK = .5 * (nrmCjr * JacInfoJ.AbsJacobian + nrmCkr * JacInfoK.AbsJacobian); // Test Rajout Viscosité type Rusanov v2 si chgt signe vp .. + bool anySignChgJ = EvaluateChangingSignVpAlongNormal( // rho_n[j], u_n[j], // E_n[j], gamma, c_n[j], // p_n[j], rho_n[K], @@ -1452,7 +1454,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 bool anySignDif = false; if (anySignChgJ or anySignChgK) - anySignDif = true; + anySignDif = true; // false; // true; if (anySignDif) { Rpxp AddedViscosity(identity); @@ -1479,6 +1481,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 Gjr[j][l] *= 1. / node_to_cell.size(); } }); + synchronize(Gjr); if (checkLocalConservation) { auto is_boundary_node = p_mesh->connectivity().isBoundaryNode(); @@ -1493,13 +1496,15 @@ class RoeViscousFormEulerianCompositeSolver_v2 if (not is_boundary_node[l]) { Rp SumGjr(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjr[K][R])); } // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr)); - MaxErrorConservationNode[l] = l2Norm(SumGjr); + MaxErrorConservationNode[l] = l2Norm(SumGjr) / maxGjrAbs; } }); std::cout << " Max Error Node " << max(MaxErrorConservationNode) << " Max Error RoeMatrice " << max(MaxErrNode) @@ -1572,7 +1577,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 bool anySignDif = false; if (anySignChgJ or anySignChgK) - anySignDif = true; + anySignDif = true; // false; // true; if (anySignDif) { Rpxp AddedViscosity(identity); @@ -1599,6 +1604,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 Gjf[j][l] *= 1. / face_to_cell.size(); } }); + synchronize(Gjf); if (checkLocalConservation) { auto is_boundary_face = p_mesh->connectivity().isBoundaryFace(); @@ -1613,12 +1619,14 @@ class RoeViscousFormEulerianCompositeSolver_v2 if (not is_boundary_face[l]) { Rp SumGjf(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjf[K][R])); } - MaxErrorConservationFace[l] = l2Norm(SumGjf); + MaxErrorConservationFace[l] = l2Norm(SumGjf) / maxGjrAbs; // MaxErrorConservationFace = std::max(MaxErrorConservationFace, l2Norm(SumGjf)); } }); @@ -1702,7 +1710,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 bool anySignDif = false; if (anySignChgJ or anySignChgK) - anySignDif = true; + anySignDif = true; // false; // true; if (anySignDif) { Rpxp AddedViscosity(identity); @@ -1729,6 +1737,7 @@ class RoeViscousFormEulerianCompositeSolver_v2 Gje[j][l] *= 1. / edge_to_cell.size(); } }); + synchronize(Gje); if (checkLocalConservation) { auto is_boundary_edge = p_mesh->connectivity().isBoundaryEdge(); @@ -1743,13 +1752,15 @@ class RoeViscousFormEulerianCompositeSolver_v2 if (not is_boundary_edge[l]) { Rp SumGje(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gje[K][R])); } // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2Norm(SumGje)); - MaxErrorConservationEdge[l] = l2Norm(SumGje); + MaxErrorConservationEdge[l] = l2Norm(SumGje) / maxGjrAbs; } }); std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << " Max Error RoeMatrice " << max(MaxErrEdge) @@ -1761,8 +1772,8 @@ class RoeViscousFormEulerianCompositeSolver_v2 double theta = 2. / 3.; //.5; double eta = 1. / 6.; //.2; if constexpr (Dimension == 2) { - eta = 0; - theta = 1; // pour schema aux aretes + eta = 0; + // theta = 1; // pour schema aux aretes // theta=0; //pour schema aux noeuds } diff --git a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp index 76251d822d80e812cd70ad89923f48e243190267..d22e940a9c3d53575747cf3d09bb1160377f209a 100644 --- a/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp +++ b/src/scheme/RoeViscousFormEulerianCompositeSolver_v2_o2.cpp @@ -1336,6 +1336,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 _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); _applyWallBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); // @@ -1626,7 +1627,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 bool anySignDif = false; if (anySignChgJ or anySignChgK) - anySignDif = true; + anySignDif = true; // false; // true; if (anySignDif) { Rpxp AddedViscosity(identity); @@ -1667,13 +1668,15 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 if (not is_boundary_node[l]) { Rp SumGjr(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjr[K][R])); } // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr)); - MaxErrorConservationNode[l] = l2Norm(SumGjr); + MaxErrorConservationNode[l] = l2Norm(SumGjr) / maxGjrAbs; } }); std::cout << " Max Error Node " << max(MaxErrorConservationNode) << " Max Error RoeMatrice " << max(MaxErrNode) @@ -1746,7 +1749,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 bool anySignDif = false; if (anySignChgJ or anySignChgK) - anySignDif = true; + anySignDif = true; // false; // true; if (anySignDif) { Rpxp AddedViscosity(identity); @@ -1788,12 +1791,14 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 if (not is_boundary_face[l]) { Rp SumGjf(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjf[K][R])); } - MaxErrorConservationFace[l] = l2Norm(SumGjf); + MaxErrorConservationFace[l] = l2Norm(SumGjf) / maxGjrAbs; // MaxErrorConservationFace = std::max(MaxErrorConservationFace, l2Norm(SumGjf)); } }); @@ -1877,7 +1882,7 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 bool anySignDif = false; if (anySignChgJ or anySignChgK) - anySignDif = true; + anySignDif = true; // false; // true; if (anySignDif) { Rpxp AddedViscosity(identity); @@ -1919,13 +1924,15 @@ class RoeViscousFormEulerianCompositeSolver_v2_o2 if (not is_boundary_edge[l]) { Rp SumGje(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gje[K][R])); } // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2Norm(SumGje)); - MaxErrorConservationEdge[l] = l2Norm(SumGje); + MaxErrorConservationEdge[l] = l2Norm(SumGje) / maxGjrAbs; } }); std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << " Max Error RoeMatrice " << max(MaxErrEdge) diff --git a/src/scheme/RusanovEulerianCompositeSolver.cpp b/src/scheme/RusanovEulerianCompositeSolver.cpp index afd27e7dc6374dd30d603c82f856eb7b0797372d..a13b34fb94582505b86683af9eed5efdca6f8af3 100644 --- a/src/scheme/RusanovEulerianCompositeSolver.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver.cpp @@ -12,9 +12,9 @@ #include <mesh/MeshNodeBoundary.hpp> #include <mesh/MeshTraits.hpp> #include <mesh/MeshVariant.hpp> +#include <mesh/SubItemValuePerItemUtils.hpp> #include <scheme/DiscreteFunctionUtils.hpp> #include <scheme/InflowListBoundaryConditionDescriptor.hpp> - #include <variant> template <MeshConcept MeshTypeT> @@ -786,6 +786,7 @@ class RusanovEulerianCompositeSolver _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); _applyWallBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); // @@ -1122,6 +1123,7 @@ class RusanovEulerianCompositeSolver Gjr[j][l] *= 1. / node_to_cell.size(); } }); + synchronize(Gjr); if (checkLocalConservation) { auto is_boundary_node = p_mesh->connectivity().isBoundaryNode(); @@ -1136,13 +1138,15 @@ class RusanovEulerianCompositeSolver if (not is_boundary_node[l]) { Rp SumGjr(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjr[K][R])); } // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr)); - MaxErrorConservationNode[l] = l2Norm(SumGjr); + MaxErrorConservationNode[l] = l2Norm(SumGjr) / maxGjrAbs; } }); std::cout << " Max Error Node " << max(MaxErrorConservationNode) << "\n"; @@ -1184,6 +1188,7 @@ class RusanovEulerianCompositeSolver Gjf[j][l] *= 1. / face_to_cell.size(); } }); + synchronize(Gjf); if (checkLocalConservation) { auto is_boundary_face = p_mesh->connectivity().isBoundaryFace(); @@ -1201,10 +1206,12 @@ class RusanovEulerianCompositeSolver if (not is_boundary_face[l]) { Rp SumGjf(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjf[K][R])); } /* std::cout << " maxErr face " << xl[l] << " err " << l2Norm(SumGjf) << " size face cell" @@ -1214,7 +1221,7 @@ class RusanovEulerianCompositeSolver if (l2Norm(SumGjf) > 1e-3) throw UnexpectedError("Pb consevvation face loc"); */ - MaxErrorConservationFace[l] = l2Norm(SumGjf); + MaxErrorConservationFace[l] = l2Norm(SumGjf) / maxGjrAbs; // MaxErrorConservationFace = std::max(MaxErrorConservationFace, l2Norm(SumGjf)); } }); @@ -1307,6 +1314,7 @@ class RusanovEulerianCompositeSolver Gje[j][l] *= 1. / edge_to_cell.size(); } }); + synchronize(Gje); if (checkLocalConservation) { auto is_boundary_edge = p_mesh->connectivity().isBoundaryEdge(); @@ -1321,13 +1329,15 @@ class RusanovEulerianCompositeSolver if (not is_boundary_edge[l]) { Rp SumGje(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gje[K][R])); } // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2Norm(SumGje)); - MaxErrorConservationEdge[l] = l2Norm(SumGje); + MaxErrorConservationEdge[l] = l2Norm(SumGje) / maxGjrAbs; } }); std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << "\n"; diff --git a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp index a27926f9cba7176aec6d40f649f76f0f690a28f8..867e6c4e5f7c2f9865820e954c26c219800f2cb9 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_o2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_o2.cpp @@ -12,6 +12,7 @@ #include <mesh/MeshNodeBoundary.hpp> #include <mesh/MeshTraits.hpp> #include <mesh/MeshVariant.hpp> +#include <mesh/SubItemValuePerItemUtils.hpp> #include <scheme/DiscreteFunctionDPk.hpp> #include <scheme/DiscreteFunctionDPkVariant.hpp> #include <scheme/DiscreteFunctionDPkVector.hpp> @@ -955,6 +956,8 @@ class RusanovEulerianCompositeSolver_o2 _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); + _applyWallBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); // // Construction du flux .. ok pour l'ordre 1 @@ -1290,6 +1293,7 @@ class RusanovEulerianCompositeSolver_o2 Gjr[j][l] *= 1. / node_to_cell.size(); } }); + synchronize(Gjr); if (checkLocalConservation) { auto is_boundary_node = p_mesh->connectivity().isBoundaryNode(); @@ -1304,13 +1308,15 @@ class RusanovEulerianCompositeSolver_o2 if (not is_boundary_node[l]) { Rp SumGjr(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjr[K][R])); } // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr)); - MaxErrorConservationNode[l] = l2Norm(SumGjr); + MaxErrorConservationNode[l] = l2Norm(SumGjr) / maxGjrAbs; } }); std::cout << " Max Error Node " << max(MaxErrorConservationNode) << "\n"; @@ -1352,6 +1358,7 @@ class RusanovEulerianCompositeSolver_o2 Gjf[j][l] *= 1. / face_to_cell.size(); } }); + synchronize(Gjf); if (checkLocalConservation) { auto is_boundary_face = p_mesh->connectivity().isBoundaryFace(); @@ -1369,10 +1376,12 @@ class RusanovEulerianCompositeSolver_o2 if (not is_boundary_face[l]) { Rp SumGjf(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjf[K][R])); } /* std::cout << " maxErr face " << xl[l] << " err " << l2Norm(SumGjf) << " size face cell" @@ -1382,7 +1391,7 @@ class RusanovEulerianCompositeSolver_o2 if (l2Norm(SumGjf) > 1e-3) throw UnexpectedError("Pb consevvation face loc"); */ - MaxErrorConservationFace[l] = l2Norm(SumGjf); + MaxErrorConservationFace[l] = l2Norm(SumGjf) / maxGjrAbs; // MaxErrorConservationFace = std::max(MaxErrorConservationFace, l2Norm(SumGjf)); } }); @@ -1475,6 +1484,7 @@ class RusanovEulerianCompositeSolver_o2 Gje[j][l] *= 1. / edge_to_cell.size(); } }); + synchronize(Gje); if (checkLocalConservation) { auto is_boundary_edge = p_mesh->connectivity().isBoundaryEdge(); @@ -1489,13 +1499,15 @@ class RusanovEulerianCompositeSolver_o2 if (not is_boundary_edge[l]) { Rp SumGje(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gje[K][R])); } // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2Norm(SumGje)); - MaxErrorConservationEdge[l] = l2Norm(SumGje); + MaxErrorConservationEdge[l] = l2Norm(SumGje) / maxGjrAbs; } }); std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << "\n"; diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp index 1f6ca8f3cdbafeefc3a6d070ae0d5fe0be16b8eb..cde234af1fcfe35fa66938cbc30b5a37942f3392 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_v2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_v2.cpp @@ -211,6 +211,7 @@ class RusanovEulerianCompositeSolver_v2 size_t node_local_number_in_cell = node_local_number_in_its_cells[i_cell]; normal += Cjr(node_cell_id, node_local_number_in_cell); } + normal *= 1. / node_cell_list.size(); normal *= 1. / l2Norm(normal); @@ -226,7 +227,7 @@ class RusanovEulerianCompositeSolver_v2 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]; + // stateNode[node_cell_id][node_local_number_in_cell][dim + 1] = 0; } } @@ -251,6 +252,7 @@ class RusanovEulerianCompositeSolver_v2 for (size_t dim = 0; dim < Dimension; ++dim) stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = vectorSym[dim]; + // stateFace[face_cell_id][face_local_number_in_cell][dim + 1] = 0; } if constexpr (Dimension == 3) { @@ -291,6 +293,7 @@ class RusanovEulerianCompositeSolver_v2 for (size_t dim = 0; dim < Dimension; ++dim) stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = vectorSym[dim]; + // stateEdge[edge_cell_id][edge_local_number_in_cell][dim + 1] = 0; } } } @@ -799,6 +802,7 @@ class RusanovEulerianCompositeSolver_v2 //_applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); _applySymmetricBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); _applyNeumannflatBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + _applyWallBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); // // Construction du flux .. ok pour l'ordre 1 @@ -1071,13 +1075,15 @@ class RusanovEulerianCompositeSolver_v2 if (not is_boundary_node[l]) { Rp SumGjr(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjr[K][R])); } // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr)); - MaxErrorConservationNode[l] = l2Norm(SumGjr); + MaxErrorConservationNode[l] = l2Norm(SumGjr) / maxGjrAbs; } }); std::cout << " Max Error Node " << max(MaxErrorConservationNode) << "\n"; @@ -1145,12 +1151,14 @@ class RusanovEulerianCompositeSolver_v2 if (not is_boundary_face[l]) { Rp SumGjf(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjf[K][R])); } - MaxErrorConservationFace[l] = l2Norm(SumGjf); + MaxErrorConservationFace[l] = l2Norm(SumGjf) / maxGjrAbs; // MaxErrorConservationFace = std::max(MaxErrorConservationFace, l2Norm(SumGjf)); } }); @@ -1228,13 +1236,15 @@ class RusanovEulerianCompositeSolver_v2 if (not is_boundary_edge[l]) { Rp SumGje(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gje[K][R])); } // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2Norm(SumGje)); - MaxErrorConservationEdge[l] = l2Norm(SumGje); + MaxErrorConservationEdge[l] = l2Norm(SumGje) / maxGjrAbs; } }); std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << "\n"; diff --git a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp index b05e33c3b8c5b7e8eb45dd13ef571aeda066c616..7386f4e84451e2d1b99d284df97eb5f9f44693d1 100644 --- a/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp +++ b/src/scheme/RusanovEulerianCompositeSolver_v2_o2.cpp @@ -965,6 +965,7 @@ class RusanovEulerianCompositeSolver_v2_o2 //_applyOutflowBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); _applySymmetricBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); _applyNeumannflatBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); + _applyWallBoundaryCondition(bc_list, *p_mesh, StateAtNode, StateAtEdge, StateAtFace); // // Construction du flux .. ok pour l'ordre 1 @@ -1236,13 +1237,15 @@ class RusanovEulerianCompositeSolver_v2_o2 if (not is_boundary_node[l]) { Rp SumGjr(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjr[K][R])); } // MaxErrorConservationNode = std::max(MaxErrorConservationNode, l2Norm(SumGjr)); - MaxErrorConservationNode[l] = l2Norm(SumGjr); + MaxErrorConservationNode[l] = l2Norm(SumGjr) / maxGjrAbs; } }); std::cout << " Max Error Node " << max(MaxErrorConservationNode) << "\n"; @@ -1309,12 +1312,14 @@ class RusanovEulerianCompositeSolver_v2_o2 if (not is_boundary_face[l]) { Rp SumGjf(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gjf[K][R])); } - MaxErrorConservationFace[l] = l2Norm(SumGjf); + MaxErrorConservationFace[l] = l2Norm(SumGjf) / maxGjrAbs; // MaxErrorConservationFace = std::max(MaxErrorConservationFace, l2Norm(SumGjf)); } }); @@ -1392,13 +1397,15 @@ class RusanovEulerianCompositeSolver_v2_o2 if (not is_boundary_edge[l]) { Rp SumGje(zero); + double maxGjrAbs = 0; 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]; + maxGjrAbs = std::max(maxGjrAbs, l2Norm(Gje[K][R])); } // MaxErrorConservationEdge = std::max(MaxErrorConservationEdge, l2Norm(SumGje)); - MaxErrorConservationEdge[l] = l2Norm(SumGje); + MaxErrorConservationEdge[l] = l2Norm(SumGje) / maxGjrAbs; } }); std::cout << " Max Error Edge " << max(MaxErrorConservationEdge) << "\n";