mise a jour Rusanov v1 et rajout Rusanov v2 correction Cjf en 2d (vu pour certains maillages)

src/language/modules/SchemeModule.cpp

@@ -44,6 +44,7 @@
 #include <scheme/NeumannBoundaryConditionDescriptor.hpp>
 #include <scheme/OutflowBoundaryConditionDescriptor.hpp>
 #include <scheme/RusanovEulerianCompositeSolver.hpp>
+#include <scheme/RusanovEulerianCompositeSolver_v2.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <scheme/VariableBCDescriptor.hpp>
 #include <utils/Socket.hpp>
@@ -507,7 +508,7 @@ SchemeModule::SchemeModule()
 
                               [](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>& E, const double& gamma,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& c,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
@@ -515,7 +516,7 @@ SchemeModule::SchemeModule()
                                  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);
+                                return rusanovEulerianCompositeSolver(rho, u, E, gamma, c, p, bc_descriptor_list, dt);
                               }
 
                               ));
@@ -525,7 +526,43 @@ SchemeModule::SchemeModule()
 
                               [](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>& E, const double& gamma,
+                                 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, gamma, c, p, bc_descriptor_list, dt,
+                                                                      true);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("rusanov_eulerian_composite_solver_version2",
+                            std::function(
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 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::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);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("rusanov_eulerian_composite_solver_version2_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 double& gamma,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& c,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& p,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
@@ -533,7 +570,8 @@ SchemeModule::SchemeModule()
                                  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);
+                                return rusanovEulerianCompositeSolver_v2(rho, u, E, gamma, c, p, bc_descriptor_list, dt,
+                                                                         true);
                               }
 
                               ));
@@ -797,7 +835,7 @@ SchemeModule::SchemeModule()
 
                                             [](const std::shared_ptr<const DiscreteFunctionVariant>& u,
                                                const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
-                                              return compute_dt(u, c);
+                                              return toolsCompositeSolver::compute_dt(u, c);
                                             }));
 
   this->_addBuiltinFunction("cell_volume",

src/mesh/MeshData.hpp

@@ -516,22 +516,7 @@ class MeshData<Mesh<Dimension>>
           Cjf(j, 1) = Rd{1};
         });
       m_Cjf = Cjf;
-    } else if constexpr (Dimension == 2) {
-      const NodeValue<const Rd>& xr   = m_mesh.xr();
-      const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
-
-      FaceValuePerCell<Rd> Cjf(m_mesh.connectivity());
-      parallel_for(
-        m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-          const auto& cell_nodes = cell_to_node_matrix[j];
-          for (size_t R = 0; R < cell_nodes.size(); ++R) {
-            int Rp           = (R + 1) % cell_nodes.size();
-            const Rd& xrp_xr = (xr[cell_nodes[Rp]] - xr[cell_nodes[R]]);
-            Cjf(j, R)        = Rd{xrp_xr[1], -xrp_xr[0]};
-          }
-        });
-      m_Cjf = Cjf;
-    } else if (Dimension == 3) {
+    } else if (Dimension >= 2) {
       auto Nl = this->Nl();
 
       const auto& cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
@@ -609,20 +594,54 @@ class MeshData<Mesh<Dimension>>
         });
       m_Cje = Cje;
     } else if constexpr (Dimension == 2) {
+      // EdgeValuePerCell<Rd> Cje(m_mesh.connectivity());
+      // Cje = this->Cjf();
+      // auto Cjf = this->Cjf();
+      // m_xl                            = FaceValue<const Rd>{m_mesh.connectivity(), m_mesh.xr().arrayView()};
+      m_Cje = EdgeValuePerCell<const Rd>{m_mesh.connectivity(), this->Cjf().arrayView()};
+      /*
       const NodeValue<const Rd>& xr   = m_mesh.xr();
+      const auto& cell_to_edge_matrix = m_mesh.connectivity().cellToEdgeMatrix();
+      const auto& edge_to_node_matrix = m_mesh.connectivity().edgeToNodeMatrix();
       const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
 
+      // const auto& cell_edge_is_reversed = m_mesh.connectivity().cellEgdeIsReversed();
+
       EdgeValuePerCell<Rd> Cje(m_mesh.connectivity());
       parallel_for(
         m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-          const auto& cell_nodes = cell_to_node_matrix[j];
-          for (size_t R = 0; R < cell_nodes.size(); ++R) {
-            int Rp           = (R + 1) % cell_nodes.size();
-            const Rd& xrp_xr = (xr[cell_nodes[Rp]] - xr[cell_nodes[R]]);
+          const auto& cell_edge = cell_to_edge_matrix[j];
+          const auto& cell_node = cell_to_node_matrix[j];
+          assert(cell_edge.size() == cell_node.size());
+          // size_t node_locl = 0;
+          for (size_t R = 0; R < cell_edge.size(); ++R) {
+            {
+              const EdgeId& edge_id   = cell_edge[R];
+              const auto& edges_nodes = edge_to_node_matrix[edge_id];
+
+              const NodeId& node0_id = edges_nodes[0];
+              const NodeId& node1_id = edges_nodes[1];
+              assert((node0_id == cell_node[R]) or (node0_id == cell_node[R]));
+              assert((node0_id == cell_node[(R + 1) % cell_node.size()]) or
+                     (node0_id == cell_node[(R + 1) % cell_node.size()]));
+
+              bool reversed = true;
+              if ((node0_id == cell_node[R]) and (node1_id == cell_node[(R + 1) % cell_node.size()]))
+                reversed = false;
+
+              const Rd& xrp_xr = (xr[node1_id] - xr[node0_id]);
+
+              if (reversed)
+                Cje(j, R) = -Rd{xrp_xr[1], -xrp_xr[0]};
+              else
                 Cje(j, R) = Rd{xrp_xr[1], -xrp_xr[0]};
+
+              // node_locl++;
+            }
           }
         });
       m_Cje = Cje;
+      */
     } else if (Dimension == 3) {
       auto Nlr = this->Nlr();
 

src/scheme/CMakeLists.txt

@@ -21,7 +21,9 @@ add_library(
   DiscreteFunctionVectorIntegrator.cpp
   DiscreteFunctionVectorInterpoler.cpp
   FluxingAdvectionSolver.cpp
+  RusanovEulerianCompositeSolverTools.cpp
   RusanovEulerianCompositeSolver.cpp
+  RusanovEulerianCompositeSolver_v2.cpp
 )
 
 target_link_libraries(

src/scheme/RusanovEulerianCompositeSolver.hpp

@@ -4,12 +4,13 @@
 #include <mesh/MeshVariant.hpp>
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
 
 #include <memory>
 #include <vector>
 
-double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
-                  const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+// 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>,
@@ -18,6 +19,7 @@ rusanovEulerianCompositeSolver(
   const std::shared_ptr<const DiscreteFunctionVariant>& rho,
   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::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,

src/scheme/RusanovEulerianCompositeSolverTools.cpp

+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+
+template <class Rd>
+double
+toolsCompositeSolver::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
+toolsCompositeSolver::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());
+}

src/scheme/RusanovEulerianCompositeSolverTools.hpp

+#ifndef RUSANOV_EULERIAN_COMPOSITE_SOLVER_TOOLS_HPP
+#define RUSANOV_EULERIAN_COMPOSITE_SOLVER_TOOLS_HPP
+
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+#include <memory>
+#include <vector>
+
+namespace toolsCompositeSolver
+{
+template <class Rd>
+double EvaluateMaxEigenValueTimesNormalLengthInGivenDirection(   // const double& rho_mean,
+  const Rd& U_mean,
+  // const double& E_mean,
+  const double& c_mean,
+  const Rd& normal);
+
+double compute_dt(const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+                  const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
+}   // namespace toolsCompositeSolver
+#endif

src/scheme/RusanovEulerianCompositeSolver_v2.hpp

+#ifndef RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+#define RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP
+
+#include <mesh/MeshVariant.hpp>
+#include <scheme/DiscreteFunctionVariant.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/RusanovEulerianCompositeSolverTools.hpp>
+
+#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>>
+rusanovEulerianCompositeSolver_v2(
+  const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+  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::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const double& dt,
+  const bool check = false);
+
+#endif   // RUSANOV_EULERIAN_COMPOSITE_SOLVER_V2_HPP