diff --git a/src/main.cpp b/src/main.cpp
index e351d6471a50bb133286a1b53bbd9bc7ee12f1dc..e31ab068ed1f41f064cec799bb0de0f222bdd6fb 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -28,356 +28,349 @@
 int main(int argc, char *argv[])
 {
   std::string filename = initialize(argc, argv);
-
   std::map<std::string, double> method_cost_map;
 
-  try  {
-    if (filename != "") {
-      pout() << "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->meshDimension()) {
-        case 1: {
-          std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
-          std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
-          for (const auto& sym_boundary_name : sym_boundary_name_list){
-            std::shared_ptr<BoundaryDescriptor> boudary_descriptor
-                = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
-            SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
-                = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
-            bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
-          }
+  if (filename != "") {
+    pout() << "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->meshDimension()) {
+      case 1: {
+        std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
+        std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
+        for (const auto& sym_boundary_name : sym_boundary_name_list){
+          std::shared_ptr<BoundaryDescriptor> boudary_descriptor
+              = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+          SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
+              = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+          bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(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 BoundaryConditionDescriptor::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().numberOfRefNodeList();
-                       ++i_ref_node_list) {
-                    const RefNodeList& ref_node_list = mesh.connectivity().refNodeList(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));
-                    }
+        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 BoundaryConditionDescriptor::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().numberOfRefNodeList();
+                     ++i_ref_node_list) {
+                  const RefNodeList& ref_node_list = mesh.connectivity().refNodeList(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: {
-                  perr() << "Unknown BCDescription\n";
-                  std::exit(1);
                 }
+                break;
+              }
+              default: {
+                perr() << "Unknown BCDescription\n";
+                std::exit(1);
               }
             }
           }
+        }
 
-          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
-
-          pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
-                    << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+                                  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);
 
-          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
+
+        pout() << "* " << 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<BoundaryConditionDescriptor>> bc_descriptor_list;
+        for (const auto& sym_boundary_name : sym_boundary_name_list){
+          std::shared_ptr<BoundaryDescriptor> boudary_descriptor
+              = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+          SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
+              = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+          bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(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<BoundaryConditionDescriptor>> bc_descriptor_list;
-          for (const auto& sym_boundary_name : sym_boundary_name_list){
-            std::shared_ptr<BoundaryDescriptor> boudary_descriptor
-                = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
-            SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
-                = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
-            bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(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 BoundaryConditionDescriptor::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().numberOfRefFaceList();
-                       ++i_ref_face_list) {
-                    const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(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 BoundaryConditionDescriptor::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().numberOfRefFaceList();
+                     ++i_ref_face_list) {
+                  const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(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: {
-                  perr() << "Unknown BCDescription\n";
-                  std::exit(1);
                 }
+                break;
+              }
+              default: {
+                perr() << "Unknown BCDescription\n";
+                std::exit(1);
               }
             }
           }
+        }
 
-          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);
+        VTKWriter vtk_writer("mesh", 0.01);
 
-            block_eos.updatePandCFromRhoE();
-
-            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
+                                  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);
 
-          pout() << "* " << 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;
         }
-        case 3: {
-          std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX", "YMIN", "YMAX", "ZMIN", "ZMAX"};
-          std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
-          for (const auto& sym_boundary_name : sym_boundary_name_list){
-            std::shared_ptr<BoundaryDescriptor> boudary_descriptor
-                = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
-            SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
-                = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
-
-            bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
-          }
+        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
+
+        pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
+               << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
-          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 BoundaryConditionDescriptor::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().numberOfRefFaceList();
-                       ++i_ref_face_list) {
-                    const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(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));
-                    }
+        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<BoundaryConditionDescriptor>> bc_descriptor_list;
+        for (const auto& sym_boundary_name : sym_boundary_name_list){
+          std::shared_ptr<BoundaryDescriptor> boudary_descriptor
+              = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+          SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
+              = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+          bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(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 BoundaryConditionDescriptor::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().numberOfRefFaceList();
+                     ++i_ref_face_list) {
+                  const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(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: {
-                  perr() << "Unknown BCDescription\n";
-                  std::exit(1);
                 }
+                break;
+              }
+              default: {
+                perr() << "Unknown BCDescription\n";
+                std::exit(1);
               }
             }
           }
+        }
 
-          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
-
-          pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
-                    << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
+                                  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();
 
-          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
 
-      pout() << "* "  << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n";
+        pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
+               << ":  " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
 
-    } else {
-      perr() << "Connectivity1D defined by number of nodes no more implemented\n";
-      std::exit(0);
+        method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
+        break;
+      }
     }
-  }
-  catch (const AssertError& error) {
-    perr() << error << '\n';
-    std::exit(1);
+
+    pout() << "* "  << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n";
+
+  } else {
+    perr() << "Connectivity1D defined by number of nodes no more implemented\n";
+    std::exit(0);
   }
 
   finalize();
@@ -389,14 +382,14 @@ int main(int argc, char *argv[])
 
   for (const auto& method_cost : method_cost_map) {
     pout() << "* ["
-              << 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';
+           << 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/PastisAssert.hpp b/src/utils/PastisAssert.hpp
index 6bcab320db2daad8624407f8a4741395df1d9600..c2f2ee440cbb5ec701a6ef8f1759611b9d3f9387 100644
--- a/src/utils/PastisAssert.hpp
+++ b/src/utils/PastisAssert.hpp
@@ -7,12 +7,19 @@
 #include <iostream>
 #include <string>
 
+template <typename ErrorType>
+void printAndThrow(const ErrorType& error)
+{
+  throw error;
+}
+
 class AssertError
 {
  private:
   const std::string m_file;
   const int m_line;
   const std::string m_function;
+  const std::string m_test;
   const std::string m_message;
 
  public:
@@ -22,11 +29,16 @@ class AssertError
   {
     os << '\n'
        << rang::style::bold
-       << "*** Assertion error ***\n"
+       << "---------- Assertion error -----------\n"
        << " at " << assert_error.m_file << ':' << assert_error.m_line << '\n'
        << " in " << assert_error.m_function << '\n'
-       << "*** " << rang::fgB::red << assert_error.m_message << rang::fg::reset
-       << rang::style::reset << '\n';
+       << " assertion (" << rang::fgB::red << assert_error.m_test << rang::fg::reset
+       << ") failed!\n";
+    if (not assert_error.m_message.empty()) {
+      os << ' ' << rang::fgB::yellow << assert_error.m_message
+         << rang::fg::reset << '\n';
+    }
+    os << "--------------------------------------" << rang::style::reset << '\n';
 
     return os;
   }
@@ -35,10 +47,12 @@ class AssertError
   AssertError(std::string file,
               int line,
               std::string function,
-              std::string message)
+              std::string test,
+              std::string message="")
       : m_file(file),
         m_line(line),
         m_function(function),
+        m_test(test),
         m_message(message)
   {
     ;
@@ -59,17 +73,32 @@ PRAGMA_DIAGNOSTIC_POP
 #ifdef NDEBUG
 
 // Useless test is there to check syntax even in optimized mode. Costs nothing.
-#define Assert(assertion)                       \
-  if (not _pastis_assert(assertion)) {}
+#define Assert(assertion,...)                                   \
+  if (not _pastis_assert(assertion)) {                          \
+    using vargs_t = decltype(std::make_tuple(__VA_ARGS__));     \
+    static_assert(std::tuple_size_v<vargs_t> <= 1,              \
+                  "too many arguments");                        \
+  }
 
 #else // NDEBUG
 
-#define Assert(assertion)                       \
-  if (not _pastis_assert(assertion)) {          \
-    throw AssertError(__FILE__,                 \
-                      __LINE__,                 \
-                      __PRETTY_FUNCTION__,      \
-                      (#assertion));            \
+#define Assert(assertion,...)                                   \
+  if (not _pastis_assert(assertion)) {                          \
+    using vargs_t = decltype(std::make_tuple(__VA_ARGS__));     \
+    static_assert(std::tuple_size_v<vargs_t> <= 1,              \
+                  "too many arguments");                        \
+    if constexpr(std::tuple_size_v<vargs_t> == 0) {             \
+        printAndThrow(AssertError(__FILE__,                     \
+                                  __LINE__,                     \
+                                  __PRETTY_FUNCTION__,          \
+                                  #assertion));                 \
+      } else {                                                  \
+      printAndThrow(AssertError(__FILE__,                       \
+                                __LINE__,                       \
+                                __PRETTY_FUNCTION__,            \
+                                #assertion,                     \
+                                #__VA_ARGS__));                 \
+    }                                                           \
   }
 
 #endif // NDEBUG
diff --git a/src/utils/SignalManager.cpp b/src/utils/SignalManager.cpp
index 09a1c5ab8f9a6b56f58e61cd645ba3da86db8a5d..cd3c4251029b7dcc3c9aaee81627b96deea4e40f 100644
--- a/src/utils/SignalManager.cpp
+++ b/src/utils/SignalManager.cpp
@@ -1,5 +1,6 @@
 #include <SignalManager.hpp>
 
+#include <PastisAssert.hpp>
 #include <PastisOStream.hpp>
 
 #include <BacktraceManager.hpp>
@@ -60,6 +61,19 @@ void SignalManager::handler(int signal)
   std::signal(SIGINT,  SIG_DFL);
   std::signal(SIGABRT, SIG_DFL);
 
+  std::exception_ptr eptr = std::current_exception();
+  try {
+    if (eptr) {
+      std::rethrow_exception(eptr);
+    }
+  }
+  catch(const AssertError& assert_error) {
+    perr() << assert_error << '\n';
+  }
+  catch(...) {
+    perr() << "Unknown exception!\n";
+  }
+
   perr() << "\n *** "
 	    << rang::style::reset
 	    << rang::fg::reset
diff --git a/tests/test_PastisAssert.cpp b/tests/test_PastisAssert.cpp
index d55ad8d826d75ac10e81b5b110495e74d131e921..90995c30f9782baa8f2c0468909983ac166e3a9c 100644
--- a/tests/test_PastisAssert.cpp
+++ b/tests/test_PastisAssert.cpp
@@ -8,10 +8,22 @@ TEST_CASE("PastisAssert", "[utils]") {
     const std::string filename = "filename";
     const int line = 10;
     const std::string function = "function";
+    const std::string test = "test";
+
+    AssertError assert_error(filename, line, function, test);
+
+    REQUIRE(Catch::Detail::stringify(assert_error) == "\n---------- Assertion error -----------\n at filename:10\n in function\n assertion (test) failed!\n--------------------------------------\n");
+  }
+
+  SECTION("checking for assert error with message") {
+    const std::string filename = "filename";
+    const int line = 10;
+    const std::string function = "function";
+    const std::string test = "test";
     const std::string message = "message";
 
-    AssertError assert_error(filename, line, function, message);
+    AssertError assert_error(filename, line, function, test, message);
 
-    REQUIRE(Catch::Detail::stringify(assert_error) == "\n*** Assertion error ***\n at filename:10\n in function\n*** message\n");
+    REQUIRE(Catch::Detail::stringify(assert_error) == "\n---------- Assertion error -----------\n at filename:10\n in function\n assertion (test) failed!\n message\n--------------------------------------\n");
   }
 }