diff --git a/CMakeLists.txt b/CMakeLists.txt
index 58fc97c458dbb7fc65fd0ebbf8d519a19f254f5e..8b5a8d153925bbe9c57dd04c817ad25ad67e01c6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -425,6 +425,7 @@ target_link_libraries(
   PugsLanguage
   PugsLanguageAST
   PugsLanguageModules
+  PugsLanguageAlgorithms
   PugsMesh
   PugsUtils
   PugsLanguageUtils
diff --git a/src/language/CMakeLists.txt b/src/language/CMakeLists.txt
index b9c6094e4403cb1d2289d8cfcae58de9e2e3c4e5..3565082161700da019088f7bfca70b0d9dfe6fd8 100644
--- a/src/language/CMakeLists.txt
+++ b/src/language/CMakeLists.txt
@@ -1,5 +1,6 @@
 # ------------------- Source files --------------------
 
+add_subdirectory(algorithms)
 add_subdirectory(ast)
 add_subdirectory(modules)
 add_subdirectory(node_processor)
@@ -12,6 +13,7 @@ add_library(
 # Additional dependencies
 add_dependencies(PugsLanguage
   PugsLanguage
+  PugsLanguageAlgorithms
   PugsLanguageAST
   PugsLanguageModules
   PugsLanguageUtils
diff --git a/src/language/algorithms/AcousticSolverAlgorithm.cpp b/src/language/algorithms/AcousticSolverAlgorithm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1a11de171d8172f9f276718d0624cd31b38312f4
--- /dev/null
+++ b/src/language/algorithms/AcousticSolverAlgorithm.cpp
@@ -0,0 +1,231 @@
+#include <language/algorithms/AcousticSolverAlgorithm.hpp>
+
+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)
+{
+  std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+  std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
+
+  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) {
+      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)});
+          }
+        }
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::velocity: {
+        using VelocityBoundaryCondition = typename AcousticSolver<MeshType>::VelocityBoundaryCondition;
+
+        const VelocityBoundaryConditionDescriptor& velocity_bc_descriptor =
+          dynamic_cast<const VelocityBoundaryConditionDescriptor&>(*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 == velocity_bc_descriptor.boundaryDescriptor()) {
+            MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
+
+            const FunctionSymbolId velocity_id = velocity_bc_descriptor.functionSymbolId();
+
+            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});
+          }
+        }
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::pressure: {
+        using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition;
+
+        const PressureBoundaryConditionDescriptor& pressure_bc_descriptor =
+          dynamic_cast<const PressureBoundaryConditionDescriptor&>(*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 == pressure_bc_descriptor.boundaryDescriptor()) {
+            const auto& face_list = ref_face_list.list();
+
+            const FunctionSymbolId pressure_id = pressure_bc_descriptor.functionSymbolId();
+
+            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});
+          }
+        }
+        break;
+      }
+      default: {
+        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]; });
+  }
+
+  VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
+
+  while ((t < tmax) and (iteration < itermax)) {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+    vtk_writer.write(mesh,
+                     {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                      NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+                      NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                      NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
+                     t);
+
+    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';
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+    vtk_writer.write(mesh,
+                     {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                      NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+                      NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                      NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
+                     t, true);   // forces last output
+  }
+}
+
+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
index 3fca4703f6d69c7790159519f6728479727183a5..bb9d9ef66d176274d305a02067e4f299c71d2c78 100644
--- a/src/language/algorithms/AcousticSolverAlgorithm.hpp
+++ b/src/language/algorithms/AcousticSolverAlgorithm.hpp
@@ -11,7 +11,7 @@
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 #include <scheme/VelocityBoundaryConditionDescriptor.hpp>
 
-template <size_t Dimension, AcousticSolverType solver_type>
+template <size_t Dimension>
 struct AcousticSolverAlgorithm
 {
   using ConnectivityType = Connectivity<Dimension>;
@@ -19,204 +19,12 @@ struct AcousticSolverAlgorithm
   using MeshDataType     = MeshData<Dimension>;
   using UnknownsType     = FiniteVolumesEulerUnknowns<MeshType>;
 
-  AcousticSolverAlgorithm(std::shared_ptr<const IMesh> i_mesh,
+  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)
-  {
-    std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
-
-    std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
-
-    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) {
-        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)});
-            }
-          }
-          break;
-        }
-        case IBoundaryConditionDescriptor::Type::velocity: {
-          using VelocityBoundaryCondition = typename AcousticSolver<MeshType>::VelocityBoundaryCondition;
-
-          const VelocityBoundaryConditionDescriptor& velocity_bc_descriptor =
-            dynamic_cast<const VelocityBoundaryConditionDescriptor&>(*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 == velocity_bc_descriptor.boundaryDescriptor()) {
-              MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
-
-              const FunctionSymbolId velocity_id = velocity_bc_descriptor.functionSymbolId();
-
-              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});
-            }
-          }
-          break;
-        }
-        case IBoundaryConditionDescriptor::Type::pressure: {
-          using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition;
-
-          const PressureBoundaryConditionDescriptor& pressure_bc_descriptor =
-            dynamic_cast<const PressureBoundaryConditionDescriptor&>(*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 == pressure_bc_descriptor.boundaryDescriptor()) {
-              const auto& face_list = ref_face_list.list();
-
-              const FunctionSymbolId pressure_id = pressure_bc_descriptor.functionSymbolId();
-
-              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});
-            }
-          }
-          break;
-        }
-        default: {
-          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]; });
-    }
-
-    VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
-
-    while ((t < tmax) and (iteration < itermax)) {
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-      vtk_writer.write(mesh,
-                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
-                        NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
-                        NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
-                       t);
-
-      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';
-    {
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
-
-      vtk_writer.write(mesh,
-                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
-                        NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
-                        NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
-                       t, true);   // forces last output
-    }
-  }
+                          const FunctionSymbolId& p_id);
 };
 
 #endif   // ACOUSTIC_SOLVER_ALGORITHM_HPP
diff --git a/src/language/algorithms/CMakeLists.txt b/src/language/algorithms/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a8599f575058185b5ceeed0814b242ae6ec70b0a
--- /dev/null
+++ b/src/language/algorithms/CMakeLists.txt
@@ -0,0 +1,10 @@
+# ------------------- Source files --------------------
+
+add_library(PugsLanguageAlgorithms
+  AcousticSolverAlgorithm.cpp
+)
+
+
+add_dependencies(PugsLanguageAlgorithms
+  PugsUtils
+  PugsMesh)
diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt
index 703f19857cfb99c0053aecf1af0984dab371d72e..f8cf2f0b445882a672d54e60dbc75858bbb246c7 100644
--- a/src/language/modules/CMakeLists.txt
+++ b/src/language/modules/CMakeLists.txt
@@ -11,5 +11,6 @@ add_library(PugsLanguageModules
 
 
 add_dependencies(PugsLanguageModules
+  PugsLanguageAlgorithms
   PugsUtils
   PugsMesh)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index b38373e672a736a8fe4c64d2d6125b265e1cb7a1..0bdf7de8ec290730ff10595ba1f4866978bbb333 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -75,63 +75,97 @@ SchemeModule::SchemeModule()
 
                               ));
 
-  this->_addBuiltinFunction(
-    "glace",
-    std::make_shared<BuiltinFunctionEmbedder<
-      void(std::shared_ptr<const IMesh>, const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-           const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
-
-      [](std::shared_ptr<const IMesh> p_mesh,
-         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-         const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id, const FunctionSymbolId& p_id) -> void {
-        switch (p_mesh->dimension()) {
-        case 1: {
-          AcousticSolverAlgorithm<1, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
-          break;
-        }
-        case 2: {
-          AcousticSolverAlgorithm<2, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
-          break;
-        }
-        case 3: {
-          AcousticSolverAlgorithm<3, AcousticSolverType::Glace>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
-          break;
-        }
-        default: {
-          throw UnexpectedError("invalid mesh dimension");
-        }
-        }
-      }
-
-      ));
-
-  this->_addBuiltinFunction(
-    "eucclhyd",
-    std::make_shared<BuiltinFunctionEmbedder<
-      void(std::shared_ptr<const IMesh>, const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
-           const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
-
-      [](std::shared_ptr<const IMesh> p_mesh,
-         const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-         const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id, const FunctionSymbolId& p_id) -> void {
-        switch (p_mesh->dimension()) {
-        case 1: {
-          AcousticSolverAlgorithm<1, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
-          break;
-        }
-        case 2: {
-          AcousticSolverAlgorithm<2, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
-          break;
-        }
-        case 3: {
-          AcousticSolverAlgorithm<3, AcousticSolverType::Eucclhyd>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
-          break;
-        }
-        default: {
-          throw UnexpectedError("invalid mesh dimension");
-        }
-        }
-      }
-
-      ));
+  this->_addBuiltinFunction("glace",
+                            std::make_shared<BuiltinFunctionEmbedder<
+                              void(std::shared_ptr<const IMesh>,
+                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
+
+                              [](std::shared_ptr<const IMesh> p_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id,
+                                 const FunctionSymbolId& p_id) -> void {
+                                switch (p_mesh->dimension()) {
+                                case 1: {
+                                  AcousticSolverAlgorithm<1>{AcousticSolverType::Glace,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
+                                  break;
+                                }
+                                case 2: {
+                                  AcousticSolverAlgorithm<2>{AcousticSolverType::Glace,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
+                                  break;
+                                }
+                                case 3: {
+                                  AcousticSolverAlgorithm<3>{AcousticSolverType::Glace,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
+                                  break;
+                                }
+                                default: {
+                                  throw UnexpectedError("invalid mesh dimension");
+                                }
+                                }
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("eucclhyd",
+                            std::make_shared<BuiltinFunctionEmbedder<
+                              void(std::shared_ptr<const IMesh>,
+                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
+
+                              [](std::shared_ptr<const IMesh> p_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id,
+                                 const FunctionSymbolId& p_id) -> void {
+                                switch (p_mesh->dimension()) {
+                                case 1: {
+                                  AcousticSolverAlgorithm<1>{AcousticSolverType::Eucclhyd,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
+                                  break;
+                                }
+                                case 2: {
+                                  AcousticSolverAlgorithm<2>{AcousticSolverType::Eucclhyd,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
+                                  break;
+                                }
+                                case 3: {
+                                  AcousticSolverAlgorithm<3>{AcousticSolverType::Eucclhyd,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
+                                  break;
+                                }
+                                default: {
+                                  throw UnexpectedError("invalid mesh dimension");
+                                }
+                                }
+                              }
+
+                              ));
 }
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 2dfd3234486b352e5b9f27e13e10001846008bc1..7cd1ed89023e619af478625b58441002ecae0af2 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -90,6 +90,7 @@ target_link_libraries (unit_tests
   PugsLanguage
   PugsLanguageAST
   PugsLanguageModules
+  PugsLanguageAlgorithms
   PugsLanguageUtils
   PugsMesh
   PugsUtils