diff --git a/src/language/PugsParser.cpp b/src/language/PugsParser.cpp
index 27bc4e1d0b228d2c629f44891f15f23059216d31..b13c77149975ac211a8e0f2f19794cb02248c434 100644
--- a/src/language/PugsParser.cpp
+++ b/src/language/PugsParser.cpp
@@ -16,8 +16,9 @@
 #include <language/utils/ASTDotPrinter.hpp>
 #include <language/utils/ASTPrinter.hpp>
 #include <language/utils/SymbolTable.hpp>
-#include <utils/Exceptions.hpp>
 #include <utils/PugsAssert.hpp>
+#include <utils/PugsUtils.hpp>
+#include <utils/SignalManager.hpp>
 
 #include <pegtl/contrib/analyze.hpp>
 #include <pegtl/contrib/parse_tree.hpp>
@@ -41,10 +42,8 @@ parser(const std::string& filename)
   std::cout << rang::style::bold << "Parsing file " << rang::style::reset << rang::style::underline << filename
             << rang::style::reset << " ...\n";
 
-  std::unique_ptr<ASTNode> root_node;
-  read_input input(filename);
-  try {
-    root_node = ASTBuilder::build(input);
+  auto parse_and_execute = [](auto& input) {
+    std::unique_ptr<ASTNode> root_node = ASTBuilder::build(input);
 
     ASTModulesImporter{*root_node};
     ASTNodeTypeCleaner<language::import_instruction>{*root_node};
@@ -82,18 +81,27 @@ parser(const std::string& filename)
     ExecutionPolicy exec_all;
     root_node->execute(exec_all);
     std::cout << *(root_node->m_symbol_table) << '\n';
-  }
-  catch (const parse_error& e) {
-    const auto p = e.positions.front();
-
-    std::ostringstream os;
-    os << rang::style::bold << p.source << ':' << p.line << ':' << p.byte_in_line << ": " << rang::style::reset
-       << rang::fgB::red << "error: " << rang::fg::reset << rang::style::bold << e.what() << rang::style::reset << '\n'
-       << input.line_at(p) << '\n'
-       << std::string(p.byte_in_line, ' ') << rang::fgB::yellow << '^' << rang::fg::reset;
+  };
 
-    throw RawError(os.str());
+  if (not SignalManager::pauseOnError()) {
+    read_input input(filename);
+    try {
+      parse_and_execute(input);
+    }
+    catch (const parse_error& e) {
+      const auto p = e.positions.front();
+
+      std::cerr << rang::style::bold << p.source << ':' << p.line << ':' << p.byte_in_line << ": " << rang::style::reset
+                << rang::fgB::red << "error: " << rang::fg::reset << rang::style::bold << e.what() << rang::style::reset
+                << '\n'
+                << input.line_at(p) << '\n'
+                << std::string(p.byte_in_line, ' ') << rang::fgB::yellow << '^' << rang::fg::reset << '\n';
+      finalize();
+      std::exit(1);
+    }
+  } else {
+    read_input input(filename);
+    parse_and_execute(input);
   }
-
   std::cout << "Parsed: " << filename << '\n';
 }
diff --git a/src/language/node_processor/BuiltinFunctionProcessor.hpp b/src/language/node_processor/BuiltinFunctionProcessor.hpp
index 6d688284b4ab456fb6b4b5b9bc60d803c2992584..61150be71154be2d08f4a1376df2df13681793d5 100644
--- a/src/language/node_processor/BuiltinFunctionProcessor.hpp
+++ b/src/language/node_processor/BuiltinFunctionProcessor.hpp
@@ -5,7 +5,6 @@
 #include <language/node_processor/FunctionArgumentConverter.hpp>
 #include <language/node_processor/INodeProcessor.hpp>
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
-#include <utils/Exceptions.hpp>
 
 class BuiltinFunctionExpressionProcessor final : public INodeProcessor
 {
@@ -63,7 +62,7 @@ class BuiltinFunctionProcessor : public INodeProcessor
     try {
       return m_function_expression_processor->execute(context_exec_policy);
     }
-    catch (NormalError& e) {
+    catch (std::runtime_error& e) {
       throw parse_error(e.what(), {m_argument_node.begin()});
     }
   }
diff --git a/src/main.cpp b/src/main.cpp
index 1bde1fe585ae2cc5a7498dbd0530bddf9ccd9fcc..47bbbc7258533ed8aefbeea1639e8afdd0646c93 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -39,395 +39,382 @@
 int
 main(int argc, char* argv[])
 {
-  try {
-    std::string filename = initialize(argc, argv);
-
-    std::regex gmsh_regex("(.*).msh");
-    if (not std::regex_match(filename, gmsh_regex)) {
-      SynchronizerManager::create();
-      parser(filename);
-      SynchronizerManager::destroy();
-      finalize();
-      return 0;
-    }
-
-    std::map<std::string, double> method_cost_map;
+  std::string filename = initialize(argc, argv);
 
+  std::regex gmsh_regex("(.*).msh");
+  if (not std::regex_match(filename, gmsh_regex)) {
     SynchronizerManager::create();
+    parser(filename);
+    SynchronizerManager::destroy();
+    finalize();
+    return 0;
+  }
 
-    if (filename != "") {
-      std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
-      Timer gmsh_timer;
-      gmsh_timer.reset();
-      GmshReader gmsh_reader(filename);
-      method_cost_map["Mesh building"] = gmsh_timer.seconds();
-
-      std::shared_ptr<IMesh> p_mesh = gmsh_reader.mesh();
-
-      switch (p_mesh->dimension()) {
-      case 1: {
-        std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
-        std::vector<std::shared_ptr<IBoundaryConditionDescriptor>> bc_descriptor_list;
-        for (const auto& sym_boundary_name : sym_boundary_name_list) {
-          std::shared_ptr<IBoundaryDescriptor> boudary_descriptor =
-            std::shared_ptr<IBoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
-          SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
-            new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
-          bc_descriptor_list.push_back(std::shared_ptr<IBoundaryConditionDescriptor>(sym_bc_descriptor));
-        }
+  std::map<std::string, double> method_cost_map;
+
+  SynchronizerManager::create();
+
+  if (filename != "") {
+    std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
+    Timer gmsh_timer;
+    gmsh_timer.reset();
+    GmshReader gmsh_reader(filename);
+    method_cost_map["Mesh building"] = gmsh_timer.seconds();
 
-        using ConnectivityType = Connectivity1D;
-        using MeshType         = Mesh<ConnectivityType>;
-        using MeshDataType     = MeshData<MeshType>;
-        using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
-
-        const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
-
-        Timer timer;
-        timer.reset();
-        MeshDataType mesh_data(mesh);
-
-        std::vector<BoundaryConditionHandler> bc_list;
-        {
-          for (const auto& bc_descriptor : bc_descriptor_list) {
-            switch (bc_descriptor->type()) {
-            case IBoundaryConditionDescriptor::Type::symmetry: {
-              const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
-                dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-              for (size_t i_ref_node_list = 0;
-                   i_ref_node_list < mesh.connectivity().numberOfRefItemList<ItemType::node>(); ++i_ref_node_list) {
-                const RefNodeList& ref_node_list = mesh.connectivity().refItemList<ItemType::node>(i_ref_node_list);
-                const RefId& ref                 = ref_node_list.refId();
-                if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-                  SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
-                    new SymmetryBoundaryCondition<MeshType::Dimension>(
-                      MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_node_list));
-                  std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
-                  bc_list.push_back(BoundaryConditionHandler(bc));
-                }
+    std::shared_ptr<IMesh> p_mesh = gmsh_reader.mesh();
+
+    switch (p_mesh->dimension()) {
+    case 1: {
+      std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
+      std::vector<std::shared_ptr<IBoundaryConditionDescriptor>> bc_descriptor_list;
+      for (const auto& sym_boundary_name : sym_boundary_name_list) {
+        std::shared_ptr<IBoundaryDescriptor> boudary_descriptor =
+          std::shared_ptr<IBoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+        SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
+          new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+        bc_descriptor_list.push_back(std::shared_ptr<IBoundaryConditionDescriptor>(sym_bc_descriptor));
+      }
+
+      using ConnectivityType = Connectivity1D;
+      using MeshType         = Mesh<ConnectivityType>;
+      using MeshDataType     = MeshData<MeshType>;
+      using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
+
+      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+      Timer timer;
+      timer.reset();
+      MeshDataType mesh_data(mesh);
+
+      std::vector<BoundaryConditionHandler> bc_list;
+      {
+        for (const auto& bc_descriptor : bc_descriptor_list) {
+          switch (bc_descriptor->type()) {
+          case IBoundaryConditionDescriptor::Type::symmetry: {
+            const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+              dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+            for (size_t i_ref_node_list = 0;
+                 i_ref_node_list < mesh.connectivity().numberOfRefItemList<ItemType::node>(); ++i_ref_node_list) {
+              const RefNodeList& ref_node_list = mesh.connectivity().refItemList<ItemType::node>(i_ref_node_list);
+              const RefId& ref                 = ref_node_list.refId();
+              if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+                SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
+                  new SymmetryBoundaryCondition<MeshType::Dimension>(
+                    MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_node_list));
+                std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
+                bc_list.push_back(BoundaryConditionHandler(bc));
               }
-              break;
-            }
-            default: {
-              throw UnexpectedError("Unknown BCDescription\n");
-            }
             }
+            break;
+          }
+          default: {
+            throw UnexpectedError("Unknown BCDescription\n");
+          }
           }
         }
+      }
 
-        UnknownsType unknowns(mesh_data);
-
-        unknowns.initializeSod();
-
-        AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+      UnknownsType unknowns(mesh_data);
 
-        using Rd = TinyVector<MeshType::Dimension>;
+      unknowns.initializeSod();
 
-        const CellValue<const double>& Vj = mesh_data.Vj();
+      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
 
-        const double tmax = 0.2;
-        double t          = 0;
+      using Rd = TinyVector<MeshType::Dimension>;
 
-        int itermax   = std::numeric_limits<int>::max();
-        int iteration = 0;
+      const CellValue<const double>& Vj = mesh_data.Vj();
 
-        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();
+      const double tmax = 0.2;
+      double t          = 0;
 
-        BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+      int itermax   = std::numeric_limits<int>::max();
+      int iteration = 0;
 
-        VTKWriter vtk_writer("mesh", 0.01);
+      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();
 
-        while ((t < tmax) and (iteration < itermax)) {
-          vtk_writer.write(mesh,
-                           {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                            NamedItemValue{"coords", mesh.xr()},
-                            NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                            NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                           t);
-          double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
-          if (t + dt > tmax) {
-            dt = tmax - t;
-          }
-          acoustic_solver.computeNextStep(t, dt, unknowns);
+      BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
-          block_eos.updatePandCFromRhoE();
+      VTKWriter vtk_writer("mesh", 0.01);
 
-          t += dt;
-          ++iteration;
-        }
+      while ((t < tmax) and (iteration < itermax)) {
         vtk_writer.write(mesh,
                          {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                           NamedItemValue{"coords", mesh.xr()},
                           NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                           NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                         t, true);   // forces last output
-
-        std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
-                  << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+                         t);
+        double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
+        if (t + dt > tmax) {
+          dt = tmax - t;
+        }
+        acoustic_solver.computeNextStep(t, dt, unknowns);
 
-        method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+        block_eos.updatePandCFromRhoE();
 
-        {   // gnuplot output for density
-          const CellValue<const Rd>& xj       = mesh_data.xj();
-          const CellValue<const double>& rhoj = unknowns.rhoj();
-          std::ofstream fout("rho");
-          for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
-            fout << xj[j][0] << ' ' << rhoj[j] << '\n';
-          }
+        t += dt;
+        ++iteration;
+      }
+      vtk_writer.write(mesh,
+                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                        NamedItemValue{"coords", mesh.xr()},
+                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                       t, true);   // forces last output
+
+      std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
+                << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+
+      method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+
+      {   // gnuplot output for density
+        const CellValue<const Rd>& xj       = mesh_data.xj();
+        const CellValue<const double>& rhoj = unknowns.rhoj();
+        std::ofstream fout("rho");
+        for (CellId j = 0; j < mesh.numberOfCells(); ++j) {
+          fout << xj[j][0] << ' ' << rhoj[j] << '\n';
         }
+      }
 
-        break;
+      break;
+    }
+    case 2: {
+      // test case boundary condition description
+      std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX"};
+      std::vector<std::shared_ptr<IBoundaryConditionDescriptor>> bc_descriptor_list;
+      for (const auto& sym_boundary_name : sym_boundary_name_list) {
+        std::shared_ptr<IBoundaryDescriptor> boudary_descriptor =
+          std::shared_ptr<IBoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+        SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
+          new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+        bc_descriptor_list.push_back(std::shared_ptr<IBoundaryConditionDescriptor>(sym_bc_descriptor));
       }
-      case 2: {
-        // test case boundary condition description
-        std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX"};
-        std::vector<std::shared_ptr<IBoundaryConditionDescriptor>> bc_descriptor_list;
-        for (const auto& sym_boundary_name : sym_boundary_name_list) {
-          std::shared_ptr<IBoundaryDescriptor> boudary_descriptor =
-            std::shared_ptr<IBoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
-          SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
-            new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
-          bc_descriptor_list.push_back(std::shared_ptr<IBoundaryConditionDescriptor>(sym_bc_descriptor));
-        }
 
-        using ConnectivityType = Connectivity2D;
-        using MeshType         = Mesh<ConnectivityType>;
-        using MeshDataType     = MeshData<MeshType>;
-        using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
-
-        const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
-
-        Timer timer;
-        timer.reset();
-        MeshDataType mesh_data(mesh);
-
-        std::vector<BoundaryConditionHandler> bc_list;
-        {
-          for (const auto& bc_descriptor : bc_descriptor_list) {
-            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().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
-                const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
-                const RefId& ref                 = ref_face_list.refId();
-                if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-                  SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
-                    new SymmetryBoundaryCondition<MeshType::Dimension>(
-                      MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
-                  std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
-                  bc_list.push_back(BoundaryConditionHandler(bc));
-                }
+      using ConnectivityType = Connectivity2D;
+      using MeshType         = Mesh<ConnectivityType>;
+      using MeshDataType     = MeshData<MeshType>;
+      using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
+
+      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+      Timer timer;
+      timer.reset();
+      MeshDataType mesh_data(mesh);
+
+      std::vector<BoundaryConditionHandler> bc_list;
+      {
+        for (const auto& bc_descriptor : bc_descriptor_list) {
+          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().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
+              const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
+              const RefId& ref                 = ref_face_list.refId();
+              if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+                SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
+                  new SymmetryBoundaryCondition<MeshType::Dimension>(
+                    MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
+                std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
+                bc_list.push_back(BoundaryConditionHandler(bc));
               }
-              break;
-            }
-            default: {
-              throw UnexpectedError("Unknown BCDescription\n");
-            }
             }
+            break;
+          }
+          default: {
+            throw UnexpectedError("Unknown BCDescription\n");
+          }
           }
         }
+      }
 
-        UnknownsType unknowns(mesh_data);
-
-        unknowns.initializeSod();
-
-        AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+      UnknownsType unknowns(mesh_data);
 
-        const CellValue<const double>& Vj = mesh_data.Vj();
+      unknowns.initializeSod();
 
-        const double tmax = 0.2;
-        double t          = 0;
+      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
 
-        int itermax   = std::numeric_limits<int>::max();
-        int iteration = 0;
+      const CellValue<const double>& Vj = mesh_data.Vj();
 
-        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();
+      const double tmax = 0.2;
+      double t          = 0;
 
-        BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+      int itermax   = std::numeric_limits<int>::max();
+      int iteration = 0;
 
-        VTKWriter vtk_writer("mesh", 0.01);
+      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();
 
-        while ((t < tmax) and (iteration < itermax)) {
-          vtk_writer.write(mesh,
-                           {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                            NamedItemValue{"coords", mesh.xr()},
-                            NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                            NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                           t);
-          double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
-          if (t + dt > tmax) {
-            dt = tmax - t;
-          }
-          acoustic_solver.computeNextStep(t, dt, unknowns);
+      BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
-          block_eos.updatePandCFromRhoE();
+      VTKWriter vtk_writer("mesh", 0.01);
 
-          t += dt;
-          ++iteration;
-        }
+      while ((t < tmax) and (iteration < itermax)) {
         vtk_writer.write(mesh,
                          {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                           NamedItemValue{"coords", mesh.xr()},
                           NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                           NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                         t, true);   // forces last output
+                         t);
+        double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
+        if (t + dt > tmax) {
+          dt = tmax - t;
+        }
+        acoustic_solver.computeNextStep(t, dt, unknowns);
 
-        std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
-                  << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+        block_eos.updatePandCFromRhoE();
 
-        method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
-        break;
+        t += dt;
+        ++iteration;
+      }
+      vtk_writer.write(mesh,
+                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                        NamedItemValue{"coords", mesh.xr()},
+                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                       t, true);   // forces last output
+
+      std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
+                << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+
+      method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+      break;
+    }
+    case 3: {
+      std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
+      std::vector<std::shared_ptr<IBoundaryConditionDescriptor>> bc_descriptor_list;
+      for (const auto& sym_boundary_name : sym_boundary_name_list) {
+        std::shared_ptr<IBoundaryDescriptor> boudary_descriptor =
+          std::shared_ptr<IBoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+        SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
+          new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+        bc_descriptor_list.push_back(std::shared_ptr<IBoundaryConditionDescriptor>(sym_bc_descriptor));
       }
-      case 3: {
-        std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
-        std::vector<std::shared_ptr<IBoundaryConditionDescriptor>> bc_descriptor_list;
-        for (const auto& sym_boundary_name : sym_boundary_name_list) {
-          std::shared_ptr<IBoundaryDescriptor> boudary_descriptor =
-            std::shared_ptr<IBoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
-          SymmetryBoundaryConditionDescriptor* sym_bc_descriptor =
-            new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
-          bc_descriptor_list.push_back(std::shared_ptr<IBoundaryConditionDescriptor>(sym_bc_descriptor));
-        }
 
-        using ConnectivityType = Connectivity3D;
-        using MeshType         = Mesh<ConnectivityType>;
-        using MeshDataType     = MeshData<MeshType>;
-        using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
-
-        const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
-
-        Timer timer;
-        timer.reset();
-        MeshDataType mesh_data(mesh);
-
-        std::vector<BoundaryConditionHandler> bc_list;
-        {
-          for (const auto& bc_descriptor : bc_descriptor_list) {
-            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().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
-                const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
-                const RefId& ref                 = ref_face_list.refId();
-                if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-                  SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
-                    new SymmetryBoundaryCondition<MeshType::Dimension>(
-                      MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
-                  std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
-                  bc_list.push_back(BoundaryConditionHandler(bc));
-                }
+      using ConnectivityType = Connectivity3D;
+      using MeshType         = Mesh<ConnectivityType>;
+      using MeshDataType     = MeshData<MeshType>;
+      using UnknownsType     = FiniteVolumesEulerUnknowns<MeshDataType>;
+
+      const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
+
+      Timer timer;
+      timer.reset();
+      MeshDataType mesh_data(mesh);
+
+      std::vector<BoundaryConditionHandler> bc_list;
+      {
+        for (const auto& bc_descriptor : bc_descriptor_list) {
+          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().numberOfRefItemList<ItemType::face>(); ++i_ref_face_list) {
+              const RefFaceList& ref_face_list = mesh.connectivity().refItemList<ItemType::face>(i_ref_face_list);
+              const RefId& ref                 = ref_face_list.refId();
+              if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+                SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
+                  new SymmetryBoundaryCondition<MeshType::Dimension>(
+                    MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list));
+                std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
+                bc_list.push_back(BoundaryConditionHandler(bc));
               }
-              break;
-            }
-            default: {
-              throw UnexpectedError("Unknown BCDescription\n");
-            }
             }
+            break;
+          }
+          default: {
+            throw UnexpectedError("Unknown BCDescription\n");
+          }
           }
         }
+      }
 
-        UnknownsType unknowns(mesh_data);
-
-        unknowns.initializeSod();
+      UnknownsType unknowns(mesh_data);
 
-        AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
+      unknowns.initializeSod();
 
-        const CellValue<const double>& Vj = mesh_data.Vj();
+      AcousticSolver<MeshDataType> acoustic_solver(mesh_data, bc_list);
 
-        const double tmax = 0.2;
-        double t          = 0;
+      const CellValue<const double>& Vj = mesh_data.Vj();
 
-        int itermax   = std::numeric_limits<int>::max();
-        int iteration = 0;
+      const double tmax = 0.2;
+      double t          = 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();
+      int itermax   = std::numeric_limits<int>::max();
+      int iteration = 0;
 
-        BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+      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();
 
-        VTKWriter vtk_writer("mesh", 0.01);
+      BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
 
-        while ((t < tmax) and (iteration < itermax)) {
-          vtk_writer.write(mesh,
-                           {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                            NamedItemValue{"coords", mesh.xr()},
-                            NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
-                            NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                           t);
-          double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
-          if (t + dt > tmax) {
-            dt = tmax - t;
-          }
-          acoustic_solver.computeNextStep(t, dt, unknowns);
-          block_eos.updatePandCFromRhoE();
+      VTKWriter vtk_writer("mesh", 0.01);
 
-          t += dt;
-          ++iteration;
-        }
+      while ((t < tmax) and (iteration < itermax)) {
         vtk_writer.write(mesh,
                          {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
                           NamedItemValue{"coords", mesh.xr()},
                           NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
                           NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
-                         t, true);   // forces last output
-
-        std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
-                  << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+                         t);
+        double dt = 0.4 * acoustic_solver.acoustic_dt(Vj, cj);
+        if (t + dt > tmax) {
+          dt = tmax - t;
+        }
+        acoustic_solver.computeNextStep(t, dt, unknowns);
+        block_eos.updatePandCFromRhoE();
 
-        method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
-        break;
-      }
+        t += dt;
+        ++iteration;
       }
-
-      std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':'
-                << __LINE__ << ")\n";
-
-    } else {
-      throw NormalError("Connectivity1D defined by number of nodes no more implemented\n");
+      vtk_writer.write(mesh,
+                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                        NamedItemValue{"coords", mesh.xr()},
+                        NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+                        NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},
+                       t, true);   // forces last output
+
+      std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset << ":  " << rang::fgB::green
+                << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+
+      method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+      break;
+    }
     }
 
-    SynchronizerManager::destroy();
+    std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':'
+              << __LINE__ << ")\n";
 
-    finalize();
+  } else {
+    throw NormalError("Connectivity1D defined by number of nodes no more implemented\n");
+  }
 
-    std::string::size_type size = 0;
-    for (const auto& method_cost : method_cost_map) {
-      size = std::max(size, method_cost.first.size());
-    }
+  SynchronizerManager::destroy();
 
-    for (const auto& method_cost : method_cost_map) {
-      std::cout << "* [" << rang::fgB::cyan << std::setw(size) << std::left << method_cost.first << rang::fg::reset
-                << "] Execution time: " << rang::style::bold << method_cost.second << rang::style::reset << '\n';
-    }
+  finalize();
+
+  std::string::size_type size = 0;
+  for (const auto& method_cost : method_cost_map) {
+    size = std::max(size, method_cost.first.size());
   }
-  catch (const IExitError& e) {
-    if (SignalManager::pauseOnError()) {
-      std::rethrow_exception(std::current_exception());
-    } else {
-      // Each failing process must write
-      std::cerr.setstate(std::ios::goodbit);
-      std::cerr << e.what() << '\n';
-
-      return 1;
-    }
+
+  for (const auto& method_cost : method_cost_map) {
+    std::cout << "* [" << rang::fgB::cyan << std::setw(size) << std::left << method_cost.first << rang::fg::reset
+              << "] Execution time: " << rang::style::bold << method_cost.second << rang::style::reset << '\n';
   }
 
   return 0;
diff --git a/src/utils/PugsUtils.cpp b/src/utils/PugsUtils.cpp
index ec024977f7cd21b7ebaca795f379b0eea6ad06ec..de9e8830fbbc200d724f3f49bf9d44db82e3438b 100644
--- a/src/utils/PugsUtils.cpp
+++ b/src/utils/PugsUtils.cpp
@@ -48,7 +48,7 @@ initialize(int& argc, char* argv[])
   {
     CLI::App app{"Pugs help"};
 
-    app.add_option("filename,-f,--filename", filename, "gmsh file");
+    app.add_option("filename,-f,--filename", filename, "pugs script file")->check(CLI::ExistingFile);
 
     int threads = -1;
     app.add_option("--threads", threads, "Number of Kokkos threads")