diff --git a/.gitlab-ci/clang10-mpi-debug.yml b/.gitlab-ci/clang10-mpi-debug.yml
index ceafce232dcaf5935f1c6be403c796364a3b8c16..a757e02d46b81e60e981c0d3d11aaa8058c6fcab 100644
--- a/.gitlab-ci/clang10-mpi-debug.yml
+++ b/.gitlab-ci/clang10-mpi-debug.yml
@@ -6,7 +6,7 @@ test:clang10-mpi-debug:
     - mkdir -p build/clang10-debug-mpi
     - cd build/clang10-debug-mpi
     - CXX=clang++-10 CC=clang-10 cmake ../.. -DCMAKE_BUILD_TYPE=Debug
-    - make
+    - make -j 4
     - make check
   cache:
     key: "${CI_COMMIT_REF_SLUG}-clang10-debug-mpi"
diff --git a/.gitlab-ci/clang10-mpi-release.yml b/.gitlab-ci/clang10-mpi-release.yml
index 58a2d2d2e3878ffd8672eb71e15b34a5d58062a6..3172c1c2c5f869e1e44153956cfde696adb6b05d 100644
--- a/.gitlab-ci/clang10-mpi-release.yml
+++ b/.gitlab-ci/clang10-mpi-release.yml
@@ -6,7 +6,7 @@ test:clang10-mpi-release:
     - mkdir -p build/clang10-release-mpi
     - cd build/clang10-release-mpi
     - CXX=clang++-10 CC=clang-10 cmake ../.. -DCMAKE_BUILD_TYPE=Release -DCLANG_FORMAT=/usr/bin/clang-format-10
-    - make
+    - make -j 4
     - make check
   cache:
     key: "${CI_COMMIT_REF_SLUG}-clang10-release-mpi"
diff --git a/.gitlab-ci/gcc10-mpi-coverage.yml b/.gitlab-ci/gcc10-mpi-coverage.yml
index d2cc109d71cdada957ae5ce65c72ede6ce3b21ea..5953c48fb68b696a466522e81d44e6ae9ecdcd86 100644
--- a/.gitlab-ci/gcc10-mpi-coverage.yml
+++ b/.gitlab-ci/gcc10-mpi-coverage.yml
@@ -6,7 +6,7 @@ coverage:gcc10-mpi-coverage:
     - mkdir -p build/gcc10-cov-mpi
     - cd build/gcc10-cov-mpi
     - CXX=g++-10 CC=gcc-10 cmake ../.. -DCMAKE_BUILD_TYPE=Coverage
-    - make
+    - make -j 4
   cache:
     key: "${CI_COMMIT_REF_SLUG}-gcc10-cov-mpi"
     paths:
diff --git a/.gitlab-ci/gcc10-mpi-release.yml b/.gitlab-ci/gcc10-mpi-release.yml
index 5e8277ad61a089d8baa72e5cc7e8d2b7219ad775..198fe1785140889fb65d54e3fbdc4f9fe8ba94af 100644
--- a/.gitlab-ci/gcc10-mpi-release.yml
+++ b/.gitlab-ci/gcc10-mpi-release.yml
@@ -6,7 +6,7 @@ test:gcc10-mpi-release:
     - mkdir -p build/gcc10-release-mpi
     - cd build/gcc10-release-mpi
     - CXX=g++-10 CC=gcc-10 cmake ../.. -DCMAKE_BUILD_TYPE=Release
-    - make
+    - make -j 4
     - make test
   cache:
     key: "${CI_COMMIT_REF_SLUG}-gcc10-release-mpi"
diff --git a/.gitlab-ci/gcc10-seq-coverage.yml b/.gitlab-ci/gcc10-seq-coverage.yml
index 5cacacd5d02a0838c36b852076c8cb6ad0573031..81d8e3b453496c73a20bfe4515c08c2faa0f5991 100644
--- a/.gitlab-ci/gcc10-seq-coverage.yml
+++ b/.gitlab-ci/gcc10-seq-coverage.yml
@@ -6,7 +6,7 @@ coverage:gcc10-seq-coverage:
     - mkdir -p build/gcc10-cov
     - cd build/gcc10-cov
     - CXX=g++-10 CC=gcc-10 cmake ../.. -DCMAKE_BUILD_TYPE=Coverage
-    - make
+    - make -j 4
   cache:
     key: "${CI_COMMIT_REF_SLUG}-gcc10-cov"
     paths:
diff --git a/.gitlab-ci/gcc10-seq-release.yml b/.gitlab-ci/gcc10-seq-release.yml
index b4423bbd31f947698a5fa4dc13171e0582fc5d32..8a34f11b0dea3b4df56a2e8883c8d08245ddac2c 100644
--- a/.gitlab-ci/gcc10-seq-release.yml
+++ b/.gitlab-ci/gcc10-seq-release.yml
@@ -6,7 +6,7 @@ test:gcc10-seq-release:
     - mkdir -p build/gcc10-release-seq
     - cd build/gcc10-release-seq
     - CXX=g++-10 CC=gcc-10 cmake ../.. -DCMAKE_BUILD_TYPE=Release
-    - make
+    - make -j 4
     - make check
   cache:
     key: "${CI_COMMIT_REF_SLUG}-gcc10-release-seq"
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 114ceecc5d49685a513708dced3dcffea83bff29..95f96478b331b3790fb57946a031cea15fe73d67 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -367,6 +367,20 @@ add_custom_target(run_unit_tests
 # unit tests coverage
 
 if("${CMAKE_BUILD_TYPE}" STREQUAL "Coverage")
+  file(GLOB_RECURSE GCNO_FILE_LIST RELATIVE "${PUGS_BINARY_DIR}"  "*.gcno")
+  list(FILTER GCNO_FILE_LIST EXCLUDE REGEX "^\.\./.*" )
+  file(GLOB_RECURSE GCDA_FILE_LIST RELATIVE "${PUGS_BINARY_DIR}"  "*.gcda")
+  list(FILTER GCDA_FILE_LIST EXCLUDE REGEX "^\.\./.*" )
+  foreach(COV_INFO_FILE IN LISTS GCNO_FILE_LIST GCDA_FILE_LIST)
+    string(REGEX REPLACE "/CMakeFiles/.*\.dir/" "/" COV_SRC_FILE "${PUGS_SOURCE_DIR}/${COV_INFO_FILE}")
+    string(REGEX REPLACE "\.gcda$" "" COV_SRC_FILE "${COV_SRC_FILE}")
+    string(REGEX REPLACE "\.gcno$" "" COV_SRC_FILE "${COV_SRC_FILE}")
+    if (NOT EXISTS "${COV_SRC_FILE}")
+      file(REMOVE "${PUGS_BINARY_DIR}/${COV_INFO_FILE}")
+      message(STATUS "removed file ${COV_INFO_FILE}: no longer needed")
+    endif()
+  endforeach()
+
   find_program(LCOV lcov)
   if(NOT LCOV)
     message(FATAL_ERROR "lcov not found, cannot perform coverage.")
diff --git a/src/language/algorithms/AcousticSolverAlgorithm.cpp b/src/language/algorithms/AcousticSolverAlgorithm.cpp
deleted file mode 100644
index edce9e025f0feb4cd39acaadcb537451c9251065..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 AcousticSolver<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 AcousticSolver<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 AcousticSolver<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 AcousticSolver<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);
-
-  AcousticSolver 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 c23d3794e4c26a95301599b6933fe0d643391010..3051b6bb2265754ed8f1a84b5a1952c59325705b 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -1,12 +1,13 @@
 #include <language/modules/SchemeModule.hpp>
 
-#include <language/algorithms/AcousticSolverAlgorithm.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/Mesh.hpp>
+#include <scheme/AcousticSolver.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionDescriptorP0.hpp>
 #include <scheme/DiscreteFunctionInterpoler.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/FourierBoundaryConditionDescriptor.hpp>
 #include <scheme/FreeBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
@@ -18,6 +19,8 @@
 #include <scheme/NumberedBoundaryDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 
+#include <scheme/PerfectGas.hpp>
+
 #include <memory>
 
 SchemeModule::SchemeModule()
@@ -109,97 +112,171 @@ 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&)>>(
+  this
+    ->_addBuiltinFunction("perfect_gas_epsilon_from_rho_p_gamma",
+                          std::make_shared<BuiltinFunctionEmbedder<
+                            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>& rho,
+                               const std::shared_ptr<const IDiscreteFunction>& p,
+                               const std::shared_ptr<const IDiscreteFunction>& gamma)
+                              -> std::shared_ptr<const IDiscreteFunction> {
+                              return perfect_gas::epsilonFromRhoPGamma(rho, p, gamma);
+                            }
+
+                            ));
 
-                              [](std::shared_ptr<const IMesh> p_mesh,
+  this
+    ->_addBuiltinFunction("perfect_gas_p_from_rho_epsilon_gamma",
+                          std::make_shared<BuiltinFunctionEmbedder<
+                            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>& rho,
+                               const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                               const std::shared_ptr<const IDiscreteFunction>& gamma)
+                              -> std::shared_ptr<const IDiscreteFunction> {
+                              return perfect_gas::pFromRhoEpsilonGamma(rho, epsilon, gamma);
+                            }
+
+                            ));
+
+  this
+    ->_addBuiltinFunction("perfect_gas_sound_speed",
+                          std::make_shared<BuiltinFunctionEmbedder<
+                            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>& rho,
+                               const std::shared_ptr<const IDiscreteFunction>& p,
+                               const std::shared_ptr<const IDiscreteFunction>& gamma)
+                              -> std::shared_ptr<const IDiscreteFunction> {
+                              return perfect_gas::soundSpeed(rho, p, gamma);
+                            }
+
+                            ));
+
+  this
+    ->_addBuiltinFunction("E_from_epsilon_u",
+                          std::make_shared<BuiltinFunctionEmbedder<
+                            std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+                                                                     const std::shared_ptr<const IDiscreteFunction>&)>>(
+
+                            [](const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                               const std::shared_ptr<const IDiscreteFunction>& u)
+                              -> std::shared_ptr<const IDiscreteFunction> {
+                              return perfect_gas::totalEnergyFromEpsilonU(epsilon, u);
+                            }
+
+                            ));
+
+  this
+    ->_addBuiltinFunction("epsilon_from_E_u",
+                          std::make_shared<BuiltinFunctionEmbedder<
+                            std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IDiscreteFunction>&,
+                                                                     const std::shared_ptr<const IDiscreteFunction>&)>>(
+
+                            [](const std::shared_ptr<const IDiscreteFunction>& E,
+                               const std::shared_ptr<const IDiscreteFunction>& u)
+                              -> std::shared_ptr<const IDiscreteFunction> {
+                              return perfect_gas::epsilonFromTotalEnergyU(E, u);
+                            }
+
+                            ));
+
+  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 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");
-                                }
-                                }
+                                 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",
-                            std::make_shared<BuiltinFunctionEmbedder<
-                              void(std::shared_ptr<const IMesh>,
-                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
+  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&)>>(
 
-                              [](std::shared_ptr<const IMesh> p_mesh,
+                              [](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 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");
-                                }
-                                }
+                                 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",
+                          std::make_shared<BuiltinFunctionEmbedder<
+                            std::shared_ptr<const IDiscreteFunction>(const std::shared_ptr<const IMesh>&,
+                                                                     const std::shared_ptr<const IDiscreteFunction>&)>>(
+
+                            [](const std::shared_ptr<const IMesh>& mesh,
+                               const std::shared_ptr<const IDiscreteFunction>& v)
+                              -> std::shared_ptr<const IDiscreteFunction> { return shallowCopy(mesh, v); }
+
+                            ));
+
+  this->_addBuiltinFunction("acoustic_dt",
+                            std::make_shared<
+                              BuiltinFunctionEmbedder<double(const std::shared_ptr<const IDiscreteFunction>&)>>(
+
+                              [](const std::shared_ptr<const IDiscreteFunction>& c) -> double { return acoustic_dt(c); }
+
                               ));
 }
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index c82e9ce06ec8b7d43ce1c6cc2fec282f89dc1fee..99cc8aeeb4ca127e8e60ba3e7ee75fcc93e8db76 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -41,6 +41,7 @@ class MeshData : public IMeshData
   CellValue<const Rd> m_cell_centroid;
   CellValue<const Rd> m_cell_iso_barycenter;
   CellValue<const double> m_Vj;
+  CellValue<const double> m_sum_r_ljr;
   FaceValue<const double> m_ll;
 
   PUGS_INLINE
@@ -435,6 +436,34 @@ class MeshData : public IMeshData
     }
   }
 
+  PUGS_INLINE
+  void
+  _computeSumOverRljr()
+  {
+    CellValue<double> sum_r_ljr(m_mesh.connectivity());
+
+    if constexpr (Dimension == 1) {
+      parallel_for(
+        m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { sum_r_ljr[j] = 2; });
+    } else {
+      auto ljr = this->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 sum             = 0;
+          for (size_t r = 0; r < cell_nodes.size(); ++r) {
+            sum += ljr(j, r);
+          }
+          sum_r_ljr[j] = sum;
+        });
+    }
+
+    m_sum_r_ljr = sum_r_ljr;
+  }
+
   void
   _checkCellVolume()
   {
@@ -563,6 +592,16 @@ class MeshData : public IMeshData
     return m_Vj;
   }
 
+  PUGS_INLINE
+  CellValue<const double>
+  sumOverRLjr()
+  {
+    if (not m_sum_r_ljr.isBuilt()) {
+      this->_computeSumOverRljr();
+    }
+    return m_sum_r_ljr;
+  }
+
  private:
   // MeshData **must** be constructed through MeshDataManager
   friend class MeshDataManager;
diff --git a/src/output/GnuplotWriter.cpp b/src/output/GnuplotWriter.cpp
index 436ae7ffa22d7f3aaa20cbd44e2331393958823c..f503c996c3f86a383f2d12fa95185fd1fe205578 100644
--- a/src/output/GnuplotWriter.cpp
+++ b/src/output/GnuplotWriter.cpp
@@ -102,8 +102,8 @@ GnuplotWriter::_writeCellValue(const CellValue<DataType>& cell_value, CellId cel
 template <typename DataType, ItemType item_type>
 void
 GnuplotWriter::_writeValue(const ItemValue<DataType, item_type>& item_value,
-                           CellId cell_id,
-                           NodeId node_id,
+                           [[maybe_unused]] CellId cell_id,
+                           [[maybe_unused]] NodeId node_id,
                            std::ostream& fout) const
 {
   if constexpr (item_type == ItemType::cell) {
diff --git a/src/output/WriterBase.cpp b/src/output/WriterBase.cpp
index 0359cf7b3d33c8aa9fb33f8b62c671425e5958ba..4a26605cd330301ddf2fc3346c544228c34f4f9b 100644
--- a/src/output/WriterBase.cpp
+++ b/src/output/WriterBase.cpp
@@ -172,7 +172,7 @@ WriterBase::getLastTime() const
 }
 
 WriterBase::WriterBase(const std::string& base_filename, const double& time_period)
-  : m_base_filename{base_filename}, m_time_period{time_period}
+  : m_base_filename{base_filename}, m_time_period{time_period}, m_next_time{0}
 {}
 
 void
@@ -183,7 +183,8 @@ WriterBase::writeIfNeeded(const std::vector<std::shared_ptr<const NamedDiscreteF
   if (time == last_time)
     return;   // output already performed
 
-  if (time >= last_time + m_time_period) {
+  if (time >= m_next_time) {
+    m_next_time += m_time_period;
     this->write(named_discrete_function_list, time);
     m_saved_times.push_back(time);
   }
diff --git a/src/output/WriterBase.hpp b/src/output/WriterBase.hpp
index 6e9a08ad123361178327ddb984e5a3da5bdfa909..babdec90733416cf9007e55c9652f2c380d8899d 100644
--- a/src/output/WriterBase.hpp
+++ b/src/output/WriterBase.hpp
@@ -13,6 +13,7 @@ class WriterBase : public IWriter
  protected:
   const std::string m_base_filename;
   const double m_time_period;
+  mutable double m_next_time;
 
   mutable std::vector<double> m_saved_times;
 
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3336fd716822a86f86fc335a3f02cb80882a26fc
--- /dev/null
+++ b/src/scheme/AcousticSolver.cpp
@@ -0,0 +1,787 @@
+#include <scheme/AcousticSolver.hpp>
+
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+
+#include <variant>
+#include <vector>
+
+template <size_t Dimension>
+double
+acoustic_dt(const DiscreteFunctionP0<Dimension, double>& c)
+{
+  const std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh =
+    std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(c.mesh());
+
+  const auto Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+  const auto Sj = MeshDataManager::instance().getMeshData(*mesh).sumOverRLjr();
+
+  CellValue<double> local_dt{mesh->connectivity()};
+  parallel_for(
+    mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); });
+
+  return min(local_dt);
+}
+
+double
+acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c)
+{
+  if ((c->descriptor().type() != DiscreteFunctionType::P0) or (c->dataType() != ASTNodeDataType::double_t)) {
+    throw NormalError("invalid discrete function type");
+  }
+
+  std::shared_ptr mesh = c->mesh();
+
+  switch (mesh->dimension()) {
+  case 1: {
+    return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*c));
+  }
+  case 2: {
+    return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*c));
+  }
+  case 3: {
+    return acoustic_dt(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*c));
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+template <size_t Dimension>
+class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler::IAcousticSolver
+{
+ private:
+  using Rdxd = TinyMatrix<Dimension>;
+  using Rd   = TinyVector<Dimension>;
+
+  using MeshType     = Mesh<Connectivity<Dimension>>;
+  using MeshDataType = MeshData<Dimension>;
+
+  using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, double>;
+  using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, Rd>;
+
+  class PressureBoundaryCondition;
+  class SymmetryBoundaryCondition;
+  class VelocityBoundaryCondition;
+
+  using BoundaryCondition =
+    std::variant<PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
+
+  using BoundaryConditionList = std::vector<BoundaryCondition>;
+
+  const SolverType m_solver_type;
+  BoundaryConditionList m_boundary_condition_list;
+
+  NodeValue<const Rd> m_ur;
+  NodeValuePerCell<const Rd> m_Fjr;
+
+  CellValue<const double>
+  _getRhoC(const DiscreteScalarFunction& rho, const DiscreteScalarFunction& c)
+  {
+    Assert(rho.mesh() == c.mesh());
+
+    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+
+    CellValue<double> rhoc{mesh->connectivity()};
+    parallel_for(
+      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rhoc[cell_id] = rho[cell_id] * c[cell_id]; });
+
+    return rhoc;
+  }
+
+  NodeValuePerCell<const Rdxd>
+  _computeGlaceAjr(const MeshType& mesh, const CellValue<const double>& rhoc)
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const TinyVector<Dimension>> Cjr = mesh_data.Cjr();
+    const NodeValuePerCell<const TinyVector<Dimension>> njr = mesh_data.njr();
+
+    NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const size_t& nb_nodes = Ajr.numberOfSubValues(j);
+        const double& rhoc_j   = rhoc[j];
+        for (size_t r = 0; r < nb_nodes; ++r) {
+          Ajr(j, r) = tensorProduct(rhoc_j * Cjr(j, r), njr(j, r));
+        }
+      });
+
+    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)
+  {
+    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+    NodeValue<Rdxd> Ar{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        Rdxd 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);
+        }
+        Ar[r] = sum;
+      });
+
+    return Ar;
+  }
+
+  NodeValue<Rd>
+  _computeBr(const Mesh<Connectivity<Dimension>>& mesh,
+             const NodeValuePerCell<const Rdxd>& Ajr,
+             const DiscreteVectorFunction& u,
+             const DiscreteScalarFunction& p)
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
+
+    const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+    const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+    NodeValue<Rd> b{mesh.connectivity()};
+
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+        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);
+
+        Rd br = zero;
+        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) * u[J] + p[J] * Cjr(J, R);
+        }
+
+        b[r] = br;
+      });
+
+    return b;
+  }
+
+  BoundaryConditionList
+  _getBCList(const std::shared_ptr<const MeshType>& mesh,
+             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+  {
+    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: {
+        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()) {
+            SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)};
+            bc_list.push_back(SymmetryBoundaryCondition{MeshFlatNodeBoundary<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") {
+          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") {
+          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());
+      }
+    }
+
+    return bc_list;
+  }
+
+  void _applyPressureBC(const MeshType& mesh, NodeValue<Rd>& br);
+  void _applySymmetryBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br);
+  void _applyVelocityBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br);
+
+  void
+  _applyBoundaryConditions(const MeshType& mesh, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br)
+  {
+    this->_applyPressureBC(mesh, br);
+    this->_applySymmetryBC(Ar, br);
+    this->_applyVelocityBC(Ar, br);
+  }
+
+  NodeValue<const Rd>
+  _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br)
+  {
+    NodeValue<Rd> u{mesh.connectivity()};
+    parallel_for(
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { u[r] = inverse(Ar[r]) * br[r]; });
+
+    return u;
+  }
+
+  NodeValuePerCell<Rd>
+  _computeFjr(const MeshType& mesh,
+              const NodeValuePerCell<const Rdxd>& Ajr,
+              const NodeValue<const Rd>& ur,
+              const DiscreteVectorFunction& u,
+              const DiscreteScalarFunction& p)
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+
+    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    NodeValuePerCell<Rd> F{mesh.connectivity()};
+    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) {
+          F(j, r) = Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + p[j] * Cjr(j, r);
+        }
+      });
+
+    return F;
+  }
+
+  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_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->_computeAjr(mesh, rhoc);
+
+    NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
+    NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, p);
+
+    this->_applyBoundaryConditions(mesh, Ar, br);
+
+    synchronize(Ar);
+    synchronize(br);
+
+    m_ur  = this->_computeUr(mesh, Ar, br);
+    m_Fjr = this->_computeFjr(mesh, Ajr, m_ur, u, p);
+  }
+
+ public:
+  std::tuple<std::shared_ptr<const IMesh>,
+             std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>,
+             std::shared_ptr<const DiscreteFunctionP0<Dimension, Rd>>,
+             std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>>
+  apply(const double& dt,
+        const std::shared_ptr<const MeshType>& mesh,
+        const DiscreteFunctionP0<Dimension, double>& rho,
+        const DiscreteFunctionP0<Dimension, Rd>& u,
+        const DiscreteFunctionP0<Dimension, double>& E) const
+  {
+    const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+
+    NodeValue<Rd> new_xr = copy(mesh->xr());
+    parallel_for(
+      mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * m_ur[r]; });
+
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr);
+
+    CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+
+    CellValue<double> new_rho = copy(rho.cellValues());
+    CellValue<Rd> new_u       = copy(u.cellValues());
+    CellValue<double> new_E   = copy(E.cellValues());
+
+    parallel_for(
+      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 += m_Fjr(j, R);
+          energy_fluxes += (m_Fjr(j, R), m_ur[r]);
+        }
+        const double dt_over_Mj = dt / (rho[j] * Vj[j]);
+        new_u[j] -= dt_over_Mj * momentum_fluxes;
+        new_E[j] -= dt_over_Mj * energy_fluxes;
+      });
+
+    CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
+
+    parallel_for(
+      mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+
+    return {new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_rho),
+            std::make_shared<DiscreteVectorFunction>(new_mesh, new_u),
+            std::make_shared<DiscreteScalarFunction>(new_mesh, new_E)};
+  }
+
+  std::tuple<std::shared_ptr<const IMesh>,
+             std::shared_ptr<const IDiscreteFunction>,
+             std::shared_ptr<const IDiscreteFunction>,
+             std::shared_ptr<const IDiscreteFunction>>
+  apply(const double& 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 i_mesh = getCommonMesh({rho, u, E});
+    if (not i_mesh) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
+      throw NormalError("acoustic solver expects P0 functions");
+    }
+
+    return this->apply(dt, std::dynamic_pointer_cast<const MeshType>(i_mesh),
+                       *std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho),
+                       *std::dynamic_pointer_cast<const DiscreteVectorFunction>(u),
+                       *std::dynamic_pointer_cast<const DiscreteScalarFunction>(E));
+  }
+
+  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{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),
+                     *std::dynamic_pointer_cast<const DiscreteScalarFunction>(p),
+                     bc_descriptor_list}
+  {}
+
+  AcousticSolver()                 = default;
+  AcousticSolver(AcousticSolver&&) = default;
+  ~AcousticSolver()                = default;
+};
+
+template <size_t Dimension>
+void
+AcousticSolverHandler::AcousticSolver<Dimension>::_applyPressureBC(const MeshType& mesh, NodeValue<Rd>& br)
+{
+  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(mesh);
+          if constexpr (Dimension == 1) {
+            const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+            const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
+            const auto& node_local_numbers_in_their_cells = mesh.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);
+
+                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               = mesh.connectivity().faceToCellMatrix();
+            const auto& face_to_node_matrix               = mesh.connectivity().faceToNodeMatrix();
+            const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
+            const auto& face_cell_is_reversed             = mesh.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];
+                br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
+              }
+            }
+          }
+        }
+      },
+      boundary_condition);
+  }
+}
+
+template <size_t Dimension>
+void
+AcousticSolverHandler::AcousticSolver<Dimension>::_applySymmetryBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br)
+{
+  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 Rdxd I   = identity;
+          const Rdxd nxn = tensorProduct(n, n);
+          const Rdxd 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];
+
+              Ar[r] = P * Ar[r] * P + nxn;
+              br[r] = P * br[r];
+            });
+        }
+      },
+      boundary_condition);
+  }
+}
+
+template <size_t Dimension>
+void
+AcousticSolverHandler::AcousticSolver<Dimension>::_applyVelocityBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br)
+{
+  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];
+
+              Ar[node_id] = identity;
+              br[node_id] = value;
+            });
+        }
+      },
+      boundary_condition);
+  }
+}
+
+template <size_t Dimension>
+class AcousticSolverHandler::AcousticSolver<Dimension>::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 AcousticSolverHandler::AcousticSolver<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 <size_t Dimension>
+class AcousticSolverHandler::AcousticSolver<Dimension>::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 <size_t Dimension>
+class AcousticSolverHandler::AcousticSolver<Dimension>::SymmetryBoundaryCondition
+{
+ public:
+  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;
+};
+
+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,
+  const std::shared_ptr<const IDiscreteFunction>& p,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+{
+  std::shared_ptr i_mesh = getCommonMesh({rho, c, u, p});
+  if (not i_mesh) {
+    throw NormalError("discrete functions are not defined on the same mesh");
+  }
+
+  if (not checkDiscretizationType({rho, c, u, p}, DiscreteFunctionType::P0)) {
+    throw NormalError("acoustic solver expects P0 functions");
+  }
+
+  switch (i_mesh->dimension()) {
+  case 1: {
+    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>>(solver_type, i_mesh, rho, c, u, p, bc_descriptor_list);
+    break;
+  }
+  case 3: {
+    m_acoustic_solver = std::make_unique<AcousticSolver<3>>(solver_type, i_mesh, rho, c, u, p, bc_descriptor_list);
+    break;
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 3e3151d25429593955cb281c2dd814690f67bfa5..2ed97075717f8984afb22799648b4b26fbc9327b 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -1,612 +1,62 @@
 #ifndef ACOUSTIC_SOLVER_HPP
 #define ACOUSTIC_SOLVER_HPP
 
-#include <algebra/TinyMatrix.hpp>
-#include <algebra/TinyVector.hpp>
-#include <mesh/ItemValueUtils.hpp>
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
-#include <mesh/MeshDataManager.hpp>
-#include <mesh/MeshNodeBoundary.hpp>
-#include <mesh/SubItemValuePerItem.hpp>
-#include <scheme/AcousticSolverType.hpp>
-#include <scheme/BlockPerfectGas.hpp>
-#include <scheme/FiniteVolumesEulerUnknowns.hpp>
-#include <utils/ArrayUtils.hpp>
-#include <utils/Exceptions.hpp>
-#include <utils/Messenger.hpp>
-#include <utils/PugsAssert.hpp>
+#include <memory>
+#include <tuple>
+#include <vector>
 
-#include <rang.hpp>
+class IDiscreteFunction;
+class IBoundaryConditionDescriptor;
+class IMesh;
 
-#include <iostream>
+double acoustic_dt(const std::shared_ptr<const IDiscreteFunction>& c);
 
-template <typename MeshType>
-class AcousticSolver
+class AcousticSolverHandler
 {
  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:
-  AcousticSolver(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 AcousticSolver<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
+  enum class SolverType
   {
-    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 AcousticSolver<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;
-};
+    Glace,
+    Eucclhyd
+  };
 
-template <typename MeshType>
-class AcousticSolver<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
+  struct IAcousticSolver
   {
-    return m_node_list;
-  }
+    virtual std::tuple<std::shared_ptr<const IMesh>,
+                       std::shared_ptr<const IDiscreteFunction>,
+                       std::shared_ptr<const IDiscreteFunction>,
+                       std::shared_ptr<const IDiscreteFunction>>
+    apply(const double& dt,
+          const std::shared_ptr<const IDiscreteFunction>& rho,
+          const std::shared_ptr<const IDiscreteFunction>& u,
+          const std::shared_ptr<const IDiscreteFunction>& E) const = 0;
 
-  const Array<const TinyVector<Dimension>>&
-  valueList() const
-  {
-    return m_value_list;
-  }
+    IAcousticSolver()                  = default;
+    IAcousticSolver(IAcousticSolver&&) = default;
+    IAcousticSolver& operator=(IAcousticSolver&&) = default;
 
-  VelocityBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list)
-    : m_value_list{value_list}, m_node_list{node_list}
-  {}
+    virtual ~IAcousticSolver() = default;
+  };
 
-  ~VelocityBoundaryCondition() = default;
-};
+  template <size_t Dimension>
+  class AcousticSolver;
 
-template <typename MeshType>
-class AcousticSolver<MeshType>::SymmetryBoundaryCondition
-{
- public:
-  static constexpr size_t Dimension = MeshType::Dimension;
-
-  using Rd = TinyVector<Dimension, double>;
-
- private:
-  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+  std::unique_ptr<IAcousticSolver> m_acoustic_solver;
 
  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)
+  const IAcousticSolver&
+  solver() const
   {
-    ;
+    return *m_acoustic_solver;
   }
 
-  ~SymmetryBoundaryCondition() = default;
+  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);
 };
 
 #endif   // ACOUSTIC_SOLVER_HPP
diff --git a/src/scheme/AcousticSolverType.hpp b/src/scheme/AcousticSolverType.hpp
deleted file mode 100644
index 3095dcc864d5598e985290302ad94106adbd6b2e..0000000000000000000000000000000000000000
--- a/src/scheme/AcousticSolverType.hpp
+++ /dev/null
@@ -1,10 +0,0 @@
-#ifndef ACOUSTIC_SOLVER_TYPE_HPP
-#define ACOUSTIC_SOLVER_TYPE_HPP
-
-enum class AcousticSolverType
-{
-  Eucclhyd,
-  Glace
-};
-
-#endif   // ACOUSTIC_SOLVER_TYPE_HPP
diff --git a/src/scheme/BlockPerfectGas.hpp b/src/scheme/BlockPerfectGas.hpp
deleted file mode 100644
index 38aef3262a39101fdc7607143a0bdb0e36d4cfe5..0000000000000000000000000000000000000000
--- a/src/scheme/BlockPerfectGas.hpp
+++ /dev/null
@@ -1,59 +0,0 @@
-#ifndef BLOCK_PERFECTGAS_HPP
-#define BLOCK_PERFECTGAS_HPP
-
-#include <mesh/ItemValue.hpp>
-
-struct BlockPerfectGas
-{
- private:
-  CellValue<double>& m_rhoj;
-  CellValue<double>& m_ej;
-  CellValue<double>& m_pj;
-  CellValue<double>& m_gammaj;
-  CellValue<double>& m_cj;
-
- public:
-  BlockPerfectGas(CellValue<double>& rhoj,
-                  CellValue<double>& ej,
-                  CellValue<double>& pj,
-                  CellValue<double>& gammaj,
-                  CellValue<double>& cj)
-    : m_rhoj(rhoj), m_ej(ej), m_pj(pj), m_gammaj(gammaj), m_cj(cj)
-  {
-    ;
-  }
-
-  void
-  updatePandCFromRhoE()
-  {
-    const size_t nj                      = m_ej.size();
-    const CellValue<const double>& rho   = m_rhoj;
-    const CellValue<const double>& e     = m_ej;
-    const CellValue<const double>& gamma = m_gammaj;
-
-    parallel_for(
-      nj, PUGS_LAMBDA(CellId j) {
-        const double gamma_minus_one = gamma[j] - 1;
-        m_pj[j]                      = gamma_minus_one * rho[j] * e[j];
-        m_cj[j]                      = std::sqrt(gamma[j] * gamma_minus_one * e[j]);
-      });
-  }
-
-  void
-  updateEandCFromRhoP()
-  {
-    const size_t nj                      = m_ej.size();
-    const CellValue<const double>& rho   = m_rhoj;
-    const CellValue<const double>& p     = m_pj;
-    const CellValue<const double>& gamma = m_gammaj;
-
-    parallel_for(
-      nj, PUGS_LAMBDA(CellId j) { m_ej[j] = p[j] / (rho[j] * (gamma[j] - 1)); });
-
-    const CellValue<const double>& e = m_ej;
-    parallel_for(
-      nj, PUGS_LAMBDA(CellId j) { m_cj[j] = std::sqrt(gamma[j] * (gamma[j] - 1) * e[j]); });
-  }
-};
-
-#endif   // BLOCK_PERFECTGAS_HPP
diff --git a/src/scheme/CMakeLists.txt b/src/scheme/CMakeLists.txt
index 1467393f2cbcbd9e31dd495e2905e9be478a1d78..b7522b9bbd08a90ee1bd95a56747753d902e10ed 100644
--- a/src/scheme/CMakeLists.txt
+++ b/src/scheme/CMakeLists.txt
@@ -2,4 +2,7 @@
 
 add_library(
   PugsScheme
-  DiscreteFunctionInterpoler.cpp)
+  AcousticSolver.cpp
+  DiscreteFunctionInterpoler.cpp
+  DiscreteFunctionUtils.cpp
+  PerfectGas.cpp)
diff --git a/src/scheme/DiscreteFunctionP0.hpp b/src/scheme/DiscreteFunctionP0.hpp
index bf97079b776c80fec10e26608cfc6ff5e40a664a..c51e8e09dbb78dc6f40f4fbb18eb737ddb7e1a30 100644
--- a/src/scheme/DiscreteFunctionP0.hpp
+++ b/src/scheme/DiscreteFunctionP0.hpp
@@ -28,7 +28,7 @@ class DiscreteFunctionP0 : public IDiscreteFunction
     return ast_node_data_type_from<DataType>;
   }
 
-  CellValue<DataType>
+  const CellValue<DataType>&
   cellValues() const
   {
     return m_cell_values;
@@ -63,6 +63,17 @@ class DiscreteFunctionP0 : public IDiscreteFunction
                                                                                                   mesh_data.xj());
   }
 
+  DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh) : m_mesh{mesh}, m_cell_values{mesh->connectivity()}
+  {
+    ;
+  }
+
+  DiscreteFunctionP0(const std::shared_ptr<const MeshType>& mesh, const CellValue<DataType>& cell_value)
+    : m_mesh{mesh}, m_cell_values{cell_value}
+  {
+    ;
+  }
+
   DiscreteFunctionP0(const DiscreteFunctionP0&) noexcept = default;
   DiscreteFunctionP0(DiscreteFunctionP0&&) noexcept      = default;
 
diff --git a/src/scheme/DiscreteFunctionUtils.cpp b/src/scheme/DiscreteFunctionUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..56453a746a19b9ee999fe1548d7067250e2884eb
--- /dev/null
+++ b/src/scheme/DiscreteFunctionUtils.cpp
@@ -0,0 +1,115 @@
+#include <scheme/DiscreteFunctionUtils.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/IMesh.hpp>
+#include <mesh/Mesh.hpp>
+#include <scheme/DiscreteFunctionP0.hpp>
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<const IDiscreteFunction>
+shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh,
+            const std::shared_ptr<const DiscreteFunctionP0<Dimension, DataType>>& discrete_function)
+{
+  return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(mesh, discrete_function->cellValues());
+}
+
+template <size_t Dimension>
+std::shared_ptr<const IDiscreteFunction>
+shallowCopy(const std::shared_ptr<const Mesh<Connectivity<Dimension>>>& mesh,
+            const std::shared_ptr<const IDiscreteFunction>& discrete_function)
+{
+  const std::shared_ptr function_mesh =
+    std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(discrete_function->mesh());
+
+  if (mesh->shared_connectivity() != function_mesh->shared_connectivity()) {
+    throw NormalError("incompatible connectivities");
+  }
+
+  switch (discrete_function->descriptor().type()) {
+  case DiscreteFunctionType::P0: {
+    switch (discrete_function->dataType()) {
+    case ASTNodeDataType::double_t: {
+      return shallowCopy(mesh,
+                         std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, double>>(discrete_function));
+    }
+    case ASTNodeDataType::vector_t: {
+      switch (discrete_function->dataType().dimension()) {
+      case 1: {
+        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<1>>>(
+                                   discrete_function));
+      }
+      case 2: {
+        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<2>>>(
+                                   discrete_function));
+      }
+      case 3: {
+        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyVector<3>>>(
+                                   discrete_function));
+      }
+      default: {
+        throw UnexpectedError("invalid data vector dimension: " +
+                              std::to_string(discrete_function->dataType().dimension()));
+      }
+      }
+    }
+    case ASTNodeDataType::matrix_t: {
+      if (discrete_function->dataType().nbRows() != discrete_function->dataType().nbColumns()) {
+        throw UnexpectedError(
+          "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().nbRows()) + "x" +
+          std::to_string(discrete_function->dataType().nbColumns()));
+      }
+      switch (discrete_function->dataType().nbRows()) {
+      case 1: {
+        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<1>>>(
+                                   discrete_function));
+      }
+      case 2: {
+        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<2>>>(
+                                   discrete_function));
+      }
+      case 3: {
+        return shallowCopy(mesh, std::dynamic_pointer_cast<const DiscreteFunctionP0<Dimension, TinyMatrix<3>>>(
+                                   discrete_function));
+      }
+      default: {
+        throw UnexpectedError(
+          "invalid data matrix dimensions: " + std::to_string(discrete_function->dataType().nbRows()) + "x" +
+          std::to_string(discrete_function->dataType().nbColumns()));
+      }
+      }
+    }
+    default: {
+      throw UnexpectedError("invalid kind of P0 function: invalid data type");
+    }
+    }
+  }
+  default: {
+    throw NormalError("invalid discretization type");
+  }
+  }
+}
+
+std::shared_ptr<const IDiscreteFunction>
+shallowCopy(const std::shared_ptr<const IMesh>& mesh, const std::shared_ptr<const IDiscreteFunction>& discrete_function)
+{
+  if (mesh == discrete_function->mesh()) {
+    return discrete_function;
+  } else if (mesh->dimension() != discrete_function->mesh()->dimension()) {
+    throw NormalError("incompatible mesh dimensions");
+  }
+
+  switch (mesh->dimension()) {
+  case 1: {
+    return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<1>>>(mesh), discrete_function);
+  }
+  case 2: {
+    return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<2>>>(mesh), discrete_function);
+  }
+  case 3: {
+    return shallowCopy(std::dynamic_pointer_cast<const Mesh<Connectivity<3>>>(mesh), discrete_function);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
diff --git a/src/scheme/DiscreteFunctionUtils.hpp b/src/scheme/DiscreteFunctionUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1eea93a7eae8db4c592b6bbec35341ad33e35bdc
--- /dev/null
+++ b/src/scheme/DiscreteFunctionUtils.hpp
@@ -0,0 +1,43 @@
+#ifndef DISCRETE_FUNCTION_UTILS_HPP
+#define DISCRETE_FUNCTION_UTILS_HPP
+
+#include <scheme/DiscreteFunctionType.hpp>
+#include <scheme/IDiscreteFunction.hpp>
+#include <scheme/IDiscreteFunctionDescriptor.hpp>
+
+#include <vector>
+
+PUGS_INLINE
+bool
+checkDiscretizationType(const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list,
+                        const DiscreteFunctionType& discrete_function_type)
+{
+  for (auto discrete_function : discrete_function_list) {
+    if (discrete_function->descriptor().type() != discrete_function_type) {
+      return false;
+    }
+  }
+  return true;
+}
+
+PUGS_INLINE
+std::shared_ptr<const IMesh>
+getCommonMesh(const std::vector<std::shared_ptr<const IDiscreteFunction>>& discrete_function_list)
+{
+  std::shared_ptr<const IMesh> i_mesh;
+  for (auto discrete_function : discrete_function_list) {
+    if (not i_mesh.use_count()) {
+      i_mesh = discrete_function->mesh();
+    } else {
+      if (i_mesh != discrete_function->mesh()) {
+        return {};
+      }
+    }
+  }
+  return i_mesh;
+}
+
+std::shared_ptr<const IDiscreteFunction> shallowCopy(const std::shared_ptr<const IMesh>& mesh,
+                                                     const std::shared_ptr<const IDiscreteFunction>& discrete_function);
+
+#endif   // DISCRETE_FUNCTION_UTILS_HPP
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
deleted file mode 100644
index ec7e8851ed98b64971810cfccbf2b587615e8939..0000000000000000000000000000000000000000
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ /dev/null
@@ -1,152 +0,0 @@
-#ifndef FINITE_VOLUMES_EULER_UNKNOWNS_HPP
-#define FINITE_VOLUMES_EULER_UNKNOWNS_HPP
-
-#include <algebra/TinyVector.hpp>
-#include <mesh/ItemValue.hpp>
-#include <scheme/BlockPerfectGas.hpp>
-
-template <typename MeshType>
-class FiniteVolumesEulerUnknowns
-{
- public:
-  static constexpr size_t Dimension = MeshType::Dimension;
-
-  using Rd = TinyVector<Dimension>;
-
- private:
-  CellValue<double> m_rhoj;
-  CellValue<Rd> m_uj;
-  CellValue<double> m_Ej;
-
-  CellValue<double> m_ej;
-  CellValue<double> m_pj;
-  CellValue<double> m_gammaj;
-  CellValue<double> m_cj;
-  CellValue<double> m_mj;
-  CellValue<double> m_inv_mj;
-
- public:
-  CellValue<double>&
-  rhoj()
-  {
-    return m_rhoj;
-  }
-
-  const CellValue<const double>
-  rhoj() const
-  {
-    return m_rhoj;
-  }
-
-  CellValue<Rd>&
-  uj()
-  {
-    return m_uj;
-  }
-
-  const CellValue<const Rd>
-  uj() const
-  {
-    return m_uj;
-  }
-
-  CellValue<double>&
-  Ej()
-  {
-    return m_Ej;
-  }
-
-  const CellValue<const double>
-  Ej() const
-  {
-    return m_Ej;
-  }
-
-  CellValue<double>&
-  ej()
-  {
-    return m_ej;
-  }
-
-  const CellValue<const double>
-  ej() const
-  {
-    return m_ej;
-  }
-
-  CellValue<double>&
-  pj()
-  {
-    return m_pj;
-  }
-
-  const CellValue<const double>
-  pj() const
-  {
-    return m_pj;
-  }
-
-  CellValue<double>&
-  gammaj()
-  {
-    return m_gammaj;
-  }
-
-  const CellValue<const double>
-  gammaj() const
-  {
-    return m_gammaj;
-  }
-
-  CellValue<double>&
-  cj()
-  {
-    return m_cj;
-  }
-
-  const CellValue<const double>
-  cj() const
-  {
-    return m_cj;
-  }
-
-  CellValue<double>&
-  mj()
-  {
-    return m_mj;
-  }
-
-  const CellValue<const double>
-  mj() const
-  {
-    return m_mj;
-  }
-
-  CellValue<double>&
-  invMj()
-  {
-    return m_inv_mj;
-  }
-
-  const CellValue<const double>
-  invMj() const
-  {
-    return m_inv_mj;
-  }
-
-  FiniteVolumesEulerUnknowns(const MeshType& mesh)
-    : m_rhoj(mesh.connectivity()),
-      m_uj(mesh.connectivity()),
-      m_Ej(mesh.connectivity()),
-      m_ej(mesh.connectivity()),
-      m_pj(mesh.connectivity()),
-      m_gammaj(mesh.connectivity()),
-      m_cj(mesh.connectivity()),
-      m_mj(mesh.connectivity()),
-      m_inv_mj(mesh.connectivity())
-  {
-    ;
-  }
-};
-
-#endif   // FINITE_VOLUMES_EULER_UNKNOWNS_HPP
diff --git a/src/scheme/PerfectGas.cpp b/src/scheme/PerfectGas.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b371dea3c9e4aac8ee867dd0f938b72bcd02791e
--- /dev/null
+++ b/src/scheme/PerfectGas.cpp
@@ -0,0 +1,248 @@
+#include <scheme/PerfectGas.hpp>
+
+#include <scheme/DiscreteFunctionP0.hpp>
+#include <scheme/DiscreteFunctionUtils.hpp>
+
+namespace perfect_gas
+{
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+soundSpeed(const DiscreteFunctionP0<Dimension, DataType>& rho,
+           const DiscreteFunctionP0<Dimension, DataType>& p,
+           const DiscreteFunctionP0<Dimension, DataType>& gamma)
+{
+  using MeshType       = Mesh<Connectivity<Dimension>>;
+  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+
+  DiscreteFunctionP0<Dimension, DataType> sound_speed{mesh};
+
+  parallel_for(
+    mesh->numberOfCells(),
+    PUGS_LAMBDA(CellId cell_id) { sound_speed[cell_id] = std::sqrt(gamma[cell_id] * p[cell_id] / rho[cell_id]); });
+
+  return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(sound_speed);
+}
+
+std::shared_ptr<IDiscreteFunction>
+soundSpeed(const std::shared_ptr<const IDiscreteFunction>& rho,
+           const std::shared_ptr<const IDiscreteFunction>& p,
+           const std::shared_ptr<const IDiscreteFunction>& gamma)
+{
+  std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, p, gamma});
+  if (not mesh.use_count()) {
+    throw NormalError("rho, p and gamma are not defined on the same mesh");
+  }
+
+  switch (mesh->dimension()) {
+  case 1: {
+    return soundSpeed(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho),
+                      dynamic_cast<const DiscreteFunctionP0<1, double>&>(*p),
+                      dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma));
+  }
+  case 2: {
+    return soundSpeed(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho),
+                      dynamic_cast<const DiscreteFunctionP0<2, double>&>(*p),
+                      dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma));
+  }
+  case 3: {
+    return soundSpeed(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho),
+                      dynamic_cast<const DiscreteFunctionP0<3, double>&>(*p),
+                      dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma));
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+epsilonFromRhoPGamma(const DiscreteFunctionP0<Dimension, DataType>& rho,
+                     const DiscreteFunctionP0<Dimension, DataType>& p,
+                     const DiscreteFunctionP0<Dimension, DataType>& gamma)
+{
+  using MeshType       = Mesh<Connectivity<Dimension>>;
+  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+
+  DiscreteFunctionP0<Dimension, DataType> epsilon{mesh};
+
+  parallel_for(
+    mesh->numberOfCells(),
+    PUGS_LAMBDA(CellId cell_id) { epsilon[cell_id] = p[cell_id] / ((gamma[cell_id] - 1) * rho[cell_id]); });
+
+  return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon);
+}
+
+std::shared_ptr<IDiscreteFunction>
+epsilonFromRhoPGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
+                     const std::shared_ptr<const IDiscreteFunction>& p,
+                     const std::shared_ptr<const IDiscreteFunction>& gamma)
+{
+  std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, p, gamma});
+  if (not mesh.use_count()) {
+    throw NormalError("rho, p and gamma are not defined on the same mesh");
+  }
+
+  switch (mesh->dimension()) {
+  case 1: {
+    return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho),
+                                dynamic_cast<const DiscreteFunctionP0<1, double>&>(*p),
+                                dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma));
+  }
+  case 2: {
+    return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho),
+                                dynamic_cast<const DiscreteFunctionP0<2, double>&>(*p),
+                                dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma));
+  }
+  case 3: {
+    return epsilonFromRhoPGamma(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho),
+                                dynamic_cast<const DiscreteFunctionP0<3, double>&>(*p),
+                                dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma));
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+pFromRhoEpsilonGamma(const DiscreteFunctionP0<Dimension, DataType>& rho,
+                     const DiscreteFunctionP0<Dimension, DataType>& epsilon,
+                     const DiscreteFunctionP0<Dimension, DataType>& gamma)
+{
+  using MeshType       = Mesh<Connectivity<Dimension>>;
+  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+
+  DiscreteFunctionP0<Dimension, DataType> p{mesh};
+
+  parallel_for(
+    mesh->numberOfCells(),
+    PUGS_LAMBDA(CellId cell_id) { p[cell_id] = ((gamma[cell_id] - 1) * rho[cell_id]) * epsilon[cell_id]; });
+
+  return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(p);
+}
+
+std::shared_ptr<IDiscreteFunction>
+pFromRhoEpsilonGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
+                     const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                     const std::shared_ptr<const IDiscreteFunction>& gamma)
+{
+  std::shared_ptr<const IMesh> mesh = getCommonMesh({rho, epsilon, gamma});
+  if (not mesh.use_count()) {
+    throw NormalError("rho, epsilon and gamma are not defined on the same mesh");
+  }
+
+  switch (mesh->dimension()) {
+  case 1: {
+    return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*rho),
+                                dynamic_cast<const DiscreteFunctionP0<1, double>&>(*epsilon),
+                                dynamic_cast<const DiscreteFunctionP0<1, double>&>(*gamma));
+  }
+  case 2: {
+    return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*rho),
+                                dynamic_cast<const DiscreteFunctionP0<2, double>&>(*epsilon),
+                                dynamic_cast<const DiscreteFunctionP0<2, double>&>(*gamma));
+  }
+  case 3: {
+    return pFromRhoEpsilonGamma(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*rho),
+                                dynamic_cast<const DiscreteFunctionP0<3, double>&>(*epsilon),
+                                dynamic_cast<const DiscreteFunctionP0<3, double>&>(*gamma));
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+totalEnergyFromEpsilonU(const DiscreteFunctionP0<Dimension, DataType>& epsilon,
+                        const DiscreteFunctionP0<Dimension, TinyVector<Dimension, DataType>>& u)
+{
+  using MeshType       = Mesh<Connectivity<Dimension>>;
+  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(epsilon.mesh());
+
+  DiscreteFunctionP0<Dimension, DataType> E{mesh};
+
+  parallel_for(
+    mesh->numberOfCells(),
+    PUGS_LAMBDA(CellId cell_id) { E[cell_id] = epsilon[cell_id] + 0.5 * (u[cell_id], u[cell_id]); });
+
+  return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon);
+}
+
+std::shared_ptr<IDiscreteFunction>
+totalEnergyFromEpsilonU(const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                        const std::shared_ptr<const IDiscreteFunction>& u)
+{
+  std::shared_ptr<const IMesh> mesh = getCommonMesh({epsilon, u});
+  if (not mesh.use_count()) {
+    throw NormalError("epsilon and u are not defined on the same mesh");
+  }
+
+  switch (mesh->dimension()) {
+  case 1: {
+    return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*epsilon),
+                                   dynamic_cast<const DiscreteFunctionP0<1, TinyVector<1, double>>&>(*u));
+  }
+  case 2: {
+    return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*epsilon),
+                                   dynamic_cast<const DiscreteFunctionP0<2, TinyVector<2, double>>&>(*u));
+  }
+  case 3: {
+    return totalEnergyFromEpsilonU(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*epsilon),
+                                   dynamic_cast<const DiscreteFunctionP0<3, TinyVector<3, double>>&>(*u));
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+template <size_t Dimension, typename DataType>
+std::shared_ptr<IDiscreteFunction>
+epsilonFromTotalEnergyU(const DiscreteFunctionP0<Dimension, DataType>& E,
+                        const DiscreteFunctionP0<Dimension, TinyVector<Dimension, DataType>>& u)
+{
+  using MeshType       = Mesh<Connectivity<Dimension>>;
+  std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(u.mesh());
+
+  DiscreteFunctionP0<Dimension, DataType> epsilon{mesh};
+
+  parallel_for(
+    mesh->numberOfCells(),
+    PUGS_LAMBDA(CellId cell_id) { epsilon[cell_id] = E[cell_id] - 0.5 * (u[cell_id], u[cell_id]); });
+
+  return std::make_shared<DiscreteFunctionP0<Dimension, DataType>>(epsilon);
+}
+
+std::shared_ptr<IDiscreteFunction>
+epsilonFromTotalEnergyU(const std::shared_ptr<const IDiscreteFunction>& E,
+                        const std::shared_ptr<const IDiscreteFunction>& u)
+{
+  std::shared_ptr<const IMesh> mesh = getCommonMesh({E, u});
+  if (not mesh.use_count()) {
+    throw NormalError("E and u are not defined on the same mesh");
+  }
+
+  switch (mesh->dimension()) {
+  case 1: {
+    return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*E),
+                                   dynamic_cast<const DiscreteFunctionP0<1, TinyVector<1, double>>&>(*u));
+  }
+  case 2: {
+    return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*E),
+                                   dynamic_cast<const DiscreteFunctionP0<2, TinyVector<2, double>>&>(*u));
+  }
+  case 3: {
+    return epsilonFromTotalEnergyU(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*E),
+                                   dynamic_cast<const DiscreteFunctionP0<3, TinyVector<3, double>>&>(*u));
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+}   // namespace perfect_gas
diff --git a/src/scheme/PerfectGas.hpp b/src/scheme/PerfectGas.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fb7b9d212b4681f365faa04e4565dfd78c316f66
--- /dev/null
+++ b/src/scheme/PerfectGas.hpp
@@ -0,0 +1,32 @@
+#ifndef PERFECT_GAS_HPP
+#define PERFECT_GAS_HPP
+
+#include <scheme/IDiscreteFunction.hpp>
+
+namespace perfect_gas
+{
+std::shared_ptr<IDiscreteFunction> soundSpeed(const std::shared_ptr<const IDiscreteFunction>& rho,
+                                              const std::shared_ptr<const IDiscreteFunction>& p,
+                                              const std::shared_ptr<const IDiscreteFunction>& gamma);
+
+std::shared_ptr<IDiscreteFunction> epsilonFromRhoPGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
+                                                        const std::shared_ptr<const IDiscreteFunction>& p,
+                                                        const std::shared_ptr<const IDiscreteFunction>& gamma);
+
+std::shared_ptr<IDiscreteFunction> pFromRhoEpsilonGamma(const std::shared_ptr<const IDiscreteFunction>& rho,
+                                                        const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                                                        const std::shared_ptr<const IDiscreteFunction>& gamma);
+
+// This is a temporary function
+// Todo: TO REMOVE
+std::shared_ptr<IDiscreteFunction> totalEnergyFromEpsilonU(const std::shared_ptr<const IDiscreteFunction>& epsilon,
+                                                           const std::shared_ptr<const IDiscreteFunction>& u);
+
+// This is a temporary function
+// Todo: TO REMOVE
+std::shared_ptr<IDiscreteFunction> epsilonFromTotalEnergyU(const std::shared_ptr<const IDiscreteFunction>& E,
+                                                           const std::shared_ptr<const IDiscreteFunction>& u);
+
+}   // namespace perfect_gas
+
+#endif   // PERFECT_GAS_HPP