From e416477aba5b41cd80b8655de978525e0c016937 Mon Sep 17 00:00:00 2001
From: Robin Dunn <robin@simpetus.com>
Date: Fri, 31 Jul 2020 13:29:50 -0700
Subject: [PATCH] Add verbosity flag for MPB (#1302)

* Return the prior verbosity level setting. Change the parameter name to `level`

* Add verbosity level global variable, and a Python function for setting it

* Switch the verbosity value and related code to a class

* switch parameter name to `level`, prep for reuse with other cvars

* Move verbosity module to meep

* Rename module to avoid name conflicts with names of instances

* Restore returning the former level from set()

* Use the Verbosity class in meep

* Add ability to set the initial level

* Make sure we use mpb's verbosity flag, and not colliding with meep's

* Move the creation of the verbosity obj so it's available when submodules are imported

* Add explicit rules for the meep.mpb .py files so they will be copied into place with a simple `make`

* Add verbosity checks in solver.py

* Add verbosity checks in libpympb/pympb.cpp

* Revert a Makefile change

* Handle the case where a cvar object does not yet have a verbosity attribute.
This will help with transitioning while the feature rolls out.

* Python2 doesn't have f-strings

* empty commit to trigger build

* Support building with an older MPB which doesn't have the mpb_verbosity global

* Set the verbosity default to 1

* match MPB

Co-authored-by: Steven G. Johnson <stevenj@mit.edu>
---
 libpympb/pympb.cpp                     | 165 ++++++++++++++++---------
 libpympb/pympb.hpp                     |  10 ++
 python/Makefile.am                     |   4 +-
 python/adjoint/optimization_problem.py |  56 ++++-----
 python/mpb.i                           |   6 +
 python/simulation.py                   |  27 ++--
 python/solver.py                       | 103 ++++++++-------
 python/verbosity_mgr.py                | 105 ++++++++++++++++
 python/visualization.py                |   6 +-
 9 files changed, 332 insertions(+), 150 deletions(-)
 create mode 100644 python/verbosity_mgr.py

diff --git a/libpympb/pympb.cpp b/libpympb/pympb.cpp
index 8003ba27c..6bfe04ae7 100644
--- a/libpympb/pympb.cpp
+++ b/libpympb/pympb.cpp
@@ -8,6 +8,14 @@
 #include "pympb.hpp"
 #include "meep/mympi.hpp"
 
+// If the MPB lib is not new enough to have the mpb_verbosity global then make
+// one here to give the swig wrapper something to link with.
+#if MPB_VERSION_MAJOR > 1 || (MPB_VERSION_MAJOR == 1 && MPB_VERSION_MINOR >= 11)
+// do nothing, libmpb should have the mpb_verbosity symbol
+#else
+extern "C" int mpb_verbosity = 2;
+#endif
+
 // xyz_loop.h
 #ifndef HAVE_MPI
 #define LOOP_XYZ(md)                                                                               \
@@ -55,6 +63,8 @@ typedef mpb_real real; // needed for the CASSIGN macros below
 
 namespace py_mpb {
 
+
+
 // TODO: Placeholder
 int mpb_comm;
 
@@ -684,20 +694,23 @@ void mode_solver::init(int p, bool reset_fields, geometric_object_list geometry,
   n[1] = grid_size.y;
   n[2] = grid_size.z;
 
-  if (target_freq != 0.0) { meep::master_printf("Target frequency is %g\n", target_freq); }
+  if (target_freq != 0.0 && mpb_verbosity >= 2) { meep::master_printf("Target frequency is %g\n", target_freq); }
 
   int true_rank = n[2] > 1 ? 3 : (n[1] > 1 ? 2 : 1);
   if (true_rank < dimensions) { dimensions = true_rank; }
   else if (true_rank > dimensions) {
-    meep::master_printf("WARNING: rank of grid is > dimensions.\n"
+    if (mpb_verbosity >= 2)
+      meep::master_printf("WARNING: rank of grid is > dimensions.\n"
                         "         setting extra grid dims. to 1.\n");
     // force extra dims to be 1
     if (dimensions <= 2) { n[2] = 1; }
     if (dimensions <= 1) { n[1] = 1; }
   }
 
-  meep::master_printf("Working in %d dimensions.\n", dimensions);
-  meep::master_printf("Grid size is %d x %d x %d.\n", n[0], n[1], n[2]);
+  if (mpb_verbosity >= 2) {
+    meep::master_printf("Working in %d dimensions.\n", dimensions);
+    meep::master_printf("Grid size is %d x %d x %d.\n", n[0], n[1], n[2]);
+  }
 
   int block_size;
 
@@ -708,7 +721,8 @@ void mode_solver::init(int p, bool reset_fields, geometric_object_list geometry,
       block_size = (num_bands - block_size - 1) / (-block_size);
       block_size = (num_bands + block_size - 1) / block_size;
     }
-    meep::master_printf("Solving for %d bands at a time.\n", block_size);
+    if (mpb_verbosity >= 2)
+      meep::master_printf("Solving for %d bands at a time.\n", block_size);
   }
   else {
     block_size = num_bands;
@@ -748,7 +762,8 @@ void mode_solver::init(int p, bool reset_fields, geometric_object_list geometry,
     srand(314159); // * (rank + 1));
   }
 
-  meep::master_printf("Creating Maxwell data...\n");
+  if (mpb_verbosity >= 2)
+    meep::master_printf("Creating Maxwell data...\n");
   mdata = create_maxwell_data(n[0], n[1], n[2], &local_N, &N_start, &alloc_N, block_size,
                               NUM_FFT_BANDS);
 
@@ -759,7 +774,8 @@ void mode_solver::init(int p, bool reset_fields, geometric_object_list geometry,
   if (check_maxwell_dielectric(mdata, 0)) { meep::abort("invalid dielectric function for MPB"); }
 
   if (!have_old_fields) {
-    meep::master_printf("Allocating fields...\n");
+    if (mpb_verbosity >= 2)
+      meep::master_printf("Allocating fields...\n");
 
     int N = n[0] * n[1] * n[2];
     int c = 2;
@@ -798,32 +814,39 @@ void mode_solver::init_epsilon(geometric_object_list *geometry) {
   mpb_real no_size_y = geometry_lattice.size.y == 0 ? 1 : geometry_lattice.size.y;
   mpb_real no_size_z = geometry_lattice.size.z == 0 ? 1 : geometry_lattice.size.z;
 
-  meep::master_printf("Mesh size is %d.\n", mesh_size);
+  if (mpb_verbosity >= 2)
+    meep::master_printf("Mesh size is %d.\n", mesh_size);
 
   Rm.c0 = vector3_scale(no_size_x, geometry_lattice.basis.c0);
   Rm.c1 = vector3_scale(no_size_y, geometry_lattice.basis.c1);
   Rm.c2 = vector3_scale(no_size_z, geometry_lattice.basis.c2);
 
-  meep::master_printf("Lattice vectors:\n");
-  meep::master_printf("     (%g, %g, %g)\n", Rm.c0.x, Rm.c0.y, Rm.c0.z);
-  meep::master_printf("     (%g, %g, %g)\n", Rm.c1.x, Rm.c1.y, Rm.c1.z);
-  meep::master_printf("     (%g, %g, %g)\n", Rm.c2.x, Rm.c2.y, Rm.c2.z);
+  if (mpb_verbosity >= 2) {
+    meep::master_printf("Lattice vectors:\n");
+    meep::master_printf("     (%g, %g, %g)\n", Rm.c0.x, Rm.c0.y, Rm.c0.z);
+    meep::master_printf("     (%g, %g, %g)\n", Rm.c1.x, Rm.c1.y, Rm.c1.z);
+    meep::master_printf("     (%g, %g, %g)\n", Rm.c2.x, Rm.c2.y, Rm.c2.z);
+  }
 
   vol = fabs(matrix3x3_determinant(Rm));
-  meep::master_printf("Cell volume = %g\n", vol);
+  if (mpb_verbosity >= 2)
+    meep::master_printf("Cell volume = %g\n", vol);
 
   Gm = matrix3x3_inverse(matrix3x3_transpose(Rm));
-  meep::master_printf("Reciprocal lattice vectors (/ 2 pi):\n");
-  meep::master_printf("     (%g, %g, %g)\n", Gm.c0.x, Gm.c0.y, Gm.c0.z);
-  meep::master_printf("     (%g, %g, %g)\n", Gm.c1.x, Gm.c1.y, Gm.c1.z);
-  meep::master_printf("     (%g, %g, %g)\n", Gm.c2.x, Gm.c2.y, Gm.c2.z);
+  if (mpb_verbosity >= 2) {
+    meep::master_printf("Reciprocal lattice vectors (/ 2 pi):\n");
+    meep::master_printf("     (%g, %g, %g)\n", Gm.c0.x, Gm.c0.y, Gm.c0.z);
+    meep::master_printf("     (%g, %g, %g)\n", Gm.c1.x, Gm.c1.y, Gm.c1.z);
+    meep::master_printf("     (%g, %g, %g)\n", Gm.c2.x, Gm.c2.y, Gm.c2.z);
+  }
 
   matrix3x3_to_arr(R, Rm);
   matrix3x3_to_arr(G, Gm);
 
   geom_fix_object_list(*geometry);
 
-  meep::master_printf("Geometric objects:\n");
+  if (mpb_verbosity >= 2)
+    meep::master_printf("Geometric objects:\n");
   if (meep::am_master()) {
     for (int i = 0; i < geometry->num_items; ++i) {
 
@@ -863,16 +886,18 @@ void mode_solver::init_epsilon(geometric_object_list *geometry) {
   }
 
   if (verbose && meep::am_master()) {
-    meep::master_printf("Geometry object bounding box tree:\n");
+    if (mpb_verbosity >= 2)
+      meep::master_printf("Geometry object bounding box tree:\n");
     display_geom_box_tree(5, geometry_tree);
   }
 
   int tree_depth;
   int tree_nobjects;
   geom_box_tree_stats(geometry_tree, &tree_depth, &tree_nobjects);
-  meep::master_printf("Geometric object tree has depth %d and %d object nodes"
-                      " (vs. %d actual objects)\n",
-                      tree_depth, tree_nobjects, geometry->num_items);
+  if (mpb_verbosity >= 2)
+    meep::master_printf("Geometric object tree has depth %d and %d object nodes"
+                        " (vs. %d actual objects)\n",
+                        tree_depth, tree_nobjects, geometry->num_items);
 
   // restricted_tree = geometry_tree;
 
@@ -894,12 +919,14 @@ void mode_solver::reset_epsilon(geometric_object_list *geometry) {
   // if (!mu_input_file.empty()) {
   // }
 
-  meep::master_printf("Initializing epsilon function...\n");
+  if (mpb_verbosity >= 2)
+    meep::master_printf("Initializing epsilon function...\n");
   set_maxwell_dielectric(mdata, mesh, R, G, dielectric_function, mean_epsilon_func,
                          static_cast<void *>(this));
 
   if (has_mu(geometry)) {
-    meep::master_printf("Initializing mu function...\n");
+    if (mpb_verbosity >= 2)
+      meep::master_printf("Initializing mu function...\n");
     eps = false;
     set_maxwell_mu(mdata, mesh, R, G, dielectric_function, mean_epsilon_func,
                    static_cast<void *>(this));
@@ -956,7 +983,8 @@ void mode_solver::set_parity(integer p) {
     meep::master_fprintf(stderr, "k vector incompatible with parity\n");
     exit(EXIT_FAILURE);
   }
-  meep::master_printf("Solving for band polarization: %s.\n", parity_string(mdata));
+  if (mpb_verbosity >= 2)
+    meep::master_printf("Solving for band polarization: %s.\n", parity_string(mdata));
 
   last_parity = p;
   set_kpoint_index(0); /* reset index */
@@ -974,7 +1002,8 @@ void mode_solver::set_kpoint_index(int i) { kpoint_index = i; }
 void mode_solver::randomize_fields() {
 
   if (!mdata) { return; }
-  meep::master_printf("Initializing fields to random numbers...\n");
+  if (mpb_verbosity >= 2)
+    meep::master_printf("Initializing fields to random numbers...\n");
 
   for (int i = 0; i < H.n * H.p; ++i) {
     ASSIGN_SCALAR(H.data[i], rand() * 1.0 / RAND_MAX, rand() * 1.0 / RAND_MAX);
@@ -987,12 +1016,14 @@ void mode_solver::solve_kpoint(vector3 kvector) {
   // special handling of this k
   if (vector3_norm(kvector) < 1e-10) { kvector.x = kvector.y = kvector.z = 0; }
 
-  meep::master_printf("solve_kpoint (%g,%g,%g):\n", kvector.x, kvector.y, kvector.z);
+  if (mpb_verbosity >= 2)
+    meep::master_printf("solve_kpoint (%g,%g,%g):\n", kvector.x, kvector.y, kvector.z);
 
   curfield_reset();
 
   if (num_bands == 0) {
-    meep::master_printf("  num-bands is zero, not solving for any bands\n");
+    if (mpb_verbosity >= 2)
+      meep::master_printf("  num-bands is zero, not solving for any bands\n");
     return;
   }
 
@@ -1003,14 +1034,16 @@ void mode_solver::solve_kpoint(vector3 kvector) {
 
   // If this is the first k point, print out a header line for the frequency
   // grep data.
-  if (!kpoint_index && meep::am_master()) {
-    meep::master_printf("%sfreqs:, k index, k1, k2, k3, kmag/2pi", parity_string(mdata));
+  if (mpb_verbosity >= 2) {
+    if (!kpoint_index && meep::am_master()) {
+      meep::master_printf("%sfreqs:, k index, k1, k2, k3, kmag/2pi", parity_string(mdata));
 
-    for (int i = 0; i < num_bands; ++i) {
-      meep::master_printf(", %s%sband %d", parity_string(mdata),
-                          mdata->parity == NO_PARITY ? "" : " ", i + 1);
+      for (int i = 0; i < num_bands; ++i) {
+        meep::master_printf(", %s%sband %d", parity_string(mdata),
+                            mdata->parity == NO_PARITY ? "" : " ", i + 1);
+      }
+      meep::master_printf("\n");
     }
-    meep::master_printf("\n");
   }
 
   cur_kvector = kvector;
@@ -1023,7 +1056,7 @@ void mode_solver::solve_kpoint(vector3 kvector) {
 
   // TODO: Get flags from python
   int flags = EIGS_DEFAULT_FLAGS;
-  if (verbose) { flags |= EIGS_VERBOSE; }
+  if (verbose || mpb_verbosity >= 2) { flags |= EIGS_VERBOSE; }
 
   // Constant (zero frequency) bands at k=0 are handled specially, so remove
   // them from the solutions for the eigensolver.
@@ -1068,7 +1101,8 @@ void mode_solver::solve_kpoint(vector3 kvector) {
       evectmatrix_resize(&Hblock, num_bands - ib, 0);
     }
 
-    meep::master_printf("Solving for bands %d to %d...\n", ib + 1, ib + Hblock.p);
+    if (mpb_verbosity >= 2)
+      meep::master_printf("Solving for bands %d to %d...\n", ib + 1, ib + Hblock.p);
 
     constraints = NULL;
     constraints = evect_add_constraint(constraints, maxwell_parity_constraint, (void *)mdata);
@@ -1126,13 +1160,14 @@ void mode_solver::solve_kpoint(vector3 kvector) {
 
     evect_destroy_constraints(constraints);
 
-    meep::master_printf("Finished solving for bands %d to %d after %d iterations.\n", ib + 1,
-                        ib + Hblock.p, num_iters);
+    if (mpb_verbosity >= 2)
+      meep::master_printf("Finished solving for bands %d to %d after %d iterations.\n", ib + 1,
+                          ib + Hblock.p, num_iters);
 
     total_iters += num_iters * Hblock.p;
   }
 
-  if (num_bands - ib0 > Hblock.alloc_p) {
+  if (num_bands - ib0 > Hblock.alloc_p && mpb_verbosity >= 2) {
     meep::master_printf("Finished k-point with %g mean iterations/band.\n",
                         total_iters * 1.0 / num_bands);
   }
@@ -1169,15 +1204,18 @@ void mode_solver::solve_kpoint(vector3 kvector) {
 
   set_kpoint_index(kpoint_index + 1);
 
-  meep::master_printf("%sfreqs:, %d, %g, %g, %g, %g", parity_string(mdata), kpoint_index,
-                      (double)k[0], (double)k[1], (double)k[2],
-                      vector3_norm(matrix3x3_vector3_mult(Gm, kvector)));
+  if (mpb_verbosity >= 2)
+    meep::master_printf("%sfreqs:, %d, %g, %g, %g, %g", parity_string(mdata), kpoint_index,
+                        (double)k[0], (double)k[1], (double)k[2],
+                        vector3_norm(matrix3x3_vector3_mult(Gm, kvector)));
 
   for (int i = 0; i < num_bands; ++i) {
     freqs[i] = negative_epsilon_ok ? eigvals[i] : sqrt(eigvals[i]);
-    meep::master_printf(", %g", freqs[i]);
+    if (mpb_verbosity >= 2)
+      meep::master_printf(", %g", freqs[i]);
   }
-  meep::master_printf("\n");
+  if (mpb_verbosity >= 2)
+    meep::master_printf("\n");
 
   eigensolver_flops = evectmatrix_flops;
 }
@@ -1230,11 +1268,12 @@ void mode_solver::get_epsilon() {
   eps_mean /= N;
   eps_inv_mean = N / eps_inv_mean;
 
-  meep::master_printf("epsilon: %g-%g, mean %g, harm. mean %g, %g%% > 1, %g%% \"fill\"\n", eps_low,
-                      eps_high, eps_mean, eps_inv_mean, (100.0 * fill_count) / N,
-                      eps_high == eps_low ? 100.0
-                                          : 100.0 * (eps_mean - eps_low) / (eps_high - eps_low));
-}
+  if (mpb_verbosity >= 2)
+    meep::master_printf("epsilon: %g-%g, mean %g, harm. mean %g, %g%% > 1, %g%% \"fill\"\n", eps_low,
+                        eps_high, eps_mean, eps_inv_mean, (100.0 * fill_count) / N,
+                        eps_high == eps_low ? 100.0
+                                            : 100.0 * (eps_mean - eps_low) / (eps_high - eps_low));
+  }
 
 /* get the mu function, and compute some statistics */
 void mode_solver::get_mu() {
@@ -1285,10 +1324,11 @@ void mode_solver::get_mu() {
   eps_mean /= N;
   mu_inv_mean = N / mu_inv_mean;
 
-  meep::master_printf("mu: %g-%g, mean %g, harm. mean %g, %g%% > 1, %g%% \"fill\"\n", eps_low,
-                      eps_high, eps_mean, mu_inv_mean, (100.0 * fill_count) / N,
-                      eps_high == eps_low ? 100.0
-                                          : 100.0 * (eps_mean - eps_low) / (eps_high - eps_low));
+  if (mpb_verbosity >= 2)
+    meep::master_printf("mu: %g-%g, mean %g, harm. mean %g, %g%% > 1, %g%% \"fill\"\n", eps_low,
+                        eps_high, eps_mean, mu_inv_mean, (100.0 * fill_count) / N,
+                        eps_high == eps_low ? 100.0
+                                            : 100.0 * (eps_mean - eps_low) / (eps_high - eps_low));
 }
 
 void mode_solver::curfield_reset() {
@@ -1613,12 +1653,16 @@ std::vector<mpb_real> mode_solver::compute_field_energy() {
   mpb_real comp_sum[6];
   mpb_real energy_sum = compute_field_energy_internal(comp_sum);
 
-  meep::master_printf("%c-energy-components:, %d, %d", curfield_type, kpoint_index, curfield_band);
+  if (mpb_verbosity >= 2)
+    meep::master_printf("%c-energy-components:, %d, %d", curfield_type, kpoint_index, curfield_band);
   for (int i = 0; i < 6; ++i) {
     comp_sum[i] /= (energy_sum == 0 ? 1 : energy_sum);
-    if (i % 2 == 1) { meep::master_printf(", %g", comp_sum[i] + comp_sum[i - 1]); }
+    if (i % 2 == 1 && mpb_verbosity >= 2) {
+      meep::master_printf(", %g", comp_sum[i] + comp_sum[i - 1]);
+      }
   }
-  meep::master_printf("\n");
+  if (mpb_verbosity >= 2)
+    meep::master_printf("\n");
 
   /* The return value is a list of 7 items: the total energy,
      followed by the 6 elements of the comp_sum array (the fraction
@@ -2032,9 +2076,10 @@ void mode_solver::fix_field_phase() {
 
   ASSIGN_SCALAR(phase, SCALAR_RE(phase) * maxabs_sign, SCALAR_IM(phase) * maxabs_sign);
 
-  meep::master_printf("Fixing %c-field (band %d) phase by %g + %gi; "
-                      "max ampl. = %g\n",
-                      curfield_type, curfield_band, SCALAR_RE(phase), SCALAR_IM(phase), maxabs);
+  if (mpb_verbosity >= 2)
+    meep::master_printf("Fixing %c-field (band %d) phase by %g + %gi; "
+                        "max ampl. = %g\n",
+                        curfield_type, curfield_band, SCALAR_RE(phase), SCALAR_IM(phase), maxabs);
 
   /* Now, multiply everything by this phase, *including* the
      stored "raw" eigenvector in H, so that any future fields
@@ -2614,7 +2659,7 @@ void map_data(mpb_real *d_in_re, int size_in_re, mpb_real *d_in_im, int size_in_
 #undef IN_INDEX
       }
 
-  if (verbose) {
+  if (verbose || mpb_verbosity >= 2) {
     meep::master_printf("real part range: %g .. %g\n", min_out_re, max_out_re);
     if (size_out_im > 0) meep::master_printf("imag part range: %g .. %g\n", min_out_im, max_out_im);
   }
diff --git a/libpympb/pympb.hpp b/libpympb/pympb.hpp
index f59abdf8e..76499cc3c 100644
--- a/libpympb/pympb.hpp
+++ b/libpympb/pympb.hpp
@@ -9,6 +9,16 @@
 
 namespace py_mpb {
 
+// A global "verbosity" level.
+//    0: minimal output
+//    1: a little (default)
+//    2: a lot
+//    3: debugging
+//
+// This is redeclared here so SWIG will see it.
+extern "C" int mpb_verbosity;
+
+
 #define TWOPI 6.2831853071795864769252867665590057683943388
 
 void map_data(mpb_real *d_in_re, int size_in_re, mpb_real *d_in_im, int size_in_im, int n_in[3],
diff --git a/python/Makefile.am b/python/Makefile.am
index 7e5ceadcc..2b213809d 100644
--- a/python/Makefile.am
+++ b/python/Makefile.am
@@ -197,7 +197,9 @@ HL_IFACE =                  \
     $(srcdir)/simulation.py \
     $(srcdir)/source.py     \
     $(srcdir)/visualization.py  \
-    $(srcdir)/materials.py
+    $(srcdir)/materials.py  \
+    $(srcdir)/verbosity_mgr.py
+
 
 pkgpython_PYTHON = __init__.py $(HL_IFACE)
 
diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py
index 2c578197c..0026c4a9d 100644
--- a/python/adjoint/optimization_problem.py
+++ b/python/adjoint/optimization_problem.py
@@ -24,7 +24,7 @@ def get_gradient(self,fields_a,fields_f,frequencies,geom_list,f):
         grad = np.zeros((self.num_design_params*num_freqs,)) # preallocate
 
         # compute the gradient
-        mp._get_gradient(grad,fields_a,fields_f,self.volume,np.array(frequencies),geom_list,f) 
+        mp._get_gradient(grad,fields_a,fields_f,self.volume,np.array(frequencies),geom_list,f)
 
         return np.squeeze(grad.reshape(self.num_design_params,num_freqs,order='F'))
 
@@ -42,7 +42,7 @@ class OptimizationProblem(object):
 
     """
 
-    def __init__(self, 
+    def __init__(self,
                 simulation,
                 objective_functions,
                 objective_arguments,
@@ -70,7 +70,7 @@ def __init__(self,
             self.design_regions = design_regions
         else:
             self.design_regions = [design_regions]
-        
+
         self.num_design_params = [ni.num_design_params for ni in self.design_regions]
         self.num_design_regions = len(self.design_regions)
 
@@ -83,7 +83,7 @@ def __init__(self,
                 self.nf = nf
                 self.frequencies = [fcen]
             else:
-                fmax = fcen+0.5*df 
+                fmax = fcen+0.5*df
                 fmin = fcen-0.5*df
                 dfreq = (fmax-fmin)/(nf-1)
                 self.frequencies = np.linspace(fmin, fmin+dfreq*nf, num=nf, endpoint=False)
@@ -98,7 +98,7 @@ def __init__(self,
         self.decay_fields=decay_fields
         self.decay_dt=decay_dt
         self.minimum_run_time=minimum_run_time
-        
+
         # store sources for finite difference estimations
         self.forward_sources = self.sim.sources
 
@@ -115,7 +115,7 @@ def __call__(self, rho_vector=None, need_value=True, need_gradient=True):
             self.update_design(rho_vector=rho_vector)
 
         # Run forward run if requested
-        if need_value and self.current_state == "INIT": 
+        if need_value and self.current_state == "INIT":
             print("Starting forward run...")
             self.forward_run()
 
@@ -162,7 +162,7 @@ def _df(x=None):
             return df
 
         return _f, _df
-    
+
     def prepare_forward_run(self):
         # prepare forward run
         self.sim.reset_meep()
@@ -226,7 +226,7 @@ def adjoint_run(self,objective_idx=0):
         self.sim.change_sources(self.adjoint_sources)
 
         # register design flux
-        # TODO use yee grid directly 
+        # TODO use yee grid directly
         self.design_region_monitors = [self.sim.add_dft_fields([mp.Ex,mp.Ey,mp.Ez],self.frequencies,where=dr.volume,yee_grid=False) for dr in self.design_regions]
 
         # Adjoint run
@@ -245,7 +245,7 @@ def adjoint_run(self,objective_idx=0):
     def calculate_gradient(self):
         # Iterate through all design regions and calculate gradient
         self.gradient = [[dr.get_gradient(self.a_E[ar][dri],self.d_E[dri],self.frequencies,self.sim.geometry,self.sim.fields) for dri, dr in enumerate(self.design_regions)] for ar in range(len(self.objective_functions))]
-        
+
         # Cleanup list of lists
         if len(self.gradient) == 1:
             self.gradient = self.gradient[0] # only one objective function
@@ -256,7 +256,7 @@ def calculate_gradient(self):
                 self.gradient = [g[0] for g in self.gradient] # multiple objective functions bu one design region
         # Return optimizer's state to initialization
         self.current_state = "INIT"
-    
+
     def calculate_fd_gradient(self,num_gradients=1,db=1e-4,design_variables_idx=0,filter=None):
         '''
         Estimate central difference gradients.
@@ -272,7 +272,7 @@ def calculate_fd_gradient(self,num_gradients=1,db=1e-4,design_variables_idx=0,fi
 
         Returns
         -----------
-        fd_gradient ... : lists 
+        fd_gradient ... : lists
             [number of objective functions][number of gradients]
 
         '''
@@ -292,25 +292,25 @@ def calculate_fd_gradient(self,num_gradients=1,db=1e-4,design_variables_idx=0,fi
         fd_gradient_idx = np.random.choice(self.num_design_params[design_variables_idx],num_gradients,replace=False)
 
         for k in fd_gradient_idx:
-            
+
             b0 = np.ones((self.num_design_params[design_variables_idx],))
             b0[:] = (self.design_regions[design_variables_idx].design_parameters.design_parameters)
             # -------------------------------------------- #
             # left function evaluation
             # -------------------------------------------- #
             self.sim.reset_meep()
-            
+
             # assign new design vector
             b0[k] -= db
             self.design_regions[design_variables_idx].update_design_parameters(b0)
-            
+
             # initialize design monitors
             self.forward_monitors = []
             for m in self.objective_arguments:
                 self.forward_monitors.append(m.register_monitors(self.frequencies))
-            
+
             self.sim.run(until_after_sources=stop_when_dft_decayed(self.sim, self.forward_monitors, self.decay_dt, self.decay_fields, self.fcen_idx, self.decay_by, self.minimum_run_time))
-            
+
             # record final objective function value
             results_list = []
             for m in self.objective_arguments:
@@ -330,10 +330,10 @@ def calculate_fd_gradient(self,num_gradients=1,db=1e-4,design_variables_idx=0,fi
             self.forward_monitors = []
             for m in self.objective_arguments:
                 self.forward_monitors.append(m.register_monitors(self.frequencies))
-            
+
             # add monitor used to track dft convergence
             self.sim.run(until_after_sources=stop_when_dft_decayed(self.sim, self.forward_monitors, self.decay_dt, self.decay_fields, self.fcen_idx, self.decay_by, self.minimum_run_time))
-            
+
             # record final objective function value
             results_list = []
             for m in self.objective_arguments:
@@ -344,13 +344,13 @@ def calculate_fd_gradient(self,num_gradients=1,db=1e-4,design_variables_idx=0,fi
             # estimate derivative
             # -------------------------------------------- #
             fd_gradient.append( [np.squeeze((fp[fi] - fm[fi]) / (2*db)) for fi in range(len(self.objective_functions))] )
-        
+
         # Cleanup singleton dimensions
         if len(fd_gradient) == 1:
             fd_gradient = fd_gradient[0]
 
         return fd_gradient, fd_gradient_idx
-    
+
     def update_design(self, rho_vector):
         """Update the design permittivity function.
 
@@ -360,7 +360,7 @@ def update_design(self, rho_vector):
             if np.array(rho_vector[bi]).ndim > 1:
                 raise ValueError("Each vector of design variables must contain only one dimension.")
             b.update_design_parameters(rho_vector[bi])
-        
+
         self.sim.reset_meep()
         self.current_state = "INIT"
     def get_objective_arguments(self):
@@ -368,7 +368,7 @@ def get_objective_arguments(self):
         '''
         objective_args_evaluation = [m.get_evaluation() for m in self.objective_arguments]
         return objective_args_evaluation
-        
+
     def plot2D(self,init_opt=False, **kwargs):
         """Produce a graphical visualization of the geometry and/or fields,
            as appropriately autodetermined based on the current state of
@@ -397,7 +397,7 @@ def stop_when_dft_decayed(simob, mon, dt, c, fcen_idx, decay_by, minimum_run_tim
         x.append(len(xi))
         y.append(len(yi))
         z.append(len(zi))
-    
+
     # Record data in closure so that we can persitently edit
     closure = {
         'previous_fields': [np.ones((x[mi],y[mi],z[mi],len(c)),dtype=np.complex128) for mi, m in enumerate(mon)],
@@ -405,7 +405,7 @@ def stop_when_dft_decayed(simob, mon, dt, c, fcen_idx, decay_by, minimum_run_tim
     }
 
     def _stop(sim):
-        
+
         if sim.round_time() <= dt + closure['t0']:
             return False
         else:
@@ -427,11 +427,11 @@ def _stop(sim):
             relative_change = np.mean(relative_change) # average across monitors
             closure['previous_fields'] = current_fields
             closure['t0'] = sim.round_time()
-            
-            if mp.cvar.verbosity > 0:
+
+            if mp.verbosity > 0:
                 fmt = "DFT decay(t = {0:1.1f}): {1:0.4e}"
                 print(fmt.format(sim.meep_time(), np.real(relative_change)))
-            return relative_change <= decay_by and sim.round_time() >= minimum_run_time 
+            return relative_change <= decay_by and sim.round_time() >= minimum_run_time
     return _stop
 
 
@@ -440,7 +440,7 @@ def atleast_3d(*arys):
     '''
     Modified version of numpy's `atleast_3d`
 
-    Keeps one dimensional array data in first dimension, as 
+    Keeps one dimensional array data in first dimension, as
     opposed to moving it to the second dimension as numpy's
     version does. Keeps the meep dimensionality convention.
 
diff --git a/python/mpb.i b/python/mpb.i
index 0e54b360d..f96319cd9 100644
--- a/python/mpb.i
+++ b/python/mpb.i
@@ -346,9 +346,14 @@ static mpb_real field_integral_energy_callback(mpb_real energy, mpb_real epsilon
 
 %apply double { number };
 
+%rename(verbosity) mpb_verbosity;
 %include "pympb.hpp"
 
+
 %pythoncode %{
+    from meep.verbosity_mgr import Verbosity
+    verbosity = Verbosity(_mpb.cvar, 1)
+
     from .solver import (
         MPBArray,
         ModeSolver,
@@ -390,4 +395,5 @@ static mpb_real field_integral_energy_callback(mpb_real energy, mpb_real epsilon
     from .mpb_data import (
         MPBData,
     )
+
 %}
diff --git a/python/simulation.py b/python/simulation.py
index 390dd639e..9c97bc408 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -20,7 +20,7 @@
 from meep.geom import Vector3, init_do_averaging
 from meep.source import EigenModeSource, check_positive
 import meep.visualization as vis
-
+from meep.verbosity_mgr import Verbosity
 
 try:
     basestring
@@ -34,6 +34,8 @@
 except ImportError:
     do_progress = False
 
+verbosity = Verbosity(mp.cvar, 1)
+
 
 # Send output from Meep, ctlgeom, and MPB to Python's stdout
 mp.set_meep_printf_callback(mp.py_master_printf_wrap)
@@ -1556,7 +1558,7 @@ def _compute_fragment_stats(self, gv):
         return stats
 
     def _init_structure(self, k=False):
-        if mp.cvar.verbosity > 0:
+        if verbosity > 0:
             print('-' * 11)
             print('Initializing structure...')
 
@@ -1579,7 +1581,7 @@ def _init_structure(self, k=False):
         if self.collect_stats and isinstance(self.default_material, mp.Medium):
             self.fragment_stats = self._compute_fragment_stats(gv)
 
-        if self._output_stats and isinstance(self.default_material, mp.Medium) and mp.cvar.verbosity > 0:
+        if self._output_stats and isinstance(self.default_material, mp.Medium) and verbosity > 0:
             stats = self._compute_fragment_stats(gv)
             print("STATS: aniso_eps: {}".format(stats.num_anisotropic_eps_pixels))
             print("STATS: anis_mu: {}".format(stats.num_anisotropic_mu_pixels))
@@ -1851,7 +1853,7 @@ def use_real(self):
 
         if use_real(self):
             self.fields.use_real_fields()
-        elif mp.cvar.verbosity > 0:
+        elif verbosity > 0:
             print("Meep: using complex fields.")
 
         if self.k_point:
@@ -2039,7 +2041,7 @@ def use_output_directory(self, dname=''):
         closure = {'trashed': False}
 
         def hook():
-            if mp.cvar.verbosity > 0:
+            if verbosity > 0:
                 print("Meep: using output directory '{}'".format(dname))
             self.fields.set_output_directory(dname)
             if not closure['trashed']:
@@ -2100,7 +2102,7 @@ def stop_cond(sim):
             self.progress.value = t0 + stop_time
             self.progress.description = "100% done "
 
-        if mp.cvar.verbosity > 0:
+        if verbosity > 0:
             print("run {} finished at t = {} ({} timesteps)".format(self.run_index, self.meep_time(), self.fields.t))
         self.run_index += 1
 
@@ -4057,7 +4059,7 @@ def _stop(sim):
             closure['cur_max'] = 0
             closure['t0'] = sim.round_time()
             closure['max_abs'] = max(closure['max_abs'], old_cur)
-            if closure['max_abs'] != 0 and mp.cvar.verbosity > 0:
+            if closure['max_abs'] != 0 and verbosity > 0:
                 fmt = "field decay(t = {}): {} / {} = {}"
                 print(fmt.format(sim.meep_time(), old_cur, closure['max_abs'], old_cur / closure['max_abs']))
             return old_cur <= closure['max_abs'] * decay_by
@@ -4175,7 +4177,7 @@ def _disp(sim):
                 sim.progress.value = val1
                 sim.progress.description = "{}% done ".format(int(val2))
 
-            if mp.cvar.verbosity > 0:
+            if verbosity > 0:
                 print(msg_fmt.format(val1, t, val2, val3, val4))
             closure['tlast'] = t1
 
@@ -4891,15 +4893,8 @@ def quiet(quietval=True):
     output can be enabled again by passing `False`. This sets a global variable, so the
     value will persist across runs within the same script.
     """
-    mp.cvar.verbosity = int(not quietval)
+    verbosity(int(not quietval))
 
-def verbosity(v=1):
-    """
-    Given a number `v`, specify the degree of Meep's output: `0` is quiet mode, `1` (the
-    default) is ordinary output, `2` is extra debugging output, and `3` is all debugging
-    output.
-    """
-    mp.cvar.verbosity = v
 
 def get_num_groups():
     # Lazy import
diff --git a/python/solver.py b/python/solver.py
index 44de2c428..90dd0a212 100644
--- a/python/solver.py
+++ b/python/solver.py
@@ -12,7 +12,7 @@
 import h5py
 import numpy as np
 import meep as mp
-from . import mode_solver, with_hermitian_epsilon
+from . import mode_solver, with_hermitian_epsilon, verbosity
 from meep.geom import init_do_averaging
 from meep.simulation import get_num_args
 try:
@@ -510,12 +510,13 @@ def update_brd(brd, freqs, br_start):
         return update_brd(brd, freqs, [])
 
     def output_band_range_data(self, br_data):
-        for tup, band in zip(br_data, range(1, len(br_data) + 1)):
-            fmt = "Band {} range: {} at {} to {} at {}"
-            min_band, max_band = tup
-            min_freq, min_kpoint = min_band
-            max_freq, max_kpoint = max_band
-            print(fmt.format(band, min_freq, min_kpoint, max_freq, max_kpoint))
+        if verbosity >= 2:
+            for tup, band in zip(br_data, range(1, len(br_data) + 1)):
+                fmt = "Band {} range: {} at {} to {} at {}"
+                min_band, max_band = tup
+                min_freq, min_kpoint = min_band
+                max_freq, max_kpoint = max_band
+                print(fmt.format(band, min_freq, min_kpoint, max_freq, max_kpoint))
 
     # Output any gaps in the given band ranges, and return a list of the gaps as
     # a list of (percent, freq-min, freq-max) tuples.
@@ -536,8 +537,9 @@ def ogaps(br_cur, br_rest, i, gaps):
                 else:
                     gap_size = ((200 * (br_rest_min_f - br_cur_max_f)) /
                                 (br_rest_min_f + br_cur_max_f))
-                    fmt = "Gap from band {} ({}) to band {} ({}), {}%"
-                    print(fmt.format(i, br_cur_max_f, i + 1, br_rest_min_f, gap_size))
+                    if verbosity >= 2:
+                        fmt = "Gap from band {} ({}) to band {} ({}), {}%"
+                        print(fmt.format(i, br_cur_max_f, i + 1, br_rest_min_f, gap_size))
                     return ogaps(br_rest[0], br_rest[1:], i + 1,
                                  [gap_size, br_cur_max_f, br_rest_min_f] + gaps)
         if not br_data:
@@ -640,7 +642,8 @@ def _output_complex_scalar_field(self, fname_prefix, output_k):
             self.freqs[curfield_band - 1]
         )
         fname = self._create_fname(fname, fname_prefix, True)
-        print("Outputting complex scalar field to {}...".format(fname))
+        if verbosity >= 2:
+            print("Outputting complex scalar field to {}...".format(fname))
 
         with h5py.File(fname, 'w') as f:
             f['description'] = description.encode()
@@ -673,7 +676,8 @@ def _output_vector_field(self, curfield_type, fname_prefix, output_k, component)
         )
 
         fname = self._create_fname(fname, fname_prefix, True)
-        print("Outputting fields to {}...".format(fname))
+        if verbosity >= 2:
+            print("Outputting fields to {}...".format(fname))
 
         with h5py.File(fname, 'w') as f:
             f['description'] = description.encode()
@@ -717,7 +721,8 @@ def _output_scalar_field(self, curfield_type, fname_prefix):
 
         parity_suffix = False if curfield_type in 'mn' else True
         fname = self._create_fname(fname, fname_prefix, parity_suffix)
-        print("Outputting {}...".format(fname))
+        if verbosity >= 2:
+            print("Outputting {}...".format(fname))
 
         with h5py.File(fname, 'w') as f:
             f['description'] = description.encode()
@@ -805,11 +810,12 @@ def randomize_fields(self):
 
     def display_kpoint_data(self, name, data):
         k_index = self.mode_solver.get_kpoint_index()
-        print("{}{}:, {}".format(self.parity, name, k_index), end='')
+        if verbosity >= 2:
+            print("{}{}:, {}".format(self.parity, name, k_index), end='')
 
-        for d in data:
-            print(", {}".format(d), end='')
-        print()
+            for d in data:
+                print(", {}".format(d), end='')
+            print()
 
     def display_eigensolver_stats(self):
         num_runs = len(self.eigensolver_iters)
@@ -820,20 +826,24 @@ def display_eigensolver_stats(self):
         min_iters = min(self.eigensolver_iters)
         max_iters = max(self.eigensolver_iters)
         mean_iters = np.mean(self.eigensolver_iters)
-        fmt = "eigensolver iterations for {} kpoints: {}-{}, mean = {}"
-        print(fmt.format(num_runs, min_iters, max_iters, mean_iters), end='')
+        if verbosity >= 2:
+            fmt = "eigensolver iterations for {} kpoints: {}-{}, mean = {}"
+            print(fmt.format(num_runs, min_iters, max_iters, mean_iters), end='')
 
         sorted_iters = sorted(self.eigensolver_iters)
         idx1 = num_runs // 2
         idx2 = ((num_runs + 1) // 2) - 1
         median_iters = 0.5 * (sorted_iters[idx1] + sorted_iters[idx2])
-        print(", median = {}".format(median_iters))
+        if verbosity >= 2:
+            print(", median = {}".format(median_iters))
 
         mean_flops = self.eigensolver_flops / (num_runs * mean_iters)
-        print("mean flops per iteration = {}".format(mean_flops))
+        if verbosity >= 2:
+            print("mean flops per iteration = {}".format(mean_flops))
 
         mean_time = self.total_run_time / (mean_iters * num_runs)
-        print("mean time per iteration = {} s".format(mean_time))
+        if verbosity >= 2:
+            print("mean time per iteration = {} s".format(mean_time))
 
     def _get_grid_size(self):
         grid_size = mp.Vector3(self.resolution[0] * self.geometry_lattice.size.x,
@@ -886,20 +896,22 @@ def run_parity(self, p, reset_fields, *band_functions):
 
         init_time = time.time()
 
-        print("Initializing eigensolver data")
-        print("Computing {} bands with {} tolerance".format(self.num_bands, self.tolerance))
+        if verbosity >= 2:
+            print("Initializing eigensolver data")
+            print("Computing {} bands with {} tolerance".format(self.num_bands, self.tolerance))
 
         self.init_params(p, reset_fields)
 
         if isinstance(reset_fields, basestring):
             self.load_eigenvectors(reset_fields)
 
-        print("{} k-points".format(len(self.k_points)))
+        if verbosity >= 2:
+            print("{} k-points".format(len(self.k_points)))
 
-        for kp in self.k_points:
-            print("  {}".format(kp))
+            for kp in self.k_points:
+                print("  {}".format(kp))
 
-        print("elapsed time for initialization: {}".format(time.time() - init_time))
+            print("elapsed time for initialization: {}".format(time.time() - init_time))
 
         # TODO: Split over multiple processes
         # k_split = list_split(self.k_points, self.k_split_num, self.k_split_index)
@@ -912,7 +924,8 @@ def run_parity(self, p, reset_fields, *band_functions):
                 solve_kpoint_time = time.time()
                 self.mode_solver.solve_kpoint(k)
                 self.iterations = self.mode_solver.get_iterations()
-                print("elapsed time for k point: {}".format(time.time() - solve_kpoint_time))
+                if verbosity >= 2:
+                    print("elapsed time for k point: {}".format(time.time() - solve_kpoint_time))
                 self.freqs = self.get_freqs()
                 self.all_freqs[i, :] = np.array(self.freqs)
                 self.band_range_data = self.update_band_range_data(self.band_range_data,
@@ -939,11 +952,13 @@ def run_parity(self, p, reset_fields, *band_functions):
                 self.gap_list = []
 
         end = time.time() - start
-        print("total elapsed time for run: {}".format(end))
+        if verbosity >= 2:
+            print("total elapsed time for run: {}".format(end))
         self.total_run_time += end
         self.eigensolver_flops = self.mode_solver.get_eigensolver_flops()
         self.parity = self.mode_solver.get_parity_string()
-        print("done")
+        if verbosity >= 2:
+            print("done")
 
     def run(self, *band_functions):
         self.run_parity(mp.NO_PARITY, True, *band_functions)
@@ -1009,7 +1024,8 @@ def _rootfun(k):
                 # First, look in the cached table
                 tab_val = bktab.get((b, k), None)
                 if tab_val:
-                    print("find-k {} at {}: {} (cached)".format(b, k, tab_val[0]))
+                    if verbosity >= 2:
+                        print("find-k {} at {}: {} (cached)".format(b, k, tab_val[0]))
                     return tab_val
                 # Otherwise, compute bands and cache results
                 else:
@@ -1030,7 +1046,8 @@ def _rootfun(k):
                         bktab[(_b, k)] = (_f - omega, _v)
 
                     fun = self.freqs[-1] - omega
-                    print("find-k {} at {}: {}".format(b, k, fun))
+                    if verbosity >= 2:
+                        print("find-k {} at {}: {}".format(b, k, fun))
                     return (fun, v[-1])
 
             return _rootfun
@@ -1058,14 +1075,15 @@ def bfunc(ms, b_prime):
         self.num_bands = num_bands_save
         self.k_points = kpoints_save
         ks = list(reversed(ks))
-        print("{}kvals:, {}, {}, {}".format(self.parity, omega, band_min, band_max), end='')
-        for k in korig:
-            print(", {}".format(k), end='')
-        for k in kdir1:
-            print(", {}".format(k), end='')
-        for k in ks:
-            print(", {}".format(k), end='')
-        print()
+        if verbosity >= 2:
+            print("{}kvals:, {}, {}, {}".format(self.parity, omega, band_min, band_max), end='')
+            for k in korig:
+                print(", {}".format(k), end='')
+            for k in kdir1:
+                print(", {}".format(k), end='')
+            for k in ks:
+                print(", {}".format(k), end='')
+            print()
 
         return ks
 
@@ -1221,8 +1239,9 @@ def _output(ms, which_band):
         ms.get_dfield(which_band, False)
         ms.compute_field_energy()
         energy = ms.compute_energy_in_objects(objects)
-        fmt = "dpwr:, {}, {}, {} "
-        print(fmt.format(which_band, ms.freqs[which_band - 1], energy))
+        if verbosity >= 2:
+            fmt = "dpwr:, {}, {}, {} "
+            print(fmt.format(which_band, ms.freqs[which_band - 1], energy))
         if energy >= min_energy:
             apply_band_func(ms, output_func, which_band)
 
diff --git a/python/verbosity_mgr.py b/python/verbosity_mgr.py
new file mode 100644
index 000000000..de55ee281
--- /dev/null
+++ b/python/verbosity_mgr.py
@@ -0,0 +1,105 @@
+# -*- coding: utf-8 -*-
+from __future__ import absolute_import, print_function
+
+
+class Verbosity(object):
+    """
+    A class to help make accessing and setting the global verbosity level a bit
+    more pythonic.
+
+    Verbosity levels are:
+        0: minimal output
+        1: a little
+        2: a lot (default)
+        3: debugging
+    """
+    def __init__(self, cvar=None, initial_level=None):
+        """
+        Initialize the Verbosity manager. `cvar` should be some object that has
+        a `verbosity` attribute, such as meep.cvar or mpb.cvar.
+        """
+        super(Verbosity, self).__init__()
+        self.cvar = cvar
+        if cvar is None or not hasattr(cvar, 'verbosity'):
+            # If we're not given a module.cvar (e.g., while testing) or if the
+            # cvar does not have a verbosity member (the lib hasn't been updated
+            # yet) then use a dummy object so things can still run without it.
+            class _dummy():
+                verbosity = 1
+            self.cvar = _dummy()
+        if initial_level is not None:
+            self.set(initial_level)
+
+    def get(self):
+        """
+        Returns the current verbosity level.
+        """
+        return self.cvar.verbosity
+
+    def set(self, level):
+        """
+        Validates the range, and sets the global verbosity level. Returns the
+        former value.
+        """
+        if level < 0 or level > 3:
+            raise ValueError('Only verbosity levels 0-3 are supported')
+        old = self.cvar.verbosity
+        self.cvar.verbosity = level
+        return old
+
+    def __call__(self, level):
+        """
+        Convenience for setting the verbosity level. This lets you set the level
+        by calling the instance like a function. For example, if `verbosity` is
+        an instance of this class, then it's value can be changed like this:
+
+        ```
+        verbosity(0)
+        ```
+        """
+        return self.set(level)
+
+    def __int__(self):
+        """
+        A convenience for getting the verbosity level anywhere an integer is
+        expected.
+        """
+        return self.get()
+
+    def __repr__(self):
+        return "Verbosity: level={}".format(self.get())
+
+    # Some comparison operators
+    def __gt__(self, o): return self.get() > o
+    def __lt__(self, o): return self.get() < o
+    def __eq__(self, o): return self.get() == o
+    def __ne__(self, o): return self.get() != o
+    def __le__(self, o): return self.get() <= o
+    def __ge__(self, o): return self.get() >= o
+
+
+
+
+def main():
+    # simple test code
+    verbosity = Verbosity()
+    print('initial value: {}'.format(verbosity.get()))
+
+    print('v==1: {}'.format(verbosity == 1))
+    print('v==2: {}'.format(verbosity == 2))
+
+    print('v>1: {}'.format(verbosity > 1))
+    print('v<3: {}'.format(verbosity < 3))
+    print('v>=2: {}'.format(verbosity >= 2))
+    print('v<=1: {}'.format(verbosity <= 1))
+
+    verbosity(1)
+    print('v==1: {}'.format(verbosity == 1))
+    print('v==2: {}'.format(verbosity == 2))
+
+    # should raise ValueError
+    verbosity(5)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/python/visualization.py b/python/visualization.py
index 60d7c01d6..616b27eea 100644
--- a/python/visualization.py
+++ b/python/visualization.py
@@ -846,7 +846,7 @@ def __call__(self,sim,todo):
 
             # Normalize the frames, if requested, and export
             if self.normalize and mp.am_master():
-                if mp.cvar.verbosity > 0:
+                if mp.verbosity > 0:
                     print("Normalizing field data...")
                 fields = np.array(self.cumulative_fields) / np.max(np.abs(self.cumulative_fields),axis=(0,1,2))
                 for k in range(len(self.cumulative_fields)):
@@ -951,7 +951,7 @@ def to_gif(self,fps,filename):
                         '-vf', 'pad=width=ceil(iw/2)*2:height=ceil(ih/2)*2',
                         '-an', filename  # output filename
             ]
-            if mp.cvar.verbosity > 0:
+            if mp.verbosity > 0:
                 print("Generating GIF...")
             proc = Popen(command, stdin=PIPE, stdout=PIPE, stderr=PIPE)
             for i in range(len(self._saved_frames)):
@@ -988,7 +988,7 @@ def to_mp4(self,fps,filename):
                         '-vf', 'pad=width=ceil(iw/2)*2:height=ceil(ih/2)*2',
                         '-an', filename  # output filename
             ]
-            if mp.cvar.verbosity > 0:
+            if mp.verbosity > 0:
                 print("Generating MP4...")
             proc = Popen(command, stdin=PIPE, stdout=PIPE, stderr=PIPE)
             for i in range(len(self._saved_frames)):