diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index e3c6f1b96bf5ae660f7e3648ab70858edd4c20a3..f324138fb8a27cbab81a6ac38992c46d621cbe98 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -288,6 +288,43 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("acoustic_rho_dt",
+                            std::function(
+
+                              [](const std::shared_ptr<const IDiscreteFunction>& rho,
+                                 const std::shared_ptr<const ItemValueVariant>& ur,
+                                 const double& dt_max) -> double { return acoustic_rho_dt(rho, ur, dt_max); }
+
+                              ));
+
+  this->_addBuiltinFunction("acoustic_epsilon_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_epsilon_dt(rho, u, E, p, Ajr, ur, dt_max);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("acoustic_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("glace_fluxes", std::function(
 
                                               [](const std::shared_ptr<const IDiscreteFunction>& rho,
@@ -296,7 +333,8 @@ SchemeModule::SchemeModule()
                                                  const std::shared_ptr<const IDiscreteFunction>& p,
                                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                                    bc_descriptor_list)
-                                                -> std::tuple<std::shared_ptr<const ItemValueVariant>,
+                                                -> std::tuple<std::shared_ptr<const SubItemValuePerItemVariant>,
+                                                              std::shared_ptr<const ItemValueVariant>,
                                                               std::shared_ptr<const SubItemValuePerItemVariant>> {
                                                 return AcousticSolverHandler{getCommonMesh({rho, c, u, p})}
                                                   .solver()
@@ -337,7 +375,8 @@ SchemeModule::SchemeModule()
                                  const std::shared_ptr<const IDiscreteFunction>& p,
                                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
                                    bc_descriptor_list)
-                                -> std::tuple<std::shared_ptr<const ItemValueVariant>,
+                                -> std::tuple<std::shared_ptr<const SubItemValuePerItemVariant>,
+                                              std::shared_ptr<const ItemValueVariant>,
                                               std::shared_ptr<const SubItemValuePerItemVariant>> {
                                 return AcousticSolverHandler{getCommonMesh({rho, c, u, p})}
                                   .solver()
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 04d0545bbaec039d6f5d8395ec05de45fd15347f..5cfc32a8832e49666fa61e5a7ce143500fc9a355 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -64,6 +64,446 @@ acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c)
   }
 }
 
+template <size_t Dimension>
+double
+acoustic_rho_dt(const DiscreteFunctionP0<Dimension, double>& rho,
+                const std::shared_ptr<const ItemValueVariant>& ur_variant,
+                const double& dt_max)
+{
+  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 MeshType& mesh         = dynamic_cast<const MeshType&>(*rho.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;
+    }
+  }
+
+  double dt = dt_max;
+
+  for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+    Polynomial1D p_tau({1. / rho[j], dVj[j] / (rho[j] * Vj[j])});
+    if constexpr (Dimension == 2) {
+      p_tau = p_tau + Polynomial1D({0, 0, Gj[j] / (rho[j] * Vj[j])});
+    }
+    Polynomial1DRootComputer computer{p_tau, 0, dt};
+    std::optional dt_tau = computer.getFirstRoot();
+    if (dt_tau.has_value()) {
+      Assert(dt_tau.value() <= dt);
+      dt = dt_tau.value();
+    }
+  }
+
+  // if (dt != dt_max) {
+  //   std::cout << "volume variation imposes time step\n";
+  // }
+
+  // if (dt < dt_max) {
+  //   dt *= 0.95 * dt;
+  // }
+
+  return dt;
+}
+
+template <size_t Dimension>
+double
+acoustic_epsilon_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)
+{
+  using Rd           = TinyVector<Dimension>;
+  using Rdxd         = TinyMatrix<Dimension>;
+  using MeshType     = Mesh<Connectivity<Dimension>>;
+  using MeshDataType = MeshData<Dimension>;
+
+  double dt = dt_max;
+
+  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> 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];
+
+        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];
+
+        const Rd du = u[j] - ur[r];
+        dP += Ajr(j, R) * du;
+      }
+
+      dPj[j] = dP;
+    });
+
+  for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+    const double inv_mj = 1 / (rho[j] * Vj[j]);
+
+    Polynomial1D p_e(
+      {E[j] - 0.5 * dot(u[j], u[j]), inv_mj * (dSj[j] - p[j] * dVj[j]), -0.5 * inv_mj * inv_mj * dot(dPj[j], dPj[j])});
+    Polynomial1DRootComputer computer{p_e, 0, dt};
+    std::optional dt_e = computer.getFirstRoot();
+    if (dt_e.has_value()) {
+      Assert(dt_e.value() <= dt);
+      dt = dt_e.value();
+    }
+  }
+
+  // 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
+acoustic_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>;
+
+#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>>();
+  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] /*, 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);
+}
+
+double
+acoustic_rho_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                const std::shared_ptr<const ItemValueVariant>& ur,
+                const double& dt_max)
+{
+  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_rho_dt(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho), ur, dt_max);
+  }
+  case 2: {
+    return acoustic_rho_dt(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho), ur, dt_max);
+  }
+  case 3: {
+    return acoustic_rho_dt(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho), ur, dt_max);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+double
+acoustic_epsilon_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 acoustic_epsilon_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 acoustic_epsilon_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 acoustic_epsilon_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,
+                    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 acoustic_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 acoustic_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 acoustic_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");
+  }
+  }
+}
+
 template <size_t Dimension>
 class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler::IAcousticSolver
 {
@@ -380,296 +820,135 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
     return F;
   }
 
-  double
-  _computeMaxDensityDt(const MeshType& mesh,
-                       const DiscreteFunctionP0<Dimension, double>& rho,
-                       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;
-    //   }
-    // }
-
-    // for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
-    //   Polynomial1D p_tau({1. / rho[j], dVj[j] / (rho[j] * Vj[j])});
-    //   if constexpr (Dimension == 2) {
-    //     p_tau = p_tau + Polynomial1D({0, 0, Gj[j] / (rho[j] * Vj[j])});
-    //   }
-    //   Polynomial1DRootComputer computer{p_tau, 0, dt};
-    //   std::optional dt_tau = computer.getFirstRoot();
-    //   if (dt_tau.has_value()) {
-    //     Assert(dt_tau.value() <= dt);
-    //     dt = dt_tau.value();
-    //   }
-    // }
-
-    // if (dt != dt_max) {
-    //   std::cout << "volume variation imposes time step\n";
-    // }
-
-    // if (dt < dt_max) {
-    //   dt *= 0.95 * dt;
-    // }
-
-    return parallel::allReduceMin(dt);
-  }
-
-  double
-  _computeMaxEnergyDt(const MeshType& mesh,
-                      const DiscreteScalarFunction& rho,
-                      const DiscreteVectorFunction& u,
-                      const DiscreteScalarFunction& E,
-                      const DiscreteScalarFunction& p,
-                      const double dt_max) const
-  {
-    double dt = dt_max;
-
-    // MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+  /* double */
+  /* _computeMaxEntropyDt(const MeshType& mesh, */
+  /*                      const DiscreteScalarFunction& rho, */
+  /*                      const DiscreteVectorFunction& u, */
+  /*                      const DiscreteScalarFunction& p, */
+  /*                      const double dt_max) const */
+  /* { */
+  /*   double dt = dt_max; */
 
-    // const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
-    // const CellValue<const double> Vj     = mesh_data.Vj();
+  /*   //     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh); */
 
-    // const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+  /*   //     const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr(); */
+  /*   //     const CellValue<const double> Vj     = mesh_data.Vj(); */
 
-    // CellValue<double> dVj{mesh.connectivity()};
+  /*   //     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix(); */
 
-    // parallel_for(
-    //   mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-    //     const auto& cell_nodes = cell_to_node_matrix[j];
+  /*   //     CellValue<double> dVj{mesh.connectivity()}; */
 
-    //     double dV = 0;
+  /*   //     parallel_for( */
+  /*   //       mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { */
+  /*   //         const auto& cell_nodes = cell_to_node_matrix[j]; */
 
-    //     for (size_t R = 0; R < cell_nodes.size(); ++R) {
-    //       const NodeId r = cell_nodes[R];
-    //       dV += dot(Cjr(j, R), m_ur[r]);
-    //     }
+  /*   //         double dV = 0; */
 
-    //     dVj[j] = dV;
-    //   });
+  /*   //         for (size_t R = 0; R < cell_nodes.size(); ++R) { */
+  /*   //           const NodeId r = cell_nodes[R]; */
+  /*   //           dV += dot(Cjr(j, R), m_ur[r]); */
+  /*   //         } */
 
-    // CellValue<double> dSj{mesh.connectivity()};
+  /*   //         dVj[j] = dV; */
+  /*   //       }); */
 
-    // parallel_for(
-    //   mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-    //     const auto& cell_nodes = cell_to_node_matrix[j];
+  /*   //     CellValue<double> Gj{mesh.connectivity()}; */
 
-    //     double dS = 0;
+  /*   //     if constexpr (Dimension == 2) { */
+  /*   //       for (CellId j = 0; j < mesh.numberOfCells(); ++j) { */
+  /*   //         const auto& cell_nodes = cell_to_node_matrix[j]; */
 
-    //     for (size_t R = 0; R < cell_nodes.size(); ++R) {
-    //       const NodeId r = cell_nodes[R];
+  /*   //         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; */
+  /*   //       } */
+  /*   //     } */
 
-    //       Rd du = u[j] - m_ur[r];
-    //       dS += dot(m_Ajr(j, R) * du, du);
-    //     }
+  /*   //     CellValue<double> dSj{mesh.connectivity()}; */
 
-    //     dSj[j] = dS;
-    //   });
+  /*   //     parallel_for( */
+  /*   //       mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { */
+  /*   //         const auto& cell_nodes = cell_to_node_matrix[j]; */
 
-    // CellValue<Rd> dPj{mesh.connectivity()};
+  /*   //         double dS = 0; */
 
-    // parallel_for(
-    //   mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-    //     const auto& cell_nodes = cell_to_node_matrix[j];
+  /*   //         for (size_t R = 0; R < cell_nodes.size(); ++R) { */
+  /*   //           const NodeId r = cell_nodes[R]; */
 
-    //     Rd dP = zero;
+  /*   //           const Rd du = u[j] - m_ur[r]; */
+  /*   //           dS += dot(m_Ajr(j, R) * du, du); */
+  /*   //         } */
 
-    //     for (size_t R = 0; R < cell_nodes.size(); ++R) {
-    //       const NodeId r = cell_nodes[R];
+  /*   //         dSj[j] = dS; */
+  /*   //       }); */
 
-    //       Rd du = u[j] - m_ur[r];
-    //       dP += m_Ajr(j, R) * du;
-    //     }
+  /*   //     CellValue<Rd> dPj{mesh.connectivity()}; */
 
-    //     dPj[j] = dP;
-    //   });
+  /*   //     parallel_for( */
+  /*   //       mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { */
+  /*   //         const auto& cell_nodes = cell_to_node_matrix[j]; */
 
-    // for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
-    //   const double inv_mj = 1 / (rho[j] * Vj[j]);
+  /*   //         Rd dP = zero; */
 
-    //   Polynomial1D p_e(
-    //     {E[j] - 0.5 * dot(u[j], u[j]), inv_mj * (dSj[j] - p[j] * dVj[j]), -inv_mj * inv_mj * dot(dPj[j], dPj[j])});
-    //   Polynomial1DRootComputer computer{p_e, 0, dt};
-    //   std::optional dt_e = computer.getFirstRoot();
-    //   if (dt_e.has_value()) {
-    //     Assert(dt_e.value() <= dt);
-    //     dt = dt_e.value();
-    //   }
-    // }
+  /*   //         for (size_t R = 0; R < cell_nodes.size(); ++R) { */
+  /*   //           const NodeId r = cell_nodes[R]; */
 
-    // if (dt < dt_max) {
-    //   dt *= 0.95 * dt;
-    // }
+  /*   //           Rd du = u[j] - m_ur[r]; */
+  /*   //           dP += m_Ajr(j, R) * du; */
+  /*   //         } */
 
-    // if (dt != dt_max) {
-    //   std::cout << "internal energy variation imposes time step\n";
-    // }
+  /*   //         dPj[j] = dP; */
+  /*   //       }); */
 
-    return parallel::allReduceMin(dt);
-  }
+  /*   // #warning fixed gamma value */
+  /*   //     const double gamma = 1.4; */
 
-  double
-  _computeMaxEntropyDt(const MeshType& mesh,
-                       const DiscreteScalarFunction& rho,
-                       const DiscreteVectorFunction& u,
-                       const DiscreteScalarFunction& p,
-                       const double dt_max) const
-  {
-    double dt = dt_max;
+  /*   //     for (CellId j = 0; j < mesh.numberOfCells(); ++j) { */
+  /*   //       const double inv_mj = 1 / (rho[j] * Vj[j]); */
+  /*   //       const double inv_Vj = 1 / Vj[j]; */
 
-    //     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+  /*   //       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); */
 
-    //     const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
-    //     const CellValue<const double> Vj     = mesh_data.Vj();
+  /*   //       std::vector<double> delta_V{dVj[j] /\*, dG *\/}; */
+  /*   //       if constexpr (Dimension == 2) { */
+  /*   //         delta_V.push_back(Gj[j]); */
+  /*   //       } */
+  /*   //       Polynomial1D p_deltaV(delta_V); */
 
-    //     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+  /*   //       Polynomial1D p_S1({dSj[j] - p[j] * dVj[j], -0.5 * inv_mj * dot(dPj[j], dPj[j])}); */
 
-    //     CellValue<double> dVj{mesh.connectivity()};
+  /*   //       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; */
 
-    //     parallel_for(
-    //       mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-    //         const auto& cell_nodes = cell_to_node_matrix[j];
+  /*   //       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(); */
+  /*   //       } */
+  /*   //     } */
 
-    //         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);
-  }
+  /*   //     if (dt != dt_max) { */
+  /*   //       std::cout << "entropy variation imposes time step\n"; */
+  /*   //     } */
+  /*   return parallel::allReduceMin(dt); */
+  /* } */
 
  public:
-  std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
+  std::tuple<const std::shared_ptr<const SubItemValuePerItemVariant>,
+             const std::shared_ptr<const ItemValueVariant>,
+             const std::shared_ptr<const SubItemValuePerItemVariant>>
   compute_fluxes(const SolverType& solver_type,
                  const std::shared_ptr<const IDiscreteFunction>& i_rho,
                  const std::shared_ptr<const IDiscreteFunction>& i_c,
@@ -707,7 +986,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
     NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
     NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, p);
 
-    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
+    return std::make_tuple(std::make_shared<const SubItemValuePerItemVariant>(Ajr),
+                           std::make_shared<const ItemValueVariant>(ur),
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));
   }
 
@@ -826,7 +1106,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
         const std::shared_ptr<const IDiscreteFunction>& p,
         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
-    auto [ur, Fjr] = compute_fluxes(solver_type, rho, c, u, p, bc_descriptor_list);
+    auto [Ajr, ur, Fjr] = compute_fluxes(solver_type, rho, c, u, p, bc_descriptor_list);
     return apply_fluxes(dt, rho, u, E, ur, Fjr);
   }
 
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 489b8e795d843f9677daf562214bd0266f803baa..573b25672c74a194762e751e4f9de32dab7a2148 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -13,6 +13,26 @@ class SubItemValuePerItemVariant;
 
 double acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c);
 
+double acoustic_rho_dt(const std::shared_ptr<const IDiscreteFunction>& rho,
+                       const std::shared_ptr<const ItemValueVariant>& ur,
+                       const double& dt_max);
+
+double acoustic_epsilon_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:
@@ -25,7 +45,8 @@ class AcousticSolverHandler
  private:
   struct IAcousticSolver
   {
-    virtual std::tuple<const std::shared_ptr<const ItemValueVariant>,
+    virtual std::tuple<const std::shared_ptr<const SubItemValuePerItemVariant>,
+                       const std::shared_ptr<const ItemValueVariant>,
                        const std::shared_ptr<const SubItemValuePerItemVariant>>
     compute_fluxes(
       const SolverType& solver_type,