From b59181f737978785bffb472ee37439828ea940aa Mon Sep 17 00:00:00 2001
From: Clovis <clovis.schoeck@etudiant.univ-rennes.fr>
Date: Mon, 17 Jun 2024 11:52:40 +0200
Subject: [PATCH] =?UTF-8?q?Premier=20commit,=20cr=C3=A9ation=20des=20fichi?=
 =?UTF-8?q?ers=20li=C3=A9s=20au=20solveur=20GKS?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/language/modules/SchemeModule.cpp | 15 +++++
 src/scheme/CMakeLists.txt             |  4 +-
 src/scheme/GKS.cpp                    | 96 +++++++++++++++++++++++++++
 src/scheme/GKS.hpp                    | 16 +++++
 4 files changed, 130 insertions(+), 1 deletion(-)
 create mode 100644 src/scheme/GKS.cpp
 create mode 100644 src/scheme/GKS.hpp

diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 7b404fcf6..0be440945 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/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(
 
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index dadbe73bb..f3aa4dc1b 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -10,7 +10,9 @@ add_library(
   DiscreteFunctionVectorIntegrator.cpp
   DiscreteFunctionVectorInterpoler.cpp
   FluxingAdvectionSolver.cpp
-)
+  HyperelasticSolver.cpp
+  GKS.cpp
+  )
 
 target_link_libraries(
   PugsScheme
diff --git a/src/scheme/GKS.cpp b/src/scheme/GKS.cpp
new file mode 100644
index 000000000..09fa5f4ff
--- /dev/null
+++ b/src/scheme/GKS.cpp
@@ -0,0 +1,96 @@
+#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());
+}
diff --git a/src/scheme/GKS.hpp b/src/scheme/GKS.hpp
new file mode 100644
index 000000000..cd3643003
--- /dev/null
+++ b/src/scheme/GKS.hpp
@@ -0,0 +1,16 @@
+#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
-- 
GitLab