Premier commit, création des fichiers liés au solveur GKS

src/language/modules/SchemeModule.cpp

@@ -35,6 +35,7 @@
 #include <scheme/FluxingAdvectionSolver.hpp>
 #include <scheme/FourierBoundaryConditionDescriptor.hpp>
 #include <scheme/FreeBoundaryConditionDescriptor.hpp>
+#include <scheme/GKS.hpp>
 #include <scheme/HyperelasticSolver.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
@@ -420,6 +421,20 @@ SchemeModule::SchemeModule()
 
                                               ));
 
+  this->_addBuiltinFunction("gks",
+                            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& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>,
+                                                                 std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return gks(rho, U, E, dt);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("glace_solver",
                             std::function(
 

src/scheme/CMakeLists.txt

@@ -10,6 +10,8 @@ add_library(
   DiscreteFunctionVectorIntegrator.cpp
   DiscreteFunctionVectorInterpoler.cpp
   FluxingAdvectionSolver.cpp
+  HyperelasticSolver.cpp
+  GKS.cpp
   )
 
 target_link_libraries(

src/scheme/GKS.cpp

+#include <scheme/GKS.hpp>
+
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshTraits.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+
+template <MeshConcept MeshType>
+class GKS
+{
+ public:
+  std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  solve(std::shared_ptr<const MeshType> p_mesh,
+        std::shared_ptr<const DiscreteFunctionVariant> rho_v,
+        std::shared_ptr<const DiscreteFunctionVariant> U_v,
+        std::shared_ptr<const DiscreteFunctionVariant> E_v,
+        double dt)
+  {
+    using Rd = TinyVector<MeshType::Dimension>;
+
+    const MeshType& mesh = *p_mesh;
+
+    DiscreteFunctionP0<const double> rho_n = rho_v->get<DiscreteFunctionP0<const double>>();
+    DiscreteFunctionP0<const Rd> U_n       = U_v->get<DiscreteFunctionP0<const Rd>>();
+    DiscreteFunctionP0<const double> E_n   = E_v->get<DiscreteFunctionP0<const double>>();
+
+    DiscreteFunctionP0<double> rho = copy(rho_n);
+    DiscreteFunctionP0<Rd> U       = copy(U_n);
+    DiscreteFunctionP0<double> E   = copy(E_n);
+
+    auto& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+    auto Vj         = mesh_data.Vj();
+
+    auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    NodeValue<double> rho_flux(mesh.connectivity());
+    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+      auto node_list         = cell_to_node_matrix[cell_id];
+      rho_flux[node_list[1]] = rho_n[cell_id];
+    }
+    rho_flux[NodeId{0}] = rho_flux[NodeId{static_cast<typename CellId::base_type>(mesh.numberOfNodes() - 1)}];
+
+    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+      auto node_list            = cell_to_node_matrix[cell_id];
+      const double rho_flux_sum = (rho_flux[node_list[1]] - rho_flux[node_list[0]]) / Vj[cell_id];
+
+      rho[cell_id] -= dt * 0.5 * rho_flux_sum;
+    }
+
+    return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(rho), std::make_shared<DiscreteFunctionVariant>(U),
+                           std::make_shared<DiscreteFunctionVariant>(E));
+  }
+
+  GKS() = default;
+};
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,   // rho
+           std::shared_ptr<const DiscreteFunctionVariant>,   // U
+           std::shared_ptr<const DiscreteFunctionVariant>>   // E
+gks(std::shared_ptr<const DiscreteFunctionVariant> rho_v,
+    std::shared_ptr<const DiscreteFunctionVariant> U_v,
+    std::shared_ptr<const DiscreteFunctionVariant> E_v,
+    const double dt)
+{
+  std::shared_ptr mesh_v = getCommonMesh({rho_v, U_v, E_v});
+  if (not mesh_v) {
+    throw NormalError("discrete functions are not defined on the same mesh");
+  }
+
+  if (not checkDiscretizationType({rho_v, U_v, E_v}, DiscreteFunctionType::P0)) {
+    throw NormalError("GKS solver expects P0 functions");
+  }
+
+  return std::visit(
+    [&](auto&& p_mesh)
+      -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                    std::shared_ptr<const DiscreteFunctionVariant>> {
+      using MeshType = std::decay_t<decltype(*p_mesh)>;
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+        if constexpr (MeshType::Dimension == 1) {
+          GKS<MeshType> gks;
+          return gks.solve(p_mesh, rho_v, U_v, E_v, dt);
+
+        } else {
+          throw NormalError("dimension not treated");
+        }
+
+      } else {
+        throw NormalError("unexpected mesh type");
+      }
+    },
+    mesh_v->variant());
+}

src/scheme/GKS.hpp

+#ifndef GKS_HPP
+#define GKS_HPP
+
+#include <scheme/DiscreteFunctionVariant.hpp>
+
+#include <tuple>
+
+std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,   // rho
+           std::shared_ptr<const DiscreteFunctionVariant>,   // U
+           std::shared_ptr<const DiscreteFunctionVariant>>   // E
+gks(std::shared_ptr<const DiscreteFunctionVariant> rho,
+    std::shared_ptr<const DiscreteFunctionVariant> U,
+    std::shared_ptr<const DiscreteFunctionVariant> E,
+    const double dt);
+
+#endif   // GKS_HPP