diff --git a/Applications/Utils/MeshGeoTools/AssignRasterDataToMesh.cpp b/Applications/Utils/MeshGeoTools/AssignRasterDataToMesh.cpp
index 5bcb11317435f8fa5080e6986e90f81efe566b57..d92e9956b7eb906f6c78a8a0cf812390a127735c 100644
--- a/Applications/Utils/MeshGeoTools/AssignRasterDataToMesh.cpp
+++ b/Applications/Utils/MeshGeoTools/AssignRasterDataToMesh.cpp
@@ -13,6 +13,10 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
+#ifdef USE_PETSC
+#include <mpi.h>
+#endif
+
 #include "Applications/FileIO/AsciiRasterInterface.h"
 #include "GeoLib/Raster.h"
 #include "InfoLib/GitInfo.h"
@@ -70,6 +74,10 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
+#ifdef USE_PETSC
+    MPI_Init(&argc, &argv);
+#endif
+
     bool const create_cell_array(set_cells_arg.isSet());
     bool const create_node_array =
         (create_cell_array) ? set_nodes_arg.isSet() : true;
@@ -83,6 +91,9 @@ int main(int argc, char* argv[])
     if (mesh->getDimension() > 2)
     {
         ERR("Method can not be applied to 3D meshes.");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
 
@@ -97,6 +108,9 @@ int main(int argc, char* argv[])
         if (!assigned)
         {
             ERR("Error assigning raster data to scalar node array");
+#ifdef USE_PETSC
+            MPI_Finalize();
+#endif
             return EXIT_FAILURE;
         }
         INFO("Created node array {:s}", array_name_arg.getValue());
@@ -110,6 +124,9 @@ int main(int argc, char* argv[])
         if (!assigned)
         {
             ERR("Error assigning raster data to scalar cell array");
+#ifdef USE_PETSC
+            MPI_Finalize();
+#endif
             return EXIT_FAILURE;
         }
         INFO("Created cell array {:s}", array_name_arg.getValue());
@@ -117,5 +134,8 @@ int main(int argc, char* argv[])
 
     MeshLib::IO::VtuInterface vtu(mesh.get());
     vtu.writeToFile(output_name);
+#ifdef USE_PETSC
+    MPI_Finalize();
+#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp b/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp
index 0ece0942cd4cb5d072dd6863e42c1d23b6d67bb6..8a229c25b0dc62958f50279035830cb740c25b14 100644
--- a/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp
+++ b/Applications/Utils/MeshGeoTools/IntegrateBoreholesIntoMesh.cpp
@@ -16,6 +16,10 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
+#ifdef USE_PETSC
+#include <mpi.h>
+#endif
+
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
 #include "InfoLib/GitInfo.h"
@@ -118,6 +122,10 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
+#ifdef USE_PETSC
+    MPI_Init(&argc, &argv);
+#endif
+
     std::pair<int, int> mat_limits(0, std::numeric_limits<int>::max());
     std::pair<double, double> elevation_limits(
         std::numeric_limits<double>::lowest(), dmax);
@@ -126,6 +134,9 @@ int main(int argc, char* argv[])
     {
         ERR("If minimum MaterialID is set, maximum ID must be set, too (and "
             "vice versa).");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
     if (min_id_arg.isSet() && max_id_arg.isSet())
@@ -140,12 +151,18 @@ int main(int argc, char* argv[])
     if (min_id_arg.isSet() && (mat_limits.first < 0 || mat_limits.second < 0))
     {
         ERR("Specified MaterialIDs must have non-negative values.");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
     if (min_elevation_arg.isSet() != max_elevation_arg.isSet())
     {
         ERR("If minimum elevation is set, maximum elevation must be set, too "
             "(and vice versa).");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
     if (min_elevation_arg.isSet() && max_elevation_arg.isSet())
@@ -167,6 +184,9 @@ int main(int argc, char* argv[])
     if (!xml_io.readFile(geo_name))
     {
         ERR("Failed to read geometry file `{:s}'.", geo_name);
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
     std::vector<GeoLib::Point*> const& points =
@@ -177,11 +197,17 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         ERR("Failed to read input mesh file `{:s}'.", mesh_name);
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
     if (mesh->getDimension() != 3)
     {
         ERR("Method can only be applied to 3D meshes.");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
 
@@ -190,6 +216,9 @@ int main(int argc, char* argv[])
     if (mat_ids == nullptr)
     {
         ERR("Mesh is required to have MaterialIDs");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
 
@@ -226,5 +255,8 @@ int main(int argc, char* argv[])
     MeshLib::Mesh const result("result", new_nodes, new_elems, props);
     MeshLib::IO::VtuInterface vtu(&result);
     vtu.writeToFile(output_name);
+#ifdef USE_PETSC
+    MPI_Finalize();
+#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/Raster2Mesh.cpp b/Applications/Utils/MeshGeoTools/Raster2Mesh.cpp
index 150ac100bf035355e51fba5cca1b2b36ecb54a96..5d7ec4094273a2fc9b92732d26c5af7dc7d4e78a 100644
--- a/Applications/Utils/MeshGeoTools/Raster2Mesh.cpp
+++ b/Applications/Utils/MeshGeoTools/Raster2Mesh.cpp
@@ -13,6 +13,10 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
+#ifdef USE_PETSC
+#include <mpi.h>
+#endif
+
 #include "Applications/FileIO/AsciiRasterInterface.h"
 #include "GeoLib/Raster.h"
 #include "InfoLib/GitInfo.h"
@@ -68,6 +72,10 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
+#ifdef USE_PETSC
+    MPI_Init(&argc, &argv);
+#endif
+
     std::string const input_name = input_arg.getValue().c_str();
     std::string const output_name = output_arg.getValue().c_str();
 
@@ -98,10 +106,16 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         ERR("Conversion failed.");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
 
     MeshLib::IO::VtuInterface vtu(mesh.get());
     vtu.writeToFile(output_name);
+#ifdef USE_PETSC
+    MPI_Finalize();
+#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp b/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
index 7e602968ef739cf812c9aa27240b5ff8088756bd..5cf029d440d912f4e3bf5166843c87de5a2ad12d 100644
--- a/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
+++ b/Applications/Utils/MeshGeoTools/VerticalSliceFromLayers.cpp
@@ -16,6 +16,10 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
+#ifdef USE_PETSC
+#include <mpi.h>
+#endif
+
 #include <QCoreApplication>
 
 #include "Applications/FileIO/Gmsh/GMSHInterface.h"
@@ -409,6 +413,10 @@ int main(int argc, char* argv[])
     cmd.add(input_arg);
     cmd.parse(argc, argv);
 
+#ifdef USE_PETSC
+    MPI_Init(&argc, &argv);
+#endif
+
     std::string const input_name = input_arg.getValue();
     std::string const output_name = output_arg.getValue();
 
@@ -426,6 +434,9 @@ int main(int argc, char* argv[])
     if (layer_names.size() < 2)
     {
         ERR("At least two layers are required to extract a slice.");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
 
@@ -437,6 +448,9 @@ int main(int argc, char* argv[])
     {
         ERR("Less than two geometries could be created from layers. Aborting "
             "extraction...");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
 
@@ -456,6 +470,9 @@ int main(int argc, char* argv[])
     if (mesh == nullptr)
     {
         ERR("Error generating mesh... (GMSH was unable to output mesh)");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
     rotateMesh(*mesh, rot_mat, z_shift);
@@ -463,6 +480,9 @@ int main(int argc, char* argv[])
     if (new_mesh == nullptr)
     {
         ERR("Error generating mesh... (GMSH created line mesh)");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
 
@@ -478,5 +498,8 @@ int main(int argc, char* argv[])
     }
     MeshLib::IO::VtuInterface vtu(revised_mesh.get());
     vtu.writeToFile(output_name + ".vtu");
+#ifdef USE_PETSC
+    MPI_Finalize();
+#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp b/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp
index 5435cf3997a5757ea9e141e92d4deeaaa3bba33e..118377c8afa72022d9d6cac8309af17185f4d997 100644
--- a/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp
+++ b/Applications/Utils/MeshGeoTools/computeSurfaceNodeIDsInPolygonalRegion.cpp
@@ -12,6 +12,10 @@
 
 #include <tclap/CmdLine.h>
 
+#ifdef USE_PETSC
+#include <mpi.h>
+#endif
+
 #include <algorithm>
 #include <memory>
 #include <string>
@@ -96,6 +100,10 @@ int main(int argc, char* argv[])
     cmd.add(gmsh_path_arg);
     cmd.parse(argc, argv);
 
+#ifdef USE_PETSC
+    MPI_Init(&argc, &argv);
+#endif
+
     std::unique_ptr<MeshLib::Mesh const> mesh(
         MeshLib::IO::readMeshFromFile(mesh_in.getValue()));
     INFO("Mesh read: {:d} nodes, {:d} elements.", mesh->getNumberOfNodes(),
@@ -183,5 +191,8 @@ int main(int argc, char* argv[])
     std::for_each(all_sfc_nodes.begin(), all_sfc_nodes.end(),
                   std::default_delete<MeshLib::Node>());
 
-    return 0;
+#ifdef USE_PETSC
+    MPI_Finalize();
+#endif
+    return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/createIntermediateRasters.cpp b/Applications/Utils/MeshGeoTools/createIntermediateRasters.cpp
index bba4b5aa590b09eb67e5b95eb22ff117a7584e21..1ecc2200a5ce7e8555166ab4fe7d5674c7aaccdd 100644
--- a/Applications/Utils/MeshGeoTools/createIntermediateRasters.cpp
+++ b/Applications/Utils/MeshGeoTools/createIntermediateRasters.cpp
@@ -9,6 +9,10 @@
 
 #include <tclap/CmdLine.h>
 
+#ifdef USE_PETSC
+#include <mpi.h>
+#endif
+
 #include <memory>
 #include <string>
 
@@ -46,6 +50,10 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
+#ifdef USE_PETSC
+    MPI_Init(&argc, &argv);
+#endif
+
     std::unique_ptr<GeoLib::Raster> dem1(
         FileIO::AsciiRasterInterface::readRaster(input1_arg.getValue()));
     std::unique_ptr<GeoLib::Raster> dem2(
@@ -53,6 +61,9 @@ int main(int argc, char* argv[])
 
     if (dem1 == nullptr || dem2 == nullptr)
     {
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return 1;
     }
 
@@ -88,6 +99,9 @@ int main(int argc, char* argv[])
 
     if (errors_found)
     {
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return 2;
     }
 
@@ -106,6 +120,9 @@ int main(int argc, char* argv[])
         if (it2 == dem2->end())
         {
             ERR("Error: File 2 is shorter than File 1.");
+#ifdef USE_PETSC
+            MPI_Finalize();
+#endif
             return 3;
         }
         if (*it1 == h1.no_data || *it2 == h2.no_data)
@@ -130,6 +147,9 @@ int main(int argc, char* argv[])
     if (it2 != dem2->end())
     {
         ERR("Error: File 1 is shorter than File 2.");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return 3;
     }
 
@@ -144,5 +164,8 @@ int main(int argc, char* argv[])
             r, basename + std::to_string(i) + ext);
         INFO("Layer {:d} written.", i + 1);
     }
+#ifdef USE_PETSC
+    MPI_Finalize();
+#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp b/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp
index f25b3d7241912af4d9807797254875c1a9dfe1f3..8062d4be8a56414370bd5c102674b1aab30c363b 100644
--- a/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp
+++ b/Applications/Utils/MeshGeoTools/geometryToGmshGeo.cpp
@@ -12,6 +12,10 @@
 // ThirdParty
 #include <tclap/CmdLine.h>
 
+#ifdef USE_PETSC
+#include <mpi.h>
+#endif
+
 #include "Applications/FileIO/Gmsh/GMSHInterface.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
@@ -66,6 +70,10 @@ int main(int argc, char* argv[])
 
     cmd.parse(argc, argv);
 
+#ifdef USE_PETSC
+    MPI_Init(&argc, &argv);
+#endif
+
     GeoLib::GEOObjects geo_objects;
     for (auto const& geometry_name : geo_input_arg.getValue())
     {
@@ -74,6 +82,9 @@ int main(int argc, char* argv[])
         {
             if (!xml.readFile(geometry_name))
             {
+#ifdef USE_PETSC
+                MPI_Finalize();
+#endif
                 return EXIT_FAILURE;
             }
         }
@@ -81,6 +92,9 @@ int main(int argc, char* argv[])
         {
             ERR("Failed to read file '{:s}'.", geometry_name);
             ERR("{:s}", err.what());
+#ifdef USE_PETSC
+            MPI_Finalize();
+#endif
             return EXIT_FAILURE;
         }
         INFO("Successfully read file '{:s}'.", geometry_name);
@@ -141,7 +155,13 @@ int main(int argc, char* argv[])
                 "for better analysis of the problem.");
         }
         ERR("An error has occurred - programme will be terminated.");
+#ifdef USE_PETSC
+        MPI_Finalize();
+#endif
         return EXIT_FAILURE;
     }
+#ifdef USE_PETSC
+    MPI_Finalize();
+#endif
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/MeshGeoTools/identifySubdomains.cpp b/Applications/Utils/MeshGeoTools/identifySubdomains.cpp
index 19aeb30c9a385e1565330f5d20f289f71470eed1..3191231e08f43f4382e82261c38c785f4e4d1f44 100644
--- a/Applications/Utils/MeshGeoTools/identifySubdomains.cpp
+++ b/Applications/Utils/MeshGeoTools/identifySubdomains.cpp
@@ -9,6 +9,10 @@
 
 #include <tclap/CmdLine.h>
 
+#ifdef USE_PETSC
+#include <mpi.h>
+#endif
+
 #include "InfoLib/GitInfo.h"
 #include "MeshGeoToolsLib/IdentifySubdomainMesh.h"
 #include "MeshGeoToolsLib/MeshNodeSearcher.h"
@@ -89,6 +93,10 @@ int main(int argc, char* argv[])
     cmd.add(subdomain_meshes_filenames_arg);
     cmd.parse(argc, argv);
 
+#ifdef USE_PETSC
+    MPI_Init(&argc, &argv);
+#endif
+
     //
     // The bulk mesh.
     //
@@ -139,5 +147,8 @@ int main(int argc, char* argv[])
             output_prefix_arg.getValue() + mesh_ptr->getName() + ".vtu");
     }
 
+#ifdef USE_PETSC
+    MPI_Finalize();
+#endif
     return EXIT_SUCCESS;
 }