diff --git a/WENO.F90 b/WENO.F90
index 712b8b773cbce21fb5f00be9b3fc365330225ebb..8013b2cbf01c89fa816dd332602dc813f5154ada 100644
--- a/WENO.F90
+++ b/WENO.F90
@@ -982,18 +982,18 @@ end subroutine
 ! prg = lp.parse_fortran(lp.c_preprocess(SOURCE), FILENAME)
 ! prg = lp.fix_parameters(prg, ndim=3, nvars=5, _remove=False)
 !
-! cfd = prg["compute_flux_derivatives"]
+! knl = prg["compute_flux_derivatives"]
 !
-! cfd = lp.assume(cfd, "nx > 0 and ny > 0 and nz > 0")
+! knl = lp.assume(knl, "nx > 0 and ny > 0 and nz > 0")
 !
-! cfd = lp.set_temporary_scope(cfd, "flux_derivatives_generalized",
+! knl = lp.set_temporary_scope(knl, "flux_derivatives_generalized",
 !         lp.AddressSpace.GLOBAL)
-! cfd = lp.set_temporary_scope(cfd, "generalized_fluxes",
+! knl = lp.set_temporary_scope(knl, "generalized_fluxes",
 !         lp.AddressSpace.GLOBAL)
-! cfd = lp.set_temporary_scope(cfd, "weno_flux_tmp",
+! knl = lp.set_temporary_scope(knl, "weno_flux_tmp",
 !         lp.AddressSpace.GLOBAL)
 !
-! prg = prg.with_kernel(cfd)
+! prg = prg.with_kernel(knl)
 !
 ! RESULT = prg
 !
diff --git a/build-and-test-py-project.sh b/build-and-test-py-project.sh
index 8bfbc06940aa87368d22708756076e22b9e615a8..22bf66974d1b41338b9a8eb336cb91142cf9a1ea 100644
--- a/build-and-test-py-project.sh
+++ b/build-and-test-py-project.sh
@@ -85,6 +85,6 @@ if test -f $REQUIREMENTS_TXT; then
   $PIP install -r $REQUIREMENTS_TXT
 fi
 
-${PY_EXE} -m pytest -rw --durations=10 --tb=native --junitxml=pytest.xml -rxsw --runslow test.py
+${PY_EXE} -m pytest -rw --durations=10 --tb=native --junitxml=pytest.xml -rxsw --runslow test
 
 # vim: foldmethod=marker
diff --git a/benchmark.py b/test/benchmark.py
similarity index 100%
rename from benchmark.py
rename to test/benchmark.py
diff --git a/conftest.py b/test/conftest.py
similarity index 100%
rename from conftest.py
rename to test/conftest.py
diff --git a/data_for_test.py b/test/data_for_test.py
similarity index 100%
rename from data_for_test.py
rename to test/data_for_test.py
diff --git a/ideal_gas.py b/test/ideal_gas.py
similarity index 100%
rename from ideal_gas.py
rename to test/ideal_gas.py
diff --git a/pytest.ini b/test/pytest.ini
similarity index 100%
rename from pytest.ini
rename to test/pytest.ini
diff --git a/reference_implementation.py b/test/reference_implementation.py
similarity index 100%
rename from reference_implementation.py
rename to test/reference_implementation.py
diff --git a/test/test_eigensystem.py b/test/test_eigensystem.py
new file mode 100644
index 0000000000000000000000000000000000000000..4b0a7a6b14af212778076242e3e0fbef81e7f2e7
--- /dev/null
+++ b/test/test_eigensystem.py
@@ -0,0 +1,117 @@
+import numpy as np
+import numpy.linalg as la  # noqa: F401
+import pyopencl as cl  # noqa: F401
+import pyopencl.array  # noqa
+import pyopencl.tools  # noqa
+import pyopencl.clrandom  # noqa
+import loopy as lp  # noqa
+
+import sys
+import logging
+
+import pytest
+from pyopencl.tools import (  # noqa
+        pytest_generate_tests_for_pyopencl
+        as pytest_generate_tests)
+
+import utilities as u
+from data_for_test import (  # noqa: F401
+        flux_test_data_fixture,
+        single_data as sd
+        )
+
+
+def test_pointwise_eigenvalues_ideal_gas(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
+    prg = u.get_weno_program_with_root_kernel("pointwise_eigenvalues")
+    queue = u.get_queue(ctx_factory)
+
+    lam_dev = u.empty_array_on_device(queue, data.nvars, 6)
+
+    prg(queue, nvars=data.nvars, d=data.direction,
+            states=data.states, lambda_pointwise=lam_dev)
+
+    u.compare_arrays(lam_dev.get(), data.lam_pointwise)
+
+
+def test_roe_uniform_grid_ideal_gas(ctx_factory, flux_test_data_fixture):
+    data = flux_test_data_fixture
+
+    def check_roe_identity(states, R, R_inv):
+        d_state = states[:, 1] - states[:, 0]
+        u.compare_arrays(R@(R_inv@d_state), d_state)
+
+    def check_roe_property(states, fluxes, R, R_inv, lam):
+        d_state = states[:, 1] - states[:, 0]
+        d_flux = fluxes[:, 1] - fluxes[:, 0]
+
+        temp = R_inv@d_state
+        temp = np.multiply(lam, temp)
+        u.compare_arrays(R@temp, d_flux)
+
+    prg = u.get_weno_program_with_root_kernel("roe_eigensystem")
+    queue = u.get_queue(ctx_factory)
+
+    R_dev = u.empty_array_on_device(queue, data.nvars, data.nvars)
+    R_inv_dev = u.empty_array_on_device(queue, data.nvars, data.nvars)
+    lam_dev = u.empty_array_on_device(queue, data.nvars)
+
+    prg(queue, nvars=data.nvars, ndim=data.ndim, d=data.direction,
+            states=data.state_pair, metrics_frozen=data.frozen_metrics,
+            R=R_dev, R_inv=R_inv_dev, lambda_roe=lam_dev)
+
+    R = R_dev.get()
+    R_inv = R_inv_dev.get()
+    lam = lam_dev.get()
+
+    check_roe_identity(data.state_pair, R, R_inv)
+    check_roe_property(data.state_pair, data.flux_pair, R, R_inv, lam)
+
+
+@pytest.mark.parametrize("lam_pointwise_str,lam_roe_str,lam_expected_str", [
+    ("1 2 3 4 5,2 4 6 8 10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
+    ("1 2 3 4 5,-2 -4 -6 -8 -10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
+    ("1 2 3 4 5,-2 -4 -6 -8 -10", "3 6 9 12 15", "3.3 6.6 9.9 13.2 16.5"),
+    ("1 2 3 4 5,2 4 6 8 10", "-3 -6 -9 -12 -15", "3.3 6.6 9.9 13.2 16.5"),
+    ("3 2 9 4 5,2 6 6 12 10", "-1 -4 -3 -8 -15", "3.3 6.6 9.9 13.2 16.5")
+    ])
+def test_lax_wavespeeds(
+        ctx_factory, lam_pointwise_str, lam_roe_str, lam_expected_str):
+    prg = u.get_weno_program_with_root_kernel("lax_wavespeeds")
+    queue = u.get_queue(ctx_factory)
+
+    nvars = 5
+
+    lam_pointwise = u.expand_to_6(u.transposed_array_from_string(lam_pointwise_str))
+    lam_roe = u.array_from_string(lam_roe_str)
+    lam_dev = u.empty_array_on_device(queue, nvars)
+
+    prg(queue, nvars=nvars, lambda_pointwise=lam_pointwise,
+            lambda_roe=lam_roe, wavespeeds=lam_dev)
+
+    lam_expected = u.array_from_string(lam_expected_str)
+    u.compare_arrays(lam_dev.get(), lam_expected)
+
+
+def test_matvec(ctx_factory):
+    prg = u.get_weno_program_with_root_kernel("mult_mat_vec")
+    queue = u.get_queue(ctx_factory)
+
+    a = u.random_array_on_device(queue, 10, 10)
+    b = u.random_array_on_device(queue, 10)
+
+    c = u.empty_array_on_device(queue, 10)
+
+    prg(queue, alpha=1.0, a=a, b=b, c=c)
+
+    u.compare_arrays(a.get()@b.get(), c.get())
+
+
+# This lets you run 'python test.py test_case(cl._csc)' without pytest.
+if __name__ == "__main__":
+    if len(sys.argv) > 1:
+        logging.basicConfig(level="INFO")
+        exec(sys.argv[1])
+    else:
+        pytest.main([__file__])
diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py
new file mode 100644
index 0000000000000000000000000000000000000000..4fb26d42deea0311a0cc89d5315149c97ff47d67
--- /dev/null
+++ b/test/test_flux_derivatives.py
@@ -0,0 +1,53 @@
+import numpy as np
+import numpy.linalg as la  # noqa: F401
+import pyopencl as cl  # noqa: F401
+import pyopencl.array  # noqa
+import pyopencl.tools  # noqa
+import pyopencl.clrandom  # noqa
+import loopy as lp  # noqa
+
+import sys
+import logging
+
+import pytest
+from pyopencl.tools import (  # noqa
+        pytest_generate_tests_for_pyopencl
+        as pytest_generate_tests)
+
+import utilities as u
+
+
+@pytest.mark.slow
+def test_compute_flux_derivatives(ctx_factory):
+    prg = u.get_weno_program()
+
+    queue = u.get_queue(ctx_factory)
+    prg = prg.copy(target=lp.PyOpenCLTarget(queue.device))
+
+    lp.auto_test_vs_ref(prg, ctx_factory(), warmup_rounds=1,
+            parameters=dict(ndim=3, nvars=5, nx=16, ny=16, nz=16))
+
+
+@pytest.mark.slow
+def test_compute_flux_derivatives_gpu(ctx_factory, write_code=False):
+    prg = u.get_weno_program()
+    prg = u.transform_weno_for_gpu(prg)
+
+    queue = u.get_queue(ctx_factory)
+    prg = prg.copy(target=lp.PyOpenCLTarget(queue.device))
+    prg = lp.set_options(prg, no_numpy=True)
+
+    if write_code:
+        u.write_target_device_code(prg)
+
+    lp.auto_test_vs_ref(prg, ctx_factory(), warmup_rounds=1,
+            parameters=dict(ndim=3, nvars=5, nx=16, ny=16, nz=16))
+
+
+# This lets you run 'python test.py test_case(cl._csc)' without pytest.
+if __name__ == "__main__":
+    if len(sys.argv) > 1:
+        logging.basicConfig(level="INFO")
+        exec(sys.argv[1])
+    else:
+        pytest.main([__file__])
diff --git a/test.py b/test/test_weno_flux.py
similarity index 59%
rename from test.py
rename to test/test_weno_flux.py
index da3abad7faeaf936b440bde78d3c1da10fbb2a84..76976cd4a19965cfb72528fa01d369383b5a1d9f 100644
--- a/test.py
+++ b/test/test_weno_flux.py
@@ -17,7 +17,7 @@ from pyopencl.tools import (  # noqa
 import utilities as u
 from data_for_test import (  # noqa: F401
         flux_test_data_fixture,
-        single_data as std  # "single_test_data", sorry
+        single_data as sd
         )
 
 
@@ -178,125 +178,6 @@ def test_flux_splitting_uniform_grid(ctx_factory, flux_test_data_fixture):
     u.compare_arrays(fluxes_neg_dev.get(), data.char_fluxes_neg)
 
 
-def test_pointwise_eigenvalues_ideal_gas(ctx_factory, flux_test_data_fixture):
-    data = flux_test_data_fixture
-
-    prg = u.get_weno_program_with_root_kernel("pointwise_eigenvalues")
-    queue = u.get_queue(ctx_factory)
-
-    lam_dev = u.empty_array_on_device(queue, data.nvars, 6)
-
-    prg(queue, nvars=data.nvars, d=data.direction,
-            states=data.states, lambda_pointwise=lam_dev)
-
-    u.compare_arrays(lam_dev.get(), data.lam_pointwise)
-
-
-def test_roe_uniform_grid_ideal_gas(ctx_factory, flux_test_data_fixture):
-    data = flux_test_data_fixture
-
-    def identity_matrix(n):
-        return np.identity(n).astype(np.float64).copy(order="F")
-
-    def check_roe_identity(states, R, R_inv):
-        d_state = states[:, 1] - states[:, 0]
-        u.compare_arrays(R@(R_inv@d_state), d_state)
-
-    def check_roe_property(states, fluxes, R, R_inv, lam):
-        d_state = states[:, 1] - states[:, 0]
-        d_flux = fluxes[:, 1] - fluxes[:, 0]
-
-        temp = R_inv@d_state
-        temp = np.multiply(lam, temp)
-        u.compare_arrays(R@temp, d_flux)
-
-    prg = u.get_weno_program_with_root_kernel("roe_eigensystem")
-    queue = u.get_queue(ctx_factory)
-
-    metrics_frozen = identity_matrix(data.ndim)
-
-    R_dev = u.empty_array_on_device(queue, data.nvars, data.nvars)
-    R_inv_dev = u.empty_array_on_device(queue, data.nvars, data.nvars)
-    lam_dev = u.empty_array_on_device(queue, data.nvars)
-
-    prg(queue, nvars=data.nvars, ndim=data.ndim, d=data.direction,
-            states=data.state_pair, metrics_frozen=metrics_frozen,
-            R=R_dev, R_inv=R_inv_dev, lambda_roe=lam_dev)
-
-    R = R_dev.get()
-    R_inv = R_inv_dev.get()
-    lam = lam_dev.get()
-
-    check_roe_identity(data.state_pair, R, R_inv)
-    check_roe_property(data.state_pair, data.flux_pair, R, R_inv, lam)
-
-
-@pytest.mark.parametrize("lam_pointwise_str,lam_roe_str,lam_expected_str", [
-    ("1 2 3 4 5,2 4 6 8 10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
-    ("1 2 3 4 5,-2 -4 -6 -8 -10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
-    ("1 2 3 4 5,-2 -4 -6 -8 -10", "3 6 9 12 15", "3.3 6.6 9.9 13.2 16.5"),
-    ("1 2 3 4 5,2 4 6 8 10", "-3 -6 -9 -12 -15", "3.3 6.6 9.9 13.2 16.5"),
-    ("3 2 9 4 5,2 6 6 12 10", "-1 -4 -3 -8 -15", "3.3 6.6 9.9 13.2 16.5")
-    ])
-def test_lax_wavespeeds(
-        ctx_factory, lam_pointwise_str, lam_roe_str, lam_expected_str):
-    prg = u.get_weno_program_with_root_kernel("lax_wavespeeds")
-    queue = u.get_queue(ctx_factory)
-
-    nvars = 5
-
-    lam_pointwise = u.expand_to_6(u.transposed_array_from_string(lam_pointwise_str))
-    lam_roe = u.array_from_string(lam_roe_str)
-    lam_dev = u.empty_array_on_device(queue, nvars)
-
-    prg(queue, nvars=nvars, lambda_pointwise=lam_pointwise,
-            lambda_roe=lam_roe, wavespeeds=lam_dev)
-
-    lam_expected = u.array_from_string(lam_expected_str)
-    u.compare_arrays(lam_dev.get(), lam_expected)
-
-
-def test_matvec(ctx_factory):
-    prg = u.get_weno_program_with_root_kernel("mult_mat_vec")
-    queue = u.get_queue(ctx_factory)
-
-    a = u.random_array_on_device(queue, 10, 10)
-    b = u.random_array_on_device(queue, 10)
-
-    c = u.empty_array_on_device(queue, 10)
-
-    prg(queue, alpha=1.0, a=a, b=b, c=c)
-
-    u.compare_arrays(a.get()@b.get(), c.get())
-
-
-@pytest.mark.slow
-def test_compute_flux_derivatives(ctx_factory):
-    prg = u.get_weno_program()
-
-    queue = u.get_queue(ctx_factory)
-    prg = prg.copy(target=lp.PyOpenCLTarget(queue.device))
-
-    lp.auto_test_vs_ref(prg, ctx_factory(), warmup_rounds=1,
-            parameters=dict(ndim=3, nvars=5, nx=16, ny=16, nz=16))
-
-
-@pytest.mark.slow
-def test_compute_flux_derivatives_gpu(ctx_factory, write_code=False):
-    prg = u.get_weno_program()
-    prg = u.transform_weno_for_gpu(prg)
-
-    queue = u.get_queue(ctx_factory)
-    prg = prg.copy(target=lp.PyOpenCLTarget(queue.device))
-    prg = lp.set_options(prg, no_numpy=True)
-
-    if write_code:
-        u.write_target_device_code(prg)
-
-    lp.auto_test_vs_ref(prg, ctx_factory(), warmup_rounds=1,
-            parameters=dict(ndim=3, nvars=5, nx=16, ny=16, nz=16))
-
-
 # This lets you run 'python test.py test_case(cl._csc)' without pytest.
 if __name__ == "__main__":
     if len(sys.argv) > 1:
diff --git a/utilities.py b/test/utilities.py
similarity index 84%
rename from utilities.py
rename to test/utilities.py
index 139500153ca15f809015cd43033f1d2fa4a1897c..35ea205de6bcec427579a775a3126d9108fe2fc2 100644
--- a/utilities.py
+++ b/test/utilities.py
@@ -7,6 +7,8 @@ import pyopencl.clrandom  # noqa
 import loopy as lp  # noqa
 from pytest import approx
 
+import os.path
+
 
 # {{{ arrays
 
@@ -91,9 +93,11 @@ def get_queue(ctx_factory):
 
 _WENO_PRG = {}
 
+project_root = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
+
 
 def parse_weno():
-    fn = "WENO.F90"
+    fn = os.path.join(project_root, "WENO.F90")
 
     with open(fn, "r") as infile:
         infile_content = infile.read()
@@ -133,26 +137,26 @@ def get_weno_program_with_root_kernel(knl):
 
 
 def transform_weno_for_gpu(prg, print_kernel=False):
-    cfd = prg["compute_flux_derivatives"]
+    knl = prg["compute_flux_derivatives"]
 
     for suffix in ["", "_1", "_2", "_3", "_4", "_5", "_6", "_7"]:
-        cfd = lp.split_iname(cfd, "i"+suffix, 16,
+        knl = lp.split_iname(knl, "i"+suffix, 16,
                 outer_tag="g.0", inner_tag="l.0")
-        cfd = lp.split_iname(cfd, "j"+suffix, 16,
+        knl = lp.split_iname(knl, "j"+suffix, 16,
                 outer_tag="g.1", inner_tag="l.1")
 
     for var_name in ["delta_xi", "delta_eta", "delta_zeta"]:
-        cfd = lp.assignment_to_subst(cfd, var_name)
+        knl = lp.assignment_to_subst(knl, var_name)
 
-    cfd = lp.add_barrier(cfd, "tag:to_generalized", "tag:flux_x_compute")
-    cfd = lp.add_barrier(cfd, "tag:flux_x_compute", "tag:flux_x_diff")
-    cfd = lp.add_barrier(cfd, "tag:flux_x_diff", "tag:flux_y_compute")
-    cfd = lp.add_barrier(cfd, "tag:flux_y_compute", "tag:flux_y_diff")
-    cfd = lp.add_barrier(cfd, "tag:flux_y_diff", "tag:flux_z_compute")
-    cfd = lp.add_barrier(cfd, "tag:flux_z_compute", "tag:flux_z_diff")
-    cfd = lp.add_barrier(cfd, "tag:flux_z_diff", "tag:from_generalized")
+    knl = lp.add_barrier(knl, "tag:to_generalized", "tag:flux_x_compute")
+    knl = lp.add_barrier(knl, "tag:flux_x_compute", "tag:flux_x_diff")
+    knl = lp.add_barrier(knl, "tag:flux_x_diff", "tag:flux_y_compute")
+    knl = lp.add_barrier(knl, "tag:flux_y_compute", "tag:flux_y_diff")
+    knl = lp.add_barrier(knl, "tag:flux_y_diff", "tag:flux_z_compute")
+    knl = lp.add_barrier(knl, "tag:flux_z_compute", "tag:flux_z_diff")
+    knl = lp.add_barrier(knl, "tag:flux_z_diff", "tag:from_generalized")
 
-    prg = prg.with_kernel(cfd)
+    prg = prg.with_kernel(knl)
 
     # FIXME: These should work, but don't
     # FIXME: Undo the hand-inlining in WENO.F90