From 74b62f638e6de24a22b8d49ce8155565c559d03a Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Sun, 21 Aug 2022 09:23:12 +0200
Subject: [PATCH] Continue implementation of entropic time step  [ci skip]

---
 src/analysis/Polynomial1D.hpp         |  60 +++--
 src/language/modules/SchemeModule.cpp |  32 ++-
 src/scheme/AcousticSolver.cpp         | 310 +++++++++++++++++++++++---
 src/scheme/AcousticSolver.hpp         |  18 +-
 4 files changed, 368 insertions(+), 52 deletions(-)

diff --git a/src/analysis/Polynomial1D.hpp b/src/analysis/Polynomial1D.hpp
index 1e0ba6ead..1ffa4d96c 100644
--- a/src/analysis/Polynomial1D.hpp
+++ b/src/analysis/Polynomial1D.hpp
@@ -12,17 +12,25 @@ class [[nodiscard]] Polynomial1D
   std::vector<double> m_coefficients;
 
   PUGS_INLINE
-  size_t _getRealDegree() const
+  size_t
+  _getRealDegree() const
   {
+    double abs_max_coef = 0;
+    for (size_t i = 0; i < m_coefficients.size(); ++i) {
+      if (std::abs(m_coefficients[i]) > abs_max_coef) {
+        abs_max_coef = std::abs(m_coefficients[i]);
+      }
+    }
     size_t real_degree = this->degree();
-    while (real_degree > 0 and std::abs(m_coefficients[real_degree]) < 1E-14) {
+    while (real_degree > 0 and std::abs(m_coefficients[real_degree]) < 1E-15 * abs_max_coef) {
       --real_degree;
     }
     return real_degree;
   }
 
   PUGS_INLINE
-  friend Polynomial1D _simplify(Polynomial1D && p)
+  friend Polynomial1D
+  _simplify(Polynomial1D&& p)
   {
     // return std::move(p);
     size_t real_degree = p._getRealDegree();
@@ -39,7 +47,8 @@ class [[nodiscard]] Polynomial1D
   }
 
  public:
-  friend std::ostream& operator<<(std::ostream& os, const Polynomial1D& p)
+  friend std::ostream&
+  operator<<(std::ostream& os, const Polynomial1D& p)
   {
     bool has_written = false;
     for (size_t i = 0; i < p.m_coefficients.size(); ++i) {
@@ -73,7 +82,8 @@ class [[nodiscard]] Polynomial1D
   }
 
   PUGS_INLINE
-  Polynomial1D operator-() const
+  Polynomial1D
+  operator-() const
   {
     Polynomial1D opposite(this->degree());
     for (size_t i = 0; i < m_coefficients.size(); ++i) {
@@ -83,7 +93,8 @@ class [[nodiscard]] Polynomial1D
   }
 
   PUGS_INLINE
-  friend Polynomial1D operator*(const double a, const Polynomial1D& p)
+  friend Polynomial1D
+  operator*(const double a, const Polynomial1D& p)
   {
     Polynomial1D product(p.degree());
     for (size_t i = 0; i < p.m_coefficients.size(); ++i) {
@@ -93,7 +104,8 @@ class [[nodiscard]] Polynomial1D
   }
 
   PUGS_INLINE
-  Polynomial1D operator*(const Polynomial1D& p) const
+  Polynomial1D
+  operator*(const Polynomial1D& p) const
   {
     Polynomial1D product(this->degree() + p.degree());
     std::fill(begin(product.m_coefficients), end(product.m_coefficients), 0);
@@ -106,7 +118,8 @@ class [[nodiscard]] Polynomial1D
   }
 
   PUGS_INLINE
-  Polynomial1D operator+(const Polynomial1D& p) const
+  Polynomial1D
+  operator+(const Polynomial1D& p) const
   {
     auto sum_left_is_bigger = [](const Polynomial1D& greater, const Polynomial1D& smaller) {
       Polynomial1D result(greater.degree());
@@ -128,7 +141,8 @@ class [[nodiscard]] Polynomial1D
   }
 
   PUGS_INLINE
-  Polynomial1D operator-(const Polynomial1D& p) const
+  Polynomial1D
+  operator-(const Polynomial1D& p) const
   {
     if (m_coefficients.size() >= p.m_coefficients.size()) {
       Polynomial1D result(this->degree());
@@ -152,7 +166,8 @@ class [[nodiscard]] Polynomial1D
   }
 
   PUGS_INLINE
-  Polynomial1D operator%(const Polynomial1D& q) const
+  Polynomial1D
+  operator%(const Polynomial1D& q) const
   {
     const Polynomial1D& p = *this;
     Polynomial1D ratio(this->degree());
@@ -186,7 +201,8 @@ class [[nodiscard]] Polynomial1D
   }
 
   PUGS_INLINE
-  Polynomial1D operator/(const Polynomial1D& q) const
+  Polynomial1D
+  operator/(const Polynomial1D& q) const
   {
     const Polynomial1D& p = *this;
     Polynomial1D ratio(this->degree());
@@ -209,28 +225,32 @@ class [[nodiscard]] Polynomial1D
     return _simplify(std::move(ratio));
   }
 
-  PUGS_INLINE size_t degree() const
+  PUGS_INLINE size_t
+  degree() const
   {
     Assert(m_coefficients.size() > 0);
     return m_coefficients.size() - 1;
   }
 
   PUGS_INLINE
-  double& coefficient(const size_t i)
+  double&
+  coefficient(const size_t i)
   {
     Assert(i < m_coefficients.size(), "invalid coefficient number");
     return m_coefficients[i];
   }
 
   PUGS_INLINE
-  const double& coefficient(const size_t i) const
+  const double&
+  coefficient(const size_t i) const
   {
     Assert(i < m_coefficients.size(), "invalid coefficient number");
     return m_coefficients[i];
   }
 
   PUGS_INLINE
-  friend Polynomial1D derive(const Polynomial1D& p)
+  friend Polynomial1D
+  derive(const Polynomial1D& p)
   {
     if (p.degree() > 0) {
       Polynomial1D derivative{p.degree() - 1};
@@ -244,7 +264,8 @@ class [[nodiscard]] Polynomial1D
   }
 
   PUGS_INLINE
-  friend Polynomial1D primitive(const Polynomial1D& p)
+  friend Polynomial1D
+  primitive(const Polynomial1D& p)
   {
     if (p.degree() > 0) {
       Polynomial1D primitive{p.degree() + 1};
@@ -259,7 +280,8 @@ class [[nodiscard]] Polynomial1D
   }
 
   PUGS_INLINE
-  double operator()(const double x) const
+  double
+  operator()(const double x) const
   {
     double value = m_coefficients[this->degree()];
     for (ssize_t i = this->degree() - 1; i >= 0; --i) {
@@ -269,7 +291,7 @@ class [[nodiscard]] Polynomial1D
     return value;
   }
 
-  Polynomial1D& operator=(Polynomial1D&&) = default;
+  Polynomial1D& operator=(Polynomial1D&&)      = default;
   Polynomial1D& operator=(const Polynomial1D&) = default;
 
   PUGS_INLINE
@@ -297,7 +319,7 @@ class [[nodiscard]] Polynomial1D
   }
 
   Polynomial1D(const Polynomial1D&) = default;
-  Polynomial1D(Polynomial1D &&)     = default;
+  Polynomial1D(Polynomial1D&&)      = default;
 
   ~Polynomial1D() = default;
 };
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 6e234910d..ec2f9c3b7 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -297,6 +297,28 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("acoustic_tau_min_max",
+                            std::function(
+
+                              [](const std::shared_ptr<const IDiscreteFunction>& rho,
+                                 const std::shared_ptr<const ItemValueVariant>& ur,
+                                 const double& dt) -> std::tuple<std::shared_ptr<const IDiscreteFunction>,
+                                                                 std::shared_ptr<const IDiscreteFunction>> {
+                                return acoustic_tau_min_max(rho, ur, dt);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("acoustic_dtau_max", std::function(
+
+                                                   [](const std::shared_ptr<const IDiscreteFunction>& rho,
+                                                      const std::shared_ptr<const ItemValueVariant>& ur,
+                                                      const double& dt) -> std::shared_ptr<const IDiscreteFunction> {
+                                                     return acoustic_dtau_max(rho, ur, dt);
+                                                   }
+
+                                                   ));
+
   this->_addBuiltinFunction("acoustic_pg_epsilon_dt",
                             std::function(
 
@@ -350,15 +372,21 @@ SchemeModule::SchemeModule()
                               [](const std::shared_ptr<const IDiscreteFunction>& rho,
                                  const std::shared_ptr<const IDiscreteFunction>& u,
                                  const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                                 const std::shared_ptr<const IDiscreteFunction>& epsilon_k,
                                  const std::shared_ptr<const IDiscreteFunction>& p,
                                  const std::shared_ptr<const IDiscreteFunction>& p_k,
                                  const std::shared_ptr<const IDiscreteFunction>& Gamma_tau,
                                  const std::shared_ptr<const IDiscreteFunction>& Gamma_inf,
+                                 const std::shared_ptr<const IDiscreteFunction>& rho0k0_expN0p1,
+                                 const std::shared_ptr<const IDiscreteFunction>& C_eps,
+                                 const std::shared_ptr<const IDiscreteFunction>& C_V,
+                                 const std::shared_ptr<const IDiscreteFunction>& C_Phi,
                                  const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
                                  const std::shared_ptr<const ItemValueVariant>& ur_variant,
                                  const double& dt_max) -> double {
-                                return acoustic_mg_entropy_dt(rho, u, epsilon, p, p_k, Gamma_tau, Gamma_inf,
-                                                              Ajr_variant, ur_variant, dt_max);
+                                return acoustic_mg_entropy_dt(rho, u, epsilon, epsilon_k, p, p_k, Gamma_tau, Gamma_inf,
+                                                              rho0k0_expN0p1, C_eps, C_V, C_Phi, Ajr_variant,
+                                                              ur_variant, dt_max);
                               }
 
                               ));
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index afeabedb9..a7b0885c0 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -131,16 +131,168 @@ acoustic_rho_dt(const DiscreteFunctionP0<Dimension, double>& rho,
     if (dt_tau.has_value()) {
       Assert(dt_tau.value() <= dt);
       dt = dt_tau.value();
+      if (dt == 0) {
+        std::cout << "Pb: q_tau=" << q_tau << '\n';
+        std::exit(0);
+      }
     }
   }
 
   if (dt < dt_max) {
     dt *= 0.95;
+    std::cout << "tau variation imposes time step: dt=" << dt << '\n';
   }
 
   return parallel::allReduceMin(dt);
 }
 
+template <size_t Dimension>
+std::tuple<const std::shared_ptr<const IDiscreteFunction>, const std::shared_ptr<const IDiscreteFunction>>
+acoustic_tau_min_max(const DiscreteFunctionP0<Dimension, double>& rho,
+                     const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                     const double& dt)
+{
+  using Rd           = TinyVector<Dimension>;
+  using MeshType     = Mesh<Connectivity<Dimension>>;
+  using MeshDataType = MeshData<Dimension>;
+
+  const NodeValue<const Rd> ur                 = ur_variant->get<NodeValue<const Rd>>();
+  const std::shared_ptr<const MeshType> p_mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+  const MeshType& mesh                         = *p_mesh;
+  if (mesh.shared_connectivity() != ur.connectivity_ptr()) {
+    throw "cannot use different connectivities";
+  }
+
+  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+  const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+  const CellValue<const double> Vj     = mesh_data.Vj();
+
+  const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+  CellValue<double> dVj{mesh.connectivity()};
+
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      double dV = 0;
+
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r = cell_nodes[R];
+        dV += dot(Cjr(j, R), ur[r]);
+      }
+
+      dVj[j] = dV;
+    });
+
+  CellValue<double> Gj{mesh.connectivity()};
+
+  if constexpr (Dimension == 2) {
+    for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      double G = 0;
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r   = cell_nodes[R];
+        const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()];
+        G += dot(ur[rp1], Rd{-ur[r][1], ur[r][0]});
+      }
+      Gj[j] = 0.5 * G;
+    }
+  }
+
+  DiscreteFunctionP0<Dimension, double> tau_min{p_mesh};
+  DiscreteFunctionP0<Dimension, double> tau_max{p_mesh};
+
+  for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+    Polynomial1D q_tau({1. / rho[j], dVj[j] / (rho[j] * Vj[j])});
+    if constexpr (Dimension == 2) {
+      q_tau = q_tau + Polynomial1D({0, 0, Gj[j] / (rho[j] * Vj[j])});
+    }
+
+    double q_tau_0  = q_tau(0);
+    double q_tau_dt = q_tau(dt);
+#warning invalid for dimension > 1
+    tau_min[j] = std::min(q_tau_0, q_tau_dt);
+    tau_max[j] = std::max(q_tau_0, q_tau_dt);
+  }
+
+  return {std::make_shared<DiscreteFunctionP0<Dimension, double>>(tau_min),
+          std::make_shared<DiscreteFunctionP0<Dimension, double>>(tau_max)};
+}
+
+template <size_t Dimension>
+std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>
+acoustic_dtau_max(const DiscreteFunctionP0<Dimension, double>& rho,
+                  const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                  const double& dt)
+{
+  using Rd           = TinyVector<Dimension>;
+  using MeshType     = Mesh<Connectivity<Dimension>>;
+  using MeshDataType = MeshData<Dimension>;
+
+  const NodeValue<const Rd> ur                 = ur_variant->get<NodeValue<const Rd>>();
+  const std::shared_ptr<const MeshType> p_mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+  const MeshType& mesh                         = *p_mesh;
+  if (mesh.shared_connectivity() != ur.connectivity_ptr()) {
+    throw "cannot use different connectivities";
+  }
+
+  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+  const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+  const CellValue<const double> Vj     = mesh_data.Vj();
+
+  const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+  CellValue<double> dVj{mesh.connectivity()};
+
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      double dV = 0;
+
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r = cell_nodes[R];
+        dV += dot(Cjr(j, R), ur[r]);
+      }
+
+      dVj[j] = dV;
+    });
+
+  CellValue<double> Gj{mesh.connectivity()};
+
+  if constexpr (Dimension == 2) {
+    for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      double G = 0;
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r   = cell_nodes[R];
+        const NodeId rp1 = cell_nodes[(R + 1) % cell_nodes.size()];
+        G += dot(ur[rp1], Rd{-ur[r][1], ur[r][0]});
+      }
+      Gj[j] = 0.5 * G;
+    }
+  }
+
+  DiscreteFunctionP0<Dimension, double> dtau_max{p_mesh};
+
+  for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+    Polynomial1D q_dtau(std::vector<double>{dVj[j] / (rho[j] * Vj[j])});
+    if constexpr (Dimension == 2) {
+      q_dtau = q_dtau + Polynomial1D({0, Gj[j] / (rho[j] * Vj[j])});
+    }
+
+#warning invalid for dimension > 1
+    dtau_max[j] = std::max(q_dtau(0), q_dtau(dt));
+  }
+
+  return std::make_shared<DiscreteFunctionP0<Dimension, double>>(dtau_max);
+}
+
 template <size_t Dimension>
 double
 acoustic_pg_epsilon_dt(const DiscreteFunctionP0<Dimension, double>& rho,
@@ -233,7 +385,7 @@ acoustic_pg_epsilon_dt(const DiscreteFunctionP0<Dimension, double>& rho,
   }
 
   if (dt < dt_max) {
-    std::cout << "energy variation imposes time step\n";
+    std::cout << "## energy variation imposes time step\n";
     dt *= 0.95;
   }
 
@@ -354,8 +506,7 @@ acoustic_pg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
       } else if constexpr (Dimension == 2) {
         return mj * p[j] * Gj[j];
       } else {
-#warning 3d
-        // throw NotImplementedError("3d");
+        throw NotImplementedError("3d");
       }
     }();
 
@@ -391,6 +542,9 @@ acoustic_pg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
     std::optional dt_S = computer.getFirstRoot();
     if (dt_S.has_value()) {
       Assert(dt_S.value() <= dt);
+      std::cout << "## entropy pg imposes time step\n";
+      std::cout << "   cell_id=" << j << '\n';
+      std::cout << "   q_S=" << p_S << '\n';
       dt = dt_S.value();
     }
   }
@@ -511,7 +665,7 @@ acoustic_mg_epsilon_dt(const DiscreteFunctionP0<Dimension, double>& rho,
   }
 
   if (dt < dt_max) {
-    std::cout << "energy variation imposes time step\n";
+    std::cout << "energy variation imposes time step: dt=" << dt << '\n';
     dt *= 0.95;
   }
 
@@ -523,10 +677,15 @@ double
 acoustic_mg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
                        const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>& u,
                        const DiscreteFunctionP0<Dimension, double>& epsilon,
+                       const DiscreteFunctionP0<Dimension, double>& epsilon_k,
                        const DiscreteFunctionP0<Dimension, double>& p,
                        const DiscreteFunctionP0<Dimension, double>& p_k,
                        const DiscreteFunctionP0<Dimension, double>& Gamma_tau,
                        const DiscreteFunctionP0<Dimension, double>& Gamma_inf,
+                       const DiscreteFunctionP0<Dimension, double>& rho0k0_expN0p1,
+                       const DiscreteFunctionP0<Dimension, double>& C_eps,
+                       const DiscreteFunctionP0<Dimension, double>& C_V,
+                       const DiscreteFunctionP0<Dimension, double>& C_Phi,
                        const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
                        const std::shared_ptr<const ItemValueVariant>& ur_variant,
                        const double& dt_max)
@@ -629,8 +788,7 @@ acoustic_mg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
       } else if constexpr (Dimension == 2) {
         return mj * p[j] * Gj[j];
       } else {
-#warning 3D
-        //        throw NotImplementedError("3d");
+        throw NotImplementedError("3d");
       }
     }();
 
@@ -649,27 +807,47 @@ acoustic_mg_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
       q_tau = q_tau + Polynomial1D({0, 0, inv_mj * Gj[j]});
     }
 
-    if (Gamma_inf[j] >= 1) {
-      Polynomial1D S1 = Polynomial1D(std::vector<double>{(mj_pi_Gj - 0.5 * dPj2)});
+    Polynomial1D S1 = Polynomial1D(std::vector<double>{(mj_pi_Gj - 0.5 * dPj2)});
 
-      // if (gamma[j] < 2) {
-      //   const double tauj_min = std::min(q_tau(0), q_tau(dt_max));
+    // if (gamma[j] < 2) {
+    //   const double tauj_min = std::min(q_tau(0), q_tau(dt_max));
 
-      S1 = S1 + Gamma_tau[j] * rho[j] * dV * dS1;
+    S1 = S1 + Gamma_tau[j] * rho[j] * dV * dS1;
+    S1 = S1 - 0.5 * C_eps[j] * inv_mj * dV * dV * DT *
+                Polynomial1D(std::vector<double>{1, inv_mj * C_V[j] * Gamma_tau[j] * rho[j]});
 
-#warning APPROX
-      S1 = S1 - 0.5 * 5e+17 * inv_mj * dV * dV * DT;
+    Polynomial1D p_S = Polynomial1D(std::vector<double>{dSj[j]}) + (inv_mj * DT) * S1;
 
-      Polynomial1D p_S = Polynomial1D(std::vector<double>{dSj[j]}) + (inv_mj * DT) * S1;
+    if ((j == 89) or (j == 90)) {
+      std::cout << "0: cell_id=" << j << '\n';
+      std::cout << "0: q_S=" << p_S << '\n';
+    }
 
-      Polynomial1DRootComputer computer{p_S, 0, dt};
-      std::optional dt_S = computer.getFirstRoot();
-      if (dt_S.has_value()) {
-        Assert(dt_S.value() <= dt);
-        dt = dt_S.value();
-      }
-    } else {
-      throw NotImplementedError("Gamma_inf<1");
+    if (Gamma_inf[j] < 1) {
+      const double tau_j_Gamma_inf         = std::pow(rho[j], -Gamma_inf[j]);
+      const double exp_Gamma_tau_Gamma_inf = std::exp(Gamma_tau[j] - Gamma_inf[j]);
+
+      p_S = (tau_j_Gamma_inf * exp_Gamma_tau_Gamma_inf) * p_S;
+
+      Polynomial1D dS2 = Polynomial1D(std::vector<double>{epsilon[j] - epsilon_k[j]}) + inv_mj * dS1 * DT -
+                         0.5 * inv_mj * inv_mj * rho0k0_expN0p1[j] * (dV * dV) * DT * DT;
+
+      p_S = p_S + (inv_mj * DT) * dS2 * (((Gamma_inf[j] - 1) * Gamma_inf[j]) * C_Phi[j] * 0.5 * (dV * dV));
+    }
+
+    Polynomial1DRootComputer computer{p_S, 0, dt};
+    std::optional dt_S = computer.getFirstRoot();
+    if (dt_S.has_value()) {
+      Assert(dt_S.value() <= dt);
+      std::cout << "## entropy mg imposes time step\n";
+      std::cout << "   cell_id=" << j << '\n';
+      std::cout << "   q_S=" << p_S << '\n';
+      dt = dt_S.value();
+    }
+
+    if ((j == 89) or (j == 90)) {
+      std::cout << "-  cell_id=" << j << '\n';
+      std::cout << "-  q_S=" << p_S << '\n';
     }
   }
 
@@ -703,6 +881,60 @@ acoustic_rho_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
   }
 }
 
+std::tuple<const std::shared_ptr<const IDiscreteFunction>, const std::shared_ptr<const IDiscreteFunction>>
+acoustic_tau_min_max(const std::shared_ptr<const IDiscreteFunction>& rho,
+                     const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                     const double& dt)
+{
+  if ((rho->descriptor().type() != DiscreteFunctionType::P0)) {
+    throw NormalError("invalid discrete function type");
+  }
+
+  std::shared_ptr mesh = rho->mesh();
+
+  switch (mesh->dimension()) {
+  case 1: {
+    return acoustic_tau_min_max(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho), ur_variant, dt);
+  }
+  case 2: {
+    return acoustic_tau_min_max(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho), ur_variant, dt);
+  }
+  case 3: {
+    return acoustic_tau_min_max(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho), ur_variant, dt);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+acoustic_dtau_max(const std::shared_ptr<const IDiscreteFunction>& rho,
+                  const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                  const double& dt)
+{
+  if ((rho->descriptor().type() != DiscreteFunctionType::P0)) {
+    throw NormalError("invalid discrete function type");
+  }
+
+  std::shared_ptr mesh = rho->mesh();
+
+  switch (mesh->dimension()) {
+  case 1: {
+    return acoustic_dtau_max(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho), ur_variant, dt);
+  }
+  case 2: {
+    return acoustic_dtau_max(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho), ur_variant, dt);
+  }
+  case 3: {
+    return acoustic_dtau_max(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho), ur_variant, dt);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
 double
 acoustic_pg_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
                        const std::shared_ptr<const IDiscreteFunction>& u,
@@ -857,19 +1089,25 @@ double
 acoustic_mg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
                        const std::shared_ptr<const IDiscreteFunction>& u,
                        const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                       const std::shared_ptr<const IDiscreteFunction>& epsilon_k,
                        const std::shared_ptr<const IDiscreteFunction>& p,
                        const std::shared_ptr<const IDiscreteFunction>& p_k,
                        const std::shared_ptr<const IDiscreteFunction>& Gamma_tau,
                        const std::shared_ptr<const IDiscreteFunction>& Gamma_inf,
+                       const std::shared_ptr<const IDiscreteFunction>& rho0k0_expN0p1,
+                       const std::shared_ptr<const IDiscreteFunction>& C_eps,
+                       const std::shared_ptr<const IDiscreteFunction>& C_V,
+                       const std::shared_ptr<const IDiscreteFunction>& C_Phi,
                        const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
                        const std::shared_ptr<const ItemValueVariant>& ur_variant,
                        const double& dt_max)
 {
-  if (not checkDiscretizationType({rho, u, epsilon, p, p_k, Gamma_tau, Gamma_inf}, DiscreteFunctionType::P0)) {
+  if (not checkDiscretizationType({rho, u, epsilon, p, p_k, Gamma_tau, Gamma_inf, C_eps, C_V},
+                                  DiscreteFunctionType::P0)) {
     throw NormalError("acoustic solver expects P0 functions");
   }
 
-  std::shared_ptr mesh = getCommonMesh({rho, u, epsilon, p, p_k, Gamma_tau, Gamma_inf});
+  std::shared_ptr mesh = getCommonMesh({rho, u, epsilon, p, p_k, Gamma_tau, Gamma_inf, C_eps, C_V});
 
   switch (mesh->dimension()) {
   case 1: {
@@ -877,10 +1115,15 @@ acoustic_mg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
     return acoustic_mg_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon_k),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_k),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_tau),
-                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_inf), Ajr_variant,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_inf),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho0k0_expN0p1),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*C_eps),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*C_V),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*C_Phi), Ajr_variant,
                                   ur_variant, dt_max);
   }
   case 2: {
@@ -888,10 +1131,15 @@ acoustic_mg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
     return acoustic_mg_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon_k),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_k),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_tau),
-                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_inf), Ajr_variant,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_inf),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho0k0_expN0p1),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*C_eps),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*C_V),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*C_Phi), Ajr_variant,
                                   ur_variant, dt_max);
   }
   case 3: {
@@ -899,10 +1147,15 @@ acoustic_mg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
     return acoustic_mg_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*epsilon_k),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p_k),
                                   dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_tau),
-                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_inf), Ajr_variant,
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*Gamma_inf),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho0k0_expN0p1),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*C_eps),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*C_V),
+                                  dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*C_Phi), Ajr_variant,
                                   ur_variant, dt_max);
   }
   default: {
@@ -1348,8 +1601,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
 
           new_rho[j] = 1. / (1. / rho[j] + dt / mj * (dV + dt * G));
         } else {
-#warning 3d
-          // throw NotImplementedError("3D with GCL error");
+          throw NotImplementedError("3D with GCL error");
         }
       });
 
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index e7bd26fa8..0d6f05bad 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -17,6 +17,15 @@ double acoustic_rho_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
                        const std::shared_ptr<const ItemValueVariant>& ur,
                        const double& dt_max);
 
+std::tuple<const std::shared_ptr<const IDiscreteFunction>, const std::shared_ptr<const IDiscreteFunction>>
+acoustic_tau_min_max(const std::shared_ptr<const IDiscreteFunction>& rho,
+                     const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                     const double& dt);
+
+std::shared_ptr<const IDiscreteFunction> acoustic_dtau_max(const std::shared_ptr<const IDiscreteFunction>& rho,
+                                                           const std::shared_ptr<const ItemValueVariant>& ur,
+                                                           const double& dt_max);
+
 double acoustic_pg_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
                               const std::shared_ptr<const IDiscreteFunction>& u,
                               const std::shared_ptr<const IDiscreteFunction>& E,
@@ -55,10 +64,15 @@ double acoustic_mg_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rh
 double acoustic_mg_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
                               const std::shared_ptr<const IDiscreteFunction>& u,
                               const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                              const std::shared_ptr<const IDiscreteFunction>& epsilon_k,
                               const std::shared_ptr<const IDiscreteFunction>& p,
                               const std::shared_ptr<const IDiscreteFunction>& p_k,
                               const std::shared_ptr<const IDiscreteFunction>& Gamma_tau,
                               const std::shared_ptr<const IDiscreteFunction>& Gamma_inf,
+                              const std::shared_ptr<const IDiscreteFunction>& rho0k0_expN0p1,
+                              const std::shared_ptr<const IDiscreteFunction>& C_eps,
+                              const std::shared_ptr<const IDiscreteFunction>& C_V,
+                              const std::shared_ptr<const IDiscreteFunction>& C_Phi,
                               const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
                               const std::shared_ptr<const ItemValueVariant>& ur_variant,
                               const double& dt_max);
@@ -110,8 +124,8 @@ class AcousticSolverHandler
           const std::shared_ptr<const IDiscreteFunction>& p,
           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
 
-    IAcousticSolver()                  = default;
-    IAcousticSolver(IAcousticSolver&&) = default;
+    IAcousticSolver()                             = default;
+    IAcousticSolver(IAcousticSolver&&)            = default;
     IAcousticSolver& operator=(IAcousticSolver&&) = default;
 
     virtual ~IAcousticSolver() = default;
-- 
GitLab