From 45243ef22eae1552e93585775245e6ea6aea9809 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 22 Dec 2023 16:39:05 +0100
Subject: [PATCH] Add a discrete function that stores local acoustic time step

---
 src/language/modules/SchemeModule.cpp |  8 ++++++
 src/scheme/ImplicitAcousticSolver.cpp | 41 +++++++++++++++++++++++++++
 src/scheme/ImplicitAcousticSolver.hpp |  3 ++
 3 files changed, 52 insertions(+)

diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index b950437fa..d308ba9e6 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -812,6 +812,14 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("local_acoustic_dt",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& c)
+                                -> std::shared_ptr<const DiscreteFunctionVariant> { return local_acoustic_dt(c); }
+
+                              ));
+
   this->_addBuiltinFunction("linear_euler",
                             std::function(
 
diff --git a/src/scheme/ImplicitAcousticSolver.cpp b/src/scheme/ImplicitAcousticSolver.cpp
index df37a7591..f823c1b8f 100644
--- a/src/scheme/ImplicitAcousticSolver.cpp
+++ b/src/scheme/ImplicitAcousticSolver.cpp
@@ -38,6 +38,47 @@ Timer HJ_A_t;
 int count_getA     = 0;
 int count_getGradF = 0;
 
+template <size_t Dimension>
+std::shared_ptr<const DiscreteFunctionVariant>
+local_acoustic_dt(const DiscreteFunctionP0<Dimension, const double>& c)
+{
+  const Mesh<Connectivity<Dimension>>& mesh = dynamic_cast<const Mesh<Connectivity<Dimension>>&>(*c.mesh());
+  std::shared_ptr<const Mesh<Connectivity<Dimension>>> p_mesh =
+    std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(c.mesh());
+
+  const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+  const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
+
+  CellValue<double> local_dt{mesh.connectivity()};
+  local_dt.fill(std::numeric_limits<double>::max());
+  parallel_for(
+    mesh.numberOfCells(),
+    PUGS_LAMBDA(const CellId cell_id) { local_dt[cell_id] = 2 * Vj[cell_id] / (Sj[cell_id] * c[cell_id]); });
+
+  return std::make_shared<DiscreteFunctionVariant>(DiscreteFunctionP0<Dimension, const double>(p_mesh, local_dt));
+}
+
+std::shared_ptr<const DiscreteFunctionVariant>
+local_acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c)
+{
+  std::shared_ptr mesh = getCommonMesh({c});
+
+  switch (mesh->dimension()) {
+  case 1: {
+    return local_acoustic_dt(c->get<DiscreteFunctionP0<1, const double>>());
+  }
+  case 2: {
+    return local_acoustic_dt(c->get<DiscreteFunctionP0<2, const double>>());
+  }
+  case 3: {
+    return local_acoustic_dt(c->get<DiscreteFunctionP0<3, const double>>());
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
 template <size_t Dimension>
 double
 acoustic_dt(const DiscreteFunctionP0<Dimension, const double>& c,
diff --git a/src/scheme/ImplicitAcousticSolver.hpp b/src/scheme/ImplicitAcousticSolver.hpp
index c3c7eb4fe..4dcbaa6e5 100644
--- a/src/scheme/ImplicitAcousticSolver.hpp
+++ b/src/scheme/ImplicitAcousticSolver.hpp
@@ -30,6 +30,9 @@ class IZoneDescriptor;
 double acoustic_dt(const std::shared_ptr<const DiscreteFunctionVariant>& c,
                    const std::vector<std::shared_ptr<const IZoneDescriptor>>& explicit_zone_list);
 
+std::shared_ptr<const DiscreteFunctionVariant> local_acoustic_dt(
+  const std::shared_ptr<const DiscreteFunctionVariant>& c);
+
 class ImplicitAcousticSolverHandler
 {
  public:
-- 
GitLab