diff --git a/src/language/algorithms/AcousticSolverAlgorithm.cpp b/src/language/algorithms/AcousticSolverAlgorithm.cpp
deleted file mode 100644
index 0af71d5a5eb7e07db2625bfa405309a47213fcf2..0000000000000000000000000000000000000000
--- a/src/language/algorithms/AcousticSolverAlgorithm.cpp
+++ /dev/null
@@ -1,227 +0,0 @@
-#include <language/algorithms/AcousticSolverAlgorithm.hpp>
-
-#include <language/utils/InterpolateItemValue.hpp>
-#include <scheme/AcousticSolver.hpp>
-#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
-#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-
-#include <iomanip>
-
-template <size_t Dimension>
-AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
-  const AcousticSolverType& solver_type,
-  std::shared_ptr<const IMesh> i_mesh,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-  const FunctionSymbolId& rho_id,
-  const FunctionSymbolId& u_id,
-  const FunctionSymbolId& p_id)
-{
-  using ConnectivityType = Connectivity<Dimension>;
-  using MeshType         = Mesh<ConnectivityType>;
-  using MeshDataType     = MeshData<Dimension>;
-  using UnknownsType     = FiniteVolumesEulerUnknowns<MeshType>;
-
-  std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-  typename AcousticSolverOld<MeshType>::BoundaryConditionList bc_list;
-  {
-    constexpr ItemType FaceType = [] {
-      if constexpr (Dimension > 1) {
-        return ItemType::face;
-      } else {
-        return ItemType::node;
-      }
-    }();
-
-    for (const auto& bc_descriptor : bc_descriptor_list) {
-      bool is_valid_boundary_condition = true;
-
-      switch (bc_descriptor->type()) {
-      case IBoundaryConditionDescriptor::Type::symmetry: {
-        using SymmetryBoundaryCondition = typename AcousticSolverOld<MeshType>::SymmetryBoundaryCondition;
-
-        const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
-          dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-        for (size_t i_ref_face_list = 0;
-             i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
-          const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
-          const RefId& ref          = ref_face_list.refId();
-          if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-            bc_list.push_back(
-              SymmetryBoundaryCondition{MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list)});
-          }
-        }
-        is_valid_boundary_condition = true;
-        break;
-      }
-
-      case IBoundaryConditionDescriptor::Type::dirichlet: {
-        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
-          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
-        if (dirichlet_bc_descriptor.name() == "velocity") {
-          using VelocityBoundaryCondition = typename AcousticSolverOld<MeshType>::VelocityBoundaryCondition;
-
-          for (size_t i_ref_face_list = 0;
-               i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
-            const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
-            const RefId& ref          = ref_face_list.refId();
-            if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
-              MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
-
-              const FunctionSymbolId velocity_id = dirichlet_bc_descriptor.rhsSymbolId();
-
-              const auto& node_list = mesh_node_boundary.nodeList();
-
-              Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-                TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, mesh->xr(), node_list);
-
-              bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
-            }
-          }
-        } else if (dirichlet_bc_descriptor.name() == "pressure") {
-          using PressureBoundaryCondition = typename AcousticSolverOld<MeshType>::PressureBoundaryCondition;
-
-          for (size_t i_ref_face_list = 0;
-               i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
-            const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
-            const RefId& ref          = ref_face_list.refId();
-            if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
-              const auto& face_list = ref_face_list.list();
-
-              const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
-
-              Array<const double> face_values = [&] {
-                if constexpr (Dimension == 1) {
-                  return InterpolateItemValue<double(
-                    TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh->xr(), face_list);
-                } else {
-                  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-                  return InterpolateItemValue<double(
-                    TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(), face_list);
-                }
-              }();
-
-              bc_list.push_back(PressureBoundaryCondition{face_list, face_values});
-            }
-          }
-        } else {
-          is_valid_boundary_condition = false;
-        }
-        break;
-      }
-      default: {
-        is_valid_boundary_condition = false;
-      }
-      }
-      if (not is_valid_boundary_condition) {
-        std::ostringstream error_msg;
-        error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver";
-        throw NormalError(error_msg.str());
-      }
-    }
-  }
-
-  UnknownsType unknowns(*mesh);
-
-  {
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-    unknowns.rhoj() =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id, mesh_data.xj());
-
-    unknowns.pj() =
-      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj());
-
-    unknowns.uj() =
-      InterpolateItemValue<TinyVector<Dimension>(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id,
-                                                                                                               mesh_data
-                                                                                                                 .xj());
-  }
-  unknowns.gammaj().fill(1.4);
-
-  AcousticSolverOld acoustic_solver(mesh, bc_list, solver_type);
-
-  const double tmax = 0.2;
-  double t          = 0;
-
-  int itermax   = std::numeric_limits<int>::max();
-  int iteration = 0;
-
-  CellValue<double>& rhoj              = unknowns.rhoj();
-  CellValue<double>& ej                = unknowns.ej();
-  CellValue<double>& pj                = unknowns.pj();
-  CellValue<double>& gammaj            = unknowns.gammaj();
-  CellValue<double>& cj                = unknowns.cj();
-  CellValue<TinyVector<Dimension>>& uj = unknowns.uj();
-  CellValue<double>& Ej                = unknowns.Ej();
-  CellValue<double>& mj                = unknowns.mj();
-  CellValue<double>& inv_mj            = unknowns.invMj();
-
-  BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-  block_eos.updateEandCFromRhoP();
-
-  {
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-    const CellValue<const double> Vj = mesh_data.Vj();
-
-    parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
-
-    parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; });
-
-    parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
-  }
-
-  while ((t < tmax) and (iteration < itermax)) {
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-    double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
-    if (t + dt > tmax) {
-      dt = tmax - t;
-    }
-
-    std::cout.setf(std::cout.scientific);
-    std::cout << "iteration " << rang::fg::cyan << std::setw(4) << iteration << rang::style::reset
-              << " time=" << rang::fg::green << t << rang::style::reset << " dt=" << rang::fgB::blue << dt
-              << rang::style::reset << '\n';
-
-    mesh = acoustic_solver.computeNextStep(dt, unknowns);
-
-    block_eos.updatePandCFromRhoE();
-
-    t += dt;
-    ++iteration;
-  }
-  std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after "
-            << rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations"
-            << rang::style::reset << '\n';
-}
-
-template AcousticSolverAlgorithm<1>::AcousticSolverAlgorithm(
-  const AcousticSolverType&,
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
-
-template AcousticSolverAlgorithm<2>::AcousticSolverAlgorithm(
-  const AcousticSolverType&,
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
-
-template AcousticSolverAlgorithm<3>::AcousticSolverAlgorithm(
-  const AcousticSolverType&,
-  std::shared_ptr<const IMesh>,
-  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&,
-  const FunctionSymbolId&);
diff --git a/src/language/algorithms/AcousticSolverAlgorithm.hpp b/src/language/algorithms/AcousticSolverAlgorithm.hpp
deleted file mode 100644
index a6973480c4cd913c6463d39012f6eea906540003..0000000000000000000000000000000000000000
--- a/src/language/algorithms/AcousticSolverAlgorithm.hpp
+++ /dev/null
@@ -1,20 +0,0 @@
-#ifndef ACOUSTIC_SOLVER_ALGORITHM_HPP
-#define ACOUSTIC_SOLVER_ALGORITHM_HPP
-
-#include <language/utils/FunctionSymbolId.hpp>
-#include <mesh/IMesh.hpp>
-#include <scheme/AcousticSolverType.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
-
-template <size_t Dimension>
-struct AcousticSolverAlgorithm
-{
-  AcousticSolverAlgorithm(const AcousticSolverType& solver_type,
-                          std::shared_ptr<const IMesh> i_mesh,
-                          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-                          const FunctionSymbolId& rho_id,
-                          const FunctionSymbolId& u_id,
-                          const FunctionSymbolId& p_id);
-};
-
-#endif   // ACOUSTIC_SOLVER_ALGORITHM_HPP
diff --git a/src/language/algorithms/CMakeLists.txt b/src/language/algorithms/CMakeLists.txt
index a8599f575058185b5ceeed0814b242ae6ec70b0a..63c2374987289e7e5c7352dbc773d77ef4e675fb 100644
--- a/src/language/algorithms/CMakeLists.txt
+++ b/src/language/algorithms/CMakeLists.txt
@@ -1,8 +1,8 @@
 # ------------------- Source files --------------------
 
 add_library(PugsLanguageAlgorithms
-  AcousticSolverAlgorithm.cpp
-)
+  INTERFACE # THIS SHOULD BE REMOVED IF FILES ARE FINALY PROVIDED
+  )
 
 
 add_dependencies(PugsLanguageAlgorithms
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 212464f79ce7ccdf24c0164aab84fe5c9b527c31..3051b6bb2265754ed8f1a84b5a1952c59325705b 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -1,6 +1,5 @@
 #include <language/modules/SchemeModule.hpp>
 
-#include <language/algorithms/AcousticSolverAlgorithm.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/Mesh.hpp>
@@ -113,100 +112,6 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction("glace",
-                            std::make_shared<BuiltinFunctionEmbedder<
-                              void(std::shared_ptr<const IMesh>,
-                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
-
-                              [](std::shared_ptr<const IMesh> p_mesh,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id,
-                                 const FunctionSymbolId& p_id) -> void {
-                                switch (p_mesh->dimension()) {
-                                case 1: {
-                                  AcousticSolverAlgorithm<1>{AcousticSolverType::Glace,
-                                                             p_mesh,
-                                                             bc_descriptor_list,
-                                                             rho_id,
-                                                             u_id,
-                                                             p_id};
-                                  break;
-                                }
-                                case 2: {
-                                  AcousticSolverAlgorithm<2>{AcousticSolverType::Glace,
-                                                             p_mesh,
-                                                             bc_descriptor_list,
-                                                             rho_id,
-                                                             u_id,
-                                                             p_id};
-                                  break;
-                                }
-                                case 3: {
-                                  AcousticSolverAlgorithm<3>{AcousticSolverType::Glace,
-                                                             p_mesh,
-                                                             bc_descriptor_list,
-                                                             rho_id,
-                                                             u_id,
-                                                             p_id};
-                                  break;
-                                }
-                                default: {
-                                  throw UnexpectedError("invalid mesh dimension");
-                                }
-                                }
-                              }
-
-                              ));
-
-  this->_addBuiltinFunction("eucclhyd",
-                            std::make_shared<BuiltinFunctionEmbedder<
-                              void(std::shared_ptr<const IMesh>,
-                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
-
-                              [](std::shared_ptr<const IMesh> p_mesh,
-                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
-                                   bc_descriptor_list,
-                                 const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id,
-                                 const FunctionSymbolId& p_id) -> void {
-                                switch (p_mesh->dimension()) {
-                                case 1: {
-                                  AcousticSolverAlgorithm<1>{AcousticSolverType::Eucclhyd,
-                                                             p_mesh,
-                                                             bc_descriptor_list,
-                                                             rho_id,
-                                                             u_id,
-                                                             p_id};
-                                  break;
-                                }
-                                case 2: {
-                                  AcousticSolverAlgorithm<2>{AcousticSolverType::Eucclhyd,
-                                                             p_mesh,
-                                                             bc_descriptor_list,
-                                                             rho_id,
-                                                             u_id,
-                                                             p_id};
-                                  break;
-                                }
-                                case 3: {
-                                  AcousticSolverAlgorithm<3>{AcousticSolverType::Eucclhyd,
-                                                             p_mesh,
-                                                             bc_descriptor_list,
-                                                             rho_id,
-                                                             u_id,
-                                                             p_id};
-                                  break;
-                                }
-                                default: {
-                                  throw UnexpectedError("invalid mesh dimension");
-                                }
-                                }
-                              }
-
-                              ));
-
   this
     ->_addBuiltinFunction("perfect_gas_epsilon_from_rho_p_gamma",
                           std::make_shared<BuiltinFunctionEmbedder<
@@ -283,31 +188,77 @@ SchemeModule::SchemeModule()
 
                             ));
 
-  this->_addBuiltinFunction(
-    "glace_solver",
-    std::make_shared<BuiltinFunctionEmbedder<std::tuple<
-      std::shared_ptr<const IMesh>,
-      std::shared_ptr<const IDiscreteFunction>,
-      std::shared_ptr<const IDiscreteFunction>,
-      std::shared_ptr<const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&,
-                                                const std::shared_ptr<const IDiscreteFunction>&,
-                                                const std::shared_ptr<const IDiscreteFunction>&,
-                                                const std::shared_ptr<const IDiscreteFunction>&,
-                                                const std::shared_ptr<const IDiscreteFunction>&,
-                                                const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-                                                const double&)>>(
-
-      [](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>& c,
-         const std::shared_ptr<const IDiscreteFunction>& p,
-         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, const double& dt)
-        -> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
-                      std::shared_ptr<const IDiscreteFunction>,
-                      std::shared_ptr<const IDiscreteFunction>> {
-        return AcousticSolverHandler{rho, c, u, p, bc_descriptor_list}.solver().apply(dt, rho, u, E);
-      }
+  this->_addBuiltinFunction("glace_solver",
+                            std::make_shared<BuiltinFunctionEmbedder<std::tuple<
+                              std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
+                              std::shared_ptr<const IDiscreteFunction>,
+                              std::shared_ptr<const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::vector<std::shared_ptr<
+                                                                          const IBoundaryConditionDescriptor>>&,
+                                                                        const double&)>>(
+
+                              [](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>& c,
+                                 const std::shared_ptr<const IDiscreteFunction>& p,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt)
+                                -> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
+                                              std::shared_ptr<const IDiscreteFunction>,
+                                              std::shared_ptr<const IDiscreteFunction>> {
+                                return AcousticSolverHandler{AcousticSolverHandler::SolverType::Glace,
+                                                             rho,
+                                                             c,
+                                                             u,
+                                                             p,
+                                                             bc_descriptor_list}
+                                  .solver()
+                                  .apply(dt, rho, u, E);
+                              }
 
-      ));
+                              ));
+
+  this->_addBuiltinFunction("eucclhyd_solver",
+                            std::make_shared<BuiltinFunctionEmbedder<std::tuple<
+                              std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
+                              std::shared_ptr<const IDiscreteFunction>,
+                              std::shared_ptr<const IDiscreteFunction>>(const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::shared_ptr<const IDiscreteFunction>&,
+                                                                        const std::vector<std::shared_ptr<
+                                                                          const IBoundaryConditionDescriptor>>&,
+                                                                        const double&)>>(
+
+                              [](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>& c,
+                                 const std::shared_ptr<const IDiscreteFunction>& p,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt)
+                                -> std::tuple<std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,
+                                              std::shared_ptr<const IDiscreteFunction>,
+                                              std::shared_ptr<const IDiscreteFunction>> {
+                                return AcousticSolverHandler{AcousticSolverHandler::SolverType::Eucclhyd,
+                                                             rho,
+                                                             c,
+                                                             u,
+                                                             p,
+                                                             bc_descriptor_list}
+                                  .solver()
+                                  .apply(dt, rho, u, E);
+                              }
+
+                              ));
 
   this
     ->_addBuiltinFunction("lagrangian",
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 05a8adaf5f212e548e4a286ecd36f7295ad2b13b..76017501f2c26dc690b8c815a02ea4c0f35fa0f8 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -76,6 +76,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
 
   using BoundaryConditionList = std::vector<BoundaryCondition>;
 
+  const SolverType m_solver_type;
   BoundaryConditionList m_boundary_condition_list;
 
   NodeValue<const Rd> m_ur;
@@ -117,6 +118,74 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
     return Ajr;
   }
 
+  NodeValuePerCell<const Rdxd>
+  _computeEucclhydAjr(const MeshType& mesh, const CellValue<const double>& rhoc)
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+    const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
+
+    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+    const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
+
+    NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
+
+    parallel_for(
+      Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Ajr[jr] = zero; });
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        const auto& cell_faces = cell_to_face_matrix[j];
+
+        const double& rho_c = rhoc[j];
+
+        for (size_t L = 0; L < cell_faces.size(); ++L) {
+          const FaceId& l        = cell_faces[L];
+          const auto& face_nodes = face_to_node_matrix[l];
+
+          auto local_node_number_in_cell = [&](NodeId node_number) {
+            for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+              if (node_number == cell_nodes[i_node]) {
+                return i_node;
+              }
+            }
+            return std::numeric_limits<size_t>::max();
+          };
+
+          for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+            const size_t R = local_node_number_in_cell(face_nodes[rl]);
+            Ajr(j, R) += tensorProduct(rho_c * Nlr(l, rl), nlr(l, rl));
+          }
+        }
+      });
+
+    return Ajr;
+  }
+
+  NodeValuePerCell<const Rdxd>
+  _computeAjr(const MeshType& mesh, const CellValue<const double>& rhoc)
+  {
+    if constexpr (Dimension == 1) {
+      return _computeGlaceAjr(mesh, rhoc);
+    } else {
+      switch (m_solver_type) {
+      case SolverType::Glace: {
+        return _computeGlaceAjr(mesh, rhoc);
+      }
+      case SolverType::Eucclhyd: {
+        return _computeEucclhydAjr(mesh, rhoc);
+      }
+      default: {
+        throw UnexpectedError("invalid solver type");
+      }
+      }
+    }
+  }
+
   NodeValue<Rdxd>
   _computeAr(const MeshType& mesh, const NodeValuePerCell<const Rdxd>& Ajr)
   {
@@ -322,18 +391,19 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
     return F;
   }
 
-  AcousticSolver(const std::shared_ptr<const MeshType>& p_mesh,
+  AcousticSolver(SolverType solver_type,
+                 const std::shared_ptr<const MeshType>& p_mesh,
                  const DiscreteScalarFunction& rho,
                  const DiscreteScalarFunction& c,
                  const DiscreteVectorFunction& u,
                  const DiscreteScalarFunction& p,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-    : m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)}
+    : m_solver_type{solver_type}, m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)}
   {
     const MeshType& mesh = *p_mesh;
 
     CellValue<const double> rhoc     = this->_getRhoC(rho, c);
-    NodeValuePerCell<const Rdxd> Ajr = this->_computeGlaceAjr(mesh, rhoc);
+    NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(mesh, rhoc);
 
     NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
     NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, p);
@@ -422,13 +492,15 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
                        *std::dynamic_pointer_cast<const DiscreteScalarFunction>(E));
   }
 
-  AcousticSolver(const std::shared_ptr<const IMesh>& mesh,
+  AcousticSolver(SolverType solver_type,
+                 const std::shared_ptr<const IMesh>& mesh,
                  const std::shared_ptr<const IDiscreteFunction>& rho,
                  const std::shared_ptr<const IDiscreteFunction>& c,
                  const std::shared_ptr<const IDiscreteFunction>& u,
                  const std::shared_ptr<const IDiscreteFunction>& p,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-    : AcousticSolver{std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(mesh),
+    : AcousticSolver{solver_type,
+                     std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(mesh),
                      *std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho),
                      *std::dynamic_pointer_cast<const DiscreteScalarFunction>(c),
                      *std::dynamic_pointer_cast<const DiscreteVectorFunction>(u),
@@ -678,6 +750,7 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::SymmetryBoundaryConditio
 };
 
 AcousticSolverHandler::AcousticSolverHandler(
+  SolverType solver_type,
   const std::shared_ptr<const IDiscreteFunction>& rho,
   const std::shared_ptr<const IDiscreteFunction>& c,
   const std::shared_ptr<const IDiscreteFunction>& u,
@@ -695,15 +768,15 @@ AcousticSolverHandler::AcousticSolverHandler(
 
   switch (i_mesh->dimension()) {
   case 1: {
-    m_acoustic_solver = std::make_unique<AcousticSolver<1>>(i_mesh, rho, c, u, p, bc_descriptor_list);
+    m_acoustic_solver = std::make_unique<AcousticSolver<1>>(solver_type, i_mesh, rho, c, u, p, bc_descriptor_list);
     break;
   }
   case 2: {
-    m_acoustic_solver = std::make_unique<AcousticSolver<2>>(i_mesh, rho, c, u, p, bc_descriptor_list);
+    m_acoustic_solver = std::make_unique<AcousticSolver<2>>(solver_type, i_mesh, rho, c, u, p, bc_descriptor_list);
     break;
   }
   case 3: {
-    m_acoustic_solver = std::make_unique<AcousticSolver<3>>(i_mesh, rho, c, u, p, bc_descriptor_list);
+    m_acoustic_solver = std::make_unique<AcousticSolver<3>>(solver_type, i_mesh, rho, c, u, p, bc_descriptor_list);
     break;
   }
   default: {
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 6618361cc56a6efe55620fda83c8976ecd60b3f1..99e8a8617d8a7606deb56b56095b34f75b0a9629 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -28,6 +28,13 @@ double acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c);
 
 class AcousticSolverHandler
 {
+ public:
+  enum class SolverType
+  {
+    Glace,
+    Eucclhyd
+  };
+
  private:
   struct IAcousticSolver
   {
@@ -59,599 +66,12 @@ class AcousticSolverHandler
     return *m_acoustic_solver;
   }
 
-  AcousticSolverHandler(const std::shared_ptr<const IDiscreteFunction>& rho,
+  AcousticSolverHandler(SolverType solver_type,
+                        const std::shared_ptr<const IDiscreteFunction>& rho,
                         const std::shared_ptr<const IDiscreteFunction>& c,
                         const std::shared_ptr<const IDiscreteFunction>& u,
                         const std::shared_ptr<const IDiscreteFunction>& p,
                         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list);
 };
 
-template <typename MeshType>
-class AcousticSolverOld
-{
- public:
-  class PressureBoundaryCondition;
-  class SymmetryBoundaryCondition;
-  class VelocityBoundaryCondition;
-
-  using BoundaryCondition =
-    std::variant<PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
-
-  using BoundaryConditionList = std::vector<BoundaryCondition>;
-
-  constexpr static size_t Dimension = MeshType::Dimension;
-
- private:
-  using MeshDataType = MeshData<Dimension>;
-  using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>;
-
-  std::shared_ptr<const MeshType> m_mesh;
-  const typename MeshType::Connectivity& m_connectivity;
-
-  using Rd  = TinyVector<Dimension>;
-  using Rdd = TinyMatrix<Dimension>;
-
-  const BoundaryConditionList m_boundary_condition_list;
-
-  const AcousticSolverType m_solver_type;
-
-  void
-  _applyPressureBC()
-  {
-    for (const auto& boundary_condition : m_boundary_condition_list) {
-      std::visit(
-        [&](auto&& bc) {
-          using T = std::decay_t<decltype(bc)>;
-          if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
-            MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-            if constexpr (Dimension == 1) {
-              const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
-
-              const auto& node_to_cell_matrix               = m_connectivity.nodeToCellMatrix();
-              const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
-
-              const auto& node_list  = bc.faceList();
-              const auto& value_list = bc.valueList();
-              parallel_for(
-                node_list.size(), PUGS_LAMBDA(size_t i_node) {
-                  const NodeId node_id       = node_list[i_node];
-                  const auto& node_cell_list = node_to_cell_matrix[node_id];
-                  Assert(node_cell_list.size() == 1);
-
-                  CellId node_cell_id              = node_cell_list[0];
-                  size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
-
-                  m_br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
-                });
-            } else {
-              const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
-
-              const auto& face_to_cell_matrix               = m_connectivity.faceToCellMatrix();
-              const auto& face_to_node_matrix               = m_connectivity.faceToNodeMatrix();
-              const auto& face_local_numbers_in_their_cells = m_connectivity.faceLocalNumbersInTheirCells();
-              const auto& face_cell_is_reversed             = m_connectivity.cellFaceIsReversed();
-
-              const auto& face_list  = bc.faceList();
-              const auto& value_list = bc.valueList();
-              for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                const FaceId face_id       = face_list[i_face];
-                const auto& face_cell_list = face_to_cell_matrix[face_id];
-                Assert(face_cell_list.size() == 1);
-
-                CellId face_cell_id              = face_cell_list[0];
-                size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
-
-                const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
-
-                const auto& face_nodes = face_to_node_matrix[face_id];
-
-                for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
-                  NodeId node_id = face_nodes[i_node];
-                  m_br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
-                }
-              }
-            }
-          }
-        },
-        boundary_condition);
-    }
-  }
-
-  void
-  _applySymmetryBC()
-  {
-    for (const auto& boundary_condition : m_boundary_condition_list) {
-      std::visit(
-        [&](auto&& bc) {
-          using T = std::decay_t<decltype(bc)>;
-          if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
-            const Rd& n = bc.outgoingNormal();
-
-            const Rdd I   = identity;
-            const Rdd nxn = tensorProduct(n, n);
-            const Rdd P   = I - nxn;
-
-            const Array<const NodeId>& node_list = bc.nodeList();
-            parallel_for(
-              bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
-                const NodeId r = node_list[r_number];
-
-                m_Ar[r] = P * m_Ar[r] * P + nxn;
-                m_br[r] = P * m_br[r];
-              });
-          }
-        },
-        boundary_condition);
-    }
-  }
-
-  void
-  _applyVelocityBC()
-  {
-    for (const auto& boundary_condition : m_boundary_condition_list) {
-      std::visit(
-        [&](auto&& bc) {
-          using T = std::decay_t<decltype(bc)>;
-          if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
-            const auto& node_list  = bc.nodeList();
-            const auto& value_list = bc.valueList();
-
-            parallel_for(
-              node_list.size(), PUGS_LAMBDA(size_t i_node) {
-                NodeId node_id    = node_list[i_node];
-                const auto& value = value_list[i_node];
-
-                m_Ar[node_id] = identity;
-                m_br[node_id] = value;
-              });
-          }
-        },
-        boundary_condition);
-    }
-  }
-
-  PUGS_INLINE
-  const CellValue<const double>
-  computeRhoCj(const CellValue<const double>& rhoj, const CellValue<const double>& cj)
-  {
-    parallel_for(
-      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { m_rhocj[j] = rhoj[j] * cj[j]; });
-    return m_rhocj;
-  }
-
-  PUGS_INLINE
-  void
-  _computeGlaceAjr(const CellValue<const double>& rhocj)
-  {
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
-    const NodeValuePerCell<const Rd> njr = mesh_data.njr();
-
-    parallel_for(
-      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
-        const size_t& nb_nodes = m_Ajr.numberOfSubValues(j);
-        const double& rho_c    = rhocj[j];
-        for (size_t r = 0; r < nb_nodes; ++r) {
-          m_Ajr(j, r) = tensorProduct(rho_c * Cjr(j, r), njr(j, r));
-        }
-      });
-  }
-
-  PUGS_INLINE
-  void
-  _computeEucclhydAjr(const CellValue<const double>& rhocj)
-  {
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-    const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
-    const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
-
-    const auto& face_to_node_matrix = m_connectivity.faceToNodeMatrix();
-    const auto& cell_to_node_matrix = m_connectivity.cellToNodeMatrix();
-    const auto& cell_to_face_matrix = m_connectivity.cellToFaceMatrix();
-
-    parallel_for(
-      m_Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { m_Ajr[jr] = zero; });
-
-    parallel_for(
-      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
-        const auto& cell_nodes = cell_to_node_matrix[j];
-
-        const auto& cell_faces = cell_to_face_matrix[j];
-
-        const double& rho_c = rhocj[j];
-
-        for (size_t L = 0; L < cell_faces.size(); ++L) {
-          const FaceId& l        = cell_faces[L];
-          const auto& face_nodes = face_to_node_matrix[l];
-
-          auto local_node_number_in_cell = [&](NodeId node_number) {
-            for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
-              if (node_number == cell_nodes[i_node]) {
-                return i_node;
-              }
-            }
-            return std::numeric_limits<size_t>::max();
-          };
-
-          for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
-            const size_t R = local_node_number_in_cell(face_nodes[rl]);
-            m_Ajr(j, R) += tensorProduct(rho_c * Nlr(l, rl), nlr(l, rl));
-          }
-        }
-      });
-  }
-
-  PUGS_INLINE
-  void
-  computeAjr(const CellValue<const double>& rhocj)
-  {
-    switch (m_solver_type) {
-    case AcousticSolverType::Glace: {
-      this->_computeGlaceAjr(rhocj);
-      break;
-    }
-    case AcousticSolverType::Eucclhyd: {
-      if constexpr (Dimension > 1) {
-        this->_computeEucclhydAjr(rhocj);
-      } else {
-        this->_computeGlaceAjr(rhocj);
-      }
-      break;
-    }
-    }
-  }
-
-  PUGS_INLINE
-  const NodeValue<const Rdd>
-  computeAr(const NodeValuePerCell<const Rdd>& Ajr)
-  {
-    const auto& node_to_cell_matrix               = m_connectivity.nodeToCellMatrix();
-    const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
-
-    parallel_for(
-      m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-        Rdd sum                                    = zero;
-        const auto& node_to_cell                   = node_to_cell_matrix[r];
-        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemValues(r);
-
-        for (size_t j = 0; j < node_to_cell.size(); ++j) {
-          const CellId J       = node_to_cell[j];
-          const unsigned int R = node_local_number_in_its_cells[j];
-          sum += Ajr(J, R);
-        }
-        m_Ar[r] = sum;
-      });
-
-    return m_Ar;
-  }
-
-  PUGS_INLINE
-  const NodeValue<const Rd>
-  computeBr(const CellValue<const Rd>& uj, const CellValue<const double>& pj)
-  {
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-    const NodeValuePerCell<const Rd>& Cjr  = mesh_data.Cjr();
-    const NodeValuePerCell<const Rdd>& Ajr = m_Ajr;
-
-    const auto& node_to_cell_matrix               = m_connectivity.nodeToCellMatrix();
-    const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
-
-    parallel_for(
-      m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) {
-        Rd& br                                     = m_br[r];
-        br                                         = zero;
-        const auto& node_to_cell                   = node_to_cell_matrix[r];
-        const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemValues(r);
-        for (size_t j = 0; j < node_to_cell.size(); ++j) {
-          const CellId J       = node_to_cell[j];
-          const unsigned int R = node_local_number_in_its_cells[j];
-          br += Ajr(J, R) * uj[J] + pj[J] * Cjr(J, R);
-        }
-      });
-
-    return m_br;
-  }
-
-  void
-  applyBoundaryConditions()
-  {
-    this->_applyPressureBC();
-    this->_applySymmetryBC();
-    this->_applyVelocityBC();
-  }
-
-  void
-  computeUr()
-  {
-    const NodeValue<const Rdd> Ar = m_Ar;
-    const NodeValue<const Rd> br  = m_br;
-
-    inverse(Ar, m_inv_Ar);
-    const NodeValue<const Rdd> invAr = m_inv_Ar;
-    parallel_for(
-      m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { m_ur[r] = invAr[r] * br[r]; });
-  }
-
-  void
-  computeFjr(const CellValue<const Rd>& uj, const CellValue<const double>& pj)
-  {
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-    const NodeValuePerCell<const Rd> Cjr  = mesh_data.Cjr();
-    const NodeValuePerCell<const Rdd> Ajr = m_Ajr;
-
-    const NodeValue<const Rd> ur = m_ur;
-
-    const auto& cell_to_node_matrix = m_mesh->connectivity().cellToNodeMatrix();
-
-    parallel_for(
-      m_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) {
-          m_Fjr(j, r) = Ajr(j, r) * (uj[j] - ur[cell_nodes[r]]) + pj[j] * Cjr(j, r);
-        }
-      });
-  }
-
-  void
-  inverse(const NodeValue<const Rdd>& A, NodeValue<Rdd>& inv_A) const
-  {
-    parallel_for(
-      m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { inv_A[r] = ::inverse(A[r]); });
-  }
-
-  PUGS_INLINE
-  void
-  computeExplicitFluxes(const CellValue<const double>& rhoj,
-                        const CellValue<const Rd>& uj,
-                        const CellValue<const double>& pj,
-                        const CellValue<const double>& cj)
-  {
-    const CellValue<const double> rhocj = computeRhoCj(rhoj, cj);
-    computeAjr(rhocj);
-
-    NodeValuePerCell<const Rdd> Ajr = m_Ajr;
-    this->computeAr(Ajr);
-    this->computeBr(uj, pj);
-
-    this->applyBoundaryConditions();
-
-    synchronize(m_Ar);
-    synchronize(m_br);
-
-    computeUr();
-    computeFjr(uj, pj);
-  }
-
-  NodeValue<Rd> m_br;
-  NodeValuePerCell<Rdd> m_Ajr;
-  NodeValue<Rdd> m_Ar;
-  NodeValue<Rdd> m_inv_Ar;
-  NodeValuePerCell<Rd> m_Fjr;
-  NodeValue<Rd> m_ur;
-  CellValue<double> m_rhocj;
-  CellValue<double> m_Vj_over_cj;
-
- public:
-  AcousticSolverOld(std::shared_ptr<const MeshType> p_mesh,
-                    const BoundaryConditionList& bc_list,
-                    const AcousticSolverType solver_type)
-    : m_mesh(p_mesh),
-      m_connectivity(m_mesh->connectivity()),
-      m_boundary_condition_list(bc_list),
-      m_solver_type(solver_type),
-      m_br(m_connectivity),
-      m_Ajr(m_connectivity),
-      m_Ar(m_connectivity),
-      m_inv_Ar(m_connectivity),
-      m_Fjr(m_connectivity),
-      m_ur(m_connectivity),
-      m_rhocj(m_connectivity),
-      m_Vj_over_cj(m_connectivity)
-  {
-    ;
-  }
-
-  PUGS_INLINE
-  double
-  acoustic_dt(const CellValue<const double>& Vj, const CellValue<const double>& cj) const
-  {
-    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-    const NodeValuePerCell<const double>& ljr = mesh_data.ljr();
-    const auto& cell_to_node_matrix           = m_mesh->connectivity().cellToNodeMatrix();
-
-    parallel_for(
-      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
-        const auto& cell_nodes = cell_to_node_matrix[j];
-
-        double S = 0;
-        for (size_t r = 0; r < cell_nodes.size(); ++r) {
-          S += ljr(j, r);
-        }
-        m_Vj_over_cj[j] = 2 * Vj[j] / (S * cj[j]);
-      });
-
-    return min(m_Vj_over_cj);
-  }
-
-  [[nodiscard]] std::shared_ptr<const MeshType>
-  computeNextStep(double dt, UnknownsType& unknowns)
-  {
-    CellValue<double>& rhoj = unknowns.rhoj();
-    CellValue<Rd>& uj       = unknowns.uj();
-    CellValue<double>& Ej   = unknowns.Ej();
-
-    CellValue<double>& ej = unknowns.ej();
-    CellValue<double>& pj = unknowns.pj();
-    CellValue<double>& cj = unknowns.cj();
-
-    MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-    const CellValue<const double> Vj = mesh_data.Vj();
-
-    computeExplicitFluxes(rhoj, uj, pj, cj);
-
-    const NodeValuePerCell<Rd>& Fjr = m_Fjr;
-    const NodeValue<const Rd> ur    = m_ur;
-    const auto& cell_to_node_matrix = m_mesh->connectivity().cellToNodeMatrix();
-
-    const CellValue<const double> inv_mj = unknowns.invMj();
-    parallel_for(
-      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
-        const auto& cell_nodes = cell_to_node_matrix[j];
-
-        Rd momentum_fluxes   = zero;
-        double energy_fluxes = 0;
-        for (size_t R = 0; R < cell_nodes.size(); ++R) {
-          const NodeId r = cell_nodes[R];
-          momentum_fluxes += Fjr(j, R);
-          energy_fluxes += (Fjr(j, R), ur[r]);
-        }
-        uj[j] -= (dt * inv_mj[j]) * momentum_fluxes;
-        Ej[j] -= (dt * inv_mj[j]) * energy_fluxes;
-      });
-
-    parallel_for(
-      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { ej[j] = Ej[j] - 0.5 * (uj[j], uj[j]); });
-
-    NodeValue<Rd> new_xr = copy(m_mesh->xr());
-    parallel_for(
-      m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
-
-    m_mesh                         = std::make_shared<MeshType>(m_mesh->shared_connectivity(), new_xr);
-    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*m_mesh).Vj();
-
-    const CellValue<const double> mj = unknowns.mj();
-    parallel_for(
-      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { rhoj[j] = mj[j] / new_Vj[j]; });
-
-    return m_mesh;
-  }
-};
-
-template <typename MeshType>
-class AcousticSolverOld<MeshType>::PressureBoundaryCondition
-{
- private:
-  const Array<const double> m_value_list;
-  const Array<const FaceId> m_face_list;
-
- public:
-  const Array<const FaceId>&
-  faceList() const
-  {
-    return m_face_list;
-  }
-
-  const Array<const double>&
-  valueList() const
-  {
-    return m_value_list;
-  }
-
-  PressureBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
-    : m_value_list{value_list}, m_face_list{face_list}
-  {}
-
-  ~PressureBoundaryCondition() = default;
-};
-
-template <>
-class AcousticSolverOld<Mesh<Connectivity<1>>>::PressureBoundaryCondition
-{
- private:
-  const Array<const double> m_value_list;
-  const Array<const NodeId> m_face_list;
-
- public:
-  const Array<const NodeId>&
-  faceList() const
-  {
-    return m_face_list;
-  }
-
-  const Array<const double>&
-  valueList() const
-  {
-    return m_value_list;
-  }
-
-  PressureBoundaryCondition(const Array<const NodeId>& face_list, const Array<const double>& value_list)
-    : m_value_list{value_list}, m_face_list{face_list}
-  {}
-
-  ~PressureBoundaryCondition() = default;
-};
-
-template <typename MeshType>
-class AcousticSolverOld<MeshType>::VelocityBoundaryCondition
-{
- private:
-  const Array<const TinyVector<Dimension>> m_value_list;
-  const Array<const NodeId> m_node_list;
-
- public:
-  const Array<const NodeId>&
-  nodeList() const
-  {
-    return m_node_list;
-  }
-
-  const Array<const TinyVector<Dimension>>&
-  valueList() const
-  {
-    return m_value_list;
-  }
-
-  VelocityBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list)
-    : m_value_list{value_list}, m_node_list{node_list}
-  {}
-
-  ~VelocityBoundaryCondition() = default;
-};
-
-template <typename MeshType>
-class AcousticSolverOld<MeshType>::SymmetryBoundaryCondition
-{
- public:
-  static constexpr size_t Dimension = MeshType::Dimension;
-
-  using Rd = TinyVector<Dimension, double>;
-
- private:
-  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
-
- public:
-  const Rd&
-  outgoingNormal() const
-  {
-    return m_mesh_flat_node_boundary.outgoingNormal();
-  }
-
-  size_t
-  numberOfNodes() const
-  {
-    return m_mesh_flat_node_boundary.nodeList().size();
-  }
-
-  const Array<const NodeId>&
-  nodeList() const
-  {
-    return m_mesh_flat_node_boundary.nodeList();
-  }
-
-  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
-    : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
-  {
-    ;
-  }
-
-  ~SymmetryBoundaryCondition() = default;
-};
-
 #endif   // ACOUSTIC_SOLVER_HPP