From c71c4942bacb8a076b4f627e3128068f447f3b42 Mon Sep 17 00:00:00 2001
From: Alexandre Gangloff <alexandre.gangloff@cea.fr>
Date: Wed, 11 Jun 2025 11:00:02 +0200
Subject: [PATCH] Add dt computation for third and fourth orders

---
 src/language/modules/LocalFSIModule.cpp | 16 ++++++++
 src/scheme/HighOrderAcousticSolver.cpp  | 52 +++++++++++++++++++++++++
 src/scheme/HighOrderAcousticSolver.hpp  |  4 ++
 3 files changed, 72 insertions(+)

diff --git a/src/language/modules/LocalFSIModule.cpp b/src/language/modules/LocalFSIModule.cpp
index 878a93faa..4c46e950f 100644
--- a/src/language/modules/LocalFSIModule.cpp
+++ b/src/language/modules/LocalFSIModule.cpp
@@ -1077,6 +1077,22 @@ LocalFSIModule::LocalFSIModule()
                               }
 
                               ));
+                            
+  this->_addBuiltinFunction("acoustic_dt_order3", std::function(
+
+                                             [](const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
+                                               return acoustic_dt_order3(c);
+                                             }
+
+                                             ));
+
+  this->_addBuiltinFunction("acoustic_dt_order4", std::function(
+
+                                             [](const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
+                                               return acoustic_dt_order4(c);
+                                             }
+
+                                             ));
 }
 
 void
diff --git a/src/scheme/HighOrderAcousticSolver.cpp b/src/scheme/HighOrderAcousticSolver.cpp
index 8db171654..dc7c70a5c 100644
--- a/src/scheme/HighOrderAcousticSolver.cpp
+++ b/src/scheme/HighOrderAcousticSolver.cpp
@@ -30,6 +30,58 @@
 #include <variant>
 #include <vector>
 
+double 
+acoustic_dt_order3(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 = decltype(mesh);
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+        const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+        const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
+
+        CellValue<double> local_dt{mesh.connectivity()};
+        parallel_for(
+          mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * pow(Vj[j],3/2) / (Sj[j] * c[j]); });
+
+        return min(local_dt);
+      } else {
+        throw NormalError("unexpected mesh type");
+      }
+    },
+    c.meshVariant()->variant());
+}
+
+double 
+acoustic_dt_order4(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 = decltype(mesh);
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+        const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+        const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
+
+        CellValue<double> local_dt{mesh.connectivity()};
+        parallel_for(
+          mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] * Vj[j] / (Sj[j] * c[j]); });
+
+        return min(local_dt);
+      } else {
+        throw NormalError("unexpected mesh type");
+      }
+    },
+    c.meshVariant()->variant());
+}
+
 template <MeshConcept MeshType>
 class HighOrderAcousticSolverHandler::HighOrderAcousticSolver final : public HighOrderAcousticSolverHandler::IHighOrderAcousticSolver
 {
diff --git a/src/scheme/HighOrderAcousticSolver.hpp b/src/scheme/HighOrderAcousticSolver.hpp
index 8545750cc..a85522cb1 100644
--- a/src/scheme/HighOrderAcousticSolver.hpp
+++ b/src/scheme/HighOrderAcousticSolver.hpp
@@ -13,6 +13,10 @@ class MeshVariant;
 class ItemValueVariant;
 class SubItemValuePerItemVariant;
 
+double acoustic_dt_order3(const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
+double acoustic_dt_order4(const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
 class HighOrderAcousticSolverHandler
 {
  public:
-- 
GitLab