Schéma GKS écrit

Erreur oscillation fin de l'onde de détente
Déplacement de l'onde de détente, choc et discontinuité de contact lrosque dt est modifié

src/language/modules/SchemeModule.cpp

@@ -425,12 +425,12 @@ SchemeModule::SchemeModule()
                             std::function(
 
                               [](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>& rho_U,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho_E, const double gamma,
                                  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);
+                                return gks(rho, rho_U, rho_E, gamma, dt);
                               }
 
                               ));

src/scheme/GKS.cpp

@@ -15,43 +15,125 @@ class GKS
              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,
+        std::shared_ptr<const DiscreteFunctionVariant> rho_U_v,
+        std::shared_ptr<const DiscreteFunctionVariant> rho_E_v,
+        const double gamma,
         double dt)
   {
     using Rd = TinyVector<MeshType::Dimension>;
 
     const MeshType& mesh = *p_mesh;
 
+    const double pi = acos(-1);
+    const int delta = (3 - gamma) / (gamma - 1);
+
     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<const Rd> rho_U_n     = rho_U_v->get<DiscreteFunctionP0<const Rd>>();
+    DiscreteFunctionP0<const double> rho_E_n = rho_E_v->get<DiscreteFunctionP0<const double>>();
 
     DiscreteFunctionP0<double> rho   = copy(rho_n);
-    DiscreteFunctionP0<Rd> U       = copy(U_n);
-    DiscreteFunctionP0<double> E   = copy(E_n);
+    DiscreteFunctionP0<Rd> rho_U     = copy(rho_U_n);
+    DiscreteFunctionP0<double> rho_E = copy(rho_E_n);
+
+    CellValue<double> lambda{p_mesh->connectivity()};
+    CellValue<Rd> U{p_mesh->connectivity()};
 
     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_node(mesh.connectivity());
+    NodeValue<Rd> rho_U_node(mesh.connectivity());
+    NodeValue<double> rho_E_node(mesh.connectivity());
+
     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);
+
+      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_flux[node_list[1]] = rho_n[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_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) {
+      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.;
+
+      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;
+    }
+
+    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_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]]);
+    }
+
+    // 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() - 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]]) / Vj[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 * 0.5 * rho_flux_sum;
+      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];
     }
 
-    return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(rho), std::make_shared<DiscreteFunctionVariant>(U),
-                           std::make_shared<DiscreteFunctionVariant>(E));
+    return std::make_tuple(std::make_shared<DiscreteFunctionVariant>(rho),
+                           std::make_shared<DiscreteFunctionVariant>(rho_U),
+                           std::make_shared<DiscreteFunctionVariant>(rho_E));
   }
 
   GKS() = default;
@@ -61,16 +143,17 @@ 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,
+    std::shared_ptr<const DiscreteFunctionVariant> rho_U_v,
+    std::shared_ptr<const DiscreteFunctionVariant> rho_E_v,
+    const double gamma,
     const double dt)
 {
-  std::shared_ptr mesh_v = getCommonMesh({rho_v, U_v, E_v});
+  std::shared_ptr mesh_v = getCommonMesh({rho_v, rho_U_v, rho_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)) {
+  if (not checkDiscretizationType({rho_v, rho_U_v, rho_E_v}, DiscreteFunctionType::P0)) {
     throw NormalError("GKS solver expects P0 functions");
   }
 
@@ -82,7 +165,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, U_v, E_v, dt);
+          return gks.solve(p_mesh, rho_v, rho_U_v, rho_E_v, gamma, dt);
 
         } else {
           throw NormalError("dimension not treated");

src/scheme/GKS.hpp

@@ -11,6 +11,7 @@ std::tuple<std::shared_ptr<const DiscreteFunctionVariant>,   // rho
 gks(std::shared_ptr<const DiscreteFunctionVariant> rho,
     std::shared_ptr<const DiscreteFunctionVariant> U,
     std::shared_ptr<const DiscreteFunctionVariant> E,
+    const double gamma,
     const double dt);
 
 #endif   // GKS_HPP