diff --git a/include/constants.h b/include/constants.h
index dac2bfbb..25f757cb 100644
--- a/include/constants.h
+++ b/include/constants.h
@@ -45,7 +45,14 @@ enum icontrols
     optimized,
     /*! \brief Verbose level
          *
-         * Defines the level of verbose outputs.
+         * Defines the level of verbose outputs. For higher levels all output from lower levels is always given.
+         * 0: No output
+         * 1: Minimal convergence output.
+         * 2: Minimal problem setting details, iteration output and timings.
+         * 3: Print all parameters and more detailed problem settings.
+         * 4: Information on called functions.
+         * 5: Prints information on multigrid levels.
+         * 6: Print out all kind of variable values.
          */
     verbose,
     /*! \brief openmp paralellization
diff --git a/include/gyro.h b/include/gyro.h
index 69876085..1ce578dd 100644
--- a/include/gyro.h
+++ b/include/gyro.h
@@ -171,12 +171,40 @@ void trafo(double r_i, std::vector<double> theta, std::vector<double> sin_theta,
            int ntheta, std::vector<double>& x, std::vector<double>& y, int verbose);
 
 /***************************************************************************
-     * Display arrays or vectors
-     **************************************************************************/
-void disp(int na, double* a, const std::string& s_a);
-void disp(std::vector<double> a, const std::string& s_a);
-void disp(int na, int* a, const std::string& s_a);
-void disp(std::vector<int> a, const std::string& s_a);
+* Display arrays or vectors
+**************************************************************************/
+/*!
+ *  \brief Forwards all values inside an array to std::cout.
+ * 
+ * \param na: Size of the array.
+ * \param a: The array.
+ * \param s_a: Name of the array (to be printed).
+ *
+ */
+template <typename T>
+void disp(int na, T* a, const std::string& s_a)
+{
+    std::cout << s_a << "(" << na << "): ";
+    for (int i = 0; i < na; i++)
+        std::cout << a[i] << " ";
+    std::cout << "\n";
+} 
+
+/*!
+ *  \brief Forwards all values inside a vector to std::cout.
+ * 
+ * \param a: The vector.
+ * \param s_a: Name of the vector (to be printed).
+ *
+ */
+template <typename T>
+void disp(std::vector<T> a, const std::string& s_a)
+{
+    std::cout << s_a << "(" << a.size() << "): ";
+    for (std::size_t i = 0; i < a.size(); i++)
+        std::cout << a[i] << " ";
+    std::cout << "\n";
+} 
 
 /***************************************************************************
      * Matrix operations
diff --git a/src/build_bi_aniso_rdir_phiper.cpp b/src/build_bi_aniso_rdir_phiper.cpp
index 3a1f48ad..9aa78b52 100644
--- a/src/build_bi_aniso_rdir_phiper.cpp
+++ b/src/build_bi_aniso_rdir_phiper.cpp
@@ -627,8 +627,6 @@ std::vector<double> level::apply_prolongation_inj0(std::vector<double> u, int mc
         }
     }
 
-    //std::cout << "end of applyP_inj \n";
-
     return Pu;
 } /* ----- end of level::apply_prolongation_inj0 ----- */
 
diff --git a/src/create_grid_polar.cpp b/src/create_grid_polar.cpp
index 4349619f..2d348e76 100644
--- a/src/create_grid_polar.cpp
+++ b/src/create_grid_polar.cpp
@@ -37,8 +37,8 @@
  */
 void gmgpolar::create_grid_polar()
 {
-    if (gyro::icntl[Param::verbose] > 2)
-        std::cout << "Setup for the problem...\n";
+    if (gyro::icntl[Param::verbose] > 3)
+        std::cout << "Creating polar grid...\n";
 
     level* new_level = new level(0);
     v_level.push_back(new_level);
@@ -73,7 +73,7 @@ void gmgpolar::create_grid_polar()
         std::cout << "on the coordinates (r x theta): (" << gyro::dcntl[Param::R0] << ", " << gyro::dcntl[Param::R]
                   << ") x (" << gyro::dcntl[Param::THETA0] << ", " << gyro::dcntl[Param::THETA] << ")\n";
     }
-    if (gyro::icntl[Param::verbose] > 3) {
+    if (gyro::icntl[Param::verbose] > 5) {
         v_level[0]->display_r();
         v_level[0]->display_theta();
         std::cout << "\n";
@@ -243,7 +243,9 @@ void level::build_theta()
             theta[i] = fac * i;
     }
     else {
-        std::cout << "Anisotropy chosen in theta.\n";
+        if (gyro::icntl[Param::verbose] > 2) {
+            std::cout << "Anisotropy chosen in theta.\n";
+        }
         ntheta     = pow(2, ceil(log2(nr)) - 1);
         double fac = 2 * PI / ntheta;
 
@@ -280,8 +282,8 @@ void level::build_theta()
  */
 void gmgpolar::create_grid_polar_divide()
 {
-    if (gyro::icntl[Param::verbose] > 2)
-        std::cout << "dividing a coarser grid...\n";
+    if (gyro::icntl[Param::verbose] > 3)
+        std::cout << "Dividing a coarser grid...\n";
 
     int divide = gyro::icntl[Param::divideBy2];
     std::vector<double> r_tmp;
diff --git a/src/define_coarse_nodes.cpp b/src/define_coarse_nodes.cpp
index 94d23b0b..4b152031 100644
--- a/src/define_coarse_nodes.cpp
+++ b/src/define_coarse_nodes.cpp
@@ -51,13 +51,13 @@ void gmgpolar::define_coarse_nodes()
     }
     levels_orig = levels;
     for (int i = 0; i < levels; i++) {
-        if (gyro::icntl[Param::verbose] > 1)
+        if (gyro::icntl[Param::verbose] > 3)
             std::cout << "=> Defining coarse nodes on level " << i << "\n";
         v_level[i]->store_theta_n_co();
         if (i < levels - 1) {
             v_level[i]->define_coarse_nodes_onelevel(v_level[i + 1]);
         }
-        if (gyro::icntl[Param::verbose] > 1) {
+        if (gyro::icntl[Param::verbose] > 3) {
             std::cout << "level: " << i;
             if (i != levels - 1)
                 std::cout << ", coarse_nodes: " << v_level[i]->coarse_nodes;
@@ -66,7 +66,7 @@ void gmgpolar::define_coarse_nodes()
         }
         if (i != levels - 1 &&
             (v_level[i]->nr % 2 != 1 || v_level[i]->nr < 5 || v_level[i]->ntheta % 2 != 0 || v_level[i]->ntheta < 8)) {
-            std::cout << "ACHTUNG: Cannot reach the desired number of levels (" << levels << "). Replacing it with "
+            std::cout << "WARNING: Cannot reach the desired number of levels (" << levels << "). Replacing it with "
                       << i + 1 << "...\n";
             ;
             levels = i + 1;
@@ -139,7 +139,7 @@ void level::define_coarse_nodes_onelevel(level* coarser)
         itr_theta++;
     }
 
-    if (gyro::icntl[Param::verbose] > 3) {
+    if (gyro::icntl[Param::verbose] > 5) {
         display_r();
         display_theta();
         std::cout << "Coarse grid: \n";
@@ -166,7 +166,7 @@ void level::store_theta_n_co()
     if (fabs(theta[ntheta - 1] - 2 * PI) < 1e-10) // check if end is 2*Pi (for periodic boundary conditions)
         ntheta_int--;
     if (ntheta_int % 2)
-        std::cout << "ACHTUNG.... theta has an odd size\n";
+        std::cout << "WARNING: The number of thetas is odd. Please use even numbers only.\n";
 
     // To take first/last theta into account without conditions:
     // - theta_per = [theta_last, theta, 2PI]
@@ -227,7 +227,7 @@ void level::store_theta_n_co()
     for (int i = 0; i < nr_int; i++)
         hplus[i] = r[i + 1] - r[i];
 
-    if (gyro::icntl[Param::verbose] > 3) {
+    if (gyro::icntl[Param::verbose] > 5) {
         gyro::disp(theta_per, "theta_per");
         gyro::disp(cos_theta, "cos_theta");
         gyro::disp(sin_theta, "sin_theta");
diff --git a/src/direct_solve.cpp b/src/direct_solve.cpp
index 8ca92c14..160a38bb 100644
--- a/src/direct_solve.cpp
+++ b/src/direct_solve.cpp
@@ -288,7 +288,9 @@ std::vector<double> level::solve_gaussian_elimination(std::vector<int> A_row_ind
 {
     // int m_solution = f.size();
 
-    // std::cout << "Forward substitution \n";
+    if (gyro::icntl[Param::verbose] > 3) {
+        std::cout << "Starting Gaussian elimination forward substitution \n";
+    }
     // Solve Ly = b for y (Forward Substitution)
     // (assumes rows in increasing order)
     int row = A_row_indices[0];
@@ -299,7 +301,9 @@ std::vector<double> level::solve_gaussian_elimination(std::vector<int> A_row_ind
         y[A_row_indices[j]] -= y[A_col_indices[j]] * A_vals[j];
     }
 
-    // std::cout << "Backward substitution \n";
+    if (gyro::icntl[Param::verbose] > 3) {
+        std::cout << "Starting Gaussian elimination backward substitution \n";
+    }
     //Solve Ux = y for x (Backward Substitution)
     row = A_row_indices[A_row_indices.size() - 1];
     int jind;
@@ -748,8 +752,6 @@ void level::fill_in_circle(int ij, int smoother)
     int ind, ind2;
     int msc     = m_sc[smoother];
     int size_LU = 3 * msc + 2 * (msc - 3);
-    // std::cout << "ij: " << ij << ", smoother: " << smoother << ", m_sc: " << msc << ", size_LU: " << size_LU << "\n";
-    // std::cout << "A_Zebra_r_row[smoother][ij].size(): " << A_Zebra_r_row[smoother][ij].size() << "\n";
     // first line
     A_Zebra_r_LU_row[smoother][ij]    = std::vector<int>(size_LU);
     A_Zebra_c_LU_row[smoother][ij]    = std::vector<int>(size_LU);
@@ -767,9 +769,7 @@ void level::fill_in_circle(int ij, int smoother)
     for (int i = 0; i < msc - 2; i++) {
         ind  = 3 * i + 3;
         ind2 = 4 * i + 3;
-        // std::cout << "ind: " << ind << ", ind2: " << ind2 << ", ij: " << ij << "\n";
-        // std::cout << "A_Zebra_r_row[smoother][ij][ind]: " << A_Zebra_r_row[smoother][ij][ind] << "\n";
-        // std::cout << "A_Zebra_r_LU_row[smoother][ij][ind2]: " << A_Zebra_r_LU_row[smoother][ij][ind2] << "\n";
+
         A_Zebra_r_LU_row[smoother][ij][ind2]     = A_Zebra_r_row[smoother][ij][ind];
         A_Zebra_r_LU_row[smoother][ij][ind2 + 1] = A_Zebra_r_row[smoother][ij][ind + 1];
         A_Zebra_r_LU_row[smoother][ij][ind2 + 2] = A_Zebra_r_row[smoother][ij][ind + 2];
diff --git a/src/gyro.cpp b/src/gyro.cpp
index df76fef1..23bfdb57 100644
--- a/src/gyro.cpp
+++ b/src/gyro.cpp
@@ -239,14 +239,19 @@ void gyro::show_params()
  * Single
  ************************/
 /*!
- *  \brief Computes the distance of a node (x, y) to the Dirichlet boundary
+ *  \brief Computes the distance of a node (r, theta) to the Dirichlet boundary.
+ * The function assumes a boundary given by r_0=Param::r0_DB and R=Param::R.
  *
- * Computes the distance of a node (x, y) to the Dirichlet boundary, using 
- * tolerance of 1e-10 to define boundary.
- * Highly dependent on dcntl[Param::r0_DB] which is 1e6 is icntl[Param::DirBC_Interior]=0.
+ * Attention: As we compute (r-r_0)*(r-R), the output of this function is highly dependent on r_0.
+ * If r_0 is set "very close" to zero, roundoff errors can appear for function values on the innermost circles.
+ * However, as dcntl[Param::r0_DB] is advised to be chosen small if icntl[Param::DirBC_Interior]=0 is set,
+ * this function should be handled with care. Modestly small values like 1e-5 or 1e-8 should not create a problem.
  * 
- * \param x: the x coordinate of the node
- * \param y: the  coordinate of the node
+ * For more details on the Across-The-Origin heuristic, which is implemented with icntl[Param::DirBC_Interior]=0,
+ * see Kühn, Kruse, Rüde, Journal of Scientific Computing, 91 (28), (2022).
+ * 
+ * \param r_i: the r coordinate of the node
+ * \param theta_j: the theta coordinate of the node
  * \param verbose: verbose level for debug
  *
  * \return the distance
@@ -263,9 +268,10 @@ double gyro::distBoundary(double r_i, double theta_j, int verbose)
     double boundDefDim1 = fabs(r_i - dcntl[Param::r0_DB]) * fabs(r_i - dcntl[Param::R]);
     double boundDefDim2 = 1;
 
-    // if (verbose)
-    //     std::cout << "DISTBOUNDARY (" << x << ", " << y << "): " << boundDefDim1 * boundDefDim2
-    //               << " (boundDefDim1: " << boundDefDim1 << ", boundDefDim2: " << boundDefDim2 << ")\n";
+    if (verbose > 5) {
+        std::cout << "DISTBOUNDARY (" << r_i << ", " << theta_j << "): " << boundDefDim1 * boundDefDim2
+                  << " (boundDefDim1: " << boundDefDim1 << ", boundDefDim2: " << boundDefDim2 << ")\n";
+    }
     return boundDefDim1 * boundDefDim2;
 } /* ----- end of gyro::distBoundary ----- */
 
@@ -309,8 +315,9 @@ double gyro::def_solution_rt(double r_i, double theta_j, int verbose)
 
     dcntl[Param::t_sol] += TOC;
 
-    // if (verbose)
-    //     std::cout << "SOL (" << x << ", " << y << "): " << sol << "\n";
+    if (verbose > 5) {
+        std::cout << "SOL (" << r_i << ", " << theta_j << "): " << sol << "\n";
+    }
     return sol;
 } /* ----- end of gyro::eval_def_solution_vec ----- */
 
@@ -347,9 +354,10 @@ std::vector<double> gyro::def_solution_rt(double r_i, std::vector<double> theta,
 
     dcntl[Param::t_sol] += TOC;
 
-    // if (verbose)
-    //     for (int i = 0; i < ntheta; i++)
-    //         std::cout << "SOL (" << r_i << ", " << theta[i] << "): " << sol[i] << "\n";
+    if (verbose > 5) {
+        for (int i = 0; i < ntheta; i++)
+            std::cout << "SOL (" << r_i << ", " << theta[i] << "): " << sol[i] << "\n";
+    }
     return sol;
 } /* ----- end of gyro::eval_def_solution_vec ----- */
 
@@ -383,8 +391,9 @@ double gyro::coeff(double r, int verbose)
 
     dcntl[Param::t_coeff] += TOC;
 
-    // if (verbose)
-    //     std::cout << "COEFF (" << r << "): " << coeff << "\n";
+    if (verbose > 5) {
+        std::cout << "COEFF_A (" << r << "): " << coeff_a << "\n";
+    }
     return coeff_a;
 } /* ----- end of level::coeff ----- */
 
@@ -412,6 +421,9 @@ double gyro::coeff_beta(double r, int verbose)
 
     dcntl[Param::t_coeff] += TOC;
 
+    if (verbose > 5) {
+        std::cout << "COEFF_B (" << r << "): " << coeff_b << "\n";
+    }
     return coeff_b;
 } /* ----- end of level::coeff ----- */
 
@@ -440,10 +452,10 @@ std::vector<double> gyro::coeff(std::vector<double> r, int verbose)
 
     dcntl[Param::t_coeff] += TOC;
 
-    // if (verbose) {
-    //     disp(r, "r");
-    //     disp(coeff, "coeff");
-    // }
+    if (gyro::icntl[Param::verbose] > 5) {
+        disp(r, "r");
+        disp(coeff_a, "coeff_a");
+    }
     return coeff_a;
 } /* ----- end of level::coeff ----- */
 
@@ -471,7 +483,10 @@ std::vector<double> gyro::coeff_beta(std::vector<double> r, int verbose)
     functions->coeffs2(r, Rmax, coeff_b);
 
     dcntl[Param::t_coeff] += TOC;
-
+    if (gyro::icntl[Param::verbose] > 5) {
+        disp(r, "r");
+        disp(coeff_b, "coeff_b");
+    }
     return coeff_b;
 } /* ----- end of level::coeff_beta ----- */
 
@@ -509,8 +524,9 @@ double gyro::detDFinv(double r, double theta, int verbose)
 
     dcntl[Param::t_detDFinv] += TOC;
 
-    // if (verbose)
-    //     std::cout << "dETDFINV (" << r << ", " << theta << "): " << detDFinv_r << "\n";
+    if (gyro::icntl[Param::verbose] > 5) {
+        std::cout << "Value of detDFinv (" << r << ", " << theta << "): " << detDFinv_r << "\n";
+    }
     return detDFinv_r;
 } /* ----- end of level::detDFinv ----- */
 double gyro::arr(double r, double theta, int verbose)
@@ -530,8 +546,9 @@ double gyro::arr(double r, double theta, int verbose)
 
     dcntl[Param::t_arr_art_att] += TOC;
 
-    // if (verbose)
-    //     std::cout << "ARR (" << r << ", " << theta << "): " << arr_r << "\n";
+    if (gyro::icntl[Param::verbose] > 5) {
+        std::cout << "Value of arr (" << r << ", " << theta << "): " << arr_r << "\n";
+    }
     return arr_r;
 } /* ----- end of level::arr ----- */
 double gyro::art(double r, double theta, int verbose)
@@ -553,8 +570,9 @@ double gyro::art(double r, double theta, int verbose)
 
     dcntl[Param::t_arr_art_att] += TOC;
 
-    // if (verbose)
-    //     std::cout << "ART (" << r << ", " << theta << "): " << art_r << "\n";
+    if (gyro::icntl[Param::verbose] > 5) {
+        std::cout << "Value of art (" << r << ", " << theta << "): " << art_r << "\n";
+    }
     return art_r;
 } /* ----- end of level::art ----- */
 double gyro::att(double r, double theta, int verbose)
@@ -574,8 +592,9 @@ double gyro::att(double r, double theta, int verbose)
 
     dcntl[Param::t_arr_art_att] += TOC;
 
-    // if (verbose)
-    //     std::cout << "ATT (" << r << ", " << theta << "): " << att_r << "\n";
+    if (gyro::icntl[Param::verbose] > 5) {
+        std::cout << "Value of att (" << r << ", " << theta << "): " << att_r << "\n";
+    }
     return att_r;
 } /* ----- end of level::att ----- */
 void gyro::arr_att_art(double r, double theta, double& arr, double& att, double& art, int verbose)
@@ -600,11 +619,11 @@ void gyro::arr_att_art(double r, double theta, double& arr, double& att, double&
 
     dcntl[Param::t_arr_art_att] += TOC;
 
-    // if (verbose) {
-    //     std::cout << "ARR (" << r << ", " << theta << "): " << arr << "\n";
-    //     std::cout << "ATT (" << r << ", " << theta << "): " << att << "\n";
-    //     std::cout << "ART (" << r << ", " << theta << "): " << art << "\n";
-    // }
+    if (gyro::icntl[Param::verbose] > 5) {
+        std::cout << "Value of arr (" << r << ", " << theta << "): " << arr << "\n";
+        std::cout << "Value of att (" << r << ", " << theta << "): " << att << "\n";
+        std::cout << "Value of art (" << r << ", " << theta << "): " << art << "\n";
+    }
 } /* ----- end of level::arr_att_art ----- */
 
 /*!
@@ -649,11 +668,11 @@ std::vector<double> gyro::detDFinv(double r, std::vector<double> theta, std::vec
 
     dcntl[Param::t_detDFinv] += TOC;
 
-    // if (verbose) {
-    //     disp(theta, "theta");
-    //     disp(cos_theta, "cos_theta");
-    //     disp(detDFinv_r, "detDFinv_r");
-    // }
+    if (gyro::icntl[Param::verbose] > 5) {
+        disp(theta, "theta");
+        disp(cos_theta, "cos_theta");
+        disp(detDFinv_r, "detDFinv_r");
+    }
     return detDFinv_r;
 } /* ----- end of level::detDFinv ----- */
 
@@ -699,10 +718,10 @@ std::vector<double> gyro::detDFinv(std::vector<double> r, double theta, int verb
 
     dcntl[Param::t_detDFinv] += TOC;
 
-    // if (verbose) {
-    //     disp(r, "r");
-    //     disp(detDFinv_r, "detDFinv_r");
-    // }
+    if (gyro::icntl[Param::verbose] > 5) {
+        disp(r, "r");
+        disp(detDFinv_r, "detDFinv_r");
+    }
     return detDFinv_r;
 } /* ----- end of level::detDFinv ----- */
 
@@ -730,12 +749,12 @@ std::vector<double> gyro::arr(double r, std::vector<double> theta, std::vector<d
 
     dcntl[Param::t_arr_art_att] += TOC;
 
-    // if (verbose) {
-    //     disp(theta, "theta");
-    //     disp(sin_theta, "sin_theta");
-    //     disp(cos_theta, "cos_theta");
-    //     disp(arr_r, "arr_r");
-    // }
+    if (gyro::icntl[Param::verbose] > 5) {
+        disp(theta, "theta");
+        disp(sin_theta, "sin_theta");
+        disp(cos_theta, "cos_theta");
+        disp(arr_r, "arr_r");
+    }
     return arr_r;
 } /* ----- end of level::arr ----- */
 std::vector<double> gyro::art(double r, std::vector<double> theta, std::vector<double> sin_theta,
@@ -766,14 +785,12 @@ std::vector<double> gyro::art(double r, std::vector<double> theta, std::vector<d
 
     dcntl[Param::t_arr_art_att] += TOC;
 
-    // if (verbose) {
-    //     std::cout << "ART (" << r << "\n";
-    //     disp(theta, "theta");
-    //     disp(sin_theta, "sin_theta");
-    //     disp(cos_theta, "cos_theta");
-    //     disp(art_r, "art_r");
-    //     std::cout << ")\n";
-    // }
+    if (gyro::icntl[Param::verbose] > 5) {
+        disp(theta, "theta");
+        disp(sin_theta, "sin_theta");
+        disp(cos_theta, "cos_theta");
+        disp(art_r, "art_r");
+    }
     return art_r;
 } /* ----- end of level::art ----- */
 std::vector<double> gyro::att(double r, std::vector<double> theta, std::vector<double> sin_theta,
@@ -802,12 +819,12 @@ std::vector<double> gyro::att(double r, std::vector<double> theta, std::vector<d
 
     dcntl[Param::t_arr_art_att] += TOC;
 
-    // if (verbose) {
-    //     disp(theta, "theta");
-    //     disp(sin_theta, "sin_theta");
-    //     disp(cos_theta, "cos_theta");
-    //     disp(att_r, "att_r");
-    // }
+    if (gyro::icntl[Param::verbose] > 5) {
+        disp(theta, "theta");
+        disp(sin_theta, "sin_theta");
+        disp(cos_theta, "cos_theta");
+        disp(att_r, "att_r");
+    }
     return att_r;
 } /* ----- end of level::att ----- */
 void gyro::arr_att_art(double r, std::vector<double> theta, std::vector<double> sin_theta,
@@ -843,14 +860,14 @@ void gyro::arr_att_art(double r, std::vector<double> theta, std::vector<double>
 
     dcntl[Param::t_arr_art_att] += TOC;
 
-    // if (verbose) {
-    //     disp(theta, "theta");
-    //     disp(sin_theta, "sin_theta");
-    //     disp(cos_theta, "cos_theta");
-    //     disp(arr, "arr");
-    //     disp(att, "att");
-    //     disp(art, "art");
-    // }
+    if (gyro::icntl[Param::verbose] > 5) {
+        disp(theta, "theta");
+        disp(sin_theta, "sin_theta");
+        disp(cos_theta, "cos_theta");
+        disp(arr, "arr");
+        disp(att, "att");
+        disp(art, "art");
+    }
 } /* ----- end of level::arr_att_art ----- */
 void gyro::arr_att_art(std::vector<double> r, double theta, std::vector<double>& arr_r, std::vector<double>& att_r,
                        std::vector<double>& art_r, int verbose)
@@ -885,14 +902,14 @@ void gyro::arr_att_art(std::vector<double> r, double theta, std::vector<double>&
 
     dcntl[Param::t_arr_art_att] += TOC;
 
-    // if (verbose) {
-    //     disp(r, "r");
-    //     disp(detDFinv_r, "detDFinv_r");
-    //     disp(coeff_r, "coeff_r");
-    //     disp(arr_r, "arr_r");
-    //     disp(att_r, "att_r");
-    //     disp(art_r, "art_r");
-    // }
+    if (gyro::icntl[Param::verbose] > 5) {
+        disp(r, "r");
+        disp(detDFinv_r, "detDFinv_r");
+        disp(coeff_r, "coeff_r");
+        disp(arr_r, "arr_r");
+        disp(att_r, "att_r");
+        disp(art_r, "art_r");
+    }
 } /* ----- end of level::arr_att_art ----- */
 
 /***************************************************************************
@@ -934,6 +951,10 @@ void gyro::trafo(double& r_i, double& theta_j, double& x, double& y, int verbose
 /*!
  *  \brief Transform from polar to cartesian
  *
+ *  TODO: This function could be reimplemented and used for simple geometries
+ *        For more general geometries, the inverse mapping is not available 
+ *        and a check on the geometry would be needed.
+ *         
  *  Transform one couple r_i,theta_j from cartesian to polar
  * 
  * \param r_i: the r coordinate of the node (out)
@@ -1019,57 +1040,6 @@ void gyro::trafo(double r_i, std::vector<double> theta, std::vector<double> sin_
     //     }
 } /* ----- end of gyro::trafo ----- */
 
-/***************************************************************************
-     * Display arrays or vectors
-     **************************************************************************/
-/*!
- *  \brief Display arrays of int or double
- *
- *  Display arrays of int or double
- * 
- * \param na: size of the array
- * \param a: the array
- * \param s_a: name of the array (for display)
- *
- */
-void gyro::disp(int na, double* a, const std::string& s_a)
-{
-    std::cout << s_a << "(" << na << "): ";
-    for (int i = 0; i < na; i++)
-        std::cout << a[i] << " ";
-    std::cout << "\n";
-} /* ----- end of gyro::disp ----- */
-void gyro::disp(int na, int* a, const std::string& s_a)
-{
-    std::cout << s_a << "(" << na << "): ";
-    for (int i = 0; i < na; i++)
-        std::cout << a[i] << " ";
-    std::cout << "\n";
-} /* ----- end of gyro::disp ----- */
-/*!
- *  \brief Display vectors of int or double
- *
- *  Display vectors of int or double
- * 
- * \param a: the vector
- * \param s_a: name of the vector (for display)
- *
- */
-void gyro::disp(std::vector<double> a, const std::string& s_a)
-{
-    std::cout << s_a << "(" << a.size() << "): ";
-    for (std::size_t i = 0; i < a.size(); i++)
-        std::cout << a[i] << " ";
-    std::cout << "\n";
-} /* ----- end of gyro::disp ----- */
-void gyro::disp(std::vector<int> a, const std::string& s_a)
-{
-    std::cout << s_a << "(" << a.size() << "): ";
-    for (std::size_t i = 0; i < a.size(); i++)
-        std::cout << a[i] << " ";
-    std::cout << "\n";
-} /* ----- end of gyro::disp ----- */
-
 /*!
  *  \brief Sparse matrix-vector product
  *
diff --git a/src/level.cpp b/src/level.cpp
index abc27a04..aa86bef4 100644
--- a/src/level.cpp
+++ b/src/level.cpp
@@ -146,10 +146,11 @@ void level::build_bound()
         }
     }
 
-    if (gyro::icntl[Param::verbose] > 3)
+    if (gyro::icntl[Param::verbose] > 5)
+        std::cout << "Printing for (r, theta) if the node is on the boundary.\n";
         for (int j = 0; j < nr; j++)
             for (int i = 0; i < ntheta_int; i++)
-                std::cout << "DISTBOUNDARY (" << r[j] << ", " << theta[j] << "): " << is_bound[j * ntheta_int + i]
+                std::cout << "(" << r[j] << ", " << theta[j] << "): " << is_bound[j * ntheta_int + i]
                           << "\n";
 } /* ----- end of gyro::build_bound ----- */
 
diff --git a/src/main.cpp b/src/main.cpp
index b8315221..20a3d0fc 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -135,7 +135,7 @@ int main(int argc, char* argv[])
     //Set number of threads for the openmp parallelization
     omp_set_num_threads(gyro::icntl[Param::openmp]);
 
-    if (gyro::icntl[Param::verbose] > 0) {
+    if (gyro::icntl[Param::verbose] > 1) {
 #pragma omp parallel
         {
 #pragma omp master
@@ -144,7 +144,7 @@ int main(int argc, char* argv[])
             }
         }
     }
-    if (gyro::icntl[Param::verbose] > 0) {
+    if (gyro::icntl[Param::verbose] > 1) {
 
 #pragma omp parallel
         {
@@ -198,7 +198,7 @@ int main(int argc, char* argv[])
         if (gyro::dcntl[Param::kappa_eps] == 42 && gyro::dcntl[Param::delta_e] == 42) {
             gyro::get_geometry_coeffs(geom);
         }
-        if (gyro::icntl[Param::verbose] > 1)
+        if (gyro::icntl[Param::verbose] > 2)
             gyro::show_params();
         gyro::select_functions_class(gyro::icntl[Param::alpha_coeff], gyro::icntl[Param::beta_coeff],
                                      gyro::icntl[Param::mod_pk], gyro::icntl[Param::prob]);
@@ -277,7 +277,7 @@ int main(int argc, char* argv[])
                      "#######\n";
         geom = (geometry_type)gyro::icntl[Param::mod_pk];
         gyro::get_geometry_coeffs(geom);
-        if (gyro::icntl[Param::verbose] > 1)
+        if (gyro::icntl[Param::verbose] > 2)
             gyro::show_params();
         gyro::select_functions_class(gyro::icntl[Param::alpha_coeff], gyro::icntl[Param::beta_coeff],
                                      gyro::icntl[Param::mod_pk], gyro::icntl[Param::prob]);
@@ -341,7 +341,7 @@ int main(int argc, char* argv[])
 
                 geom = (geometry_type)gyro::icntl[Param::mod_pk];
                 gyro::get_geometry_coeffs(geom);
-                if (gyro::icntl[Param::verbose] > 1)
+                if (gyro::icntl[Param::verbose] > 2)
                     gyro::show_params();
                 gyro::select_functions_class(gyro::icntl[Param::alpha_coeff], gyro::icntl[Param::beta_coeff],
                                              gyro::icntl[Param::mod_pk], gyro::icntl[Param::prob]);
@@ -419,7 +419,7 @@ int main(int argc, char* argv[])
 
                             geom = (geometry_type)gyro::icntl[Param::mod_pk];
                             gyro::get_geometry_coeffs(geom);
-                            if (gyro::icntl[Param::verbose] > 1)
+                            if (gyro::icntl[Param::verbose] > 2)
                                 gyro::show_params();
                             gyro::select_functions_class(gyro::icntl[Param::alpha_coeff],
                                                          gyro::icntl[Param::beta_coeff], gyro::icntl[Param::mod_pk],
diff --git a/src/multigrid_iter.cpp b/src/multigrid_iter.cpp
index 2ff29ddf..fd889afd 100644
--- a/src/multigrid_iter.cpp
+++ b/src/multigrid_iter.cpp
@@ -48,9 +48,9 @@ void gmgpolar::multigrid_iter()
         v_level[1]->fVec_initial = v_level[1]->fVec; //Extrapolation: store the initial fVec on level 1
     if (gyro::icntl[Param::verbose] > 0) {
         if (extrapol == 1)
-            std::cout << "WITH IMPLICIT EXTRAPOLATION!!!\n";
+            std::cout << "with implicit extrapolation.\n";
         else if (extrapol == 2)
-            std::cout << "WITH ALTERNATIVE EXTRAPOLATION!!!\n";
+            std::cout << "with alternative extrapolation.\n WARNING: Alternative extrapolation option is a pure test or research setting.\n";
     }
 
     int it = 0;
@@ -111,7 +111,7 @@ void gmgpolar::multigrid_iter()
         //call the multigrid_cycle on level 0 (the finest level)
         multigrid_cycle_extrapol(0);
 
-        if (gyro::icntl[Param::verbose] > 3)
+        if (gyro::icntl[Param::verbose] > 5)
             gyro::disp(v_level[0]->u, "u");
 
         TIC;
@@ -121,7 +121,7 @@ void gmgpolar::multigrid_iter()
         if (extrapol < 2) {
             gmgpolar::compute_residual(0, extrapol); //compute residual on level 0
 
-            if (gyro::icntl[Param::verbose] > 3)
+            if (gyro::icntl[Param::verbose] > 5)
                 gyro::disp(v_level[0]->res, "res");
 
             nrm_2_res.push_back(0);
@@ -156,7 +156,7 @@ void gmgpolar::multigrid_iter()
                 convergence_criterium = nrm_inf_res[it];
             }
 
-            if (gyro::icntl[Param::verbose] > 0)
+            if (gyro::icntl[Param::verbose] > 1)
                 std::cout << "--> Iteration " << it << ": residual norm = " << nrm_2_res[it]
                           << ", relative residual = " << convergence_criterium << std::endl;
         }
@@ -169,7 +169,7 @@ void gmgpolar::multigrid_iter()
 
             std::vector<double> error = compute_error();
 
-            if (gyro::icntl[Param::verbose] > 3)
+            if (gyro::icntl[Param::verbose] > 5)
                 gyro::disp(error, "error");
 
             nrm_2_err.push_back(0);
@@ -181,21 +181,21 @@ void gmgpolar::multigrid_iter()
             double error_difference = fabs(nrm_2_err[it] - nrm_2_err[it - 1]);
             convergence_criterium   = error_difference / nrm_2_err[0];
 
-            if (gyro::icntl[Param::verbose] > 0)
+            if (gyro::icntl[Param::verbose] > 1)
                 std::cout << "--> Iteration " << it << ": error norm = " << nrm_2_err[it]
                           << ", relative error = " << convergence_criterium << std::endl;
         }
         t_fine_residual += TOC;
         TIC;
     }
-    // if (gyro::icntl[Param::verbose] > 0) {
-    if (it == gyro::icntl[Param::maxiter]) {
-        std::cout << "Multigrid reached maxiter=" << gyro::icntl[Param::maxiter] << "\n";
-    }
-    else {
-        std::cout << "Convergence after iteration " << it << std::endl;
+    if (gyro::icntl[Param::verbose] > 0) {
+        if (it == gyro::icntl[Param::maxiter]) {
+            std::cout << "Multigrid reached maxiter=" << gyro::icntl[Param::maxiter] << "\n";
+        }
+        else {
+            std::cout << "Convergence after iteration " << it << std::endl;
+        }
     }
-    // }
     //----------------------------------------------------------------------------------------------------------
     //!compute mean residual reduction factor rho
     if (extrapol < 2) {
@@ -221,17 +221,19 @@ void gmgpolar::multigrid_iter()
         }
         nrm_2 = sqrt(nrm_2) / sqrt(v_level[0]->m); //scaling by 1/sqrt(m)
 
-        std::cout << "2-norm of error = " << nrm_2 << std::endl;
-        std::cout << "inf-norm of error = " << nrm_inf_err << std::endl;
+        if (gyro::icntl[Param::verbose] > 0) {
+            std::cout << "2-norm of error = " << nrm_2 << std::endl;
+            std::cout << "inf-norm of error = " << nrm_inf_err << std::endl;
+        }
 
-        if (gyro::icntl[Param::verbose] > 3)
+        if (gyro::icntl[Param::verbose] > 5)
             gyro::disp(error, "error");
     }
 
     t_error += TOC;
     TIC;
 
-    if (gyro::icntl[Param::verbose] > 3)
+    if (gyro::icntl[Param::verbose] > 5)
         gyro::disp(v_level[0]->u, "u");
 } /* ----- end of gmgpolar::multigrid_iter ----- */
 
@@ -246,7 +248,7 @@ void gmgpolar::multigrid_iter()
 void gmgpolar::multigrid_cycle_extrapol(int l)
 {
     //l = current level we are on (from class level)
-    if (gyro::icntl[Param::verbose] > 2) {
+    if (gyro::icntl[Param::verbose] > 4) {
         std::cout << "********************************************" << std::endl;
         std::cout << "MULTIGRID ON LEVEL " << l << std::endl;
         std::cout << "nr = " << v_level[l]->nr << ", ntheta = " << v_level[l]->ntheta << std::endl;
@@ -292,7 +294,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
                 //call the smooothing function, the result is u_sc which is directly inserted into u
                 v_level[l]->multigrid_smoothing0(smoother);
 
-                if (gyro::icntl[Param::verbose] > 2)
+                if (gyro::icntl[Param::verbose] > 4)
                     std::cout << "SMOOTHER: " << smoother << " = " << TOC << "\n";
             }
         }
@@ -317,7 +319,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
                                 v_level[l]->dep_Asc[smoother], v_level[l]->dep_Asc[sm], v_level[l]->dep_Asc[1],
                                 v_level[l]->dep_Asc_ortho[smoother]);
 
-                            if (gyro::icntl[Param::verbose] > 2)
+                            if (gyro::icntl[Param::verbose] > 4)
                                 std::cout << "SMOOTHER: " << smoother << " = " << TOC << "\n";
                         }
                     }
@@ -339,7 +341,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
     }
     //std::cout << "pre-smoothing done \n";
 
-    if (gyro::icntl[Param::verbose] > 3)
+    if (gyro::icntl[Param::verbose] > 5)
         gyro::disp(v_level[l]->u, "u");
 
     t = t_smoothing_tmp;
@@ -350,7 +352,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
     //even if we have extrapolation, compute just normal residual (extrapolation-restriction follows in the next step)
     gmgpolar::compute_residual(l, 0);
 
-    if (gyro::icntl[Param::verbose] > 3)
+    if (gyro::icntl[Param::verbose] > 5)
         gyro::disp(v_level[l]->res, "res");
 
     t_residual += TOC;
@@ -379,7 +381,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
             }
         }
 
-        if (gyro::icntl[Param::verbose] > 3)
+        if (gyro::icntl[Param::verbose] > 5)
             gyro::disp(res_1, "res_1");
 
         //Au_coarse = A(l+1) * u_coarse = A(l+1) * P_inj^T * u(l)
@@ -404,7 +406,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
             }
         }
 
-        if (gyro::icntl[Param::verbose] > 3)
+        if (gyro::icntl[Param::verbose] > 5)
             gyro::disp(u_coarse, "u_coarse");
 
         if (gyro::icntl[Param::matrix_free] == 1) {
@@ -435,7 +437,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
             }
         }
 
-        if (gyro::icntl[Param::verbose] > 3)
+        if (gyro::icntl[Param::verbose] > 5)
             gyro::disp(Au_coarse, "Au_coarse");
 
         // f(l+1) = 4/3 * res_1 - 1/3 * (f(l+1) - Au_coarse)
@@ -443,7 +445,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
             v_level[l + 1]->fVec[i] = 4. / 3. * res_1[i] - 1. / 3. * (v_level[l + 1]->fVec_initial[i] - Au_coarse[i]);
         }
 
-        if (gyro::icntl[Param::verbose] > 3)
+        if (gyro::icntl[Param::verbose] > 5)
             gyro::disp(v_level[l + 1]->fVec, "fVec");
     }
     else { //no extrapolation
@@ -494,13 +496,13 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
         TIC;
     }
 
-    if (gyro::icntl[Param::verbose] > 2) {
+    if (gyro::icntl[Param::verbose] > 4) {
         std::cout << "********************************************" << std::endl;
         std::cout << "BACK ON LEVEL " << l << std::endl;
         std::cout << "********************************************" << std::endl;
     }
 
-    if (gyro::icntl[Param::verbose] > 3)
+    if (gyro::icntl[Param::verbose] > 5)
         gyro::disp(error_coarse, "error_coarse");
 
     TIC;
@@ -544,7 +546,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
     }
     //std::cout << "error prolongated \n";
 
-    if (gyro::icntl[Param::verbose] > 3)
+    if (gyro::icntl[Param::verbose] > 5)
         gyro::disp(error_fine, "error_fine");
 
     //! correction of solution (on level l)
@@ -585,7 +587,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
                 //call the smooothing function, the result is u_sc which is directly inserted into u
                 v_level[l]->multigrid_smoothing0(smoother);
 
-                if (gyro::icntl[Param::verbose] > 2)
+                if (gyro::icntl[Param::verbose] > 4)
                     std::cout << "SMOOTHER: " << smoother << " = " << TOC << "\n";
             }
         }
@@ -610,7 +612,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
                                 v_level[l]->dep_Asc[smoother], v_level[l]->dep_Asc[sm], v_level[l]->dep_Asc[1],
                                 v_level[l]->dep_Asc_ortho[smoother]);
 
-                            if (gyro::icntl[Param::verbose] > 2)
+                            if (gyro::icntl[Param::verbose] > 4)
                                 std::cout << "SMOOTHER: " << smoother << " = " << TOC << "\n";
                         }
                     }
@@ -626,7 +628,7 @@ void gmgpolar::multigrid_cycle_extrapol(int l)
     t = t_total_tmp;
     t_total += TOC;
 
-    if (gyro::icntl[Param::verbose] > 3)
+    if (gyro::icntl[Param::verbose] > 5)
         gyro::disp(v_level[l]->u, "u");
 } /* ----- end of gmgpolar::multigrid_cycle_extrapol ----- */
 
diff --git a/src/polar_multigrid.cpp b/src/polar_multigrid.cpp
index ca7963aa..d2070004 100644
--- a/src/polar_multigrid.cpp
+++ b/src/polar_multigrid.cpp
@@ -45,8 +45,8 @@ void getrands(std::vector<double>& x, Generator& gen, unsigned num)
  */
 void gmgpolar::polar_multigrid()
 {
-    if (gyro::icntl[Param::verbose] > 1)
-        std::cout << "Define the coarse grids...\n";
+    if (gyro::icntl[Param::verbose] > 3)
+        std::cout << "Define the coarse grid nodes...\n";
     check_geom();
     define_coarse_nodes();
 
@@ -131,7 +131,7 @@ void gmgpolar::polar_multigrid()
         std::cout << "t_applyA: " << t_applyA << std::endl;
     }
     else {
-        if (gyro::icntl[Param::verbose] > 1)
+        if (gyro::icntl[Param::verbose] > 3)
             std::cout
                 << "Building discretized system, restriction and interpolation operators, and defining splittings...\n";
         prepare_op_levels();
@@ -143,15 +143,17 @@ void gmgpolar::polar_multigrid()
             cycle_str = "V";
         else if (gyro::icntl[Param::cycle] == 2)
             cycle_str = "W";
-        // if (gyro::icntl[Param::verbose] > 1)
-        std::cout << "\nProb: " << gyro::icntl[Param::prob] << ", alpha_coeff: " << gyro::icntl[Param::alpha_coeff]
-                  << ", beta_coeff: " << gyro::icntl[Param::beta_coeff] << " ***** Problem size " << m << " ("
-                  << v_level[0]->nr << ", " << v_level[0]->ntheta << "), " << levels
-                  << " grids, nr_exp=" << gyro::icntl[Param::nr_exp] << ", aniso=" << gyro::icntl[Param::fac_ani]
-                  << " ***** smoother=" << gyro::icntl[Param::smoother] << ", r0=" << gyro::dcntl[Param::R0]
-                  << ", extrapolation=" << gyro::icntl[Param::extrapolation]
-                  << ", mod_pk=" << gyro::icntl[Param::mod_pk] << ", DirBC=" << gyro::icntl[Param::DirBC_Interior]
-                  << ", divide=" << gyro::icntl[Param::divideBy2] << " *****\n";
+        if (gyro::icntl[Param::verbose] > 2)
+        {
+            std::cout << "\nProb: " << gyro::icntl[Param::prob] << ", alpha_coeff: " << gyro::icntl[Param::alpha_coeff]
+                      << ", beta_coeff: " << gyro::icntl[Param::beta_coeff] << " ***** Problem size " << m << " ("
+                      << v_level[0]->nr << ", " << v_level[0]->ntheta << "), " << levels
+                      << " grids, nr_exp=" << gyro::icntl[Param::nr_exp] << ", aniso=" << gyro::icntl[Param::fac_ani]
+                      << " ***** smoother=" << gyro::icntl[Param::smoother] << ", r0=" << gyro::dcntl[Param::R0]
+                      << ", extrapolation=" << gyro::icntl[Param::extrapolation]
+                      << ", mod_pk=" << gyro::icntl[Param::mod_pk] << ", DirBC=" << gyro::icntl[Param::DirBC_Interior]
+                      << ", divide=" << gyro::icntl[Param::divideBy2] << " *****\n";
+        }
 
         double scaling = 1.0;
         if (gyro::icntl[Param::compute_rho])
@@ -163,39 +165,42 @@ void gmgpolar::polar_multigrid()
             debug();
         }
         else if (gyro::icntl[Param::write_radii_angles] == 0) {
-            if (gyro::icntl[Param::verbose] > 1)
-                std::cout << "Multigrid iteration....!\n\n";
+            if (gyro::icntl[Param::verbose] > 3)
+                std::cout << "Executing multigrid iteration....\n\n";
             multigrid_iter();
 
-            //        if (gyro::icntl[Param::verbose] > 1)
-            for (int l = 0; l < levels; l++) {
-                std::cout << "LEVEL " << l << "\n";
-                std::cout << "\nt_smoothing: " << v_level[l]->t_smoothing << ", t_f_sc: " << v_level[l]->t_f_sc
-                          << ", t_Asc_ortho: " << v_level[l]->t_Asc_ortho << ", t_Asc: " << v_level[l]->t_Asc << "\n";
-                std::cout << "\nt_get_ptr: " << v_level[l]->t_get_ptr
-                          << ", t_get_stencil: " << v_level[l]->t_get_stencil
-                          << ", t_get_smoother: " << v_level[l]->t_get_smoother
-                          << ", t_get_row: " << v_level[l]->t_get_row << "\n";
-                std::cout << "\n";
+            if (gyro::icntl[Param::verbose] > 1) {
+                for (int l = 0; l < levels; l++) {
+                    std::cout << "LEVEL " << l << "\n";
+                    std::cout << "\nt_smoothing: " << v_level[l]->t_smoothing << ", t_f_sc: " << v_level[l]->t_f_sc
+                              << ", t_Asc_ortho: " << v_level[l]->t_Asc_ortho << ", t_Asc: " << v_level[l]->t_Asc
+                              << "\n";
+                    std::cout << "\nt_get_ptr: " << v_level[l]->t_get_ptr
+                              << ", t_get_stencil: " << v_level[l]->t_get_stencil
+                              << ", t_get_smoother: " << v_level[l]->t_get_smoother
+                              << ", t_get_row: " << v_level[l]->t_get_row << "\n";
+                    std::cout << "\n";
+                }
+            }
+
+            if (gyro::icntl[Param::verbose] > 1) {
+                std::cout << "\nt_setup: " << t_setup << ", t_build: " << t_build << ", t_facto_Ac: " << t_facto_Ac
+                          << ", t_build_P: " << t_build_P << ", t_build_Asc: " << t_build_Asc
+                          << ", t_facto_Asc: " << t_facto_Asc << "\n";
+                std::cout << "t_total (fine): " << t_total << ", t_smoothing: " << t_smoothing
+                          << ", t_residual: " << t_residual << ", t_restriction: " << t_restriction
+                          << ", t_Ac: " << t_Ac << ", t_prolongation: " << t_prolongation
+                          << ", t_fine_residual: " << t_fine_residual << ", t_error: " << t_error << "\n";
+                std::cout << "t_applyA: " << t_applyA << std::endl;
             }
 
-            //        if (gyro::icntl[Param::verbose] > 0) {
-            std::cout << "\nt_setup: " << t_setup << ", t_build: " << t_build << ", t_facto_Ac: " << t_facto_Ac
-                      << ", t_build_P: " << t_build_P << ", t_build_Asc: " << t_build_Asc
-                      << ", t_facto_Asc: " << t_facto_Asc << "\n";
-            std::cout << "t_total (fine): " << t_total << ", t_smoothing: " << t_smoothing
-                      << ", t_residual: " << t_residual << ", t_restriction: " << t_restriction << ", t_Ac: " << t_Ac
-                      << ", t_prolongation: " << t_prolongation << ", t_fine_residual: " << t_fine_residual
-                      << ", t_error: " << t_error << "\n";
-            std::cout << "t_applyA: " << t_applyA << std::endl;
-            //        }
-
-            //        if (gyro::icntl[Param::verbose] > 0) {
-            std::cout << "\nt_coeff: " << gyro::dcntl[Param::t_coeff]
-                      << ", t_arr_art_att: " << gyro::dcntl[Param::t_arr_art_att]
-                      << ", t_sol: " << gyro::dcntl[Param::t_sol] << ", t_detDFinv: " << gyro::dcntl[Param::t_detDFinv]
-                      << ", t_trafo: " << gyro::dcntl[Param::t_trafo] << "\n";
-            //        }
+            if (gyro::icntl[Param::verbose] > 1) {
+                std::cout << "\nt_coeff: " << gyro::dcntl[Param::t_coeff]
+                          << ", t_arr_art_att: " << gyro::dcntl[Param::t_arr_art_att]
+                          << ", t_sol: " << gyro::dcntl[Param::t_sol]
+                          << ", t_detDFinv: " << gyro::dcntl[Param::t_detDFinv]
+                          << ", t_trafo: " << gyro::dcntl[Param::t_trafo] << "\n";
+            }
         }
     }
 } /* ----- end of gmgpolar::polar_multigrid ----- */
@@ -260,7 +265,7 @@ void gmgpolar::prepare_op_levels()
 
         if (l < levels - 1)
             v_level[l]->mc = v_level[l + 1]->m;
-        if (gyro::icntl[Param::verbose] > 2)
+        if (gyro::icntl[Param::verbose] > 3)
             std::cout << "Create boundary array\n";
         v_level[l]->build_bound();
 
@@ -268,7 +273,7 @@ void gmgpolar::prepare_op_levels()
             v_level[0]->read_sol();
         }
 
-        if (gyro::icntl[Param::verbose] > 1)
+        if (gyro::icntl[Param::verbose] > 3)
             std::cout << "Create operator on level " << l << "\n";
         if (l == levels - 1 || gyro::icntl[Param::matrix_free] == 0) {
             TIC;
@@ -292,7 +297,7 @@ void gmgpolar::prepare_op_levels()
             TIC;
 
             if (l == levels - 1) {
-                if (gyro::icntl[Param::verbose] > 1)
+                if (gyro::icntl[Param::verbose] > 3)
                     std::cout << "Factorizing coarse operator...\n";
                 TIC;
 #ifdef USE_MUMPS
@@ -336,7 +341,7 @@ void gmgpolar::prepare_op_levels()
         // Prolongation defined on all except the coarsest level
         if (l < levels - 1) {
             v_level[l]->define_line_splitting();
-            if (gyro::icntl[Param::verbose] > 1)
+            if (gyro::icntl[Param::verbose] > 3)
                 std::cout << "delete_circles: " << v_level[l]->delete_circles << "\n";
 
             // Number of blocks per smoother
@@ -372,7 +377,7 @@ void gmgpolar::prepare_op_levels()
             // Smoother matrices Asc
             TIC;
             //build matrices A_sc
-            if (gyro::icntl[Param::verbose] > 1)
+            if (gyro::icntl[Param::verbose] > 3)
                 std::cout << "build Asc on level " << l << ", m: " << v_level[l]->m << ", mc: " << v_level[l]->mc
                           << "\n";
             v_level[l]->define_m_nz_Asc();
@@ -448,7 +453,7 @@ void gmgpolar::prepare_op_levels()
                         v_level[l]->A_Zebra_r[smoother].shrink_to_fit();
                         v_level[l]->A_Zebra_c[smoother].shrink_to_fit();
                         v_level[l]->A_Zebra_v[smoother].shrink_to_fit();
-                        if (gyro::icntl[Param::verbose] > 1)
+                        if (gyro::icntl[Param::verbose] > 3)
                             std::cout << "Factorizing smoother " << smoother << "...\n";
                         v_level[l]->facto_gaussian_elimination(
                             v_level[l]->A_Zebra_r_LU[smoother], v_level[l]->A_Zebra_c_LU[smoother],
@@ -488,7 +493,7 @@ void gmgpolar::prepare_op_levels()
                             v_level[l]->m_sc[smoother] == 1 ||
                             (gyro::icntl[Param::extrapolation] && l == 0 && smoother % 2 == 0 &&
                              !(smoother == 0 && !gyro::icntl[Param::DirBC_Interior] && ij == 0))) {
-                            if (gyro::icntl[Param::verbose] > 2)
+                            if (gyro::icntl[Param::verbose] > 3)
                                 std::cout << "No diagonal factorization\n";
                             v_level[l]->A_Zebra_v_LU_row[smoother][ij] =
                                 std::vector<double>(v_level[l]->A_Zebra_v_row[smoother][ij]);
@@ -498,7 +503,7 @@ void gmgpolar::prepare_op_levels()
                         // - and not across
                         else if (smoother < 2 && !(smoother == 0 && !gyro::icntl[Param::DirBC_Interior] && ij == 0)) {
                             TIC;
-                            if (gyro::icntl[Param::verbose] > 2)
+                            if (gyro::icntl[Param::verbose] > 3)
                                 std::cout << "Circle factorization\n";
                             // Initialize the LU factors
                             v_level[l]->fill_in_circle(ij, smoother);
@@ -519,7 +524,7 @@ void gmgpolar::prepare_op_levels()
                                 std::vector<int>(v_level[l]->A_Zebra_c_row[smoother][ij]);
                             v_level[l]->A_Zebra_v_LU_row[smoother][ij] =
                                 std::vector<double>(v_level[l]->A_Zebra_v_row[smoother][ij]);
-                            if (gyro::icntl[Param::verbose] > 2)
+                            if (gyro::icntl[Param::verbose] > 3)
                                 std::cout << "Radial factorization\n";
                             // Factorization
                             v_level[l]->facto_radial(
@@ -537,12 +542,12 @@ void gmgpolar::prepare_op_levels()
                                 std::vector<int>(v_level[l]->A_Zebra_c_row[smoother][ij]);
                             v_level[l]->A_Zebra_v_LU_row[smoother][ij] =
                                 std::vector<double>(v_level[l]->A_Zebra_v_row[smoother][ij]);
-                            if (gyro::icntl[Param::verbose] > 2)
+                            if (gyro::icntl[Param::verbose] > 3)
                                 std::cout << "Across factorization...";
 #ifdef USE_MUMPS
                             if (gyro::icntl[Param::optimized] == 0) {
 #endif
-                                if (gyro::icntl[Param::verbose] > 2)
+                                if (gyro::icntl[Param::verbose] > 3)
                                     std::cout << "using in-house direct solver.\n";
                                 // Factorization
                                 v_level[l]->facto_gaussian_elimination(v_level[l]->A_Zebra_r_LU_row[smoother][ij],
@@ -586,7 +591,7 @@ void gmgpolar::prepare_op_levels()
                 v_level[l]->define_nz_P();
 
                 // // Number of nonzeros in the matrix
-                if (gyro::icntl[Param::verbose] > 2)
+                if (gyro::icntl[Param::verbose] > 3)
                     std::cout << "nz_P: " << v_level[l]->nz_P << "\n";
                 v_level[l]->ri_prol = std::vector<int>(v_level[l]->nz_P);
                 v_level[l]->ci_prol = std::vector<int>(v_level[l]->nz_P);
@@ -594,7 +599,7 @@ void gmgpolar::prepare_op_levels()
                 v_level[l]->build_prolongation_bi();
 
                 // // Number of nonzeros in the matrix
-                if (gyro::icntl[Param::verbose] > 2)
+                if (gyro::icntl[Param::verbose] > 3)
                     std::cout << "nz_P_inj: " << v_level[l]->nz_P_inj << "\n";
                 v_level[l]->ri_prol_inj = std::vector<int>(v_level[l]->nz_P_inj);
                 v_level[l]->ci_prol_inj = std::vector<int>(v_level[l]->nz_P_inj);
@@ -602,7 +607,7 @@ void gmgpolar::prepare_op_levels()
                 v_level[l]->build_prolongation_inj();
 
                 // // Number of nonzeros in the matrix
-                if (gyro::icntl[Param::verbose] > 2)
+                if (gyro::icntl[Param::verbose] > 3)
                     std::cout << "nz_P_ex: " << v_level[l]->nz_P_ex << "\n";
                 v_level[l]->ri_prol_ex = std::vector<int>(v_level[l]->nz_P_ex);
                 v_level[l]->ci_prol_ex = std::vector<int>(v_level[l]->nz_P_ex);
@@ -615,7 +620,7 @@ void gmgpolar::prepare_op_levels()
         }
     }
 
-    if (gyro::icntl[Param::verbose] > 2)
+    if (gyro::icntl[Param::verbose] > 3)
         std::cout << "Operators created.\n";
 
     t       = t_setup;
diff --git a/src/smoother.cpp b/src/smoother.cpp
index 3e269857..126c9940 100644
--- a/src/smoother.cpp
+++ b/src/smoother.cpp
@@ -82,7 +82,7 @@ void level::define_line_splitting()
     }
     delete_circles = i; // delete_circles = delete_circles(end);
 
-    if (gyro::icntl[Param::verbose] > 1)
+    if (gyro::icntl[Param::verbose] > 2)
         std::cout << "Shifting from circle to radial at radius " << delete_circles << "\n";
 } /* ----- end of destructor level::define_line_splitting ----- */
 
diff --git a/src/smoother0.cpp b/src/smoother0.cpp
index 8932cec8..e6c21258 100644
--- a/src/smoother0.cpp
+++ b/src/smoother0.cpp
@@ -81,7 +81,7 @@ void level::multigrid_smoothing0(int smoother)
         u_sc = solve_mumps(mumps_A_Zebra[smoother], f_total);
 #endif
 
-    if (gyro::icntl[Param::verbose] > 3)
+    if (gyro::icntl[Param::verbose] > 5)
         gyro::disp(u_sc, "u_sc");
 
     t_Asc += TOC;