diff --git a/python/simulation.py b/python/simulation.py
index 68c13c4d0..0e383244e 100644
--- a/python/simulation.py
+++ b/python/simulation.py
@@ -1008,7 +1008,7 @@ def _compute_fragment_stats(self, gv):
         return stats
 
     def _init_structure(self, k=False):
-        if not mp.cvar.quiet:
+        if mp.cvar.verbosity > 0:
             print('-' * 11)
             print('Initializing structure...')
 
@@ -1031,7 +1031,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 not mp.cvar.quiet:
+        if self._output_stats and isinstance(self.default_material, mp.Medium) and mp.cvar.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))
@@ -1158,7 +1158,7 @@ def init_sim(self):
             self.fields.require_component(mp.Hz)
 
         if self.verbose:
-            self.fields.verbose()
+            verbosity(2)
 
         def use_real(self):
             cond1 = self.is_cylindrical and self.m != 0
@@ -1170,7 +1170,7 @@ def use_real(self):
 
         if use_real(self):
             self.fields.use_real_fields()
-        elif not mp.cvar.quiet:
+        elif mp.cvar.verbosity > 0:
             print("Meep: using complex fields.")
 
         if self.k_point:
@@ -1295,7 +1295,7 @@ def use_output_directory(self, dname=''):
         closure = {'trashed': False}
 
         def hook():
-            if not mp.cvar.quiet:
+            if mp.cvar.verbosity > 0:
                 print("Meep: using output directory '{}'".format(dname))
             self.fields.set_output_directory(dname)
             if not closure['trashed']:
@@ -1347,7 +1347,7 @@ def stop_cond(sim):
         for func in step_funcs:
             _eval_step_func(self, func, 'finish')
 
-        if not mp.cvar.quiet:
+        if mp.cvar.verbosity > 0:
             print("run {} finished at t = {} ({} timesteps)".format(self.run_index, self.meep_time(), self.fields.t))
         self.run_index += 1
 
@@ -2480,7 +2480,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 not mp.cvar.quiet:
+            if closure['max_abs'] != 0 and mp.cvar.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
@@ -2551,7 +2551,7 @@ def display_progress(t0, t, dt):
 
     def _disp(sim):
         t1 = mp.wall_time()
-        if t1 - closure['tlast'] >= dt and not mp.cvar.quiet:
+        if t1 - closure['tlast'] >= dt and mp.cvar.verbosity > 0:
             msg_fmt = "Meep progress: {}/{} = {:.1f}% done in {:.1f}s, {:.1f}s to go"
             val1 = sim.meep_time() - t0
             val2 = val1 / (0.01 * t)
@@ -2924,4 +2924,7 @@ def complexarray(re, im):
 
 
 def quiet(quietval=True):
-    mp.cvar.quiet = quietval
+    mp.cvar.verbose = int(not quietval)
+
+def verbosity(verbose_val):
+    mp.cvar.verbose = verbose_val
diff --git a/python/visualization.py b/python/visualization.py
index 17456c659..994dc4458 100644
--- a/python/visualization.py
+++ b/python/visualization.py
@@ -8,7 +8,7 @@
 import meep as mp
 from meep.geom import Vector3, init_do_averaging
 from meep.source import EigenModeSource, check_positive
-    
+
 
 # ------------------------------------------------------- #
 # Visualization
@@ -44,7 +44,7 @@
         }
 
 default_eps_parameters = {
-        'interpolation':'spline36', 
+        'interpolation':'spline36',
         'cmap':'binary',
         'alpha':1.0
     }
@@ -110,24 +110,24 @@ def place_label(ax,label_text,x,y,centerx,centery,label_parameters=None):
         ytext = -offset
     else:
         ytext = offset
-    
-    ax.annotate(label_text, xy=(x, y), xytext=(xtext,ytext), 
+
+    ax.annotate(label_text, xy=(x, y), xytext=(xtext,ytext),
             textcoords='offset points', ha='center', va='bottom',
             bbox=dict(boxstyle='round,pad=0.2', fc=color, alpha=alpha),
-            arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0.5', 
+            arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=0.5',
                             color=color))
     return ax
 
 # ------------------------------------------------------- #
 # Helper functions used to plot volumes on a 2D plane
 
-# Returns the intersection points of 2 Volumes. 
-# Volumes must be a line, plane, or rectangular prism 
+# Returns the intersection points of 2 Volumes.
+# Volumes must be a line, plane, or rectangular prism
 # (since they are volume objects)
 def intersect_volume_volume(volume1,volume2):
     # volume1 ............... [volume]
     # volume2 ............... [volume]
-    
+
     # Represent the volumes by an "upper" and "lower" coordinate
     U1 = [volume1.center.x+volume1.size.x/2,volume1.center.y+volume1.size.y/2,volume1.center.z+volume1.size.z/2]
     L1 = [volume1.center.x-volume1.size.x/2,volume1.center.y-volume1.size.y/2,volume1.center.z-volume1.size.z/2]
@@ -167,7 +167,7 @@ def intersect_volume_volume(volume1,volume2):
 # volume. They can also do nothing and plot the XY plane through Z=0.
 #
 # Not only do we need to check for all of these possibilities, but we also need
-# to check if the user accidentally specifies a plane that stretches beyond the 
+# to check if the user accidentally specifies a plane that stretches beyond the
 # simulation domain.
 def get_2D_dimensions(sim,output_plane):
     from meep.simulation import Volume
@@ -183,7 +183,7 @@ def get_2D_dimensions(sim,output_plane):
 
     if plane_size.x!=0 and plane_size.y!=0 and plane_size.z!=0:
         raise ValueError("Plane volume must be 2D (a plane).")
-        
+
     check_volume = Volume(center=sim.geometry_center,size=sim.cell_size)
 
     vertices = intersect_volume_volume(check_volume,plane_volume)
@@ -195,7 +195,7 @@ def get_2D_dimensions(sim,output_plane):
 
     if (intersection_vol.size != plane_volume.size) or (intersection_vol.center != plane_volume.center):
         warnings.warn('The specified user volume is larger than the simulation domain and has been truncated.')
-    
+
     sim_center, sim_size = (intersection_vol.center, intersection_vol.size)
     return sim_center, sim_size
 
@@ -215,7 +215,7 @@ def plot_volume(sim,ax,volume,output_plane=None,plotting_parameters=None,label=N
 
     # Get domain measurements
     sim_center, sim_size = get_2D_dimensions(sim,output_plane)
-    
+
     plane = Volume(center=sim_center,size=sim_size)
 
     # Pull volume parameters
@@ -263,7 +263,7 @@ def sort_points(xy):
                 return ax
             else:
                 return ax
-        
+
         # Line volume
         elif len(intersection) == 2:
             line_args = {key:value for key, value in plotting_parameters.items() if key in ['color','linestyle','linewidth','alpha']}
@@ -281,7 +281,7 @@ def sort_points(xy):
                 return ax
             else:
                 return ax
-        
+
         # Planar volume
         elif len(intersection) > 2:
             planar_args = {key:value for key, value in plotting_parameters.items() if key in ['edgecolor','linewidth','facecolor','hatch','alpha']}
@@ -306,11 +306,11 @@ def sort_points(xy):
 def plot_eps(sim,ax,output_plane=None,eps_parameters=None,omega=0):
     if sim.structure is None:
         sim.init_sim()
-    
-    
+
+
     # consolidate plotting parameters
     eps_parameters = default_eps_parameters if eps_parameters is None else dict(default_eps_parameters, **eps_parameters)
-    
+
     # Get domain measurements
     sim_center, sim_size = get_2D_dimensions(sim,output_plane)
 
@@ -341,7 +341,7 @@ def plot_eps(sim,ax,output_plane=None,eps_parameters=None,omega=0):
         ylabel = 'Y'
     else:
         raise ValueError("A 2D plane has not been specified...")
-    
+
     eps_data = np.rot90(np.real(sim.get_array(center=center, size=cell_size, component=mp.Dielectric, omega=omega)))
     if mp.am_master():
         ax.imshow(eps_data, extent=extent, **eps_parameters)
@@ -396,7 +396,7 @@ def get_boundary_volumes(thickness,direction,side):
             size=Vector3(cell_x,cell_y,thickness))
         else:
             raise ValueError("Invalid boundary type")
-    
+
     import itertools
     for boundary in sim.boundary_layers:
         # All 4 side are the same
@@ -442,7 +442,7 @@ def plot_sources(sim,ax,output_plane=None,labels=False,source_parameters=None):
 def plot_monitors(sim,ax,output_plane=None,labels=False,monitor_parameters=None):
     if not sim._is_initialized:
         sim.init_sim()
-    
+
     from meep.simulation import Volume
 
      # consolidate plotting parameters
@@ -459,17 +459,17 @@ def plot_monitors(sim,ax,output_plane=None,labels=False,monitor_parameters=None)
 def plot_fields(sim,ax=None,fields=None,output_plane=None,field_parameters=None):
     if not sim._is_initialized:
         sim.init_sim()
-    
+
     if fields is None:
         return ax
-    
+
     field_parameters = default_field_parameters if field_parameters is None else dict(default_field_parameters, **field_parameters)
 
     # user specifies a field component
     if fields in [mp.Ex, mp.Ey, mp.Ez, mp.Hx, mp.Hy, mp.Hz]:
         # Get domain measurements
         sim_center, sim_size = get_2D_dimensions(sim,output_plane)
-        
+
         xmin = sim_center.x - sim_size.x/2
         xmax = sim_center.x + sim_size.x/2
         ymin = sim_center.y - sim_size.y/2
@@ -498,8 +498,8 @@ def plot_fields(sim,ax=None,fields=None,output_plane=None,field_parameters=None)
         fields = sim.get_array(center=center, size=cell_size, component=fields)
     else:
         raise ValueError('Please specify a valid field component (mp.Ex, mp.Ey, ...')
-    
-    
+
+
     fields = field_parameters['post_process'](fields)
 
     # Either plot the field, or return the array
@@ -539,21 +539,21 @@ def plot2D(sim,ax=None, output_plane=None, fields=None, labels=False,
                 except:
                     # No sources
                     omega = 0
-        
+
     # validate the output plane to ensure proper 2D coordinates
     from meep.simulation import Volume
     sim_center, sim_size = get_2D_dimensions(sim,output_plane)
     output_plane = Volume(center=sim_center,size=sim_size)
-    
+
     # Plot geometry
     ax = plot_eps(sim,ax,output_plane=output_plane,eps_parameters=eps_parameters,omega=omega)
-    
+
     # Plot boundaries
     ax = plot_boundaries(sim,ax,output_plane=output_plane,boundary_parameters=boundary_parameters)
-            
+
     # Plot sources
     ax = plot_sources(sim,ax,output_plane=output_plane,labels=labels,source_parameters=source_parameters)
-        
+
     # Plot monitors
     ax = plot_monitors(sim,ax,output_plane=output_plane,labels=labels,monitor_parameters=monitor_parameters)
 
@@ -567,10 +567,10 @@ def plot3D(sim):
 
     if not sim._is_initialized:
         sim.init_sim()
-        
+
     if sim.dimensions < 3:
         raise ValueError("Simulation must have 3 dimensions to visualize 3D")
-    
+
     eps_data = sim.get_epsilon()
     s = mlab.contour3d(eps_data, colormap="YlGnBu")
     return s
@@ -683,7 +683,7 @@ def plot_box(box, proc, fig, ax):
 #                       after simulation ends.
 # plot_modifiers ...... [list] additional functions to
 #                       modify plot
-# customization_args .. [dict] other customization args 
+# customization_args .. [dict] other customization args
 #                       to pass to plot2D()
 
 class Animate2D(object):
@@ -716,11 +716,11 @@ def __init__(self,sim,fields,f=None,realtime=False,normalize=False,
         self.default_mode = 'loop' # html5 video control mode
 
         self.init = False
-        
+
         # Needed for step functions
         self.__code__ = namedtuple('gna_hack',['co_argcount'])
         self.__code__.co_argcount=2
-    
+
     def __call__(self,sim,todo):
         from matplotlib import pyplot as plt
 
@@ -739,20 +739,20 @@ def __call__(self,sim,todo):
                     self.ax = ax
                     self.w, self.h = self.f.get_size_inches()
                 self.init = True
-            else:               
+            else:
                 # Update the plot
                 filtered_plot_fields= filter_dict(self.customization_args, plot_fields)
                 fields = sim.plot_fields(fields=self.fields,**filtered_plot_fields)
                 if mp.am_master():
                     self.ax.images[-1].set_data(fields)
                     self.ax.images[-1].set_clim(vmin=0.8*np.min(fields), vmax=0.8*np.max(fields))
-            
+
             if self.realtime and mp.am_master():
                 # Redraw the current figure if requested
                 plt.pause(0.05)
 
             if self.normalize and mp.am_master():
-                # Save fields as a numpy array to be normalized 
+                # Save fields as a numpy array to be normalized
                 # and saved later.
                 self.cumulative_fields.append(fields)
             elif mp.am_master():
@@ -761,26 +761,26 @@ def __call__(self,sim,todo):
                 self.grab_frame()
             return
         elif todo == 'finish':
-            
+
             # Normalize the frames, if requested, and export
             if self.normalize and mp.am_master():
-                if not mp.cvar.quiet:
+                if mp.cvar.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)):
                     self.ax.images[-1].set_data(fields[k,:,:])
                     self.ax.images[-1].set_clim(vmin=-0.8, vmax=0.8)
                     self.grab_frame()
-                
+
             return
-    
+
     @property
     def frame_size(self):
         # A tuple ``(width, height)`` in pixels of a movie frame.
         # modified from matplotlib library
         w, h = self.f.get_size_inches()
         return int(w * self.f.dpi), int(h * self.f.dpi)
-    
+
     def grab_frame(self):
         # Saves the figures frame to memory.
         # modified from matplotlib library
@@ -804,7 +804,7 @@ def to_jshtml(self,fps):
         # Exports a javascript enabled html object that is
         # ready for jupyter notebook embedding.
         # modified from matplotlib/animation.py code.
-        
+
         # Only works with Python3 and matplotlib > 3.1.0
         from distutils.version import LooseVersion
         import matplotlib
@@ -826,7 +826,7 @@ def to_jshtml(self,fps):
             mode_dict[self.default_mode + '_checked'] = 'checked'
 
             interval = 1000 // fps
-            
+
             html_string = ""
             html_string += JS_INCLUDE
             html_string += STYLE_INCLUDE
@@ -845,7 +845,7 @@ def to_gif(self,fps,filename):
             from subprocess import Popen, PIPE
             from io import TextIOWrapper, BytesIO
             FFMPEG_BIN = 'ffmpeg'
-            command = [FFMPEG_BIN, 
+            command = [FFMPEG_BIN,
                         '-f', 'image2pipe', # force piping of rawvideo
                         '-vcodec', self.frame_format, # raw input codec
                         '-s', '%dx%d' % (self.frame_size),
@@ -853,14 +853,14 @@ def to_gif(self,fps,filename):
                         '-i', 'pipe:', # The input comes from a pipe
                         '-vcodec', 'gif', # output gif format
                         '-r', str(fps), # frame rate in frames per second
-                        '-y', 
+                        '-y',
                         '-vf', 'pad=width=ceil(iw/2)*2:height=ceil(ih/2)*2',
                         '-an', filename  # output filename
             ]
-            if not mp.cvar.quiet:
+            if mp.cvar.verbosity > 0:
                 print("Generating GIF...")
             proc = Popen(command, stdin=PIPE, stdout=PIPE, stderr=PIPE)
-            for i in range(len(self._saved_frames)): 
+            for i in range(len(self._saved_frames)):
                 proc.stdin.write(self._saved_frames[i])
             out, err = proc.communicate() # pipe in images
             proc.stdin.close()
@@ -875,7 +875,7 @@ def to_mp4(self,fps,filename):
             from subprocess import Popen, PIPE
             from io import TextIOWrapper, BytesIO
             FFMPEG_BIN = 'ffmpeg'
-            command = [FFMPEG_BIN, 
+            command = [FFMPEG_BIN,
                         '-f', 'image2pipe', # force piping of rawvideo
                         '-vcodec', self.frame_format, # raw input codec
                         '-s', '%dx%d' % (self.frame_size),
@@ -885,20 +885,20 @@ def to_mp4(self,fps,filename):
                         '-vcodec', self.codec, # output mp4 format
                         '-pix_fmt','yuv420p',
                         '-r', str(fps), # frame rate in frames per second
-                        '-y', 
+                        '-y',
                         '-vf', 'pad=width=ceil(iw/2)*2:height=ceil(ih/2)*2',
                         '-an', filename  # output filename
             ]
-            if not mp.cvar.quiet:
+            if mp.cvar.verbosity > 0:
                 print("Generating MP4...")
             proc = Popen(command, stdin=PIPE, stdout=PIPE, stderr=PIPE)
-            for i in range(len(self._saved_frames)): 
+            for i in range(len(self._saved_frames)):
                 proc.stdin.write(self._saved_frames[i])
             out, err = proc.communicate() # pipe in images
             proc.stdin.close()
             proc.wait()
         return
-    
+
     def reset(self):
         self.cumulative_fields = []
         self.ax = None
@@ -910,7 +910,7 @@ def set_figure(self,f):
 # ------------------------------------------------------- #
 # JS_Animation
 # ------------------------------------------------------- #
-# A helper class used to make jshtml animations embed 
+# A helper class used to make jshtml animations embed
 # seamlessly within Jupyter notebooks.
 
 class JS_Animation():
diff --git a/scheme/structure.cpp b/scheme/structure.cpp
index 786a1d96d..62c36bf2b 100644
--- a/scheme/structure.cpp
+++ b/scheme/structure.cpp
@@ -278,7 +278,7 @@ geom_epsilon::geom_epsilon(geometric_object_list g, material_type_list mlist,
   geom_fix_object_list(geometry);
   geom_box box = gv2box(v);
   geometry_tree = create_geom_box_tree0(geometry, box);
-  if (!meep::quiet && verbose && meep::am_master()) {
+  if (meep::verbosity > 0 && verbose && meep::am_master()) {
     printf("Geometric-object bounding-box tree:\n");
     display_geom_box_tree(5, geometry_tree);
 
@@ -1306,7 +1306,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
         lorentzian_susceptibility *d = p->user_s.subclass.lorentzian_susceptibility_data;
         if (d->which_subclass == lorentzian_susceptibility::NOISY_LORENTZIAN_SUSCEPTIBILITY) {
           noisy_lorentzian_susceptibility *nd = d->subclass.noisy_lorentzian_susceptibility_data;
-          if (!meep::quiet)
+          if (meep::verbosity > 0)
             master_printf("noisy lorentzian susceptibility: frequency=%g, gamma=%g, amp = %g\n",
                           d->frequency, d->gamma, nd->noise_amp);
           sus = new meep::noisy_lorentzian_susceptibility(nd->noise_amp, d->frequency, d->gamma);
@@ -1315,7 +1315,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
                  lorentzian_susceptibility::GYROTROPIC_LORENTZIAN_SUSCEPTIBILITY) {
           gyrotropic_lorentzian_susceptibility *gd =
               d->subclass.gyrotropic_lorentzian_susceptibility_data;
-          if (!meep::quiet)
+          if (meep::verbosity > 0)
             master_printf(
                 "gyrotropic lorentzian susceptibility: bias=(%g,%g,%g), frequency=%g, gamma=%g\n",
                 gd->bias.x, gd->bias.y, gd->bias.z, d->frequency, d->gamma);
@@ -1323,7 +1323,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
                                                     d->gamma, 0.0, meep::GYROTROPIC_LORENTZIAN);
         }
         else { // just a Lorentzian
-          if (!meep::quiet)
+          if (meep::verbosity > 0)
             master_printf("lorentzian susceptibility: frequency=%g, gamma=%g\n", d->frequency,
                           d->gamma);
           sus = new meep::lorentzian_susceptibility(d->frequency, d->gamma);
@@ -1334,7 +1334,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
         drude_susceptibility *d = p->user_s.subclass.drude_susceptibility_data;
         if (d->which_subclass == drude_susceptibility::NOISY_DRUDE_SUSCEPTIBILITY) {
           noisy_drude_susceptibility *nd = d->subclass.noisy_drude_susceptibility_data;
-          if (!meep::quiet)
+          if (meep::verbosity > 0)
             master_printf("noisy drude susceptibility: frequency=%g, gamma=%g, amp = %g\n",
                           d->frequency, d->gamma, nd->noise_amp);
           sus = new meep::noisy_lorentzian_susceptibility(nd->noise_amp, d->frequency, d->gamma,
@@ -1342,7 +1342,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
         }
         else if (d->which_subclass == drude_susceptibility::GYROTROPIC_DRUDE_SUSCEPTIBILITY) {
           gyrotropic_drude_susceptibility *gd = d->subclass.gyrotropic_drude_susceptibility_data;
-          if (!meep::quiet)
+          if (meep::verbosity > 0)
             master_printf(
                 "gyrotropic drude susceptibility: bias=(%g,%g,%g), frequency=%g, gamma=%g\n",
                 gd->bias.x, gd->bias.y, gd->bias.z, d->frequency, d->gamma);
@@ -1350,7 +1350,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
                                                     d->gamma, 0.0, meep::GYROTROPIC_DRUDE);
         }
         else { // just a Drude
-          if (!meep::quiet)
+          if (meep::verbosity > 0)
             master_printf("drude susceptibility: frequency=%g, gamma=%g\n", d->frequency, d->gamma);
           sus = new meep::lorentzian_susceptibility(d->frequency, d->gamma, true);
         }
@@ -1359,7 +1359,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
       case susceptibility::GYROTROPIC_SATURATED_SUSCEPTIBILITY: {
         gyrotropic_saturated_susceptibility *d =
             p->user_s.subclass.gyrotropic_saturated_susceptibility_data;
-        if (!meep::quiet)
+        if (meep::verbosity > 0)
           master_printf("gyrotropic Landau-Lifshitz-Gilbert-type susceptibility: bias=(%g,%g,%g), "
                         "frequency=%g, gamma=%g, alpha=%g\n",
                         d->bias.x, d->bias.y, d->bias.z, d->frequency, d->gamma, d->alpha);
@@ -1419,7 +1419,7 @@ meep::structure *make_structure(int dims, vector3 size, vector3 center, double r
                                 const char *eps_input_file, pml_list pml_layers,
                                 symmetry_list symmetries, int num_chunks, double Courant,
                                 double global_D_conductivity_, double global_B_conductivity_) {
-  if (!meep::quiet) master_printf("-----------\nInitializing structure...\n");
+  if (meep::verbosity > 0) master_printf("-----------\nInitializing structure...\n");
 
   // only cartesian lattices are currently allowed
   geom_initialize();
@@ -1438,7 +1438,7 @@ meep::structure *make_structure(int dims, vector3 size, vector3 center, double r
   geometry_lattice.size = size;
   geometry_edge = vector3_to_vec(size) * 0.5;
 
-  if (!meep::quiet) {
+  if (meep::verbosity > 0) {
     master_printf("Working in %s dimensions.\n", meep::dimension_name(dim));
     master_printf("Computational cell is %g x %g x %g with resolution %g\n", size.x, size.y, size.z,
                   resolution);
@@ -1592,7 +1592,7 @@ meep::structure *make_structure(int dims, vector3 size, vector3 center, double r
 
   geps.add_susceptibilities(s);
 
-  if (!meep::quiet) master_printf("-----------\n");
+  if (meep::verbosity > 0) master_printf("-----------\n");
 
   return s;
 }
diff --git a/src/anisotropic_averaging.cpp b/src/anisotropic_averaging.cpp
index 050975ea1..4ccb9ef14 100644
--- a/src/anisotropic_averaging.cpp
+++ b/src/anisotropic_averaging.cpp
@@ -250,7 +250,7 @@ void structure_chunk::set_chi1inv(component c, material_function &medium,
       trivial[2] = trivial[2] && (chi1inv[c][d2][i] == trivial_val[2]);
     }
 
-    if (!quiet && (ipixel + 1) % 1000 == 0 &&
+    if (verbosity > 0 && (ipixel + 1) % 1000 == 0 &&
         wall_time() > last_output_time + MEEP_MIN_OUTPUT_TIME) {
       master_printf("subpixel-averaging is %g%% done, %g s remaining\n", ipixel * 100.0 / npixels,
                     (npixels - ipixel) * (wall_time() - last_output_time) / ipixel);
diff --git a/src/bicgstab.cpp b/src/bicgstab.cpp
index de88353af..ed22806d2 100644
--- a/src/bicgstab.cpp
+++ b/src/bicgstab.cpp
@@ -158,7 +158,7 @@ ptrdiff_t bicgstabL(const int L, const size_t n, realnum *x, bicgstab_op A, void
   double resid;
   while ((resid = norm2(n, r[0])) > tol * bnrm) {
     ++iter;
-    if (!quiet && wall_time() > last_output_wall_time + MEEP_MIN_OUTPUT_TIME) {
+    if (verbosity > 0 && wall_time() > last_output_wall_time + MEEP_MIN_OUTPUT_TIME) {
       master_printf("residual[%d] = %g\n", iter, resid / bnrm);
       last_output_wall_time = wall_time();
     }
@@ -220,7 +220,7 @@ ptrdiff_t bicgstabL(const int L, const size_t n, realnum *x, bicgstab_op A, void
     }
   }
 
-  if (!quiet) master_printf("final residual = %g\n", norm2(n, r[0]) / bnrm);
+  if (verbosity > 0) master_printf("final residual = %g\n", norm2(n, r[0]) / bnrm);
 
 finish:
   delete[] sigma;
diff --git a/src/cw_fields.cpp b/src/cw_fields.cpp
index 6ce78638c..305815404 100644
--- a/src/cw_fields.cpp
+++ b/src/cw_fields.cpp
@@ -172,9 +172,9 @@ bool fields::solve_cw(double tol, int maxiters, complex<double> frequency, int L
   data.iters = 0;
 
   int ierr = (int)bicgstabL(L, N, reinterpret_cast<realnum *>(x), fieldop, &data,
-                            reinterpret_cast<realnum *>(b), tol, &maxiters, work, quiet);
+                            reinterpret_cast<realnum *>(b), tol, &maxiters, work, verbosity == 0);
 
-  if (!quiet) {
+  if (verbosity > 0) {
     master_printf("Finished solve_cw after %d steps and %d CG iters.\n", data.iters, maxiters);
     if (ierr) master_printf(" -- CONVERGENCE FAILURE (%d) in solve_cw!\n", ierr);
   }
diff --git a/src/fields.cpp b/src/fields.cpp
index 8193c6f3b..52d85652d 100644
--- a/src/fields.cpp
+++ b/src/fields.cpp
@@ -30,7 +30,6 @@ namespace meep {
 fields::fields(structure *s, double m, double beta, bool zero_fields_near_cylorigin)
     : S(s->S), gv(s->gv), user_volume(s->user_volume), v(s->v), m(m), beta(beta) {
   shared_chunks = s->shared_chunks;
-  verbosity = 0;
   components_allocated = false;
   synchronized_magnetic_fields = 0;
   outdir = new char[strlen(s->outdir) + 1];
@@ -92,7 +91,6 @@ fields::fields(structure *s, double m, double beta, bool zero_fields_near_cylori
 fields::fields(const fields &thef)
     : S(thef.S), gv(thef.gv), user_volume(thef.user_volume), v(thef.v) {
   shared_chunks = thef.shared_chunks;
-  verbosity = 0;
   components_allocated = thef.components_allocated;
   synchronized_magnetic_fields = thef.synchronized_magnetic_fields;
   outdir = new char[strlen(thef.outdir) + 1];
@@ -156,12 +154,6 @@ fields::~fields() {
   delete[] outdir;
 }
 
-void fields::verbose(int gv) {
-  verbosity = gv;
-  for (int i = 0; i < num_chunks; i++)
-    chunks[i]->verbose(gv);
-}
-
 void fields::use_real_fields() {
   LOOP_OVER_DIRECTIONS(gv.dim, d) {
     if (boundaries[High][d] == Periodic && k[d] != 0.0)
@@ -231,7 +223,6 @@ fields_chunk::fields_chunk(structure_chunk *the_s, const char *od, double m, dou
       beta(beta) {
   s = the_s;
   s->refcount++;
-  verbosity = 0;
   outdir = od;
   new_s = NULL;
   is_real = 0;
@@ -289,7 +280,6 @@ fields_chunk::fields_chunk(structure_chunk *the_s, const char *od, double m, dou
 fields_chunk::fields_chunk(const fields_chunk &thef) : gv(thef.gv), v(thef.v) {
   s = thef.s;
   s->refcount++;
-  verbosity = thef.verbosity;
   outdir = thef.outdir;
   m = thef.m;
   zero_fields_near_cylorigin = thef.zero_fields_near_cylorigin;
diff --git a/src/h5fields.cpp b/src/h5fields.cpp
index e473495b6..0638b7767 100644
--- a/src/h5fields.cpp
+++ b/src/h5fields.cpp
@@ -451,7 +451,7 @@ const char *fields::h5file_name(const char *name, const char *prefix, bool times
 h5file *fields::open_h5file(const char *name, h5file::access_mode mode, const char *prefix,
                             bool timestamp) {
   const char *filename = h5file_name(name, prefix, timestamp);
-  if (!quiet && mode == h5file::WRITE) master_printf("creating output file \"%s\"...\n", filename);
+  if (verbosity > 0 && mode == h5file::WRITE) master_printf("creating output file \"%s\"...\n", filename);
   return new h5file(filename, mode, true);
 }
 
diff --git a/src/meep.hpp b/src/meep.hpp
index 9fa946848..2c2ec30de 100644
--- a/src/meep.hpp
+++ b/src/meep.hpp
@@ -44,7 +44,7 @@ typedef double realnum;
 
 #define MEEP_MIN_OUTPUT_TIME 4.0 // output no more often than this many seconds
 
-extern bool quiet; // if true, suppress all non-error messages from Meep
+extern int verbosity; // if 0, suppress all non-error messages from Meep; 1 is default, 2 is debug output
 
 const double pi = 3.141592653589793238462643383276;
 
@@ -1349,7 +1349,6 @@ class fields_chunk {
   void restore_component(component c);
 
   void set_output_directory(const char *name);
-  void verbose(int gv = 1) { verbosity = gv; }
 
   double count_volume(component);
   friend class fields;
@@ -1386,7 +1385,6 @@ class fields_chunk {
   bool doing_solve_cw;                 // true when inside solve_cw
   std::complex<double> solve_cw_omega; // current omega for solve_cw
 
-  int verbosity; // Turn on verbosity for debugging purposes...
   // fields.cpp
   bool have_plus_deriv[NUM_FIELD_COMPONENTS], have_minus_deriv[NUM_FIELD_COMPONENTS];
   component plus_component[NUM_FIELD_COMPONENTS], minus_component[NUM_FIELD_COMPONENTS];
@@ -1832,7 +1830,6 @@ class fields {
   double electric_energy_weighted_integral(double (*f)(const vec &), const volume &where);
 
   void set_output_directory(const char *name);
-  void verbose(int gv = 1);
   double count_volume(component);
   // fields.cpp
   bool have_component(component);
@@ -1853,7 +1850,6 @@ class fields {
   void unset_solve_cw_omega();
 
 private:
-  int verbosity;                    // Turn on verbosity for debugging purposes...
   int synchronized_magnetic_fields; // count number of nested synchs
   double last_wall_time;
 #define MEEP_TIMING_STACK_SZ 10
diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index 65770116c..1a90fde7a 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -381,12 +381,12 @@ geom_epsilon::geom_epsilon(geometric_object_list g, material_type_list mlist,
   if (meep::am_master()) {
     for (int i = 0; i < geometry.num_items; ++i) {
 
-      if (!meep::quiet) display_geometric_object_info(5, geometry.items[i]);
+      if (meep::verbosity > 0) display_geometric_object_info(5, geometry.items[i]);
 
       medium_struct *mm;
       if (is_medium(geometry.items[i].material, &mm)) {
         check_offdiag(mm);
-        if (!meep::quiet)
+        if (meep::verbosity > 0)
           master_printf("%*sdielectric constant epsilon diagonal "
                         "= (%g,%g,%g)\n",
                         5 + 5, "", mm->epsilon_diag.x, mm->epsilon_diag.y, mm->epsilon_diag.z);
@@ -397,7 +397,7 @@ geom_epsilon::geom_epsilon(geometric_object_list g, material_type_list mlist,
   geom_fix_object_list(geometry);
   geom_box box = gv2box(v);
   geometry_tree = create_geom_box_tree0(geometry, box);
-  if (!meep::quiet && verbose && meep::am_master()) {
+  if (meep::verbosity > 0 && verbose && meep::am_master()) {
     master_printf("Geometric-object bounding-box tree:\n");
     display_geom_box_tree(5, geometry_tree);
 
@@ -1383,7 +1383,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
     if (ss->transitions.size() != 0 || ss->initial_populations.size() != 0) {
       // multilevel atom
       sus = make_multilevel_sus(ss);
-      if (!meep::quiet) master_printf("multilevel atom susceptibility\n");
+      if (meep::verbosity > 0) master_printf("multilevel atom susceptibility\n");
     }
     else {
       if (noisy) {
@@ -1401,7 +1401,7 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft, meep::structure *s)
       else {
         sus = new meep::lorentzian_susceptibility(ss->frequency, ss->gamma, ss->drude);
       }
-      if (!meep::quiet) {
+      if (meep::verbosity > 0) {
         master_printf("%s%s susceptibility: frequency=%g, gamma=%g",
                       noisy ? "noisy " : gyrotropic ? "gyrotropic " : "",
                       ss->saturated_gyrotropy ? "Landau-Lifshitz-Gilbert-type"
@@ -1523,7 +1523,7 @@ void set_materials_from_geometry(meep::structure *s, geometric_object_list g, ve
   geometry_lattice.size = size;
   geometry_edge = vector3_to_vec(size) * 0.5;
 
-  if (!meep::quiet) {
+  if (meep::verbosity > 0) {
     master_printf("Working in %s dimensions.\n", meep::dimension_name(s->gv.dim));
     master_printf("Computational cell is %g x %g x %g with resolution %g\n", size.x, size.y, size.z,
                   resolution);
@@ -1557,7 +1557,7 @@ void set_materials_from_geometry(meep::structure *s, geometric_object_list g, ve
   s->remove_susceptibilities();
   geps.add_susceptibilities(s);
 
-  if (!meep::quiet) master_printf("-----------\n");
+  if (meep::verbosity > 0) master_printf("-----------\n");
 }
 
 /***************************************************************/
@@ -1596,7 +1596,7 @@ material_type make_file_material(const char *eps_input_file) {
     meep::h5file eps_file(fname, meep::h5file::READONLY, false);
     int rank; // ignored since rank < 3 is equivalent to singleton dims
     md->epsilon_data = eps_file.read(dataname, &rank, md->epsilon_dims, 3);
-    if (!meep::quiet)
+    if (meep::verbosity > 0)
       master_printf("read in %zdx%zdx%zd epsilon-input-file \"%s\"\n", md->epsilon_dims[0],
                     md->epsilon_dims[1], md->epsilon_dims[2], eps_input_file);
     delete[] fname;
diff --git a/src/mpb.cpp b/src/mpb.cpp
index 7544a0d4d..98caeb8f4 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -327,7 +327,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c
     default: abort("unsupported dimensionality in add_eigenmode_source");
   }
 
-  if (!quiet && verbose) master_printf("KPOINT: %g, %g, %g\n", k[0], k[1], k[2]);
+  if (verbosity > 0 && verbose) master_printf("KPOINT: %g, %g, %g\n", k[0], k[1], k[2]);
 
   double kcart_len = sqrt(dot_product(kcart, kcart));
 
@@ -406,7 +406,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c
           fabs(k[d - X]) > 0.4) // ensure k is well inside the Brillouin zone
         k[d - X] = k[d - X] > 0 ? 0.4 : -0.4;
     }
-    if (!quiet && verbose) master_printf("NEW KPOINT: %g, %g, %g\n", k[0], k[1], k[2]);
+    if (verbosity > 0 && verbose) master_printf("NEW KPOINT: %g, %g, %g\n", k[0], k[1], k[2]);
   }
 
   set_maxwell_data_parity(mdata, parity);
@@ -456,8 +456,8 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c
 #endif
                   maxwell_preconditioner2, (void *)mdata, evectconstraint_chain_func,
                   (void *)constraints, W, 3, eigensolver_tol, &num_iters,
-                  EIGS_DEFAULT_FLAGS | (am_master() && verbose && !quiet ? EIGS_VERBOSE : 0));
-      if (!quiet)
+                  EIGS_DEFAULT_FLAGS | (am_master() && verbose && verbosity > 0 ? EIGS_VERBOSE : 0));
+      if (verbosity > 0)
         master_printf("MPB solved for omega_%d(%g,%g,%g) = %g after %d iters\n", band_num,
                       G[0][0] * k[0], G[1][1] * k[1], G[2][2] * k[2], sqrt(eigvals[band_num - 1]),
                       num_iters);
@@ -482,7 +482,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c
         // update k via Newton step
         double dkmatch = (sqrt(eigvals[band_num - 1]) - omega_src) / vgrp;
         kmatch = kmatch - dkmatch;
-        if (!quiet && verbose)
+        if (verbosity > 0 && verbose)
           master_printf("Newton step: group velocity v=%g, kmatch=%g\n", vgrp, kmatch);
         count_dkmatch_increase += fabs(dkmatch) > fabs(dkmatch_prev);
         if (count_dkmatch_increase > 4) {
@@ -609,7 +609,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c
   if (kdom) {
 #if MPB_VERSION_MAJOR > 1 || (MPB_VERSION_MAJOR == 1 && MPB_VERSION_MINOR >= 7)
     maxwell_dominant_planewave(mdata, H, band_num, kdom);
-    if (!quiet)
+    if (verbosity > 0)
       master_printf("Dominant planewave for band %d: (%f,%f,%f)\n", band_num, kdom[0], kdom[1],
                     kdom[2]);
 #else
diff --git a/src/mympi.cpp b/src/mympi.cpp
index f4a53862d..6383d7b64 100644
--- a/src/mympi.cpp
+++ b/src/mympi.cpp
@@ -75,7 +75,7 @@ namespace meep {
 static MPI_Comm mycomm = MPI_COMM_WORLD;
 #endif
 
-bool quiet = false; // defined in meep.h
+int verbosity = 1; // defined in meep.h
 void (*master_printf_callback)(const char *s) = NULL;
 
 initialize::initialize(int &argc, char **&argv) {
@@ -83,7 +83,7 @@ initialize::initialize(int &argc, char **&argv) {
   MPI_Init(&argc, &argv);
   int major, minor;
   MPI_Get_version(&major, &minor);
-  if (!quiet)
+  if (verbosity > 0)
     master_printf("Using MPI version %d.%d, %d processes\n", major, minor, count_processors());
 #else
   UNUSED(argc);
@@ -99,7 +99,7 @@ initialize::initialize(int &argc, char **&argv) {
 }
 
 initialize::~initialize() {
-  if (!quiet) master_printf("\nElapsed run time = %g s\n", elapsed_time());
+  if (verbosity > 0) master_printf("\nElapsed run time = %g s\n", elapsed_time());
 #ifdef HAVE_MPI
   end_divide_parallel();
   MPI_Finalize();
diff --git a/src/near2far.cpp b/src/near2far.cpp
index e9b846179..fc4063692 100644
--- a/src/near2far.cpp
+++ b/src/near2far.cpp
@@ -338,7 +338,7 @@ realnum *dft_near2far::get_farfields_array(const volume &where, int &rank, size_
       for (size_t i2 = 0; i2 < dims[2]; ++i2) {
         x.set_direction(dirs[2], where.in_direction_min(dirs[2]) + i2 * dx[2]);
         double t;
-        if (!quiet && (t = wall_time()) > start + MEEP_MIN_OUTPUT_TIME) {
+        if (verbosity > 0 && (t = wall_time()) > start + MEEP_MIN_OUTPUT_TIME) {
           size_t this_point = (dims[1] * dims[2] * i0) + (dims[2] * i1) + i2 + 1;
           master_printf("get_farfields_array working on point %zu of %zu (%d%% done), %g s/point\n",
                         this_point, total_points, (int)((double)this_point / total_points * 100),
@@ -347,7 +347,7 @@ realnum *dft_near2far::get_farfields_array(const volume &where, int &rank, size_
           last_point = this_point;
         }
         farfield_lowlevel(EH1, x);
-        if (!quiet) all_wait(); // Allow consistent progress updates from master
+        if (verbosity > 0) all_wait(); // Allow consistent progress updates from master
         ptrdiff_t idx = (i0 * dims[1] + i1) * dims[2] + i2;
         for (int i = 0; i < Nfreq; ++i)
           for (int k = 0; k < 6; ++k) {
diff --git a/src/output_directory.cpp b/src/output_directory.cpp
index 2bcc67d92..f1f89a82d 100644
--- a/src/output_directory.cpp
+++ b/src/output_directory.cpp
@@ -35,7 +35,7 @@ const char symlink_name[] = "latest_output";
 void structure::set_output_directory(const char *name) {
   char buf[300];
   outdir = name;
-  if (!quiet) master_printf("Using output directory %s/\n", name);
+  if (verbosity > 0) master_printf("Using output directory %s/\n", name);
   if (readlink(symlink_name, buf, 300) > 0) {
     // Link already exists.
     unlink(symlink_name);
@@ -126,7 +126,7 @@ const char *make_output_directory(const char *exename, const char *jobname) {
   {
     int i = 0;
     while (!is_ok_dir(outdirname)) {
-      if (!quiet) master_printf("Output directory %s already exists!\n", outdirname);
+      if (verbosity > 0) master_printf("Output directory %s already exists!\n", outdirname);
       snprintf(outdirname, buflen, "%s-out-%d", basename, i++);
     }
   }
diff --git a/src/sources.cpp b/src/sources.cpp
index 121569dfb..0de3123ef 100644
--- a/src/sources.cpp
+++ b/src/sources.cpp
@@ -378,13 +378,13 @@ void fields::add_volume_source(component c, const src_time &src, const volume &w
 
   size_t re_dims[] = {1, 1, 1};
   double *real_data = eps_file.read(dataset_re.c_str(), &rank, re_dims, 3);
-  if (!quiet)
+  if (verbosity > 0)
     master_printf("read in %zdx%zdx%zd amplitude function file \"%s:%s\"\n", re_dims[0], re_dims[1],
                   re_dims[2], filename, dataset_re.c_str());
 
   size_t im_dims[] = {1, 1, 1};
   double *imag_data = eps_file.read(dataset_im.c_str(), &rank, im_dims, 3);
-  if (!quiet)
+  if (verbosity > 0)
     master_printf("read in %zdx%zdx%zd amplitude function file \"%s:%s\"\n", im_dims[0], im_dims[1],
                   im_dims[2], filename, dataset_im.c_str());
 
diff --git a/src/step.cpp b/src/step.cpp
index 3d283990b..ef66ad7c2 100644
--- a/src/step.cpp
+++ b/src/step.cpp
@@ -44,7 +44,7 @@ void fields::step() {
     last_step_output_wall_time = wall_time();
     last_step_output_t = t;
   }
-  if (!quiet && wall_time() > last_step_output_wall_time + MEEP_MIN_OUTPUT_TIME) {
+  if (verbosity > 0 && wall_time() > last_step_output_wall_time + MEEP_MIN_OUTPUT_TIME) {
     master_printf("on time step %d (time=%g), %g s/step\n", t, time(),
                   (wall_time() - last_step_output_wall_time) / (t - last_step_output_t));
     if (save_synchronized_magnetic_fields)
diff --git a/src/structure.cpp b/src/structure.cpp
index 468b4f3f7..e39037e33 100644
--- a/src/structure.cpp
+++ b/src/structure.cpp
@@ -54,7 +54,7 @@ structure::structure(const grid_volume &thegv, material_function &eps, const bou
   if (!br.check_ok(thegv)) abort("invalid boundary absorbers for this grid_volume");
   double tstart = wall_time();
   choose_chunkdivision(thegv, num, br, s);
-  if (!quiet) master_printf("time for choose_chunkdivision = %g s\n", wall_time() - tstart);
+  if (verbosity > 0) master_printf("time for choose_chunkdivision = %g s\n", wall_time() - tstart);
   set_materials(eps, use_anisotropic_averaging, tol, maxeval);
 }
 
@@ -68,7 +68,7 @@ structure::structure(const grid_volume &thegv, double eps(const vec &), const bo
   if (!br.check_ok(thegv)) abort("invalid boundary absorbers for this grid_volume");
   double tstart = wall_time();
   choose_chunkdivision(thegv, num, br, s);
-  if (!quiet) master_printf("time for choose_chunkdivision = %g s\n", wall_time() - tstart);
+  if (verbosity > 0) master_printf("time for choose_chunkdivision = %g s\n", wall_time() - tstart);
   if (eps) {
     simple_material_function epsilon(eps);
     set_materials(epsilon, use_anisotropic_averaging, tol, maxeval);
@@ -140,7 +140,7 @@ void structure::choose_chunkdivision(const grid_volume &thegv, int desired_num_c
       for (int n = 0; n < S.multiplicity(); n++)
         if (has_direction(thegv.dim, d) &&
             (S.transform(d, n).d != d || S.transform(d, n).flipped)) {
-          if (thegv.num_direction(d) & 1 && !break_this[d] && !quiet)
+          if (thegv.num_direction(d) & 1 && !break_this[d] && verbosity > 0)
             master_printf("Padding %s to even number of grid points.\n", direction_name(d));
           break_this[dd] = true;
         }
@@ -150,7 +150,7 @@ void structure::choose_chunkdivision(const grid_volume &thegv, int desired_num_c
       if (break_mult == S.multiplicity()) break_this[d] = false;
       if (break_this[d]) {
         break_mult *= 2;
-        if (!quiet)
+        if (verbosity > 0)
           master_printf("Halving computational cell along direction %s\n",
                         direction_name(direction(d)));
         gv = gv.halve((direction)d);
@@ -189,7 +189,7 @@ void structure::choose_chunkdivision(const grid_volume &thegv, int desired_num_c
   if (meep_geom::fragment_stats::resolution == 0 ||
       meep_geom::fragment_stats::has_non_medium_material() ||
       meep_geom::fragment_stats::split_chunks_evenly) {
-    if (!quiet && adjusted_num_chunks > 1)
+    if (verbosity > 0 && adjusted_num_chunks > 1)
       master_printf("Splitting into %d chunks evenly\n", adjusted_num_chunks);
     for (int i = 0; i < adjusted_num_chunks; i++) {
       grid_volume vi =
@@ -198,7 +198,7 @@ void structure::choose_chunkdivision(const grid_volume &thegv, int desired_num_c
     }
   }
   else {
-    if (!quiet && adjusted_num_chunks > 1)
+    if (verbosity > 0 && adjusted_num_chunks > 1)
       master_printf("Splitting into %d chunks by cost\n", adjusted_num_chunks);
     split_by_cost(prime_factors, gv, chunk_volumes);
     by_cost = true;
@@ -217,7 +217,7 @@ void structure::choose_chunkdivision(const grid_volume &thegv, int desired_num_c
   }
   check_chunks();
 
-  if (!quiet && by_cost) {
+  if (verbosity > 1 && by_cost) {
     double *costs = new double[count_processors()];
     for (int i = 0; i < count_processors(); i++)
       costs[i] = 0;
@@ -460,7 +460,7 @@ void structure::set_epsilon(material_function &eps, bool use_anisotropic_averagi
                             int maxeval) {
   double tstart = wall_time();
   FOR_ELECTRIC_COMPONENTS(c) { set_chi1inv(c, eps, use_anisotropic_averaging, tol, maxeval); }
-  if (!quiet) master_printf("time for set_epsilon = %g s\n", wall_time() - tstart);
+  if (verbosity > 0) master_printf("time for set_epsilon = %g s\n", wall_time() - tstart);
 }
 
 void structure::set_epsilon(double eps(const vec &), bool use_anisotropic_averaging, double tol,
@@ -473,7 +473,7 @@ void structure::set_mu(material_function &m, bool use_anisotropic_averaging, dou
                        int maxeval) {
   double tstart = wall_time();
   FOR_MAGNETIC_COMPONENTS(c) { set_chi1inv(c, m, use_anisotropic_averaging, tol, maxeval); }
-  if (!quiet) master_printf("time for set_mu = %g s\n", wall_time() - tstart);
+  if (verbosity > 0) master_printf("time for set_mu = %g s\n", wall_time() - tstart);
 }
 
 void structure::set_mu(double mufunc(const vec &), bool use_anisotropic_averaging, double tol,
@@ -488,7 +488,7 @@ void structure::set_conductivity(component c, material_function &C) {
   changing_chunks();
   for (int i = 0; i < num_chunks; i++)
     if (chunks[i]->is_mine()) chunks[i]->set_conductivity(c, C);
-  if (!quiet) master_printf("time for set_conductivity = %g s\n", wall_time() - tstart);
+  if (verbosity > 0) master_printf("time for set_conductivity = %g s\n", wall_time() - tstart);
 }
 
 void structure::set_conductivity(component c, double Cfunc(const vec &)) {
diff --git a/src/structure_dump.cpp b/src/structure_dump.cpp
index f062155f0..ed209dca2 100644
--- a/src/structure_dump.cpp
+++ b/src/structure_dump.cpp
@@ -84,7 +84,7 @@ void structure::dump_chunk_layout(const char *filename) {
 }
 
 void structure::dump(const char *filename) {
-  if (!quiet) master_printf("creating epsilon output file \"%s\"...\n", filename);
+  if (verbosity > 0) master_printf("creating epsilon output file \"%s\"...\n", filename);
 
   // make/save a num_chunks x NUM_FIELD_COMPONENTS x 5 array counting
   // the number of entries in the chi1inv array for each chunk.
@@ -479,7 +479,7 @@ void structure::load_chunk_layout(const std::vector<grid_volume> &gvs, boundary_
 void structure::load(const char *filename) {
   h5file file(filename, h5file::READONLY, true);
 
-  if (!quiet) master_printf("reading epsilon from file \"%s\"...\n", filename);
+  if (verbosity > 0) master_printf("reading epsilon from file \"%s\"...\n", filename);
 
   // make/save a num_chunks x NUM_FIELD_COMPONENTS x 5 array counting
   // the number of entries in the chi1inv array for each chunk.
diff --git a/src/time.cpp b/src/time.cpp
index 2c8bda5c3..125a96b2d 100644
--- a/src/time.cpp
+++ b/src/time.cpp
@@ -92,7 +92,7 @@ void fields::print_times() {
     pt(mean, stddev, (time_sink)i);
   master_printf("\n");
 
-  if (verbosity > 0) {
+  if (verbosity > 1) {
     master_printf("\nField time usage for all processes:\n");
     double *alltimes_tmp = new double[n * (Other + 1)];
     double *alltimes = new double[n * (Other + 1)];
diff --git a/tests/2D_convergence.cpp b/tests/2D_convergence.cpp
index 66ed7ba77..9e4f108be 100644
--- a/tests/2D_convergence.cpp
+++ b/tests/2D_convergence.cpp
@@ -102,7 +102,7 @@ void check_convergence(component c, double best_guess, double beta) {
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
 #ifdef HAVE_HARMINV
   master_printf("Running holes square-lattice resolution convergence test.\n");
   check_convergence(Ey, 0.179944, 0);  // from MPB; correct to >= 4 dec. places
diff --git a/tests/aniso_disp.cpp b/tests/aniso_disp.cpp
index 1436aa63b..acfa4c223 100644
--- a/tests/aniso_disp.cpp
+++ b/tests/aniso_disp.cpp
@@ -97,7 +97,7 @@ int main(int argc, char **argv) {
   bool ok = true;
   // we can only use one process for this 1-pixel simulation
   if (0 == divide_parallel_processes(count_processors())) {
-    quiet = true;
+    verbosity = 0;
     const double res = 200;
     grid_volume gv = vol3d(0, 0, 0, res);
     gv.center_origin();
diff --git a/tests/bench.cpp b/tests/bench.cpp
index 2905cfcf6..cf7658da7 100644
--- a/tests/bench.cpp
+++ b/tests/bench.cpp
@@ -269,7 +269,7 @@ bench bench_3d(const double xmax, const double ymax, const double zmax, double e
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   master_printf("Benchmarking with %d processor%s...\n", count_processors(),
                 count_processors() > 1 ? "s" : "");
 
diff --git a/tests/bragg_transmission.cpp b/tests/bragg_transmission.cpp
index 115c0e862..248f7c21b 100644
--- a/tests/bragg_transmission.cpp
+++ b/tests/bragg_transmission.cpp
@@ -240,7 +240,7 @@ void doit(bool use_hdf5) {
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
 
 #ifdef HAVE_HDF5
   doit(true);
diff --git a/tests/convergence_cyl_waveguide.cpp b/tests/convergence_cyl_waveguide.cpp
index 9ba32d501..8d98bdd0b 100644
--- a/tests/convergence_cyl_waveguide.cpp
+++ b/tests/convergence_cyl_waveguide.cpp
@@ -169,7 +169,7 @@ void test_convergence_with_averaging() {
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
 #ifdef HAVE_HARMINV
   master_printf("Testing convergence of a waveguide mode frequency...\n");
   test_convergence_without_averaging();
diff --git a/tests/cylindrical.cpp b/tests/cylindrical.cpp
index beaf04ad8..83a438523 100644
--- a/tests/cylindrical.cpp
+++ b/tests/cylindrical.cpp
@@ -299,7 +299,7 @@ int test_pattern(double eps(const vec &), int splitting, const char *mydirname)
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   const char *mydirname = "cylindrical-out";
   trash_output_directory(mydirname);
   master_printf("Testing cylindrical coords under different splittings...\n");
diff --git a/tests/flux.cpp b/tests/flux.cpp
index 9fc1b9d5a..b21567cf9 100644
--- a/tests/flux.cpp
+++ b/tests/flux.cpp
@@ -300,7 +300,7 @@ void attempt(const char *name, int allright) {
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   master_printf("Trying out the fluxes...\n");
 
   attempt("Split flux plane split by 7...", split_1d(cavity, 7));
diff --git a/tests/h5test.cpp b/tests/h5test.cpp
index 45c2102ed..281688323 100644
--- a/tests/h5test.cpp
+++ b/tests/h5test.cpp
@@ -354,7 +354,7 @@ int main(int argc, char **argv) {
   const double a = 10.0;
   initialize mpi(argc, argv);
   int chances;
-  quiet = true;
+  verbosity = 0;
 #ifdef HAVE_HDF5
   const double pad1 = 0.314159, pad2 = 0.27183, pad3 = 0.14142;
 
diff --git a/tests/harmonics.cpp b/tests/harmonics.cpp
index 1e9b22cde..b579605e2 100644
--- a/tests/harmonics.cpp
+++ b/tests/harmonics.cpp
@@ -93,7 +93,7 @@ int different(double a, double a0, double thresh, const char *msg) {
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   const double freq = 1.0 / 3.0;
 
   double a2, a3, a2_2, a3_2;
diff --git a/tests/integrate.cpp b/tests/integrate.cpp
index da5144ed4..9d1425977 100644
--- a/tests/integrate.cpp
+++ b/tests/integrate.cpp
@@ -315,7 +315,7 @@ void check_loop_vol(const grid_volume &gv, component c) {
 int main(int argc, char **argv) {
   const double a = 10.0;
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   const grid_volume v3d = vol3d(sz[0], sz[1], sz[2], a);
   const grid_volume v3d0 = vol3d(sz[0], sz[1], 0, a);
   const grid_volume v3d00 = vol3d(sz[0], 0, 0, a);
diff --git a/tests/known_results.cpp b/tests/known_results.cpp
index 40b1b45d2..4087d7b4a 100644
--- a/tests/known_results.cpp
+++ b/tests/known_results.cpp
@@ -150,7 +150,7 @@ double polariton_energy(const grid_volume &gv, double eps(const vec &)) {
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   const char *mydirname = "known_results-out";
   trash_output_directory(mydirname);
   master_printf("Testing with some known results...\n");
diff --git a/tests/one_dimensional.cpp b/tests/one_dimensional.cpp
index 46d64ef5a..2e198d445 100644
--- a/tests/one_dimensional.cpp
+++ b/tests/one_dimensional.cpp
@@ -151,7 +151,7 @@ int test_pattern(double eps(const vec &), int splitting, const char *mydirname)
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   const char *mydirname = "one_dimensional-out";
   master_printf("Testing one dimension under different splittings...\n");
 
diff --git a/tests/physical.cpp b/tests/physical.cpp
index 025f95927..d9956ad71 100644
--- a/tests/physical.cpp
+++ b/tests/physical.cpp
@@ -104,7 +104,7 @@ void attempt(const char *name, int allright) {
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  // quiet = true;
+  // verbosity = 0;
   master_printf("Trying out some physical tests...\n");
 
   attempt("radiating source should decay spatially as 1/sqrt(r) in 2D.", radiating_2D(8.0));
diff --git a/tests/pml.cpp b/tests/pml.cpp
index ea03f30e0..cde320351 100644
--- a/tests/pml.cpp
+++ b/tests/pml.cpp
@@ -318,7 +318,7 @@ int pmlcyl_scaling(double eps(const vec &), int m) {
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   master_printf("Running PML tests...\n");
   if (check_pml1d(one, 0)) abort("not a pml in 1d.");
   if (check_pml1d(one, 10.0)) abort("not a pml in 1d + conductivity.");
diff --git a/tests/stress_tensor.cpp b/tests/stress_tensor.cpp
index fb5ea545a..bccc13377 100644
--- a/tests/stress_tensor.cpp
+++ b/tests/stress_tensor.cpp
@@ -18,7 +18,7 @@ double two_waveguides(const vec &p) {
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   grid_volume gv = vol3d(sx + 2 * dpml, sy + 2 * dpml, 0, res);
   gv.center_origin();
   const symmetry S = mirror(X, gv) - mirror(Y, gv);
@@ -50,7 +50,7 @@ int main(int argc, char **argv) {
   double freqs_re[8], freqs_im[8];
 
   master_printf("done with timestepping, running harminv...\n");
-  int num = do_harminv(vals, iT2, f.dt, 0.19, 0.25, 8, 
+  int num = do_harminv(vals, iT2, f.dt, 0.19, 0.25, 8,
   		       amps, freqs_re, freqs_im);
   master_printf("harminv found %d modes\n",num);
   for (int i=0;i<num;i++)
diff --git a/tests/symmetry.cpp b/tests/symmetry.cpp
index ea92764aa..5e45cd305 100644
--- a/tests/symmetry.cpp
+++ b/tests/symmetry.cpp
@@ -920,7 +920,7 @@ double nonlinear_ex(const grid_volume &gv, double eps(const vec &)) {
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   trash_output_directory(mydirname);
   master_printf("Testing with various kinds of symmetry...\n");
 
diff --git a/tests/three_d.cpp b/tests/three_d.cpp
index dd7fdbf24..64a228754 100644
--- a/tests/three_d.cpp
+++ b/tests/three_d.cpp
@@ -234,7 +234,7 @@ int test_pml_splitting(double eps(const vec &), int splitting, const char *mydir
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   const char *mydirname = "three_d-out";
   trash_output_directory(mydirname);
   master_printf("Testing 3D...\n");
diff --git a/tests/two_dimensional.cpp b/tests/two_dimensional.cpp
index 2d7958411..463e742af 100644
--- a/tests/two_dimensional.cpp
+++ b/tests/two_dimensional.cpp
@@ -336,7 +336,7 @@ int test_pml_te(double eps(const vec &), int splitting, const char *mydirname) {
 
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
-  quiet = true;
+  verbosity = 0;
   const char *mydirname = "two_dimensional-out";
   trash_output_directory(mydirname);
   master_printf("Testing 2D...\n");