From 5eb27297943322a2a2b7c0fc6d7aeeb269fffab7 Mon Sep 17 00:00:00 2001
From: Clovis <clovis.schoeck@etudiant.univ-rennes.fr>
Date: Tue, 9 Jul 2024 11:04:07 +0200
Subject: [PATCH] =?UTF-8?q?Programme=20GKS=20=C3=A0=20l'ordre=201=20tourne?=
 =?UTF-8?q?=20et=20r=C3=A9sout=20les=20=C3=A9quations=20d'Euler?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/language/modules/SchemeModule.cpp |   5 +-
 src/scheme/GKS.cpp                    | 168 ++++++++++++++++++--------
 src/scheme/GKS.hpp                    |  11 +-
 3 files changed, 129 insertions(+), 55 deletions(-)

diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 76dfefd4a..4b9060790 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -426,11 +426,12 @@ SchemeModule::SchemeModule()
 
                               [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
                                  const std::shared_ptr<const DiscreteFunctionVariant>& rho_U,
-                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_E, const double gamma,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& tau, const double& delta,
                                  const double& dt) -> std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>,
                                                                  std::shared_ptr<const DiscreteFunctionVariant>> {
-                                return gks(rho, rho_U, rho_E, gamma, dt);
+                                return gks(rho, rho_U, rho_E, tau, delta, dt);
                               }
 
                               ));
diff --git a/src/scheme/GKS.cpp b/src/scheme/GKS.cpp
index c3f296b73..fed7b1ae0 100644
--- a/src/scheme/GKS.cpp
+++ b/src/scheme/GKS.cpp
@@ -17,15 +17,24 @@ class GKS
         std::shared_ptr<const DiscreteFunctionVariant> rho_v,
         std::shared_ptr<const DiscreteFunctionVariant> rho_U_v,
         std::shared_ptr<const DiscreteFunctionVariant> rho_E_v,
-        const double gamma,
+        std::shared_ptr<const DiscreteFunctionVariant> tau,
+        const double delta,
         double dt)
   {
     using Rd = TinyVector<MeshType::Dimension>;
 
     const MeshType& mesh = *p_mesh;
 
-    const double pi = acos(-1);
-    const int delta = (3 - gamma) / (gamma - 1);
+    const double pi = std::acos(-1);
+
+    DiscreteFunctionP0<const double> tau_n = tau->get<DiscreteFunctionP0<const double>>();
+    CellValue<double> eta(mesh.connectivity());
+    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+      if (tau_n[cell_id] == 0)
+        eta[cell_id] = 0;
+      else
+        eta[cell_id] = (tau_n[cell_id] / dt) * (1 - std::exp(-dt / tau_n[cell_id]));
+    }
 
     DiscreteFunctionP0<const double> rho_n   = rho_v->get<DiscreteFunctionP0<const double>>();
     DiscreteFunctionP0<const Rd> rho_U_n     = rho_U_v->get<DiscreteFunctionP0<const Rd>>();
@@ -36,101 +45,163 @@ class GKS
     DiscreteFunctionP0<double> rho_E = copy(rho_E_n);
 
     CellValue<double> lambda{p_mesh->connectivity()};
+    // lambda.fill(0);
     CellValue<Rd> U{p_mesh->connectivity()};
-
+    // U.fill(0);
     auto& mesh_data = MeshDataManager::instance().getMeshData(mesh);
     auto Vj         = mesh_data.Vj();
 
     auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
-    // const NodeId first_node_id = NodeId{0};
-    // const NodeId last_node_id  = NodeId{static_cast<typename CellId::base_type>(mesh.numberOfNodes() - 1)};
 
-    NodeValue<double> rho_flux(mesh.connectivity());
-    NodeValue<Rd> rho_U_flux(mesh.connectivity());
-    NodeValue<double> rho_E_flux(mesh.connectivity());
+    NodeValue<double> rho_flux_Euler(mesh.connectivity());
+    NodeValue<Rd> rho_U_flux_Euler(mesh.connectivity());
+    NodeValue<double> rho_E_flux_Euler(mesh.connectivity());
+    rho_flux_Euler.fill(0);
+    rho_U_flux_Euler.fill(TinyVector<1>(0));
+    rho_E_flux_Euler.fill(0);
+
+    NodeValue<double> rho_flux_Navier(mesh.connectivity());
+    NodeValue<Rd> rho_U_flux_Navier(mesh.connectivity());
+    NodeValue<double> rho_E_flux_Navier(mesh.connectivity());
+    rho_flux_Navier.fill(0);
+    rho_U_flux_Navier.fill(TinyVector<1>(0));
+    rho_E_flux_Navier.fill(0);
 
     NodeValue<double> rho_node(mesh.connectivity());
     NodeValue<Rd> rho_U_node(mesh.connectivity());
     NodeValue<double> rho_E_node(mesh.connectivity());
+    rho_node.fill(0);
+    rho_U_node.fill(TinyVector<1>(0));
+    rho_E_node.fill(0);
 
     for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
       U[cell_id][0]   = rho_U_n[cell_id][0] / rho_n[cell_id];
-      double U_2      = U[cell_id][0] * U[cell_id][0];
       double rho_U_2  = rho_U_n[cell_id][0] * U[cell_id][0];
-      lambda[cell_id] = (1. + delta) * (rho_n[cell_id]) / (4 * rho_E_n[cell_id] - 2 * rho_U_2);
+      lambda[cell_id] = 0.5 * (1. + delta) * rho_n[cell_id] / (2 * rho_E_n[cell_id] - rho_U_2);
+    }
+
+    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+      double U_2 = U[cell_id][0] * U[cell_id][0];
 
       double rho_cell_left = rho_n[cell_id] * (1 + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0]));
-      rho_cell_left /= 2.;
 
       Rd rho_U_cell_left;
       rho_U_cell_left[0] = rho_U_n[cell_id][0] * (1. + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) +
                            rho_n[cell_id] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
-      rho_U_cell_left[0] /= 2.;
 
       double rho_E_cell_left =
         rho_E_n[cell_id] * (1. + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) +
         0.5 * rho_U_n[cell_id][0] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
-      rho_E_cell_left /= 2.;
 
       auto node_list = cell_to_node_matrix[cell_id];
 
-      rho_node[node_list[1]]   = rho_cell_left;
-      rho_U_node[node_list[1]] = rho_U_cell_left;
-      rho_E_node[node_list[1]] = rho_E_cell_left;
+      rho_node[node_list[1]]      = 0.5 * rho_cell_left;
+      rho_U_node[node_list[1]][0] = 0.5 * rho_U_cell_left[0];
+      rho_E_node[node_list[1]]    = 0.5 * rho_E_cell_left;
     }
 
-    for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+    for (CellId cell_id = 1; cell_id < mesh.numberOfCells(); ++cell_id) {
       double U_2 = U[cell_id][0] * U[cell_id][0];
 
-      double rho_cell_right = rho_n[cell_id] * (1 - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0]));
-      rho_cell_right /= 2.;
+      double rho_cell_right = rho_n[cell_id] * (1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0]));
 
       Rd rho_U_cell_right;
       rho_U_cell_right[0] = rho_U_n[cell_id][0] * (1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) -
                             rho_n[cell_id] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
-      rho_U_cell_right[0] /= 2.;
 
       double rho_E_cell_right =
         rho_E_n[cell_id] * (1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) -
         0.5 * rho_U_n[cell_id][0] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
-      rho_E_cell_right /= 2.;
 
       auto node_list = cell_to_node_matrix[cell_id];
 
-      rho_node[node_list[0]] += rho_cell_right;
-      rho_U_node[node_list[0]] += rho_U_cell_right;
-      rho_E_node[node_list[0]] += rho_E_cell_right;
+      rho_node[node_list[0]] += 0.5 * rho_cell_right;
+      rho_U_node[node_list[0]][0] += 0.5 * rho_U_cell_right[0];
+      rho_E_node[node_list[0]] += 0.5 * rho_E_cell_right;
     }
 
+    for (CellId cell_id = 1; cell_id < mesh.numberOfCells(); ++cell_id) {
+      auto node_list      = cell_to_node_matrix[cell_id];
+      double rho_U_2_node = rho_U_node[node_list[0]][0] * rho_U_node[node_list[0]][0] / rho_node[node_list[0]];
+
+      rho_flux_Euler[node_list[0]] = rho_U_node[node_list[0]][0];
+      rho_U_flux_Euler[node_list[0]][0] =
+        delta * rho_U_2_node / (1. + delta) + 2 * rho_E_node[node_list[0]] / (1. + delta);
+      rho_E_flux_Euler[node_list[0]] =
+        rho_U_node[node_list[0]][0] / rho_node[node_list[0]] *
+        ((3. + delta) * rho_E_node[node_list[0]] / (1. + delta) - rho_U_2_node / (1. + delta));
+    }
+
+    //##
+    //%%%%%%%%%%%%%%%
+    //##
+
     for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+      double U_2     = U[cell_id][0] * U[cell_id][0];
+      double rho_U_2 = rho_U_n[cell_id][0] * U[cell_id][0];
+
+      Rd F2_fn_left;
+      F2_fn_left[0] = (rho_U_2 + 0.5 * rho_n[cell_id] / lambda[cell_id]) *
+                        (1. + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) +
+                      rho_U_n[cell_id][0] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
+
+      double F3_fn_left = 0.5 * rho_U_n[cell_id][0] * (U_2 + 0.5 * (delta + 3) / lambda[cell_id]) *
+                            (1. + std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) +
+                          0.5 * rho_n[cell_id] * (U_2 + 0.5 * (delta + 2) / lambda[cell_id]) *
+                            std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
+
       auto node_list = cell_to_node_matrix[cell_id];
 
-      rho_flux[node_list[1]] = rho_U_node[node_list[1]][0];
-      rho_U_flux[node_list[1]][0] =
-        delta / (1. + delta) * rho_U_node[node_list[1]][0] * rho_U_node[node_list[1]][0] / rho_node[node_list[1]] +
-        2. / (1. + delta) * rho_E_node[node_list[1]];
-      rho_E_flux[node_list[1]] =
-        rho_U_node[node_list[1]][0] / rho_node[node_list[1]] *
-        ((3. + delta) / (1. + delta) * rho_E_node[node_list[1]] -
-         1. / (1. + delta) * rho_U_node[node_list[1]][0] * rho_U_node[node_list[1]][0] / rho_node[node_list[1]]);
+      rho_U_flux_Navier[node_list[1]][0] += 0.5 * F2_fn_left[0];
+      rho_E_flux_Navier[node_list[1]] += 0.5 * F3_fn_left;
     }
 
-    // Peut être un problème ici avec les bords mais je sais pas encore
-    //    rho_flux[first_node_id]   = rho_flux[last_node_id];
-    //    rho_U_flux[first_node_id] = rho_U_flux[last_node_id];
-    //    rho_E_flux[first_node_id] = rho_E_flux[last_node_id];
+    for (CellId cell_id = 1; cell_id < mesh.numberOfCells(); ++cell_id) {
+      double U_2     = U[cell_id][0] * U[cell_id][0];
+      double rho_U_2 = rho_U_n[cell_id][0] * U[cell_id][0];
 
-    for (CellId cell_id = 1; cell_id < mesh.numberOfCells() - 1; ++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]]);
-      const Rd rho_U_flux_sum     = (rho_U_flux[node_list[1]] - rho_U_flux[node_list[0]]);
-      const double rho_E_flux_sum = (rho_E_flux[node_list[1]] - rho_E_flux[node_list[0]]);
-
-      rho[cell_id] -= dt * rho_flux_sum / Vj[cell_id];
-      rho_U[cell_id][0] -= dt * rho_U_flux_sum[0] / Vj[cell_id];
-      rho_E[cell_id] -= dt * rho_E_flux_sum / Vj[cell_id];
+      Rd F2_fn_right;
+      F2_fn_right[0] = (rho_U_2 + 0.5 * rho_n[cell_id] / lambda[cell_id]) *
+                         (1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) -
+                       rho_U_n[cell_id][0] * std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
+
+      double F3_fn_right = 0.5 * rho_U_n[cell_id][0] * (U_2 + 0.5 * (delta + 3) / lambda[cell_id]) *
+                             (1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) -
+                           0.5 * rho_n[cell_id] * (U_2 + 0.5 * (delta + 2) / lambda[cell_id]) *
+                             std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
+
+      // double F3_fn_right = U[cell_id][0] * (rho_E_n[cell_id] + 0.5 * rho_n[cell_id] / lambda[cell_id]) *
+      //                        (1. - std::erf(std::sqrt(lambda[cell_id]) * U[cell_id][0])) -
+      //                      0.5 * (rho_E_n[cell_id] + 0.25 * rho_n[cell_id] / lambda[cell_id]) *
+      //                        std::exp(-lambda[cell_id] * U_2) / std::sqrt(pi * lambda[cell_id]);
+
+      auto node_list = cell_to_node_matrix[cell_id];
+
+      rho_U_flux_Navier[node_list[0]][0] += 0.5 * F2_fn_right[0];
+      rho_E_flux_Navier[node_list[0]] += 0.5 * F3_fn_right;
     }
+    // std::cout << "lambda " << lambda << std::endl;
+    // std::cout << "rho flux " << rho_flux_Euler << std::endl;
+    // std::cout << "rhoU flux " << rho_U_flux_Euler << std::endl;
+    // std::cout << "rhoE flux " << rho_E_flux_Euler << std::endl;
+    // std::exit(0);
+    for (CellId cell_id = 1; cell_id < mesh.numberOfCells() - 1; ++cell_id) {
+      auto node_list = cell_to_node_matrix[cell_id];
 
+      const double rho_flux_Euler_sum   = (rho_flux_Euler[node_list[1]] - rho_flux_Euler[node_list[0]]);
+      const Rd rho_U_flux_Euler_sum     = (rho_U_flux_Euler[node_list[1]] - rho_U_flux_Euler[node_list[0]]);
+      const double rho_E_flux_Euler_sum = (rho_E_flux_Euler[node_list[1]] - rho_E_flux_Euler[node_list[0]]);
+
+      const Rd rho_U_flux_Navier_sum     = (rho_U_flux_Navier[node_list[1]] - rho_U_flux_Navier[node_list[0]]);
+      const double rho_E_flux_Navier_sum = (rho_E_flux_Navier[node_list[1]] - rho_E_flux_Navier[node_list[0]]);
+
+      rho[cell_id] -= dt / Vj[cell_id] * (rho_flux_Euler_sum);
+      rho_U[cell_id][0] -=
+        dt / Vj[cell_id] *
+        (rho_U_flux_Euler_sum[0] + eta[cell_id] * (rho_U_flux_Navier_sum[0] - rho_U_flux_Euler_sum[0]));
+      rho_E[cell_id] -=
+        dt / Vj[cell_id] * (rho_E_flux_Euler_sum + 0 * eta[cell_id] * (rho_E_flux_Navier_sum - rho_E_flux_Euler_sum));
+    }
     return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(rho),
                            std::make_shared<DiscreteFunctionVariant>(rho_U),
                            std::make_shared<DiscreteFunctionVariant>(rho_E));
@@ -145,7 +216,8 @@ std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,   // rho
 gks(std::shared_ptr<const DiscreteFunctionVariant> rho_v,
     std::shared_ptr<const DiscreteFunctionVariant> rho_U_v,
     std::shared_ptr<const DiscreteFunctionVariant> rho_E_v,
-    const double gamma,
+    std::shared_ptr<const DiscreteFunctionVariant> tau,
+    const double delta,
     const double dt)
 {
   std::shared_ptr mesh_v = getCommonMesh({rho_v, rho_U_v, rho_E_v});
@@ -165,7 +237,7 @@ gks(std::shared_ptr<const DiscreteFunctionVariant> rho_v,
       if constexpr (is_polygonal_mesh_v<MeshType>) {
         if constexpr (MeshType::Dimension == 1) {
           GKS<MeshType> gks;
-          return gks.solve(p_mesh, rho_v, rho_U_v, rho_E_v, gamma, dt);
+          return gks.solve(p_mesh, rho_v, rho_U_v, rho_E_v, tau, delta, dt);
 
         } else {
           throw NormalError("dimension not treated");
diff --git a/src/scheme/GKS.hpp b/src/scheme/GKS.hpp
index e29ad2aaa..4bd32fc47 100644
--- a/src/scheme/GKS.hpp
+++ b/src/scheme/GKS.hpp
@@ -6,12 +6,13 @@
 #include <tuple>
 
 std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,   // rho
-           std::shared_ptr<const DiscreteFunctionVariant>,   // U
-           std::shared_ptr<const DiscreteFunctionVariant>>   // E
+           std::shared_ptr<const DiscreteFunctionVariant>,   // rhoU
+           std::shared_ptr<const DiscreteFunctionVariant>>   // rhoE
 gks(std::shared_ptr<const DiscreteFunctionVariant> rho,
-    std::shared_ptr<const DiscreteFunctionVariant> U,
-    std::shared_ptr<const DiscreteFunctionVariant> E,
-    const double gamma,
+    std::shared_ptr<const DiscreteFunctionVariant> rhoU,
+    std::shared_ptr<const DiscreteFunctionVariant> rhoE,
+    std::shared_ptr<const DiscreteFunctionVariant> tau,
+    const double delta,
     const double dt);
 
 #endif   // GKS_HPP
-- 
GitLab