diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index f324138fb8a27cbab81a6ac38992c46d621cbe98..ef4df609f50d032d360242e2d1020d5bbe2ed518 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -311,6 +311,20 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("legacy_entropy_dt",
+                            std::function(
+
+                              [](const std::shared_ptr<const IDiscreteFunction>& rho,
+                                 const std::shared_ptr<const IDiscreteFunction>& u,
+                                 const std::shared_ptr<const IDiscreteFunction>& E,
+                                 const std::shared_ptr<const IDiscreteFunction>& p,
+                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr,
+                                 const std::shared_ptr<const ItemValueVariant>& ur, const double& dt_max) -> double {
+                                return acoustic_entropy_dt(rho, u, E, p, Ajr, ur, dt_max);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("acoustic_entropy_dt",
                             std::function(
 
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 5cfc32a8832e49666fa61e5a7ce143500fc9a355..09698e685833781c802b777396ea6a80f72f6145 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -237,16 +237,143 @@ acoustic_epsilon_dt(const DiscreteFunctionP0<Dimension, double>& rho,
     }
   }
 
-  // if (dt < dt_max) {
-  //   dt *= 0.95 * dt;
-  // }
-
-  // if (dt != dt_max) {
-  //   std::cout << "internal energy variation imposes time step\n";
-  // }
   return dt;
 }
 
+template <size_t Dimension>
+double
+legacy_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
+                  const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>& u,
+                  const DiscreteFunctionP0<Dimension, double>& E,
+                  const DiscreteFunctionP0<Dimension, double>& p,
+                  const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                  const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                  const double& dt_max)
+{
+  double dt = dt_max;
+
+  using Rd           = TinyVector<Dimension>;
+  using Rdxd         = TinyMatrix<Dimension>;
+  using MeshType     = Mesh<Connectivity<Dimension>>;
+  using MeshDataType = MeshData<Dimension>;
+
+  const MeshType& mesh = dynamic_cast<const MeshType&>(*rho.mesh());
+
+  const NodeValue<const Rd> ur           = ur_variant->get<NodeValue<const Rd>>();
+  const NodeValuePerCell<const Rdxd> Ajr = Ajr_variant->get<NodeValuePerCell<const Rdxd>>();
+  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;
+    }
+  }
+
+  CellValue<double> dSj{mesh.connectivity()};
+
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      double dS = 0;
+
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r = cell_nodes[R];
+
+        const Rd du = u[j] - ur[r];
+        dS += dot(Ajr(j, R) * du, du);
+      }
+
+      dSj[j] = dS;
+    });
+
+  CellValue<Rd> dPj{mesh.connectivity()};
+
+  parallel_for(
+    mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+      const auto& cell_nodes = cell_to_node_matrix[j];
+
+      Rd dP = zero;
+
+      for (size_t R = 0; R < cell_nodes.size(); ++R) {
+        const NodeId r = cell_nodes[R];
+
+        Rd du = u[j] - ur[r];
+        dP += Ajr(j, R) * du;
+      }
+
+      dPj[j] = dP;
+    });
+
+#warning fixed gamma value
+  const double gamma = 1.4;
+
+  for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+    const double inv_mj = 1 / (rho[j] * Vj[j]);
+    const double inv_Vj = 1 / Vj[j];
+
+    Polynomial1D DT({0, 1});
+    std::vector<double> delta_S = {dSj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])};
+    if constexpr (Dimension == 2) {
+      delta_S[1] += p[j] * Gj[j];
+    }
+    Polynomial1D p_S0(delta_S);
+
+    std::vector<double> delta_V{dVj[j]};
+    if constexpr (Dimension == 2) {
+      delta_V.push_back(Gj[j]);
+    }
+    Polynomial1D p_deltaV(delta_V);
+
+    Polynomial1D p_S1({dSj[j] - p[j] * dVj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])});
+
+    Polynomial1D p_S =                                                                                 //
+      p_S0                                                                                             //
+      + (inv_Vj * DT) * (((gamma - 1) * p_S1) * p_deltaV + (gamma - 2) * p[j] * p_deltaV * p_deltaV)   //
+      + (inv_Vj * DT) * (inv_Vj * DT) * (((gamma - 1) * (gamma - 2)) * p_S1) * p_deltaV * p_deltaV;
+
+    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();
+    }
+  }
+
+  return parallel::allReduceMin(dt);
+}
+
 template <size_t Dimension>
 double
 acoustic_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
@@ -264,7 +391,6 @@ acoustic_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
   using MeshType     = Mesh<Connectivity<Dimension>>;
   using MeshDataType = MeshData<Dimension>;
 
-#warning get mesh using all discrete functions
   const MeshType& mesh = dynamic_cast<const MeshType&>(*rho.mesh());
 
   const NodeValue<const Rd> ur           = ur_variant->get<NodeValue<const Rd>>();
@@ -358,7 +484,7 @@ acoustic_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
     }
     Polynomial1D p_S0(delta_S);
 
-    std::vector<double> delta_V{dVj[j] /*, dG */};
+    std::vector<double> delta_V{dVj[j]};
     if constexpr (Dimension == 2) {
       delta_V.push_back(Gj[j]);
     }
@@ -379,9 +505,6 @@ acoustic_entropy_dt(const DiscreteFunctionP0<Dimension, double>& rho,
     }
   }
 
-  //     if (dt != dt_max) {
-  //       std::cout << "entropy variation imposes time step\n";
-  //     }
   return parallel::allReduceMin(dt);
 }
 
@@ -458,6 +581,52 @@ acoustic_epsilon_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
   }
 }
 
+double
+legacy_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                  const std::shared_ptr<const IDiscreteFunction>& u,
+                  const std::shared_ptr<const IDiscreteFunction>& E,
+                  const std::shared_ptr<const IDiscreteFunction>& p,
+                  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, E, p}, DiscreteFunctionType::P0)) {
+    throw NormalError("acoustic solver expects P0 functions");
+  }
+
+  std::shared_ptr mesh = getCommonMesh({rho, u, E, p});
+
+  switch (mesh->dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    return legacy_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                             dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                             dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
+                             dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
+                             dt_max);
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    return legacy_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                             dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                             dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
+                             dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
+                             dt_max);
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    return legacy_entropy_dt(dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*rho),
+                             dynamic_cast<const DiscreteFunctionP0<Dimension, TinyVector<Dimension>>&>(*u),
+                             dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*E),
+                             dynamic_cast<const DiscreteFunctionP0<Dimension, double>&>(*p), Ajr_variant, ur_variant,
+                             dt_max);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
 double
 acoustic_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
                     const std::shared_ptr<const IDiscreteFunction>& u,
@@ -820,131 +989,6 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
     return F;
   }
 
-  /* double */
-  /* _computeMaxEntropyDt(const MeshType& mesh, */
-  /*                      const DiscreteScalarFunction& rho, */
-  /*                      const DiscreteVectorFunction& u, */
-  /*                      const DiscreteScalarFunction& p, */
-  /*                      const double dt_max) const */
-  /* { */
-  /*   double dt = dt_max; */
-
-  /*   //     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), m_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(m_ur[rp1], Rd{m_ur[r][1], -m_ur[r][0]}); */
-  /*   //         } */
-  /*   //         Gj[j] = 0.5 * G; */
-  /*   //       } */
-  /*   //     } */
-
-  /*   //     CellValue<double> dSj{mesh.connectivity()}; */
-
-  /*   //     parallel_for( */
-  /*   //       mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { */
-  /*   //         const auto& cell_nodes = cell_to_node_matrix[j]; */
-
-  /*   //         double dS = 0; */
-
-  /*   //         for (size_t R = 0; R < cell_nodes.size(); ++R) { */
-  /*   //           const NodeId r = cell_nodes[R]; */
-
-  /*   //           const Rd du = u[j] - m_ur[r]; */
-  /*   //           dS += dot(m_Ajr(j, R) * du, du); */
-  /*   //         } */
-
-  /*   //         dSj[j] = dS; */
-  /*   //       }); */
-
-  /*   //     CellValue<Rd> dPj{mesh.connectivity()}; */
-
-  /*   //     parallel_for( */
-  /*   //       mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { */
-  /*   //         const auto& cell_nodes = cell_to_node_matrix[j]; */
-
-  /*   //         Rd dP = zero; */
-
-  /*   //         for (size_t R = 0; R < cell_nodes.size(); ++R) { */
-  /*   //           const NodeId r = cell_nodes[R]; */
-
-  /*   //           Rd du = u[j] - m_ur[r]; */
-  /*   //           dP += m_Ajr(j, R) * du; */
-  /*   //         } */
-
-  /*   //         dPj[j] = dP; */
-  /*   //       }); */
-
-  /*   // #warning fixed gamma value */
-  /*   //     const double gamma = 1.4; */
-
-  /*   //     for (CellId j = 0; j < mesh.numberOfCells(); ++j) { */
-  /*   //       const double inv_mj = 1 / (rho[j] * Vj[j]); */
-  /*   //       const double inv_Vj = 1 / Vj[j]; */
-
-  /*   //       Polynomial1D DT({0, 1}); */
-  /*   //       std::vector<double> delta_S = {dSj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])}; */
-  /*   //       if constexpr (Dimension == 2) { */
-  /*   //         delta_S[1] += p[j] * Gj[j]; */
-  /*   //       } */
-  /*   //       Polynomial1D p_S0(delta_S); */
-
-  /*   //       std::vector<double> delta_V{dVj[j] /\*, dG *\/}; */
-  /*   //       if constexpr (Dimension == 2) { */
-  /*   //         delta_V.push_back(Gj[j]); */
-  /*   //       } */
-  /*   //       Polynomial1D p_deltaV(delta_V); */
-
-  /*   //       Polynomial1D p_S1({dSj[j] - p[j] * dVj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])}); */
-
-  /*   //       Polynomial1D p_S =                                                                                 // */
-  /*   //         p_S0                                                                                             // */
-  /*   //         + (inv_Vj * DT) * (((gamma - 1) * p_S1) * p_deltaV + (gamma - 2) * p[j] * p_deltaV * p_deltaV)   // */
-  /*   //         + (inv_Vj * DT) * (inv_Vj * DT) * (((gamma - 1) * (gamma - 2)) * p_S1) * p_deltaV * p_deltaV; */
-
-  /*   //       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(); */
-  /*   //       } */
-  /*   //     } */
-
-  /*   //     if (dt != dt_max) { */
-  /*   //       std::cout << "entropy variation imposes time step\n"; */
-  /*   //     } */
-  /*   return parallel::allReduceMin(dt); */
-  /* } */
-
  public:
   std::tuple<const std::shared_ptr<const SubItemValuePerItemVariant>,
              const std::shared_ptr<const ItemValueVariant>,
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 573b25672c74a194762e751e4f9de32dab7a2148..a00ce4e143f1d66fe8d851b9b2c2cc5a8422d39e 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -33,6 +33,22 @@ double acoustic_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
                            const std::shared_ptr<const ItemValueVariant>& ur_variant,
                            const double& dt_max);
 
+double legacy_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                         const std::shared_ptr<const IDiscreteFunction>& u,
+                         const std::shared_ptr<const IDiscreteFunction>& E,
+                         const std::shared_ptr<const IDiscreteFunction>& p,
+                         const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                         const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                         const double& dt_max);
+
+double acoustic_entropy_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                           const std::shared_ptr<const IDiscreteFunction>& u,
+                           const std::shared_ptr<const IDiscreteFunction>& E,
+                           const std::shared_ptr<const IDiscreteFunction>& p,
+                           const std::shared_ptr<const SubItemValuePerItemVariant>& Ajr_variant,
+                           const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                           const double& dt_max);
+
 class AcousticSolverHandler
 {
  public: