diff --git a/doc/docs/Python_Tutorials/Eigenmode_Source.md b/doc/docs/Python_Tutorials/Eigenmode_Source.md
index 8a2d2b82f..4f7fa99bc 100644
--- a/doc/docs/Python_Tutorials/Eigenmode_Source.md
+++ b/doc/docs/Python_Tutorials/Eigenmode_Source.md
@@ -17,7 +17,7 @@ The first structure, shown in the schematic below, is a 2d ridge waveguide with
 
 The simulation script is in [examples/oblique-source.py](https://github.com/NanoComp/meep/blob/master/python/examples/oblique-source.py). The notebook is [examples/oblique-source.ipynb](https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/oblique-source.ipynb).
 
-The simulation consists of two separate parts: (1) computing the Poynting flux and (2) plotting the field profile. The field profile is generated by setting the flag `compute_flux=False`. For the single-mode case, a constant-amplitude current source (`eig_src=False`) excites both the waveguide mode and radiating fields in *both* directions (i.e., forwards and backwards). This is shown in the main inset of the first of two figures above. The `EigenModeSource` excites only the *forward-going* waveguide mode **A** as shown in the smaller inset. Exciting this mode requires setting `eig_src=True`, `fsrc=0.15`, `kx=0.4`, and `bnum=1`. Note that `EigenModeSource` is a line centered at the origin extending the length of the entire cell. The constant-amplitude source is a line that is slightly larger than the waveguide width. The parameter `rot_angle` specifies the rotation angle of the waveguide axis and is initially 0° (i.e., straight or horizontal orientation). This enables `eig_parity` to include `EVEN_Y` in addition to `ODD_Z` and the simulation to include an overall mirror symmetry plane in the y-direction.
+The simulation consists of two separate parts: (1) computing the Poynting flux and (2) plotting the field profile. The field profile is generated by setting the flag `compute_flux=False`. For the single-mode case, a constant-amplitude current source (`eig_src=False`) excites both the waveguide mode and radiating fields in *both* directions (i.e., forwards and backwards). This is shown in the main inset of the first of two figures above. The `EigenModeSource` excites only the *forward-going* waveguide mode **A** as shown in the smaller inset. Launching this mode requires setting `eig_src=True`, `fsrc=0.15`, `kx=0.4`, and `bnum=1`. Note that the length of the `EigenModeSource` must be large enough that it captures the entire guided mode (i.e., the fields should decay to roughly zero at the edges). The line source does *not* need to span the entire length of the cell but should be slightly larger than the waveguide width. The parameter `rot_angle` specifies the rotation angle of the waveguide axis and is initially 0° (i.e., straight or horizontal orientation). This enables `eig_parity` to include `EVEN_Y` in addition to `ODD_Z` and the simulation to include an overall mirror symmetry plane in the y-direction.
 
 For the multi-mode case, a constant-amplitude current source excites a superposition of the two waveguide modes in addition to the radiating field. This is shown in the main inset of the second figure above. The `EigenModeSource` excites only a given mode: **A** (`fsrc=0.35`, `kx=0.4`, `bnum=2`) or **B** (`fsrc=0.35`, `kx=1.2`, `bnum=1`) as shown in the smaller insets.
 
@@ -35,8 +35,10 @@ pml_layers = [mp.PML(thickness=2)]
 # rotation angle (in degrees) of waveguide, counter clockwise (CCW) around z-axis
 rot_angle = np.radians(20)
 
+w = 1.0 # width of waveguide
+
 geometry = [mp.Block(center=mp.Vector3(),
-                     size=mp.Vector3(mp.inf,1,mp.inf),
+                     size=mp.Vector3(mp.inf,w,mp.inf),
                      e1=mp.Vector3(1).rotate(mp.Vector3(z=1), rot_angle),
                      e2=mp.Vector3(y=1).rotate(mp.Vector3(z=1), rot_angle),
                      material=mp.Medium(epsilon=12))]
@@ -54,7 +56,7 @@ eig_src = True # eigenmode (True) or constant-amplitude (False) source
 if eig_src:
     sources = [mp.EigenModeSource(src=mp.GaussianSource(fsrc,fwidth=0.2*fsrc) if compute_flux else mp.ContinuousSource(fsrc),
                                   center=mp.Vector3(),
-                                  size=mp.Vector3(y=14),
+                                  size=mp.Vector3(y=3*w),
                                   direction=mp.NO_DIRECTION,
                                   eig_kpoint=kpoint,
                                   eig_band=bnum,
@@ -63,7 +65,7 @@ if eig_src:
 else:
     sources = [mp.Source(src=mp.GaussianSource(fsrc,fwidth=0.2*fsrc) if compute_flux else mp.ContinuousSource(fsrc),
                          center=mp.Vector3(),
-                         size=mp.Vector3(y=2),
+                         size=mp.Vector3(y=3*w),
                          component=mp.Ez)]
 
 sim = mp.Simulation(cell_size=cell_size,
@@ -109,7 +111,7 @@ Finally, we demonstrate that the total power in a waveguide with *arbitrary* ori
 Planewaves in Homogeneous Media
 -------------------------------
 
-The eigenmode source can also be used to launch [planewaves](https://en.wikipedia.org/wiki/Plane_wave) in homogeneous media. The dispersion relation for a planewave is ω=|$\vec{k}$|/$n$ where ω is the angular frequency of the planewave and $\vec{k}$ its wavevector; $n$ is the refractive index of the homogeneous medium. This example demonstrates launching planewaves in a uniform medium with $n$ of 1.5 at three rotation angles: 0°, 20°, and 40°. Bloch-periodic boundaries via the `k_point` are used and specified by the wavevector $\vec{k}$. PML boundaries are used only along the x-direction.
+The eigenmode source can also be used to launch [planewaves](https://en.wikipedia.org/wiki/Plane_wave) in homogeneous media. The dispersion relation for a planewave is ω=|$\vec{k}$|/$n$ where ω is the angular frequency of the planewave and $\vec{k}$ its wavevector; $n$ is the refractive index of the homogeneous medium. This example demonstrates launching planewaves in a uniform medium with $n$ of 1.5 at three rotation angles: 0°, 20°, and 40°. Bloch-periodic boundaries via the `k_point` are used and specified by the wavevector $\vec{k}$. PML boundaries are used only along the $x$-direction.
 
 The simulation script is in [examples/oblique-planewave.py](https://github.com/NanoComp/meep/blob/master/python/examples/oblique-planewave.py). The notebook is in [examples/oblique-planewave.ipynb](https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/oblique-planewave.ipynb).
 
@@ -162,7 +164,7 @@ plt.axis('off')
 plt.show()
 ```
 
-Note that this example involves a `ContinuousSource` for the time profile. For a pulsed source, the oblique planewave is incident at a given angle for only a *single* frequency component of the source; see [Tutorial/Basics/Angular Reflectance Spectrum of a Planar Interface](../Python_Tutorials/Basics.md#angular-reflectance-spectrum-of-a-planar-interface). This is a fundamental feature of FDTD simulations and not of Meep per se. To simulate an incident planewave at multiple angles for a given frequency ω, you will need to do separate simulations involving different values of $\vec{k}$ (`k_point`) since each set of ($\vec{k}$,ω) specifying the Bloch-periodic boundaries and the frequency of the source will produce a different angle of the planewave. For more details, refer to Section 4.5 ("Efficient Frequency-Angle Coverage") in [Chapter 4](https://arxiv.org/abs/1301.5366) ("Electromagnetic Wave Source Conditions") of [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707).
+Note that the line source spans the *entire* length of the cell. (Planewave sources extending into the PML region must include `is_integrated=True`.) This example involves a continuous-wave (CW) time profile. For a pulse profile, the oblique planewave is incident at a given angle for only a *single* frequency component of the source; see [Tutorial/Basics/Angular Reflectance Spectrum of a Planar Interface](../Python_Tutorials/Basics.md#angular-reflectance-spectrum-of-a-planar-interface). This is a fundamental feature of FDTD simulations and not of Meep per se. To simulate an incident planewave at multiple angles for a given frequency ω, you will need to do separate simulations involving different values of $\vec{k}$ (`k_point`) since each set of ($\vec{k}$,ω) specifying the Bloch-periodic boundaries and the frequency of the source will produce a different angle of the planewave. For more details, refer to Section 4.5 ("Efficient Frequency-Angle Coverage") in [Chapter 4](https://arxiv.org/abs/1301.5366) ("Electromagnetic Wave Source Conditions") of [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707).
 
 Shown below are the steady-state field profiles generated by the planewave for the three rotation angles. Residues of the backward-propagating waves due to the discretization are slightly visible.
 
diff --git a/doc/docs/Scheme_Tutorials/Eigenmode_Source.md b/doc/docs/Scheme_Tutorials/Eigenmode_Source.md
index 622be5141..92d0f8122 100644
--- a/doc/docs/Scheme_Tutorials/Eigenmode_Source.md
+++ b/doc/docs/Scheme_Tutorials/Eigenmode_Source.md
@@ -17,7 +17,7 @@ The first structure, shown in the schematic below, is a 2d ridge waveguide with
 
 The simulation script is in [examples/oblique-source.ctl](https://github.com/NanoComp/meep/blob/master/scheme/examples/oblique-source.ctl).
 
-The simulation consists of two parts: (1) computing the Poynting flux and (2) plotting the field profile. The field profile is generated by setting the flag `compute-flux=false`. For the single-mode case, a constant-amplitude current source (`eig-src=false`) excites both the waveguide mode and radiating fields in *both* directions (i.e., forwards and backwards). This is shown in the main inset of the first of two figures above. The `eigenmode-source` excites only the *forward-going* waveguide mode **A** as shown in the smaller inset. Exciting this mode requires setting `eig-src=true`, `fsrc=0.15`, `kx=0.4`, and `bnum=1`. Note that `eigenmode-source` is a line centered at the origin extending the length of the entire cell. The constant-amplitude source is a line that is slightly larger than the waveguide width. The parameter `rot-angle` specifies the rotation angle of the waveguide axis and is initially 0° (i.e., straight or horizontal orientation). This enables `eig-parity` to include `EVEN-Y` in addition to `ODD-Z` and the cell to include an overall mirror symmetry plane in the y-direction.
+The simulation consists of two parts: (1) computing the Poynting flux and (2) plotting the field profile. The field profile is generated by setting the flag `compute-flux=false`. For the single-mode case, a constant-amplitude current source (`eig-src=false`) excites both the waveguide mode and radiating fields in *both* directions (i.e., forwards and backwards). This is shown in the main inset of the first of two figures above. The `eigenmode-source` excites only the *forward-going* waveguide mode **A** as shown in the smaller inset. Launching this mode requires setting `eig-src=true`, `fsrc=0.15`, `kx=0.4`, and `bnum=1`. Note that the length of the `EigenModeSource` must be large enough that it captures the entire guided mode (i.e., the fields should decay to roughly zero at the edges). The line source does *not* need to span the entire length of the cell but should be slightly larger than the waveguide width. The parameter `rot-angle` specifies the rotation angle of the waveguide axis and is initially 0° (i.e., straight or horizontal orientation). This enables `eig-parity` to include `EVEN-Y` in addition to `ODD-Z` and the cell to include an overall mirror symmetry plane in the y-direction.
 
 For the multi-mode case, a constant-amplitude current source excites a superposition of the two waveguide modes in addition to the radiating field. This is shown in the main inset of the second figure above. The `eigenmode-source` excites only a given mode: **A** (`fsrc=0.35`, `kx=0.4`, `bnum=2`) or **B** (`fsrc=0.35`, `kx=1.2`, `bnum=1`) as shown in the smaller insets.
 
@@ -32,9 +32,12 @@ For the multi-mode case, a constant-amplitude current source excites a superposi
 (define-param rot-angle 20)
 (set! rot-angle (deg->rad rot-angle))
 
+; width of waveguide
+(define-param w 1.0)
+
 (set! geometry (list (make block
                        (center 0 0 0)
-                       (size infinity 1 infinity)
+                       (size infinity w infinity)
                        (e1 (rotate-vector3 (vector3 0 0 1) rot-angle (vector3 1 0 0)))
                        (e2 (rotate-vector3 (vector3 0 0 1) rot-angle (vector3 0 1 0)))
                        (material (make medium (epsilon 12))))))
@@ -54,7 +57,7 @@ For the multi-mode case, a constant-amplitude current source excites a superposi
                    (make eigenmode-source
                      (src (if compute-flux? (make gaussian-src (frequency fsrc) (fwidth (* 0.2 fsrc))) (make continuous-src (frequency fsrc))))
                      (center 0 0 0)
-                     (size 0 14 0)
+                     (size 0 (* 3 w) 0)
                      (direction (if (= rot-angle 0) AUTOMATIC NO-DIRECTION))
                      (eig-kpoint kpoint)
                      (eig-band bnum)
@@ -63,7 +66,7 @@ For the multi-mode case, a constant-amplitude current source excites a superposi
                    (make source
                      (src (if compute-flux? (make gaussian-src (frequency fsrc) (fwidth (* 0.2 fsrc))) (make continuous-src (frequency fsrc))))
                      (center 0 0 0)
-                     (size 0 2 0)
+                     (size 0 (* 3 w) 0)
                      (component Ez)))))
 
 (if (= rot-angle 0)
@@ -103,7 +106,7 @@ Finally, we demonstrate that the total power in a waveguide with *arbitrary* ori
 Planewaves in Homogeneous Media
 -------------------------------
 
-The eigenmode source can also be used to launch [planewaves](https://en.wikipedia.org/wiki/Plane_wave) in homogeneous media. The dispersion relation for a planewave is ω=|$\vec{k}$|/$n$ where ω is the angular frequency of the planewave and $\vec{k}$ its wavevector; $n$ is the refractive index of the homogeneous medium. This example demonstrates launching planewaves in a uniform medium with $n$ of 1.5 at three rotation angles: 0°, 20°, and 40°. Bloch-periodic boundaries via the `k-point` are used and specified by the wavevector $\vec{k}$. PML boundaries are used only along the x-direction.
+The eigenmode source can also be used to launch [planewaves](https://en.wikipedia.org/wiki/Plane_wave) in homogeneous media. The dispersion relation for a planewave is ω=|$\vec{k}$|/$n$ where ω is the angular frequency of the planewave and $\vec{k}$ its wavevector; $n$ is the refractive index of the homogeneous medium. This example demonstrates launching planewaves in a uniform medium with $n$ of 1.5 at three rotation angles: 0°, 20°, and 40°. Bloch-periodic boundaries via the `k-point` are used and specified by the wavevector $\vec{k}$. PML boundaries are used only along the $x$-direction.
 
 The simulation script is in [examples/oblique-planewave.ctl](https://github.com/NanoComp/meep/blob/master/scheme/examples/oblique-planewave.ctl).
 
@@ -144,7 +147,7 @@ The simulation script is in [examples/oblique-planewave.ctl](https://github.com/
                           (at-end output-efield-z)))
 ```
 
-Note that this example involves a `continuous-source` for the time profile. For a pulsed source, the oblique planewave is incident at a given angle for only a *single* frequency component of the source; see [Tutorial/Basics/Angular Reflectance Spectrum of a Planar Interface](../Scheme_Tutorials/Basics.md#angular-reflectance-spectrum-of-a-planar-interface). This is a fundamental feature of FDTD simulations and not of Meep per se. To simulate an incident planewave at multiple angles for a given frequency ω, you will need to do separate simulations involving different values of $\vec{k}$ (`k-point`) since each set of ($\vec{k}$,ω) specifying the Bloch-periodic boundaries and the frequency of the source will produce a different angle of the planewave. For more details, refer to Section 4.5 ("Efficient Frequency-Angle Coverage") in [Chapter 4](https://arxiv.org/abs/1301.5366) ("Electromagnetic Wave Source Conditions") of [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707).
+Note that the line source spans the *entire* length of the cell. (Planewave sources extending into the PML region must include `(is-integrated true)`.) This example involves a continuous-wave (CW) time profile. For a pulse profile, the oblique planewave is incident at a given angle for only a *single* frequency component of the source; see [Tutorial/Basics/Angular Reflectance Spectrum of a Planar Interface](../Scheme_Tutorials/Basics.md#angular-reflectance-spectrum-of-a-planar-interface). This is a fundamental feature of FDTD simulations and not of Meep per se. To simulate an incident planewave at multiple angles for a given frequency ω, you will need to do separate simulations involving different values of $\vec{k}$ (`k-point`) since each set of ($\vec{k}$,ω) specifying the Bloch-periodic boundaries and the frequency of the source will produce a different angle of the planewave. For more details, refer to Section 4.5 ("Efficient Frequency-Angle Coverage") in [Chapter 4](https://arxiv.org/abs/1301.5366) ("Electromagnetic Wave Source Conditions") of [Advances in FDTD Computational Electrodynamics: Photonics and Nanotechnology](https://www.amazon.com/Advances-FDTD-Computational-Electrodynamics-Nanotechnology/dp/1608071707).
 
 Shown below are the steady-state field profiles generated by the planewave for the three rotation angles. Residues of the backward-propagating waves due to the discretization are slightly visible.
 
diff --git a/doc/docs/images/eigenmode_source.png b/doc/docs/images/eigenmode_source.png
index cb66527d4..d4d848eb5 100644
Binary files a/doc/docs/images/eigenmode_source.png and b/doc/docs/images/eigenmode_source.png differ
diff --git a/doc/docs/images/oblique_source_multimode.png b/doc/docs/images/oblique_source_multimode.png
index d158b7d52..cc1de2deb 100644
Binary files a/doc/docs/images/oblique_source_multimode.png and b/doc/docs/images/oblique_source_multimode.png differ
diff --git a/doc/docs/images/oblique_source_singlemode.png b/doc/docs/images/oblique_source_singlemode.png
index 7d70024f7..17a6b39eb 100644
Binary files a/doc/docs/images/oblique_source_singlemode.png and b/doc/docs/images/oblique_source_singlemode.png differ
diff --git a/python/examples/oblique-source.py b/python/examples/oblique-source.py
index 7b183b861..36de5fdc5 100644
--- a/python/examples/oblique-source.py
+++ b/python/examples/oblique-source.py
@@ -11,8 +11,10 @@
 # rotation angle (in degrees) of waveguide, counter clockwise (CCW) around z-axis
 rot_angle = np.radians(20)
 
+w = 1.0 # width of waveguide
+
 geometry = [mp.Block(center=mp.Vector3(),
-                     size=mp.Vector3(mp.inf,1,mp.inf),
+                     size=mp.Vector3(mp.inf,w,mp.inf),
                      e1=mp.Vector3(1).rotate(mp.Vector3(z=1), rot_angle),
                      e2=mp.Vector3(y=1).rotate(mp.Vector3(z=1), rot_angle),
                      material=mp.Medium(epsilon=12))]
@@ -30,7 +32,7 @@
 if eig_src:
     sources = [mp.EigenModeSource(src=mp.GaussianSource(fsrc,fwidth=0.2*fsrc) if compute_flux else mp.ContinuousSource(fsrc),
                                   center=mp.Vector3(),
-                                  size=mp.Vector3(y=14),
+                                  size=mp.Vector3(y=3*w),
                                   direction=mp.NO_DIRECTION,
                                   eig_kpoint=kpoint,
                                   eig_band=bnum,
@@ -39,7 +41,7 @@
 else:
     sources = [mp.Source(src=mp.GaussianSource(fsrc,fwidth=0.2*fsrc) if compute_flux else mp.ContinuousSource(fsrc),
                          center=mp.Vector3(),
-                         size=mp.Vector3(y=2),
+                         size=mp.Vector3(y=3*w),
                          component=mp.Ez)]
 
 sim = mp.Simulation(cell_size=cell_size,
diff --git a/scheme/examples/oblique-source.ctl b/scheme/examples/oblique-source.ctl
index 07c01109f..c5ad6f198 100644
--- a/scheme/examples/oblique-source.ctl
+++ b/scheme/examples/oblique-source.ctl
@@ -8,6 +8,9 @@
 (define-param rot-angle 20)
 (set! rot-angle (deg->rad rot-angle))
 
+; width of waveguide
+(define-param w 1.0)
+
 (set! geometry (list (make block
                        (center 0 0 0)
                        (size infinity 1 infinity)
@@ -30,7 +33,7 @@
                    (make eigenmode-source
                      (src (if compute-flux? (make gaussian-src (frequency fsrc) (fwidth (* 0.2 fsrc))) (make continuous-src (frequency fsrc))))
                      (center 0 0 0)
-                     (size 0 14 0)
+                     (size 0 (* 3 w) 0)
                      (direction (if (= rot-angle 0) AUTOMATIC NO-DIRECTION))
                      (eig-kpoint kpoint)
                      (eig-band bnum)
@@ -39,7 +42,7 @@
                    (make source
                      (src (if compute-flux? (make gaussian-src (frequency fsrc) (fwidth (* 0.2 fsrc))) (make continuous-src (frequency fsrc))))
                      (center 0 0 0)
-                     (size 0 2 0)
+                     (size 0 (* 3 w) 0)
                      (component Ez)))))
 
 (if (= rot-angle 0)
diff --git a/src/meep.hpp b/src/meep.hpp
index bd5769444..4d95e8289 100644
--- a/src/meep.hpp
+++ b/src/meep.hpp
@@ -1933,6 +1933,7 @@ class flux_vol {
 const char *make_output_directory(const char *exename, const char *jobname = NULL);
 char *make_output_directory(); // make temporary directory
 void trash_output_directory(const char *dirname);
+void delete_directory(const char *path);
 FILE *create_output_file(const char *dirname, const char *fname);
 
 // The following allows you to hit ctrl-C to tell your calculation to stop
diff --git a/src/output_directory.cpp b/src/output_directory.cpp
index 8dcdb84ac..aa50a0489 100644
--- a/src/output_directory.cpp
+++ b/src/output_directory.cpp
@@ -23,6 +23,7 @@
 #include <dirent.h>
 #include <string.h>
 #include <unistd.h>
+#include <ftw.h>
 
 #include "meep.hpp"
 
@@ -160,4 +161,15 @@ void trash_output_directory(const char *dirname) {
   if (am_master()) mkdir(dirname, 00777);
 }
 
+static int rmpath(const char *path, const struct stat *s, int t, struct FTW *ftw) {
+  (void) s; (void) t; (void) ftw; // unused
+  return remove(path);
+}
+
+// equivalent to rm -rf path
+void delete_directory(const char *path) {
+  if (am_master())
+    nftw(path, rmpath, 10, FTW_DEPTH|FTW_MOUNT|FTW_PHYS);
+}
+
 } // namespace meep
diff --git a/tests/cylindrical.cpp b/tests/cylindrical.cpp
index 55c3a29d3..de202a921 100644
--- a/tests/cylindrical.cpp
+++ b/tests/cylindrical.cpp
@@ -59,15 +59,13 @@ int compare_point(fields &f1, fields &f2, const vec &p, double eps = 4e-8) {
   return 1;
 }
 
-int test_simple_periodic(double eps(const vec &), int splitting, const char *mydirname) {
+int test_simple_periodic(double eps(const vec &), int splitting) {
   double a = 10.0;
   double ttot = 30.0;
 
   grid_volume gv = volcyl(1.5, 0.8, a);
   structure s1(gv, eps, no_pml(), identity(), 0, 0.4);
   structure s(gv, eps, no_pml(), identity(), splitting, 0.4);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
   for (int m = 0; m < 3; m++) {
     char m_str[10];
     snprintf(m_str, 10, "%d", m);
@@ -107,15 +105,13 @@ int test_simple_periodic(double eps(const vec &), int splitting, const char *myd
   return 1;
 }
 
-int test_simple_metallic(double eps(const vec &), int splitting, const char *mydirname) {
+int test_simple_metallic(double eps(const vec &), int splitting) {
   double a = 10.0;
   double ttot = 30.0;
 
   grid_volume gv = volcyl(1.5, 0.8, a);
   structure s1(gv, eps, no_pml(), identity(), 0, 0.4);
   structure s(gv, eps, no_pml(), identity(), splitting, 0.4);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
   for (int m = 0; m < 3; m++) {
     char m_str[10];
     snprintf(m_str, 10, "%d", m);
@@ -155,12 +151,11 @@ int test_simple_metallic(double eps(const vec &), int splitting, const char *myd
 
 static bool issmall(std::complex<double> x) { return abs(x) < 1e-16; }
 
-int test_r_equals_zero(double eps(const vec &), const char *mydirname) {
+int test_r_equals_zero(double eps(const vec &)) {
   double a = 10.0;
   double ttot = 3.0;
   grid_volume gv = volcyl(1.5, 0.8, a);
   structure s(gv, eps, no_pml(), identity(), 0, 0.4);
-  s.set_output_directory(mydirname);
   for (int m = 0; m < 3; m++) {
     char m_str[10];
     snprintf(m_str, 10, "%d", m);
@@ -201,15 +196,13 @@ int test_r_equals_zero(double eps(const vec &), const char *mydirname) {
   return 1;
 }
 
-int test_pml(double eps(const vec &), int splitting, const char *mydirname) {
+int test_pml(double eps(const vec &), int splitting) {
   double a = 8;
   double ttot = 25.0;
 
   grid_volume gv = volcyl(3.5, 10.0, a);
   structure s1(gv, eps, pml(2.0), identity(), 0, 0.4);
   structure s(gv, eps, pml(2.0), identity(), splitting, 0.4);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
   for (int m = 0; m < 3; m++) {
     char m_str[10];
     snprintf(m_str, 10, "%d", m);
@@ -261,13 +254,11 @@ complex<double> checkers(const vec &pt) {
   return 0.0;
 }
 
-int test_pattern(double eps(const vec &), int splitting, const char *mydirname) {
+int test_pattern(double eps(const vec &), int splitting) {
   double a = 10.0;
   grid_volume gv = volcyl(1.5, 0.8, a);
   structure s1(gv, eps);
   structure s(gv, eps, no_pml(), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
   for (int m = 0; m < 1; m++) {
     char m_str[10];
     snprintf(m_str, 10, "%d", m);
@@ -300,33 +291,31 @@ int test_pattern(double eps(const vec &), int splitting, const char *mydirname)
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
   verbosity = 0;
-  const char *mydirname = "cylindrical-out";
-  trash_output_directory(mydirname);
   master_printf("Testing cylindrical coords under different splittings...\n");
 
-  if (!test_r_equals_zero(one, mydirname)) abort("error in test_r_equals_zero");
+  if (!test_r_equals_zero(one)) abort("error in test_r_equals_zero");
 
   for (int s = 2; s < 6; s++)
-    if (!test_pattern(one, s, mydirname)) abort("error in test_pattern\n");
-  // if (!test_pattern(one, 8, mydirname)) abort("error in crazy test_pattern\n");
-  // if (!test_pattern(one, 120, mydirname)) abort("error in crazy test_pattern\n");
+    if (!test_pattern(one, s)) abort("error in test_pattern\n");
+  // if (!test_pattern(one, 8)) abort("error in crazy test_pattern\n");
+  // if (!test_pattern(one, 120)) abort("error in crazy test_pattern\n");
 
   for (int s = 2; s < 4; s++)
-    if (!test_simple_periodic(one, s, mydirname)) abort("error in test_simple_periodic\n");
-  // if (!test_simple_periodic(one, 8, mydirname))
+    if (!test_simple_periodic(one, s)) abort("error in test_simple_periodic\n");
+  // if (!test_simple_periodic(one, 8))
   //  abort("error in crazy test_simple_periodic\n");
-  // if (!test_simple_periodic(one, 120, mydirname))
+  // if (!test_simple_periodic(one, 120))
   //  abort("error in crazy test_simple_periodic\n");
 
   for (int s = 2; s < 5; s++)
-    if (!test_simple_metallic(one, s, mydirname)) abort("error in test_simple_metallic\n");
-  // if (!test_simple_metallic(one, 8, mydirname))
+    if (!test_simple_metallic(one, s)) abort("error in test_simple_metallic\n");
+  // if (!test_simple_metallic(one, 8))
   //  abort("error in crazy test_simple_metallic\n");
-  // if (!test_simple_metallic(one, 120, mydirname))
+  // if (!test_simple_metallic(one, 120))
   //  abort("error in crazy test_simple_metallic\n");
 
   for (int s = 2; s < 6; s++)
-    if (!test_pml(one, s, mydirname)) abort("error in test_pml\n");
+    if (!test_pml(one, s)) abort("error in test_pml\n");
 
   return 0;
 }
diff --git a/tests/h5test.cpp b/tests/h5test.cpp
index 281688323..bd0bad2d0 100644
--- a/tests/h5test.cpp
+++ b/tests/h5test.cpp
@@ -3,6 +3,7 @@
 #include <stdio.h>
 #include <string.h>
 #include <stdlib.h>
+#include <unistd.h>
 
 #include <meep.hpp>
 #include "meep_internals.hpp"
@@ -58,9 +59,10 @@ double get_reim(complex<double> x, int reim) { return reim ? imag(x) : real(x);
 
 bool check_2d(double eps(const vec &), double a, int splitting, symfunc Sf, double kx, double ky,
               component src_c, int file_c, volume file_gv, bool real_fields, int expected_rank,
-              const char *name) {
+              const char *name, const char *mydirname) {
   const grid_volume gv = vol2d(xsize, ysize, a);
   structure s(gv, eps, no_pml(), Sf(gv), splitting);
+  s.set_output_directory(mydirname);
   fields f(&s);
 
   f.use_bloch(X, real_fields ? 0.0 : kx);
@@ -168,9 +170,11 @@ bool check_2d(double eps(const vec &), double a, int splitting, symfunc Sf, doub
 }
 
 bool check_3d(double eps(const vec &), double a, int splitting, symfunc Sf, component src_c,
-              int file_c, volume file_gv, bool real_fields, int expected_rank, const char *name) {
+              int file_c, volume file_gv, bool real_fields, int expected_rank, const char *name,
+              const char *mydirname) {
   const grid_volume gv = vol3d(xsize, ysize, zsize, a);
   structure s(gv, eps, no_pml(), Sf(gv), splitting);
+  s.set_output_directory(mydirname);
   fields f(&s);
 
   if (real_fields) f.use_real_fields();
@@ -274,9 +278,11 @@ bool check_3d(double eps(const vec &), double a, int splitting, symfunc Sf, comp
 }
 
 bool check_2d_monitor(double eps(const vec &), double a, int splitting, symfunc Sf, component src_c,
-                      int file_c, const vec &pt, bool real_fields, const char *name) {
+                      int file_c, const vec &pt, bool real_fields, const char *name,
+                      const char *mydirname) {
   const grid_volume gv = vol2d(xsize, ysize, a);
   structure s(gv, eps, no_pml(), Sf(gv), splitting);
+  s.set_output_directory(mydirname);
   fields f(&s);
 
   if (real_fields) f.use_real_fields();
@@ -355,6 +361,7 @@ int main(int argc, char **argv) {
   initialize mpi(argc, argv);
   int chances;
   verbosity = 0;
+  const char *temp_dir = make_output_directory();
 #ifdef HAVE_HDF5
   const double pad1 = 0.314159, pad2 = 0.27183, pad3 = 0.14142;
 
@@ -377,7 +384,7 @@ int main(int argc, char **argv) {
   if (!check_2d(funky_eps_2d, a, 1,
 		Sf2[3], Sf2_kx[3], Sf2_ky[3],
 		Ez, tm_c[3], gv_2d[1],
-		1, gv_2d_rank[1], "initial check"))
+		1, gv_2d_rank[1], "initial check", temp_dir))
     return 1;
 #endif
 
@@ -397,7 +404,7 @@ int main(int argc, char **argv) {
                        gv_2d_name[igv], component_name(tm_c[ic]), use_real ? "_r" : "");
               master_printf("Checking %s...\n", name);
               if (!check_2d(funky_eps_2d, a, splitting, Sf2[iS], Sf2_kx[iS], Sf2_ky[iS], Ez,
-                            tm_c[ic], gv_2d[igv], use_real, gv_2d_rank[igv], name))
+                            tm_c[ic], gv_2d[igv], use_real, gv_2d_rank[igv], name, temp_dir))
                 return 1;
             }
 
@@ -411,7 +418,7 @@ int main(int argc, char **argv) {
                      component_name(tm_c[ic]), use_real ? "_r" : "");
             master_printf("Checking %s...\n", name);
             if (!check_2d_monitor(funky_eps_2d, a, splitting, Sf2[iS], Ez, tm_c[ic],
-                                  vec(pad1, pad2), use_real, name))
+                                  vec(pad1, pad2), use_real, name, temp_dir))
               return 1;
           }
 
@@ -438,10 +445,13 @@ int main(int argc, char **argv) {
                      gv_3d_name[igv], component_name(c3d[ic]), use_real ? "_r" : "");
             master_printf("Checking %s...\n", name);
             if (!check_3d(funky_eps_3d, a, splitting, Sf3[iS], Ez, c3d[ic], gv_3d[igv], use_real,
-                          gv_3d_rank[igv], name))
+                          gv_3d_rank[igv], name, temp_dir))
               return 1;
           }
       }
 #endif /* HAVE_HDF5 */
+
+  delete_directory(temp_dir);
+
   return 0;
 }
diff --git a/tests/known_results.cpp b/tests/known_results.cpp
index 696daeace..4fed74c08 100644
--- a/tests/known_results.cpp
+++ b/tests/known_results.cpp
@@ -151,8 +151,6 @@ double polariton_energy(const grid_volume &gv, double eps(const vec &)) {
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
   verbosity = 0;
-  const char *mydirname = "known_results-out";
-  trash_output_directory(mydirname);
   master_printf("Testing with some known results...\n");
   const double a = 10.0;
 
diff --git a/tests/one_dimensional.cpp b/tests/one_dimensional.cpp
index 4db36d98b..560afae42 100644
--- a/tests/one_dimensional.cpp
+++ b/tests/one_dimensional.cpp
@@ -62,15 +62,13 @@ int compare_point(fields &f1, fields &f2, const vec &p) {
   return 1;
 }
 
-int test_simple_periodic(double eps(const vec &), int splitting, const char *mydirname) {
+int test_simple_periodic(double eps(const vec &), int splitting) {
   double a = 10.0;
   double ttot = 170.0;
 
   grid_volume gv = volone(6.0, a);
   structure s1(gv, eps);
   structure s(gv, eps, no_pml(), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
 
   master_printf("Trying splitting into %d chunks...\n", splitting);
   fields f(&s);
@@ -115,13 +113,11 @@ complex<double> checkers(const vec &pt) {
   return 1.0;
 }
 
-int test_pattern(double eps(const vec &), int splitting, const char *mydirname) {
+int test_pattern(double eps(const vec &), int splitting) {
   double a = 10.0;
   grid_volume gv = volone(6.0, a);
   structure s1(gv, eps);
   structure s(gv, eps, no_pml(), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
 
   master_printf("Trying test pattern with %d chunks...\n", splitting);
   fields f(&s);
@@ -152,13 +148,13 @@ int test_pattern(double eps(const vec &), int splitting, const char *mydirname)
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
   verbosity = 0;
-  const char *mydirname = "one_dimensional-out";
   master_printf("Testing one dimension under different splittings...\n");
 
   for (int s = 2; s < 7; s++)
-    if (!test_pattern(one, s, mydirname)) abort("error in test_pattern\n");
+    if (!test_pattern(one, s)) abort("error in test_pattern\n");
 
   for (int s = 2; s < 7; s++)
-    if (!test_simple_periodic(one, s, mydirname)) abort("error in test_simple_periodic\n");
+    if (!test_simple_periodic(one, s)) abort("error in test_simple_periodic\n");
+
   return 0;
 }
diff --git a/tests/symmetry.cpp b/tests/symmetry.cpp
index 84f7eb757..323530cb7 100644
--- a/tests/symmetry.cpp
+++ b/tests/symmetry.cpp
@@ -23,8 +23,6 @@
 using namespace meep;
 using std::complex;
 
-const char *mydirname = "symmetry-out";
-
 double one(const vec &) { return 1.0; }
 vec the_center;
 double rods_2d(const vec &pp) {
@@ -569,8 +567,6 @@ int test_yperiodic_ymirror(double eps(const vec &)) {
   const symmetry S = mirror(Y, gv);
   structure s(gv, eps, no_pml(), S);
   structure s1(gv, eps);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
   master_printf("Testing Y periodic with mirror symmetry...\n");
 
   fields f1(&s1);
@@ -829,8 +825,6 @@ int exact_pml_rot2x_tm(double eps(const vec &)) {
 
   structure s(gv, eps, pml(1.0), S);
   structure s1(gv, eps, pml(1.0), identity());
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
   master_printf("Testing X twofold rotational symmetry with PML...\n");
 
   fields f1(&s1);
@@ -921,7 +915,6 @@ double nonlinear_ex(const grid_volume &gv, double eps(const vec &)) {
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
   verbosity = 0;
-  trash_output_directory(mydirname);
   master_printf("Testing with various kinds of symmetry...\n");
 
   if (!test_1d_periodic_mirror(one)) abort("error in test_1d_periodic_mirror vacuum\n");
diff --git a/tests/three_d.cpp b/tests/three_d.cpp
index 8293b319d..d6e9fdb19 100644
--- a/tests/three_d.cpp
+++ b/tests/three_d.cpp
@@ -92,15 +92,13 @@ int approx_point(fields &f1, fields &f2, const vec &p) {
   return 1;
 }
 
-int test_metal(double eps(const vec &), int splitting, const char *mydirname) {
+int test_metal(double eps(const vec &), int splitting) {
   double a = 10.0;
   double ttot = 17.0;
 
   grid_volume gv = vol3d(1.5, 0.5, 1.0, a);
   structure s1(gv, eps);
   structure s(gv, eps, no_pml(), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
 
   master_printf("Metal test using %d chunks...\n", splitting);
   fields f(&s);
@@ -128,15 +126,13 @@ int test_metal(double eps(const vec &), int splitting, const char *mydirname) {
   return 1;
 }
 
-int test_periodic(double eps(const vec &), int splitting, const char *mydirname) {
+int test_periodic(double eps(const vec &), int splitting) {
   double a = 10.0;
   double ttot = 17.0;
 
   grid_volume gv = vol3d(1.5, 0.5, 1.0, a);
   structure s1(gv, eps);
   structure s(gv, eps, no_pml(), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
 
   master_printf("Periodic test using %d chunks...\n", splitting);
   fields f(&s);
@@ -166,12 +162,11 @@ int test_periodic(double eps(const vec &), int splitting, const char *mydirname)
   return 1;
 }
 
-int test_pml(double eps(const vec &), const char *mydirname) {
+int test_pml(double eps(const vec &)) {
   double a = 10.0;
 
   grid_volume gv = vol3d(1.5, 1.0, 1.2, a);
   structure s(gv, eps, pml(0.401));
-  s.set_output_directory(mydirname);
 
   master_printf("Testing pml quality...\n");
   fields f(&s);
@@ -200,14 +195,12 @@ int test_pml(double eps(const vec &), const char *mydirname) {
   return 1;
 }
 
-int test_pml_splitting(double eps(const vec &), int splitting, const char *mydirname) {
+int test_pml_splitting(double eps(const vec &), int splitting) {
   double a = 10.0;
 
   grid_volume gv = vol3d(1.5, 1.0, 1.2, a);
   structure s1(gv, eps, pml(0.3));
   structure s(gv, eps, pml(0.3), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
 
   master_printf("Testing pml while splitting into %d chunks...\n", splitting);
   fields f(&s);
@@ -235,20 +228,18 @@ int test_pml_splitting(double eps(const vec &), int splitting, const char *mydir
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
   verbosity = 0;
-  const char *mydirname = "three_d-out";
-  trash_output_directory(mydirname);
   master_printf("Testing 3D...\n");
 
-  if (!test_pml(one, mydirname)) abort("error in test_pml vacuum\n");
+  if (!test_pml(one)) abort("error in test_pml vacuum\n");
 
   for (int s = 2; s < 7; s++)
-    if (!test_periodic(targets, s, mydirname)) abort("error in test_periodic targets\n");
+    if (!test_periodic(targets, s)) abort("error in test_periodic targets\n");
 
   for (int s = 2; s < 8; s++)
-    if (!test_metal(one, s, mydirname)) abort("error in test_metal vacuum\n");
+    if (!test_metal(one, s)) abort("error in test_metal vacuum\n");
 
   for (int s = 2; s < 4; s++)
-    if (!test_pml_splitting(one, s, mydirname)) abort("error in test_pml_splitting vacuum\n");
+    if (!test_pml_splitting(one, s)) abort("error in test_pml_splitting vacuum\n");
 
   return 0;
 }
diff --git a/tests/two_dimensional.cpp b/tests/two_dimensional.cpp
index 5e97764f0..632fd8cf5 100644
--- a/tests/two_dimensional.cpp
+++ b/tests/two_dimensional.cpp
@@ -70,15 +70,13 @@ int compare_point(fields &f1, fields &f2, const vec &p) {
   return 1;
 }
 
-int test_metal(double eps(const vec &), int splitting, const char *mydirname) {
+int test_metal(double eps(const vec &), int splitting) {
   double a = 10.0;
   double ttot = 17.0;
 
   grid_volume gv = voltwo(3.0, 2.0, a);
   structure s1(gv, eps);
   structure s(gv, eps, no_pml(), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
 
   s.add_susceptibility(one, E_stuff, lorentzian_susceptibility(0.3, 0.1));
   s1.add_susceptibility(one, E_stuff, lorentzian_susceptibility(0.3, 0.1));
@@ -111,15 +109,13 @@ int test_metal(double eps(const vec &), int splitting, const char *mydirname) {
   return 1;
 }
 
-int test_periodic(double eps(const vec &), int splitting, const char *mydirname) {
+int test_periodic(double eps(const vec &), int splitting) {
   double a = 10.0;
   double ttot = 17.0;
 
   grid_volume gv = voltwo(3.0, 2.0, a);
   structure s1(gv, eps);
   structure s(gv, eps, no_pml(), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
 
   master_printf("Periodic test using %d chunks...\n", splitting);
   fields f(&s);
@@ -151,15 +147,13 @@ int test_periodic(double eps(const vec &), int splitting, const char *mydirname)
   return 1;
 }
 
-int test_periodic_tm(double eps(const vec &), int splitting, const char *mydirname) {
+int test_periodic_tm(double eps(const vec &), int splitting) {
   double a = 10.0;
   double ttot = 17.0;
 
   grid_volume gv = voltwo(3.0, 2.0, a);
   structure s1(gv, eps);
   structure s(gv, eps, no_pml(), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
 
   master_printf("Periodic 2D TM test using %d chunks...\n", splitting);
   fields f(&s);
@@ -189,14 +183,12 @@ int test_periodic_tm(double eps(const vec &), int splitting, const char *mydirna
   return 1;
 }
 
-int test_pml(double eps(const vec &), int splitting, const char *mydirname) {
+int test_pml(double eps(const vec &), int splitting) {
   double a = 10.0;
 
   grid_volume gv = voltwo(3.0, 2.0, a);
   structure s1(gv, eps, pml(1.0, X) + pml(1.0, Y, High));
   structure s(gv, eps, pml(1.0, X) + pml(1.0, Y, High), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
 
   master_printf("Testing pml while splitting into %d chunks...\n", splitting);
   fields f(&s);
@@ -238,14 +230,12 @@ int test_pml(double eps(const vec &), int splitting, const char *mydirname) {
   return 1;
 }
 
-int test_pml_tm(double eps(const vec &), int splitting, const char *mydirname) {
+int test_pml_tm(double eps(const vec &), int splitting) {
   double a = 10.0;
 
   grid_volume gv = voltwo(3.0, 3.0, a);
   structure s1(gv, eps, pml(1.0));
   structure s(gv, eps, pml(1.0), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
 
   master_printf("Testing TM pml while splitting into %d chunks...\n", splitting);
   fields f(&s);
@@ -285,14 +275,12 @@ int test_pml_tm(double eps(const vec &), int splitting, const char *mydirname) {
   return 1;
 }
 
-int test_pml_te(double eps(const vec &), int splitting, const char *mydirname) {
+int test_pml_te(double eps(const vec &), int splitting) {
   double a = 10.0;
 
   grid_volume gv = voltwo(3.0, 3.0, a);
   structure s1(gv, eps, pml(1.0));
   structure s(gv, eps, pml(1.0), identity(), splitting);
-  s.set_output_directory(mydirname);
-  s1.set_output_directory(mydirname);
 
   master_printf("Testing TE pml while splitting into %d chunks...\n", splitting);
   fields f(&s);
@@ -337,34 +325,32 @@ int test_pml_te(double eps(const vec &), int splitting, const char *mydirname) {
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
   verbosity = 0;
-  const char *mydirname = "two_dimensional-out";
-  trash_output_directory(mydirname);
   master_printf("Testing 2D...\n");
 
   for (int s = 2; s < 4; s++)
-    if (!test_pml(one, s, mydirname)) abort("error in test_pml vacuum\n");
+    if (!test_pml(one, s)) abort("error in test_pml vacuum\n");
 
   for (int s = 2; s < 4; s++)
-    if (!test_pml_tm(one, s, mydirname)) abort("error in test_pml_tm vacuum\n");
+    if (!test_pml_tm(one, s)) abort("error in test_pml_tm vacuum\n");
 
   for (int s = 2; s < 4; s++)
-    if (!test_pml_te(one, s, mydirname)) abort("error in test_pml_te vacuum\n");
+    if (!test_pml_te(one, s)) abort("error in test_pml_te vacuum\n");
 
   for (int s = 2; s < 4; s++)
-    if (!test_metal(one, s, mydirname)) abort("error in test_metal vacuum\n");
-  // if (!test_metal(one, 200, mydirname)) abort("error in test_metal vacuum\n");
+    if (!test_metal(one, s)) abort("error in test_metal vacuum\n");
+  // if (!test_metal(one, 200)) abort("error in test_metal vacuum\n");
 
   for (int s = 2; s < 5; s++)
-    if (!test_metal(targets, s, mydirname)) abort("error in test_metal targets\n");
-  // if (!test_metal(targets, 60, mydirname)) abort("error in test_metal targets\n");
+    if (!test_metal(targets, s)) abort("error in test_metal targets\n");
+  // if (!test_metal(targets, 60)) abort("error in test_metal targets\n");
 
   for (int s = 2; s < 5; s++)
-    if (!test_periodic(targets, s, mydirname)) abort("error in test_periodic targets\n");
-  // if (!test_periodic(one, 200, mydirname))
+    if (!test_periodic(targets, s)) abort("error in test_periodic targets\n");
+  // if (!test_periodic(one, 200))
   //  abort("error in test_periodic targets\n");
 
   for (int s = 2; s < 4; s++)
-    if (!test_periodic_tm(one, s, mydirname)) abort("error in test_periodic_tm vacuum\n");
+    if (!test_periodic_tm(one, s)) abort("error in test_periodic_tm vacuum\n");
 
   return 0;
 }