diff --git a/config.sh b/config.sh
index d9dc4e6..6a2557f 100755
--- a/config.sh
+++ b/config.sh
@@ -11,7 +11,8 @@ MPI_PATH=/opt/homebrew
 P4EST_PATH=${HOME}/packages/p4est/p4est/build/local
 
 # PETSC_PATH : Path to PETSc install (i.e., ${PETSC_PATH}/include, ${PETSC_PATH}/lib, ...) 
-PETSC_PATH=${HOME}/packages/petsc/petsc-build-mumps
+# PETSC_PATH=${HOME}/packages/petsc/petsc-build-mumps
+PETSC_PATH=/opt/homebrew
 
 # BUILD_DIR : Build directory
 BUILD_DIR=build-$(git branch --show-current)
diff --git a/examples/elliptic-single/main.cpp b/examples/elliptic-single/main.cpp
index bca056b..6f71005 100644
--- a/examples/elliptic-single/main.cpp
+++ b/examples/elliptic-single/main.cpp
@@ -127,12 +127,6 @@ int main(int argc, char** argv) {
     int n_solves = 1;
     app.options.setOption("n-solves", n_solves);
     
-    int min_level = 0;
-    app.options.setOption("min-level", min_level);
-    
-    int max_level = 6;
-    app.options.setOption("max-level", max_level);
-
     double x_lower = -10.0;
     app.options.setOption("x-lower", x_lower);
 
@@ -145,10 +139,19 @@ int main(int argc, char** argv) {
     double y_upper = 10.0;
     app.options.setOption("y-upper", y_upper);
     
-    int nx = 8;
+    int min_level = 0;
+    int max_level = 5;
+    int nx = 16;
+    int ny = 16;
+    if (argc > 1) {
+        min_level = atoi(argv[1]);
+        max_level = atoi(argv[2]);
+        nx = atoi(argv[3]);
+        ny = atoi(argv[4]);
+    }
+    app.options.setOption("max-level", max_level);
+    app.options.setOption("min-level", min_level);
     app.options.setOption("nx", nx);
-    
-    int ny = 8;
     app.options.setOption("ny", ny);
 
     // ====================================================
diff --git a/src/EllipticForestApp.cpp b/src/EllipticForestApp.cpp
index 98e405b..d7470e4 100644
--- a/src/EllipticForestApp.cpp
+++ b/src/EllipticForestApp.cpp
@@ -2,12 +2,13 @@
 
 namespace EllipticForest {
 
-EllipticForestApp::EllipticForestApp() :
-    argc_(nullptr), argv_(nullptr)
-        {}
+EllipticForestApp::EllipticForestApp()
+    {}
     
 EllipticForestApp::EllipticForestApp(int* argc, char*** argv) :
-    argc_(argc), argv_(argv) {
+    argc_{argc},
+    argv_{argv},
+    parser{*argc, *argv} {
     
     addTimer("app-lifetime");
     timers["app-lifetime"].start();
@@ -19,14 +20,10 @@ EllipticForestApp::EllipticForestApp(int* argc, char*** argv) :
     PetscInitialize(argc_, argv_, NULL, NULL);
     PetscGetArgs(argc_, argv_);
 #endif
-
-    // Create options
-    InputParser inputParser(*argc_, *argv_);
-    inputParser.parse(options);
     
-    int myRank = -1;
-    MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
-    if (myRank == 0) {
+    int mpi_rank = -1;
+    MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
+    if (mpi_rank == 0) {
         std::cout << "[EllipticForest] Welcome to EllipticForest!" << std::endl;
     }
     this->actualClassPointer_ = this;
@@ -35,15 +32,31 @@ EllipticForestApp::EllipticForestApp(int* argc, char*** argv) :
 EllipticForestApp::~EllipticForestApp() {
     timers["app-lifetime"].stop();
 
-    int myRank = 0;
+    int mpi_rank = 0;
+    int mpi_size = 0;
     int isMPIFinalized;
     MPI_Finalized(&isMPIFinalized);
     if (!isMPIFinalized) {
-        MPI_Comm_rank(MPI_COMM_WORLD, &myRank);
+        MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
+        MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
     }
     
-    if (myRank == 0) {
+    // std::map<std::string, std::vector<double>> global_timers;
+    // for (auto& [key, value] : timers) {
+    //     if (mpi_rank == 0) {
+    //         global_timers[key] = std::vector<double>(mpi_size);
+    //     }
+    //     double time = value.time();
+    //     MPI::gather<double>(time, global_timers[key], 0, MPI_COMM_WORLD);
+    // }
+
+    if (mpi_rank == 0) {
+        // std::string filename = "ef-timers-" + getCurrentDateTimeString() + ".csv";
+        // writeMapToCSV(global_timers, filename);
+
         std::cout << "[EllipticForest] End of app life cycle, finalizing..." << std::endl;
+        std::cout << "[EllipticForest] Options:" << std::endl << options;
+        
         std::cout << "[EllipticForest] Timers: " << std::endl;
         for (auto& [key, value] : timers) {
             std::cout << "[EllipticForest]   " << key << " : " << value.time() << " [sec]" << std::endl;
diff --git a/src/EllipticForestApp.hpp b/src/EllipticForestApp.hpp
index 45a4dd5..072856c 100644
--- a/src/EllipticForestApp.hpp
+++ b/src/EllipticForestApp.hpp
@@ -17,11 +17,13 @@
 #include <petsc.h>
 #include <p4est.h>
 
+#include "MPI.hpp"
 #include "Logger.hpp"
 #include "Options.hpp"
 #include "InputParser.hpp"
 #include "Timer.hpp"
 #include "GenericSingleton.hpp"
+#include "Helpers.hpp"
 
 namespace EllipticForest {
 
@@ -52,6 +54,8 @@ class EllipticForestApp : public GenericSingleton<EllipticForestApp> {
      */
     Options options{};
 
+    InputParser parser{};
+
     /**
      * @brief Timer instance for managing timers
      * 
@@ -126,13 +130,13 @@ class EllipticForestApp : public GenericSingleton<EllipticForestApp> {
      * @brief Pointer to argc from main
      * 
      */
-    int* argc_;
+    int* argc_ = nullptr;
 
     /**
      * @brief Pointer to argv from main
      * 
      */
-    char*** argv_;
+    char*** argv_ = nullptr;
 
 };
 
diff --git a/src/HPSAlgorithm.hpp b/src/HPSAlgorithm.hpp
index c5456a7..9b0d7b7 100644
--- a/src/HPSAlgorithm.hpp
+++ b/src/HPSAlgorithm.hpp
@@ -123,10 +123,13 @@ class HPSAlgorithm : public MPI::MPIObject {
         EllipticForestApp& app = EllipticForestApp::getInstance();
         app.logHead("Begin HPS Build Stage");
         app.addTimer("build-stage");
-        app.timers["build-stage"].start();
+        app.addTimer("build-stage:leaf-callback");
+        app.addTimer("build-stage:family-callback");
 
+        app.timers["build-stage"].start();
         mesh.quadtree.merge(
             [&](Node<PatchType>* leaf_node){
+                app.timers["build-stage:leaf-callback"].start();
                 // app.log("Leaf callback: path = " + leaf_node->path);
                 // Leaf callback
                 PatchType& patch = leaf_node->data;
@@ -140,9 +143,12 @@ class HPSAlgorithm : public MPI::MPIObject {
                 else {
                     patch.matrixT() = patch_solver.buildD2N(patch.grid());
                 }
+
+                app.timers["build-stage:leaf-callback"].stop();
                 return 1;
             },
             [&](Node<PatchType>* parent_node, std::vector<Node<PatchType>*> child_nodes){
+                app.timers["build-stage:family-callback"].start();
                 // app.log("Family callback: parent path = " + parent_node->path);
                 // Family callback
                 PatchType& tau = parent_node->data;
@@ -151,11 +157,13 @@ class HPSAlgorithm : public MPI::MPIObject {
                 PatchType& gamma = child_nodes[2]->data;
                 PatchType& omega = child_nodes[3]->data;
                 merge4to1(tau, alpha, beta, gamma, omega);
+
+                app.timers["build-stage:family-callback"].stop();
                 return 1;
             }
         );
-
         app.timers["build-stage"].stop();
+
         app.logHead("End HPS Build Stage");
 
     }
@@ -181,10 +189,13 @@ class HPSAlgorithm : public MPI::MPIObject {
         EllipticForestApp& app = EllipticForestApp::getInstance();
         app.logHead("Begin HPS Upwards Stage");
         app.addTimer("upwards-stage");
-        app.timers["upwards-stage"].start();
+        app.addTimer("upwards-stage:leaf-callback");
+        app.addTimer("upwards-stage:family-callback");
 
+        app.timers["upwards-stage"].start();
         mesh.quadtree.merge(
             [&](Node<PatchType>* leaf_node){
+                app.timers["upwards-stage:leaf-callback"].start();
                 // app.log("Leaf callback: path = " + leaf_node->path);
                 // Leaf callback
                 PatchType& patch = leaf_node->data;
@@ -195,9 +206,11 @@ class HPSAlgorithm : public MPI::MPIObject {
                 // Set particular Neumann data using patch solver function
                 patch.vectorH() = patch_solver.particularNeumannData(patch.grid(), patch.vectorF());
 
+                app.timers["upwards-stage:leaf-callback"].stop();
                 return 1;
             },
             [&](Node<PatchType>* parent_node, std::vector<Node<PatchType>*> child_nodes){
+                app.timers["upwards-stage:family-callback"].start();
                 // app.log("Family callback: parent path = " + parent_node->path);
                 // Family callback
                 PatchType& tau = parent_node->data;
@@ -206,11 +219,12 @@ class HPSAlgorithm : public MPI::MPIObject {
                 PatchType& gamma = child_nodes[2]->data;
                 PatchType& omega = child_nodes[3]->data;
                 upwards4to1(tau, alpha, beta, gamma, omega);
+                app.timers["upwards-stage:family-callback"].stop();
                 return 1;
             }
         );
-
         app.timers["upwards-stage"].stop();
+
         app.logHead("End HPS Upwards Stage");
 
     }
@@ -228,10 +242,13 @@ class HPSAlgorithm : public MPI::MPIObject {
         EllipticForestApp& app = EllipticForestApp::getInstance();
         app.logHead("Begin HPS Upwards Stage");
         app.addTimer("upwards-stage");
-        app.timers["upwards-stage"].start();
+        app.addTimer("upwards-stage:leaf-callback");
+        app.addTimer("upwards-stage:family-callback");
 
+        app.timers["upwards-stage"].start();
         mesh.quadtree.merge(
             [&](Node<PatchType>* leaf_node){
+                app.timers["upwards-stage:leaf-callback"].start();
                 // app.log("Leaf callback: path = " + leaf_node->path);
                 // Leaf callback
                 PatchType& patch = leaf_node->data;
@@ -251,9 +268,11 @@ class HPSAlgorithm : public MPI::MPIObject {
                 // Set particular Neumann data using patch solver function
                 patch.vectorH() = patch_solver.particularNeumannData(grid, patch.vectorF());
 
+                app.timers["upwards-stage:leaf-callback"].stop();
                 return 1;
             },
             [&](Node<PatchType>* parent_node, std::vector<Node<PatchType>*> child_nodes){
+                app.timers["upwards-stage:family-callback"].start();
                 // app.log("Family callback: parent path = " + parent_node->path);
                 // Family callback
                 PatchType& tau = parent_node->data;
@@ -262,11 +281,12 @@ class HPSAlgorithm : public MPI::MPIObject {
                 PatchType& gamma = child_nodes[2]->data;
                 PatchType& omega = child_nodes[3]->data;
                 upwards4to1(tau, alpha, beta, gamma, omega);
+                app.timers["upwards-stage:family-callback"].stop();
                 return 1;
             }
         );
-
         app.timers["upwards-stage"].stop();
+
         app.logHead("End HPS Upwards Stage");
 
     }
@@ -293,19 +313,24 @@ class HPSAlgorithm : public MPI::MPIObject {
         EllipticForestApp& app = EllipticForestApp::getInstance();
         app.logHead("Begin HPS Solve Stage");
         app.addTimer("solve-stage");
-        app.timers["solve-stage"].start();
+        app.addTimer("solve-stage:leaf-callback");
+        app.addTimer("solve-stage:family-callback");
 
         // Set Dirichlet data on root patch
         boundary_data_function(mesh.quadtree.root());
 
+        app.timers["solve-stage"].start();
         mesh.quadtree.split(
             [&](Node<PatchType>* leaf_node){
+                app.timers["solve-stage:leaf-callback"].start();
                 // app.log("Leaf callback: path = " + leaf_node->path);
                 // Leaf callback
                 leafSolve(leaf_node->data);
+                app.timers["solve-stage:leaf-callback"].stop();
                 return 1;
             },
             [&](Node<PatchType>* parent_node, std::vector<Node<PatchType>*> child_nodes){
+                app.timers["solve-stage:family-callback"].start();
                 // app.log("Family callback: parent path = " + parent_node->path);
                 // Family callback
                 PatchType& tau = parent_node->data;
@@ -314,11 +339,12 @@ class HPSAlgorithm : public MPI::MPIObject {
                 PatchType& gamma = child_nodes[2]->data;
                 PatchType& omega = child_nodes[3]->data;
                 split1to4(tau, alpha, beta, gamma, omega);
+                app.timers["solve-stage:family-callback"].stop();
                 return 1;
             }
         );
-
         app.timers["solve-stage"].stop();
+
         app.logHead("End HPS Solve Stage");
 
     }
@@ -346,7 +372,8 @@ class HPSAlgorithm : public MPI::MPIObject {
         EllipticForestApp& app = EllipticForestApp::getInstance();
         app.logHead("Begin HPS Solve Stage");
         app.addTimer("solve-stage");
-        app.timers["solve-stage"].start();
+        app.addTimer("solve-stage:leaf-callback");
+        app.addTimer("solve-stage:family-callback");
 
         // Set up data for Dirichlet solve
         PatchType& rootPatch = mesh.quadtree.root();
@@ -419,14 +446,18 @@ class HPSAlgorithm : public MPI::MPIObject {
             g = solve(A, rhs);
         }
 
+        app.timers["solve-stage"].start();
         mesh.quadtree.split(
             [&](Node<PatchType>* leaf_node){
+                app.timers["solve-stage:leaf-callback"].start();
                 // app.log("Leaf callback: path = " + leaf_node->path);
                 // Leaf callback
                 leafSolve(leaf_node->data);
+                app.timers["solve-stage:leaf-callback"].stop();
                 return 1;
             },
             [&](Node<PatchType>* parent_node, std::vector<Node<PatchType>*> child_nodes){
+                app.timers["solve-stage:family-callback"].start();
                 // app.log("Family callback: parent path = " + parent_node->path);
                 // Family callback
                 PatchType& tau = parent_node->data;
@@ -435,11 +466,12 @@ class HPSAlgorithm : public MPI::MPIObject {
                 PatchType& gamma = child_nodes[2]->data;
                 PatchType& omega = child_nodes[3]->data;
                 split1to4(tau, alpha, beta, gamma, omega);
+                app.timers["solve-stage:family-callback"].stop();
                 return 1;
             }
         );
-
         app.timers["solve-stage"].stop();
+
         app.logHead("End HPS Solve Stage");
                 
     }
diff --git a/src/Helpers.cpp b/src/Helpers.cpp
new file mode 100644
index 0000000..27dc11f
--- /dev/null
+++ b/src/Helpers.cpp
@@ -0,0 +1,33 @@
+#include "Helpers.hpp"
+
+namespace EllipticForest {
+
+std::string getCurrentDateTimeString() {
+    std::time_t now = std::time(nullptr);
+    std::tm timeinfo = *std::localtime(&now);
+
+    std::ostringstream oss;
+    oss << std::put_time(&timeinfo, "%Y%m%d%H%M%S");
+    return oss.str();
+}
+
+void writeMapToCSV(const std::map<std::string, std::vector<double>>& data, const std::string& filename) {
+    std::ofstream file(filename);
+    if (!file.is_open()) {
+        std::cerr << "Error opening file " << filename << std::endl;
+        return;
+    }
+
+    for (const auto& entry : data) {
+        file << entry.first; // Writing the key to the file
+        const auto& vec = entry.second;
+        for (const auto& value : vec) {
+            file << "," << value; // Writing values in the vector
+        }
+        file << std::endl; // End of row
+    }
+
+    file.close();
+}
+
+} // NAMESPACE : EllipticForest
\ No newline at end of file
diff --git a/src/Helpers.hpp b/src/Helpers.hpp
new file mode 100644
index 0000000..23e5295
--- /dev/null
+++ b/src/Helpers.hpp
@@ -0,0 +1,21 @@
+#ifndef HELPERS_HPP_
+#define HELPERS_HPP_
+
+#include <iostream>
+#include <iomanip>
+#include <ctime>
+#include <sstream>
+#include <fstream>
+#include <map>
+#include <vector>
+#include <string>
+
+namespace EllipticForest {
+
+std::string getCurrentDateTimeString();
+
+void writeMapToCSV(const std::map<std::string, std::vector<double>>& data, const std::string& filename);
+
+} // NAMESPACE : EllipticForest
+
+#endif // HELPERS_HPP_
\ No newline at end of file
diff --git a/src/InputParser.cpp b/src/InputParser.cpp
index 3bad9cb..ce79d12 100644
--- a/src/InputParser.cpp
+++ b/src/InputParser.cpp
@@ -2,6 +2,9 @@
 
 namespace EllipticForest {
 
+InputParser::InputParser()
+    {}
+
 InputParser::InputParser(int& argc, char**& argv) :
     args{argv+1, argv + argc}
         {}
diff --git a/src/InputParser.hpp b/src/InputParser.hpp
index 00c7361..8ed2d0b 100644
--- a/src/InputParser.hpp
+++ b/src/InputParser.hpp
@@ -21,7 +21,9 @@ class InputParser {
      * @brief List of command-line arguments
      * 
      */
-    std::vector<std::string> args;
+    std::vector<std::string> args{};
+
+    InputParser();
 
     /**
      * @brief Construct a new Input Parser object
diff --git a/src/MPI.hpp b/src/MPI.hpp
index 28d374b..197e3aa 100644
--- a/src/MPI.hpp
+++ b/src/MPI.hpp
@@ -386,6 +386,11 @@ int broadcast(std::vector<T>& data, int root, Communicator comm) {
     return MPI_Bcast(TypeTraits<std::vector<T>>::getAddress(data), TypeTraits<std::vector<T>>::getSize(data), TypeTraits<std::vector<T>>::getType(data), root, comm);
 }
 
+template<class T>
+int gather(T& data, std::vector<T>& root_data, int root, Communicator comm) {
+    return MPI_Gather(&data, 1, TypeTraits<T>::getType(data), TypeTraits<std::vector<T>>::getAddress(root_data), 1, TypeTraits<std::vector<T>>::getType(root_data), root, comm);
+}
+
 /**
  * @brief Function overload of @sa `allgather` for std::vector<T>
  * 
diff --git a/src/Options.cpp b/src/Options.cpp
index c15d3cb..e418f78 100644
--- a/src/Options.cpp
+++ b/src/Options.cpp
@@ -13,6 +13,16 @@ Options::OptionTypes& Options::operator[](std::string const key) {
     return optionsMap[key];
 }
 
+void Options::addOption(std::string option_name, OptionTypes& value, std::string info) {
+
+    // Set option in the map
+    optionsMap[option_name] = value;
+
+    // Check for command line override of option
+    
+
+}
+
 void Options::setOption(std::string const key, OptionTypes const value) {
     optionsMap[key] = value;
 }
diff --git a/src/Options.hpp b/src/Options.hpp
index 2b821af..40f3c89 100644
--- a/src/Options.hpp
+++ b/src/Options.hpp
@@ -26,6 +26,8 @@ class Options {
      */
     using OptionTypes = std::variant<std::string, bool, int, double>;
 
+    std::vector<std::string> args;
+
     /**
      * @brief Map of options
      * 
@@ -53,6 +55,8 @@ class Options {
      */
     OptionTypes& operator[](std::string const key);
 
+    void addOption(std::string option_name, OptionTypes& value, std::string info);
+
     /**
      * @brief Set an option by key-value pair
      *