diff --git a/.github/workflows/build-ci.yml b/.github/workflows/build-ci.yml
index 246e295de..fb2a65e0f 100644
--- a/.github/workflows/build-ci.yml
+++ b/.github/workflows/build-ci.yml
@@ -140,7 +140,7 @@ jobs:
         echo "MEEP_VERSION=${MEEP_VERSION}" >> $GITHUB_ENV
 
     - name: Run configure
-      if: ${{ !(matrix.enable-mpi == false && matrix.python-version == 3.6) }}
+      if: ${{ !(matrix.enable-mpi == false && matrix.python-version == 3.6) && !(matrix.enable-mpi == true && matrix.python-version == 3.9) }}
       run: |
         mkdir -p build &&
         pushd build &&
@@ -151,6 +151,14 @@ jobs:
       if: ${{ matrix.enable-mpi == false && matrix.python-version == 3.6 }}
       run: ./configure --enable-maintainer-mode --prefix=${HOME}/local --with-libctl=${HOME}/local/share/libctl ${MPICONF}
 
+    - name: Run configure with single-precision floating point
+      if: ${{ matrix.enable-mpi == true && matrix.python-version == 3.9 }}
+      run: |
+        mkdir -p build &&
+        pushd build &&
+        ../configure --enable-maintainer-mode --prefix=${HOME}/local --with-libctl=${HOME}/local/share/libctl ${MPICONF} --enable-single &&
+        popd
+
     - name: Run make
       if: ${{ matrix.enable-mpi == false && matrix.python-version == 3.6 }}
       run: make ${MKCHECKFLAGS}
diff --git a/python/tests/test_3rd_harm_1d.py b/python/tests/test_3rd_harm_1d.py
index ce920564f..70d439d63 100644
--- a/python/tests/test_3rd_harm_1d.py
+++ b/python/tests/test_3rd_harm_1d.py
@@ -2,10 +2,9 @@
 
 import unittest
 import meep as mp
-import numpy as np
+from utils import VectorComparisonMixin
 
-
-class Test3rdHarm1d(unittest.TestCase):
+class Test3rdHarm1d(VectorComparisonMixin, unittest.TestCase):
 
     def setUp(self):
         self.sz = 100
@@ -53,7 +52,8 @@ def test_3rd_harm_1d(self):
 
         harmonics = [self.k, self.amp, mp.get_fluxes(self.trans1)[0], mp.get_fluxes(self.trans3)[0]]
 
-        np.testing.assert_allclose(expected_harmonics, harmonics)
+        tol = 3e-5 if mp.is_single_precision() else 1e-7
+        self.assertVectorsClose(expected_harmonics, harmonics, epsilon=tol)
 
 
 if __name__ == '__main__':
diff --git a/python/tests/test_adjoint_jax.py b/python/tests/test_adjoint_jax.py
index 2e1b7a92f..6a6c4ab49 100644
--- a/python/tests/test_adjoint_jax.py
+++ b/python/tests/test_adjoint_jax.py
@@ -16,7 +16,7 @@
 _FD_STEP = 1e-4
 
 # The tolerance for the adjoint and finite difference gradient comparison
-_TOL = 2e-2
+_TOL = 0.1 if mp.is_single_precision() else 0.02
 
 mp.verbosity(0)
 
diff --git a/python/tests/test_adjoint_solver.py b/python/tests/test_adjoint_solver.py
index 9fa5c11e1..e1136776b 100644
--- a/python/tests/test_adjoint_solver.py
+++ b/python/tests/test_adjoint_solver.py
@@ -8,6 +8,7 @@
 from autograd import tensor_jacobian_product
 import unittest
 from enum import Enum
+from utils import VectorComparisonMixin
 
 MonitorObject = Enum('MonitorObject', 'EIGENMODE DFT')
 
@@ -170,7 +171,7 @@ def mapping(x,filter_radius,eta,beta):
     return projected_field.flatten()
 
 
-class TestAdjointSolver(unittest.TestCase):
+class TestAdjointSolver(VectorComparisonMixin, unittest.TestCase):
 
     def test_adjoint_solver_DFT_fields(self):
         print("*** TESTING DFT ADJOINT FEATURES ***")
@@ -184,7 +185,7 @@ def test_adjoint_solver_DFT_fields(self):
             
             ## compare objective results
             print("|Ez|^2 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
-            np.testing.assert_array_almost_equal(adjsol_obj,S12_unperturbed,decimal=3)
+            self.assertVectorsClose(adjsol_obj,S12_unperturbed,epsilon=1e-3)
 
             ## compute perturbed S12
             S12_perturbed = forward_simulation(p+dp, MonitorObject.DFT, frequencies)
@@ -195,7 +196,8 @@ def test_adjoint_solver_DFT_fields(self):
             adj_scale = (dp[None,:]@adjsol_grad).flatten()
             fd_grad = S12_perturbed-S12_unperturbed
             print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
-            np.testing.assert_array_almost_equal(adj_scale,fd_grad,decimal=5)
+            tol = 0.3 if mp.is_single_precision() else 0.01
+            self.assertVectorsClose(adj_scale,fd_grad,epsilon=tol)
 
 
     def test_adjoint_solver_eigenmode(self):
@@ -207,21 +209,22 @@ def test_adjoint_solver_eigenmode(self):
 
             ## compute unperturbed S12
             S12_unperturbed = forward_simulation(p, MonitorObject.EIGENMODE, frequencies)
-            
+
             ## compare objective results
             print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
-            np.testing.assert_array_almost_equal(adjsol_obj,S12_unperturbed,decimal=3)
+            self.assertVectorsClose(adjsol_obj,S12_unperturbed,epsilon=1e-3)
 
             ## compute perturbed S12
             S12_perturbed = forward_simulation(p+dp, MonitorObject.EIGENMODE, frequencies)
-            
+
             ## compare gradients
             if adjsol_grad.ndim < 2:
                 adjsol_grad = np.expand_dims(adjsol_grad,axis=1)
             adj_scale = (dp[None,:]@adjsol_grad).flatten()
             fd_grad = S12_perturbed-S12_unperturbed
             print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
-            np.testing.assert_array_almost_equal(adj_scale,fd_grad,decimal=5)
+            tol = 0.04 if mp.is_single_precision() else 0.03
+            self.assertVectorsClose(adj_scale,fd_grad,epsilon=tol)
 
 
     def test_gradient_backpropagation(self):
@@ -250,7 +253,8 @@ def test_gradient_backpropagation(self):
 
             ## compare objective results
             print("S12 -- adjoint solver: {}, traditional simulation: {}".format(adjsol_obj,S12_unperturbed))
-            np.testing.assert_array_almost_equal(adjsol_obj,S12_unperturbed,decimal=3)
+            tol = 1e-2 if mp.is_single_precision() else 1e-3
+            self.assertVectorsClose(adjsol_obj,S12_unperturbed,epsilon=tol)
 
             ## compute perturbed S12
             S12_perturbed = forward_simulation(mapping(p+dp,filter_radius,eta,beta), MonitorObject.EIGENMODE,frequencies)
@@ -260,7 +264,8 @@ def test_gradient_backpropagation(self):
             adj_scale = (dp[None,:]@bp_adjsol_grad).flatten()
             fd_grad = S12_perturbed-S12_unperturbed
             print("Directional derivative -- adjoint solver: {}, FD: {}".format(adj_scale,fd_grad))
-            np.testing.assert_array_almost_equal(adj_scale,fd_grad,decimal=5)
+            tol = 0.04 if mp.is_single_precision() else 0.03
+            self.assertVectorsClose(adj_scale,fd_grad,epsilon=tol)
 
 
 if __name__ == '__main__':
diff --git a/python/tests/test_array_metadata.py b/python/tests/test_array_metadata.py
index 83051292b..4b433e961 100644
--- a/python/tests/test_array_metadata.py
+++ b/python/tests/test_array_metadata.py
@@ -1,5 +1,5 @@
-import unittest
 import meep as mp
+import unittest
 import numpy as np
 
 class TestArrayMetadata(unittest.TestCase):
@@ -80,7 +80,11 @@ def vec_func(r):
         vec_func_sum = np.sum(W*(xm**2 + 2*ym**2))
         pulse_modal_volume = np.sum(W*EpsE2)/np.max(EpsE2) * vec_func_sum
 
-        self.assertAlmostEqual(cw_modal_volume/pulse_modal_volume, 1.00, places=2)
+        if ((mp.count_processors() % 2 == 0) and mp.is_single_precision()):
+            ref_val = 0.94
+        else:
+            ref_val = 1.00
+        self.assertAlmostEqual(cw_modal_volume/pulse_modal_volume, ref_val, places=2)
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/python/tests/test_cavity_arrayslice.py b/python/tests/test_cavity_arrayslice.py
index 4834b614d..a8c9b5271 100644
--- a/python/tests/test_cavity_arrayslice.py
+++ b/python/tests/test_cavity_arrayslice.py
@@ -1,13 +1,11 @@
-from __future__ import division
-
+import meep as mp
+from utils import VectorComparisonMixin
 import os
 import unittest
-
-import meep as mp
 import numpy as np
 
 
-class TestCavityArraySlice(unittest.TestCase):
+class TestCavityArraySlice(VectorComparisonMixin, unittest.TestCase):
 
     data_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), 'data'))
     expected_1d = np.load(os.path.join(data_dir, 'cavity_arrayslice_1d.npy'))
@@ -60,27 +58,31 @@ def test_1d_slice(self):
         self.sim.run(until_after_sources=0)
         vol = mp.Volume(center=self.center_1d, size=self.size_1d)
         hl_slice1d = self.sim.get_array(mp.Hz, vol)
-        np.testing.assert_allclose(self.expected_1d, hl_slice1d)
+        tol = 1e-5 if mp.is_single_precision() else 1e-8
+        self.assertVectorsClose(self.expected_1d, hl_slice1d, epsilon=tol)
 
     def test_2d_slice(self):
         self.sim.run(until_after_sources=0)
         vol = mp.Volume(center=self.center_2d, size=self.size_2d)
         hl_slice2d = self.sim.get_array(mp.Hz, vol)
-        np.testing.assert_allclose(self.expected_2d, hl_slice2d)
+        tol = 1e-5 if mp.is_single_precision() else 1e-8
+        self.assertVectorsClose(self.expected_2d, hl_slice2d, epsilon=tol)
 
     def test_1d_slice_user_array(self):
         self.sim.run(until_after_sources=0)
         arr = np.zeros(126, dtype=np.float64)
         vol = mp.Volume(center=self.center_1d, size=self.size_1d)
         self.sim.get_array(mp.Hz, vol, arr=arr)
-        np.testing.assert_allclose(self.expected_1d, arr)
+        tol = 1e-5 if mp.is_single_precision() else 1e-8
+        self.assertVectorsClose(self.expected_1d, arr, epsilon=tol)
 
     def test_2d_slice_user_array(self):
         self.sim.run(until_after_sources=0)
         arr = np.zeros((126, 38), dtype=np.float64)
         vol = mp.Volume(center=self.center_2d, size=self.size_2d)
         self.sim.get_array(mp.Hz, vol, arr=arr)
-        np.testing.assert_allclose(self.expected_2d, arr)
+        tol = 1e-5 if mp.is_single_precision() else 1e-8
+        self.assertVectorsClose(self.expected_2d, arr, epsilon=tol)
 
     def test_illegal_user_array(self):
         self.sim.run(until_after_sources=0)
diff --git a/python/tests/test_cavity_farfield.py b/python/tests/test_cavity_farfield.py
index cd00d6e98..05c1e2214 100644
--- a/python/tests/test_cavity_farfield.py
+++ b/python/tests/test_cavity_farfield.py
@@ -1,11 +1,11 @@
+import meep as mp
+from utils import VectorComparisonMixin
 import os
 import unittest
 import h5py
-import numpy as np
-import meep as mp
 
 
-class TestCavityFarfield(unittest.TestCase):
+class TestCavityFarfield(VectorComparisonMixin, unittest.TestCase):
 
     data_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), 'data'))
 
@@ -50,9 +50,13 @@ def run_test(self, nfreqs):
 
         nearfield = sim.add_near2far(
             fcen, 0.1, nfreqs,
-            mp.Near2FarRegion(mp.Vector3(0, 0.5 * w + d1), size=mp.Vector3(2 * dpml - sx)),
-            mp.Near2FarRegion(mp.Vector3(-0.5 * sx + dpml, 0.5 * w + 0.5 * d1), size=mp.Vector3(0, d1), weight=-1.0),
-            mp.Near2FarRegion(mp.Vector3(0.5 * sx - dpml, 0.5 * w + 0.5 * d1), size=mp.Vector3(0, d1))
+            mp.Near2FarRegion(mp.Vector3(0, 0.5 * w + d1),
+                              size=mp.Vector3(2 * dpml - sx)),
+            mp.Near2FarRegion(mp.Vector3(-0.5 * sx + dpml, 0.5 * w + 0.5 * d1),
+                              size=mp.Vector3(0, d1),
+                              weight=-1.0),
+            mp.Near2FarRegion(mp.Vector3(0.5 * sx - dpml, 0.5 * w + 0.5 * d1),
+                              size=mp.Vector3(0, d1))
         )
         sim.run(until=200)
         d2 = 20
@@ -71,12 +75,14 @@ def run_test(self, nfreqs):
             ref_hy = mp.complexarray(f['hy.r'][()], f['hy.i'][()])
             ref_hz = mp.complexarray(f['hz.r'][()], f['hz.i'][()])
 
-            np.testing.assert_allclose(ref_ex, result['Ex'])
-            np.testing.assert_allclose(ref_ey, result['Ey'])
-            np.testing.assert_allclose(ref_ez, result['Ez'])
-            np.testing.assert_allclose(ref_hx, result['Hx'])
-            np.testing.assert_allclose(ref_hy, result['Hy'])
-            np.testing.assert_allclose(ref_hz, result['Hz'])
+            tol = 1e-5 if mp.is_single_precision() else 1e-7
+            self.assertVectorsClose(ref_ex, result['Ex'], epsilon=tol)
+            self.assertVectorsClose(ref_ey, result['Ey'], epsilon=tol)
+            self.assertVectorsClose(ref_ez, result['Ez'], epsilon=tol)
+            self.assertVectorsClose(ref_hx, result['Hx'], epsilon=tol)
+            self.assertVectorsClose(ref_hy, result['Hy'], epsilon=tol)
+            self.assertVectorsClose(ref_hz, result['Hz'], epsilon=tol)
+
 
     def test_cavity_farfield(self):
         self.run_test(nfreqs=1)
diff --git a/python/tests/test_chunks.py b/python/tests/test_chunks.py
index 36bde0fc0..886f57d4e 100644
--- a/python/tests/test_chunks.py
+++ b/python/tests/test_chunks.py
@@ -1,5 +1,6 @@
-import unittest
 import meep as mp
+import unittest
+
 
 class TestChunks(unittest.TestCase):
 
@@ -44,8 +45,12 @@ def test_chunks(self):
         sim.save_flux('tot_flux', tot_flux)
         sim1 = sim
 
-        geometry = [mp.Block(center=mp.Vector3(), size=mp.Vector3(sxy, sxy, mp.inf), material=mp.Medium(index=3.5)),
-                    mp.Block(center=mp.Vector3(), size=mp.Vector3(sxy-2*dpml, sxy-2*dpml, mp.inf), material=mp.air)]
+        geometry = [mp.Block(center=mp.Vector3(),
+                             size=mp.Vector3(sxy, sxy, mp.inf),
+                             material=mp.Medium(index=3.5)),
+                    mp.Block(center=mp.Vector3(),
+                             size=mp.Vector3(sxy-2*dpml, sxy-2*dpml, mp.inf),
+                             material=mp.air)]
 
         sim = mp.Simulation(cell_size=cell,
                             geometry=geometry,
@@ -62,7 +67,8 @@ def test_chunks(self):
 
         sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(), 1e-5))
 
-        self.assertAlmostEqual(86.90826609300862, mp.get_fluxes(tot_flux)[0])
+        places = 3 if mp.is_single_precision() else 7
+        self.assertAlmostEqual(86.90826609300862, mp.get_fluxes(tot_flux)[0], places=places)
 
 
 if __name__ == '__main__':
diff --git a/python/tests/test_dispersive_eigenmode.py b/python/tests/test_dispersive_eigenmode.py
index e12f318c7..0a56b815a 100644
--- a/python/tests/test_dispersive_eigenmode.py
+++ b/python/tests/test_dispersive_eigenmode.py
@@ -6,16 +6,14 @@
 #  * check magnetic profiles
 #  * once imaginary component is supported, check that
 
-from __future__ import division
-
-import unittest
 import meep as mp
+from utils import VectorComparisonMixin
+import unittest
 import numpy as np
-from meep import mpb
 import h5py
 import os
 
-class TestDispersiveEigenmode(unittest.TestCase):
+class TestDispersiveEigenmode(VectorComparisonMixin, unittest.TestCase):
     # ----------------------------------------- #
     # ----------- Helper Functions ------------ #
     # ----------------------------------------- #
@@ -35,8 +33,9 @@ def call_chi1(self,material,frequency):
         n = np.real(np.sqrt(np.linalg.inv(chi1inv.astype(np.complex128))))
 
         n_actual = np.real(np.sqrt(material.epsilon(frequency).astype(np.complex128)))
-        
-        np.testing.assert_allclose(n,n_actual)
+
+        tol = 1e-6 if mp.is_single_precision() else 1e-8
+        self.assertVectorsClose(n, n_actual, epsilon=tol)
 
     @classmethod
     def setUpClass(cls):
diff --git a/python/tests/test_divide_mpi_processes.py b/python/tests/test_divide_mpi_processes.py
index 88b82d837..ba7addf3c 100644
--- a/python/tests/test_divide_mpi_processes.py
+++ b/python/tests/test_divide_mpi_processes.py
@@ -47,8 +47,9 @@ def test_divide_parallel_processes(self):
 
         self.assertEqual(fcens[0],1)
         self.assertEqual(fcens[1],0.5)
-        self.assertAlmostEqual(tot_fluxes[0], 9.8628728533)
-        self.assertAlmostEqual(tot_fluxes[1], 19.6537275387)
+        places = 4 if mp.is_single_precision() else 7
+        self.assertAlmostEqual(tot_fluxes[0], 9.8628728533, places=places)
+        self.assertAlmostEqual(tot_fluxes[1], 19.6537275387, places=places)
         
 if __name__ == '__main__':
     unittest.main()
diff --git a/python/tests/test_eigfreq.py b/python/tests/test_eigfreq.py
index 5127a2ab9..f73fe73ba 100644
--- a/python/tests/test_eigfreq.py
+++ b/python/tests/test_eigfreq.py
@@ -1,18 +1,14 @@
-from __future__ import division
-
-import unittest
 import meep as mp
-import cmath
-import math
-from time import time
+import unittest
 
 class TestEigfreq(unittest.TestCase):
 
+    @unittest.skipIf(mp.is_single_precision(), "double-precision floating point specific test")
     def test_eigfreq(self):
         w = 1.2           # width of waveguide
         r = 0.36          # radius of holes
         d = 1.4           # defect spacing (ordinary spacing = 1)
-        N = 3        # number of holes on either side of defect
+        N = 3             # number of holes on either side of defect
         sy = 6            # size of cell in y direction (perpendicular to wvg.)
         pad = 2           # padding between last hole and PML edge
         dpml = 1          # PML thickness
diff --git a/python/tests/test_field_functions.py b/python/tests/test_field_functions.py
index b5ef09814..33cdc41af 100644
--- a/python/tests/test_field_functions.py
+++ b/python/tests/test_field_functions.py
@@ -1,5 +1,6 @@
-import unittest
 import meep as mp
+import unittest
+
 
 def f(r, ex, hz, eps):
     return (r.x * r.norm() + ex) - (eps * hz)
@@ -95,7 +96,8 @@ def test_integrate2_field_function(self):
         sim.run(until_after_sources=10)
 
         res1 = sim.integrate2_field_function(fields2, [mp.Ez], [mp.Ez], f2)
-        self.assertAlmostEqual(res1, 0.17158099566244897)
+        places = 6 if mp.is_single_precision() else 7
+        self.assertAlmostEqual(res1, 0.17158099566244897, places=places)
 
     def test_max_abs_field_function(self):
         sim = self.init()
diff --git a/python/tests/test_get_point.py b/python/tests/test_get_point.py
index a98fbb2d0..154336d92 100644
--- a/python/tests/test_get_point.py
+++ b/python/tests/test_get_point.py
@@ -1,10 +1,9 @@
-from __future__ import division
-
-import unittest
 import meep as mp
+import unittest
 import numpy as np
 import math
 
+
 class TestGetPoint(unittest.TestCase):
 
     def test_get_point(self):
@@ -95,9 +94,14 @@ def sinusoid(p):
                     1.0370162714 ]
 
         x = np.linspace(-0.865692,2.692867,29)
+        places = 5 if mp.is_single_precision() else 10
         for j in range(x.size):
-            self.assertAlmostEqual(np.real(sim.get_field_point(mp.Ez, mp.Vector3(x[j],-0.394862))),ez_ref[j],places=10)
-            self.assertAlmostEqual(sim.get_epsilon_point(mp.Vector3(x[j],2.967158)),eps_ref[j],places=10)
+            self.assertAlmostEqual(np.real(sim.get_field_point(mp.Ez, mp.Vector3(x[j],-0.394862))),
+                                   ez_ref[j],
+                                   places=places)
+            self.assertAlmostEqual(sim.get_epsilon_point(mp.Vector3(x[j],2.967158)),
+                                   eps_ref[j],
+                                   places=places)
 
 if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
+    unittest.main()
diff --git a/python/tests/test_holey_wvg_bands.py b/python/tests/test_holey_wvg_bands.py
index 9f0c66e6d..6d8044a28 100644
--- a/python/tests/test_holey_wvg_bands.py
+++ b/python/tests/test_holey_wvg_bands.py
@@ -1,5 +1,5 @@
-import unittest
 import meep as mp
+import unittest
 
 
 class TestHoleyWvgBands(unittest.TestCase):
@@ -65,9 +65,10 @@ def test_fields_at_kx(self):
         ]
 
         self.assertTrue(h.modes)
+        places = 5 if mp.is_single_precision() else 7
         for (r, i), m in zip(expected, h.modes):
-            self.assertAlmostEqual(m.freq, r)
-            self.assertAlmostEqual(m.decay, i)
+            self.assertAlmostEqual(m.freq, r, places=places)
+            self.assertAlmostEqual(m.decay, i, places=places)
 
 
 if __name__ == '__main__':
diff --git a/python/tests/test_holey_wvg_cavity.py b/python/tests/test_holey_wvg_cavity.py
index b9b883999..1f7cd2fc0 100644
--- a/python/tests/test_holey_wvg_cavity.py
+++ b/python/tests/test_holey_wvg_cavity.py
@@ -1,8 +1,9 @@
-import unittest
 import meep as mp
-import numpy as np
+from utils import VectorComparisonMixin
+import unittest
 
-class TestHoleyWvgCavity(unittest.TestCase):
+
+class TestHoleyWvgCavity(VectorComparisonMixin, unittest.TestCase):
 
     def setUp(self):
         eps = 13
@@ -70,7 +71,8 @@ def test_resonant_modes(self):
         m = h.modes[0]
         res = [m.freq, m.decay, m.Q, abs(m.amp), m.amp.real, m.amp.imag]
 
-        np.testing.assert_allclose(expected, res)
+        tol = 1e-6 if mp.is_single_precision() else 1e-8
+        self.assertVectorsClose(expected, res, epsilon=tol)
 
     def test_transmission_spectrum(self):
         expected = [
@@ -141,8 +143,10 @@ def test_transmission_spectrum(self):
 
         res = zip(mp.get_flux_freqs(trans), mp.get_fluxes(trans))
 
+        tol = 1e-8 if mp.is_single_precision() else 1e-10
         for e, r in zip(expected, res):
-            np.testing.assert_allclose(e, r)
+            self.assertVectorsClose(e, r, epsilon=tol)
+
 
 
 if __name__ == '__main__':
diff --git a/python/tests/test_mode_coeffs.py b/python/tests/test_mode_coeffs.py
index d3b650b1e..d60be05a4 100644
--- a/python/tests/test_mode_coeffs.py
+++ b/python/tests/test_mode_coeffs.py
@@ -1,8 +1,6 @@
-from __future__ import division
-
+import meep as mp
 import unittest
 import numpy as np
-import meep as mp
 
 
 class TestModeCoeffs(unittest.TestCase):
@@ -117,8 +115,10 @@ def test_modes(self):
         eval_point = mp.Vector3(0.7, -0.2, 0.3)
         ex_at_eval_point = emdata.amplitude(eval_point, mp.Ex)
         hz_at_eval_point = emdata.amplitude(eval_point, mp.Hz)
-        self.assertAlmostEqual(ex_at_eval_point, 0.4887779638178009+0.48424014532428294j)
-        self.assertAlmostEqual(hz_at_eval_point, 3.4249236584603495-3.455974863884166j)
+
+        places = 5 if mp.is_single_precision() else 7
+        self.assertAlmostEqual(ex_at_eval_point, 0.4887779638178009+0.484240145324284j, places=places)
+        self.assertAlmostEqual(hz_at_eval_point, 3.4249236584603495-3.455974863884166j, places=places)
 
     def test_kpoint_func(self):
 
diff --git a/python/tests/test_multilevel_atom.py b/python/tests/test_multilevel_atom.py
index 23b6d0ad4..33aa17b39 100644
--- a/python/tests/test_multilevel_atom.py
+++ b/python/tests/test_multilevel_atom.py
@@ -1,13 +1,11 @@
-from __future__ import division
-
+import meep as mp
 import math
 import unittest
 
-import meep as mp
-
 
 class TestMultiLevelAtom(unittest.TestCase):
 
+    @unittest.skipIf(mp.is_single_precision(), "double-precision floating point specific test")
     def test_multilevel_atom(self):
         resolution = 40
         ncav = 1.5
diff --git a/python/tests/test_pw_source.py b/python/tests/test_pw_source.py
index eb6a6f274..ddd1f4b15 100644
--- a/python/tests/test_pw_source.py
+++ b/python/tests/test_pw_source.py
@@ -1,10 +1,8 @@
-from __future__ import division
-
+import meep as mp
 import cmath
 import math
 import unittest
 
-import meep as mp
 
 class TestPwSource(unittest.TestCase):
 
@@ -71,7 +69,9 @@ def test_pw_source(self):
         pt1 = self.sim.get_field_point(mp.Ez, v1)
         pt2 = self.sim.get_field_point(mp.Ez, v2)
 
-        self.assertAlmostEqual(pt1 / pt2, 27.557668029008262)
+        places = 3 if mp.is_single_precision() else 7
+        self.assertAlmostEqual(pt1 / pt2, 27.557668029008262, places=places)
+
         self.assertAlmostEqual(cmath.exp(1j * self.k.dot(v1 - v2)), 0.7654030066070924 - 0.6435512702783076j)
 
 if __name__ == '__main__':
diff --git a/python/tests/test_refl_angular.py b/python/tests/test_refl_angular.py
index 240fef879..15ad269e6 100644
--- a/python/tests/test_refl_angular.py
+++ b/python/tests/test_refl_angular.py
@@ -1,12 +1,11 @@
-from __future__ import division
-
+import meep as mp
+from utils import VectorComparisonMixin
 import math
 import unittest
 import numpy as np
-import meep as mp
 
 
-class TestReflAngular(unittest.TestCase):
+class TestReflAngular(VectorComparisonMixin, unittest.TestCase):
 
     def test_refl_angular(self):
         resolution = 100
@@ -121,7 +120,8 @@ def test_refl_angular(self):
             (2.4999999999999987, -4.3484671364929225e-6),
         ]
 
-        np.testing.assert_allclose(expected, list(zip(freqs, refl_flux)), rtol=1e-6)
+        tol = 1e-7 if mp.is_single_precision() else 1e-8
+        self.assertVectorsClose(expected, list(zip(freqs, refl_flux)), epsilon=tol)
 
 
 if __name__ == '__main__':
diff --git a/python/tests/test_ring.py b/python/tests/test_ring.py
index 6d6c57284..f9f9817ff 100644
--- a/python/tests/test_ring.py
+++ b/python/tests/test_ring.py
@@ -65,8 +65,10 @@ def test_harminv(self):
         fp = self.sim.get_field_point(mp.Ez, v)
         ep = self.sim.get_epsilon_point(v)
 
-        self.assertAlmostEqual(ep, 11.559999999999999)
-        self.assertAlmostEqual(fp, -0.08185972142450348)
+        places = 5 if mp.is_single_precision() else 7
+        self.assertAlmostEqual(ep, 11.559999999999999, places=places)
+        self.assertAlmostEqual(fp, -0.08185972142450348, places=places)
+
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/python/tests/test_ring_cyl.py b/python/tests/test_ring_cyl.py
index 34691727d..862801c4d 100644
--- a/python/tests/test_ring_cyl.py
+++ b/python/tests/test_ring_cyl.py
@@ -2,10 +2,9 @@
 
 import unittest
 import meep as mp
-import numpy as np
+from utils import VectorComparisonMixin
 
-
-class TestRingCyl(unittest.TestCase):
+class TestRingCyl(VectorComparisonMixin, unittest.TestCase):
 
     def setUp(self):
         n = 3.4
@@ -66,7 +65,8 @@ def test_ring_cyl(self):
         m = h.modes[0]
         res = [m.freq, m.decay, m.Q, abs(m.amp), m.amp.real, m.amp.imag]
 
-        np.testing.assert_allclose(expected, res)
+        tol = 1e-6 if mp.is_single_precision() else 1e-7
+        self.assertVectorsClose(expected, res, epsilon=tol)
 
 
 if __name__ == '__main__':
diff --git a/python/tests/test_simulation.py b/python/tests/test_simulation.py
index 435e266af..ed751082e 100644
--- a/python/tests/test_simulation.py
+++ b/python/tests/test_simulation.py
@@ -214,14 +214,15 @@ def test_epsilon_input_file(self):
         sim.run(until=200)
         fp = sim.get_field_point(mp.Ez, mp.Vector3(x=1))
 
-        self.assertAlmostEqual(fp, -0.002989654055823199 + 0j)
+        places = 6 if mp.is_single_precision() else 7
+        self.assertAlmostEqual(fp, -0.002989654055823199, places=places)
 
         # Test unicode file name for Python 2
         if sys.version_info[0] == 2:
             sim = self.init_simple_simulation(epsilon_input_file=unicode(eps_input_path))
             sim.run(until=200)
             fp = sim.get_field_point(mp.Ez, mp.Vector3(x=1))
-            self.assertAlmostEqual(fp, -0.002989654055823199 + 0j)
+            self.assertAlmostEqual(fp, -0.002989654055823199)
 
     def test_numpy_epsilon(self):
         sim = self.init_simple_simulation()
@@ -234,7 +235,9 @@ def test_numpy_epsilon(self):
 
         sim.run(until=200)
         fp = sim.get_field_point(mp.Ez, mp.Vector3(x=1))
-        self.assertAlmostEqual(fp, -0.002989654055823199 + 0j)
+
+        places = 6 if mp.is_single_precision() else 7
+        self.assertAlmostEqual(fp, -0.002989654055823199, places=places)
 
     def test_set_materials(self):
 
@@ -507,7 +510,10 @@ def check_iterable(one, two, three, four):
         with self.assertRaises(TypeError):
             mp.vec(1, [2, 3])
 
+    @unittest.skipIf(mp.is_single_precision(), "double-precision floating point specific test")
     def test_epsilon_warning(self):
+        ## fields blow up using dispersive material
+        ## when compiled using single precision
 
         with warnings.catch_warnings(record=True) as w:
             warnings.simplefilter("always")
diff --git a/python/tests/test_user_defined_material.py b/python/tests/test_user_defined_material.py
index 2312dd428..5a363e4dc 100644
--- a/python/tests/test_user_defined_material.py
+++ b/python/tests/test_user_defined_material.py
@@ -59,7 +59,8 @@ def test_user_material_func(self):
         sim.run(until=200)
         fp = sim.get_field_point(mp.Ez, mp.Vector3(x=1))
 
-        self.assertAlmostEqual(fp, 4.816403627871773e-4 + 0j)
+        places = 3 if mp.is_single_precision() else 7
+        self.assertAlmostEqual(fp, 4.816403627871773e-4, places=places)
 
     def test_epsilon_func(self):
         sim = mp.Simulation(cell_size=self.cell,
@@ -72,7 +73,7 @@ def test_epsilon_func(self):
         sim.run(until=100)
         fp = sim.get_field_point(mp.Ez, mp.Vector3(x=1))
 
-        self.assertAlmostEqual(fp, -7.895783750440999e-4 + 0j)
+        self.assertAlmostEqual(fp, -7.895783750440999e-4)
 
     def test_geometric_obj_with_user_material(self):
         geometry = [mp.Cylinder(5, material=my_material_func)]
@@ -86,7 +87,8 @@ def test_geometric_obj_with_user_material(self):
         sim.run(until=200)
         fp = sim.get_field_point(mp.Ez, mp.Vector3(x=1))
 
-        self.assertAlmostEqual(fp, 4.816403627871773e-4 + 0j)
+        places = 3 if mp.is_single_precision() else 7
+        self.assertAlmostEqual(fp, 4.816403627871773e-4, places=places)
 
     def test_geometric_obj_with_epsilon_func(self):
         geometry = [mp.Cylinder(5, epsilon_func=my_epsilon_func)]
@@ -100,7 +102,7 @@ def test_geometric_obj_with_epsilon_func(self):
         sim.run(until=100)
         fp = sim.get_field_point(mp.Ez, mp.Vector3(x=1))
 
-        self.assertAlmostEqual(fp, -7.895783750440999e-4 + 0j)
+        self.assertAlmostEqual(fp, -7.895783750440999e-4)
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/python/tests/test_wvg_src.py b/python/tests/test_wvg_src.py
index c8ee12421..5b86a6afc 100644
--- a/python/tests/test_wvg_src.py
+++ b/python/tests/test_wvg_src.py
@@ -43,7 +43,8 @@ def test_wvg_src(self):
         flux2 = self.sim.flux_in_box(mp.X, mp.Volume(center=mp.Vector3(6.0), size=mp.Vector3(1.8, 6)))
 
         self.assertAlmostEqual(flux1, -1.775216564842667e-03)
-        self.assertAlmostEqual(flux2, 7.215785537102116e+00)
+        places = 6 if mp.is_single_precision() else 7
+        self.assertAlmostEqual(flux2, 7.215785537102116, places=places)
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/tests/aniso_disp.cpp b/tests/aniso_disp.cpp
index 4c57cad83..0b4780f10 100644
--- a/tests/aniso_disp.cpp
+++ b/tests/aniso_disp.cpp
@@ -139,8 +139,9 @@ int main(int argc, char **argv) {
     master_printf("err. real: %g\n", fabs(freqs_re[i0] - 0.41562) / 0.41562);
     master_printf("err. imag: %g\n", fabs(freqs_im[i0] + 4.8297e-07) / 4.8297e-7);
 
+    double tol = sizeof(realnum) == sizeof(float) ? 0.23 : 0.20;
     ok = fabs(freqs_re[i0] - 0.41562) / 0.41562 < 1e-4 &&
-         fabs(freqs_im[i0] + 4.8297e-07) / 4.8297e-7 < 0.2;
+         fabs(freqs_im[i0] + 4.8297e-07) / 4.8297e-7 < tol;
   }
   end_divide_parallel();
   return !and_to_all(ok);
diff --git a/tests/array-slice-ll.cpp b/tests/array-slice-ll.cpp
index 943186b2f..df7c55855 100644
--- a/tests/array-slice-ll.cpp
+++ b/tests/array-slice-ll.cpp
@@ -39,7 +39,8 @@ std::complex<double> default_field_function(const std::complex<double> *fields,
 /***************************************************************/
 /***************************************************************/
 /***************************************************************/
-#define RELTOL 1.0e-6
+const double RELTOL = sizeof(realnum) == sizeof(float) ? 1.0e-4: 1.0e-6;
+
 double Compare(realnum *d1, realnum *d2, int N, const char *Name) {
   double Norm1 = 0.0, Norm2 = 0.0, NormDelta = 0.0;
   for (int n = 0; n < N; n++) {
diff --git a/tests/cyl-ellipsoid-ll.cpp b/tests/cyl-ellipsoid-ll.cpp
index 7c7245f50..275a53a0d 100644
--- a/tests/cyl-ellipsoid-ll.cpp
+++ b/tests/cyl-ellipsoid-ll.cpp
@@ -28,7 +28,7 @@ typedef std::complex<double> cdouble;
 /***************************************************************/
 bool compare_hdf5_datasets(const char *file1, const char *name1, const char *file2,
                            const char *name2, int expected_rank = 2, double rel_tol = 1.0e-4,
-                           double abs_tol = 1.0e-8) {
+                           double abs_tol = sizeof(realnum) == sizeof(float) ? 1.0e-6 : 1.0e-8) {
   h5file f1(file1, h5file::READONLY, false);
   int rank1;
   std::vector<size_t> dims1(expected_rank);
diff --git a/tests/h5test.cpp b/tests/h5test.cpp
index bdedef0dc..75842fe81 100644
--- a/tests/h5test.cpp
+++ b/tests/h5test.cpp
@@ -46,7 +46,7 @@ symmetry make_rotate4z(const grid_volume &gv) { return rotate4(Z, gv); }
 
 typedef symmetry (*symfunc)(const grid_volume &);
 
-const double tol = sizeof(realnum) == sizeof(float) ? 1e-4 : 1e-8;
+const double tol = sizeof(realnum) == sizeof(float) ? 2e-4 : 1e-8;
 double compare(double a, double b, const char *nam, size_t i0, size_t i1, size_t i2) {
   if (fabs(a - b) > tol * tol + fabs(b) * tol || b != b) {
     master_printf("%g vs. %g differs by\t%g\n", a, b, fabs(a - b));
diff --git a/tests/harmonics.cpp b/tests/harmonics.cpp
index ac3f848b8..66dc48eef 100644
--- a/tests/harmonics.cpp
+++ b/tests/harmonics.cpp
@@ -101,7 +101,11 @@ int main(int argc, char **argv) {
   double thresh = sizeof(realnum) == sizeof(float) ? 1e-4 : 1e-5;
   harmonics(freq, 0.27e-4, 1e-4, 1.0, a2, a3);
   if (different(a2, 9.80330e-07, thresh, "2nd harmonic mismatches known val")) return 1;
-  if (different(a3, 9.97747e-07, thresh, "3rd harmonic mismatches known val")) return 1;
+  if (sizeof(realnum) == sizeof(float)) {
+    if (different(a3, 9.99349e-07, thresh, "3rd harmonic mismatches known val")) return 1;
+  } else {
+    if (different(a3, 9.97747e-07, thresh, "3rd harmonic mismatches known val")) return 1;
+  }
 
   harmonics(freq, 0.54e-4, 2e-4, 1.0, a2_2, a3_2);
   master_printf("doubling chi2, chi3 = %g x 2nd harmonic, %g x 3rd\n", a2_2 / a2, a3_2 / a3);
@@ -121,7 +125,12 @@ int main(int argc, char **argv) {
   }
 
   harmonics(freq, 0.0, 1e-4, 1.0, a2_2, a3_2);
-  if (different(a3, a3_2, 1e-3, "chi2 has too big effect on 3rd harmonic")) return 1;
+  if (sizeof(realnum) == sizeof(float)) {
+    if (different(a3, a3_2, 0.0017, "chi2 has too big effect on 3rd harmonic")) return 1;
+  } else {
+    if (different(a3, a3_2, 0.001, "chi2 has too big effect on 3rd harmonic")) return 1;
+  }
+
   if (a2_2 / a2 > 1e-5) {
     master_printf("error: too much 2nd harmonic without chi3\n");
     return 1;
diff --git a/tests/near2far.cpp b/tests/near2far.cpp
index 252a224bd..0a232a38d 100644
--- a/tests/near2far.cpp
+++ b/tests/near2far.cpp
@@ -252,7 +252,7 @@ int main(int argc, char **argv) {
 
   const double a2d = argc > 1 ? atof(argv[1]) : 20, a3d = argc > 1 ? a2d : 10;
 
-  if (!check_cyl(5.0, 10.0, 20.0)) return 1;
+  if ((sizeof(realnum) == sizeof(double)) &&  !check_cyl(5.0, 10.0, 20.0)) return 1;
 
   // NOTE: see hack above -- we require sources to be odd in Z and even in X or vice-versa
   if (!check_2d_3d(D3, 4, a3d, Ez, Hx, false)) return 1;
diff --git a/tests/pml.cpp b/tests/pml.cpp
index 4a4833264..14ce007e5 100644
--- a/tests/pml.cpp
+++ b/tests/pml.cpp
@@ -333,7 +333,9 @@ int main(int argc, char **argv) {
   // if (check_pml2d(one, Ez, 0, false, 0.5)) meep::abort("not a pml in 2d TM + offdiag.");
   // if (check_pml2d(one, Hz, 0, false, 0.5)) meep::abort("not a pml in 2d TE + offdiag.");
   // if (check_pmlcyl(one)) meep::abort("not a pml in cylincrical co-ordinates.");
-  if (pml1d_scaling(one)) meep::abort("pml doesn't scale properly with length.");
+  if (sizeof(realnum) == sizeof(double)) {
+    if (pml1d_scaling(one)) meep::abort("pml doesn't scale properly with length.");
+  }
   if (pmlcyl_scaling(one, 0)) meep::abort("m=0 cylindrical pml doesn't scale properly with length.");
   if (pmlcyl_scaling(one, 1)) meep::abort("m=1 cylindrical pml doesn't scale properly with length.");
   if (pmlcyl_scaling(one, 2)) meep::abort("m=2 cylindrical pml doesn't scale properly with length.");
diff --git a/tests/symmetry.cpp b/tests/symmetry.cpp
index b52e990e5..5a9d1d775 100644
--- a/tests/symmetry.cpp
+++ b/tests/symmetry.cpp
@@ -977,7 +977,9 @@ int main(int argc, char **argv) {
 #endif
 
   if (!nonlinear_ex(vol1d(1.0, 30.0), one)) meep::abort("error in 1D nonlinear vacuum\n");
-  if (!nonlinear_ex(vol3d(1.0, 1.2, 0.8, 10.0), one)) meep::abort("error in 3D nonlinear vacuum\n");
+  if (sizeof(realnum) == sizeof(double)) {
+    if (!nonlinear_ex(vol3d(1.0, 1.2, 0.8, 10.0), one)) meep::abort("error in 3D nonlinear vacuum\n");
+  }
 
   // disable for parallel runs due to a bug in splitting cylindrical cell
   // with z-mirror symmetry into multiple chunks