diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 334a41d0fb0b56e6582d71e9242c09d9729fdb0e..ad62d9543f8d4bb1659ee7a3afe3547316edce39 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -26,6 +26,20 @@ Python 3.6 POCL:
   except:
   - tags
 
+Python 3.6 POCL Examples:
+  script:
+  - export PY_EXE=python3.6
+  - export PYOPENCL_TEST=portable
+  - export EXTRA_INSTALL="numpy mako pyvisfile matplotlib"
+  - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-py-project-and-run-examples.sh
+  - ". ./build-py-project-and-run-examples.sh"
+  tags:
+  - python3.6
+  - pocl
+  - large-node
+  except:
+  - tags
+
 Python 3.5 Conda:
   script:
   - export SUMPY_FORCE_SYMBOLIC_BACKEND=symengine
@@ -66,7 +80,7 @@ Documentation:
 Flake8:
   script:
   - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/prepare-and-run-flake8.sh
-  - ". ./prepare-and-run-flake8.sh pytential test"
+  - ". ./prepare-and-run-flake8.sh pytential test examples"
   tags:
   - python3.5
   except:
diff --git a/examples/fmm-error.py b/examples/fmm-error.py
index fea97c99fd54398e4e6b3dd6c2092c19e60a31f5..110ec66bc166b9efc35b83e72f11584f07c7d5ee 100644
--- a/examples/fmm-error.py
+++ b/examples/fmm-error.py
@@ -6,89 +6,93 @@ from meshmode.mesh.generation import (  # noqa
 from sumpy.visualization import FieldPlotter
 from sumpy.kernel import LaplaceKernel, HelmholtzKernel
 
-import faulthandler
-faulthandler.enable()
 
-import logging
-logging.basicConfig(level=logging.INFO)
+def main():
+    import logging
+    logging.basicConfig(level=logging.WARNING)  # INFO for more progress info
 
-cl_ctx = cl.create_some_context()
-queue = cl.CommandQueue(cl_ctx)
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
 
-target_order = 16
-qbx_order = 3
-nelements = 60
-mode_nr = 0
+    target_order = 16
+    qbx_order = 3
+    nelements = 60
+    mode_nr = 0
 
-k = 0
-if k:
-    kernel = HelmholtzKernel("k")
-else:
-    kernel = LaplaceKernel()
-#kernel = OneKernel()
+    k = 0
+    if k:
+        kernel = HelmholtzKernel(2)
+    else:
+        kernel = LaplaceKernel(2)
+    #kernel = OneKernel()
 
-mesh = make_curve_mesh(
-        #lambda t: ellipse(1, t),
-        starfish,
-        np.linspace(0, 1, nelements+1),
-        target_order)
+    mesh = make_curve_mesh(
+            #lambda t: ellipse(1, t),
+            starfish,
+            np.linspace(0, 1, nelements+1),
+            target_order)
 
-from pytential.qbx import QBXLayerPotentialSource
-from meshmode.discretization import Discretization
-from meshmode.discretization.poly_element import \
-        InterpolatoryQuadratureSimplexGroupFactory
+    from pytential.qbx import QBXLayerPotentialSource
+    from meshmode.discretization import Discretization
+    from meshmode.discretization.poly_element import \
+            InterpolatoryQuadratureSimplexGroupFactory
 
-density_discr = Discretization(
-        cl_ctx, mesh,
-        InterpolatoryQuadratureSimplexGroupFactory(target_order))
+    pre_density_discr = Discretization(
+            cl_ctx, mesh,
+            InterpolatoryQuadratureSimplexGroupFactory(target_order))
 
-qbx = QBXLayerPotentialSource(
-        density_discr, fine_order=2*target_order,
-        qbx_order=qbx_order, fmm_order=qbx_order)
-slow_qbx = QBXLayerPotentialSource(
-        density_discr, fine_order=2*target_order,
-        qbx_order=qbx_order, fmm_order=False)
+    slow_qbx, _ = QBXLayerPotentialSource(
+            pre_density_discr, fine_order=2*target_order,
+            qbx_order=qbx_order, fmm_order=False,
+            target_association_tolerance=.05
+            ).with_refinement()
+    qbx = slow_qbx.copy(fmm_order=10)
+    density_discr = slow_qbx.density_discr
 
-nodes = density_discr.nodes().with_queue(queue)
+    nodes = density_discr.nodes().with_queue(queue)
 
-angle = cl.clmath.atan2(nodes[1], nodes[0])
+    angle = cl.clmath.atan2(nodes[1], nodes[0])
 
-from pytential import bind, sym
-d = sym.Derivative()
-#op = d.nabla[0] * d(sym.S(kernel, sym.var("sigma")))
-#op = sym.D(kernel, sym.var("sigma"))
-op = sym.S(kernel, sym.var("sigma"))
+    from pytential import bind, sym
+    #op = sym.d_dx(sym.S(kernel, sym.var("sigma")), qbx_forced_limit=None)
+    #op = sym.D(kernel, sym.var("sigma"), qbx_forced_limit=None)
+    op = sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None)
 
-sigma = cl.clmath.cos(mode_nr*angle)
+    sigma = cl.clmath.cos(mode_nr*angle)
 
-if isinstance(kernel, HelmholtzKernel):
-    sigma = sigma.astype(np.complex128)
+    if isinstance(kernel, HelmholtzKernel):
+        sigma = sigma.astype(np.complex128)
 
-bound_bdry_op = bind(qbx, op)
+    fplot = FieldPlotter(np.zeros(2), extent=5, npoints=600)
+    from pytential.target import PointsTarget
 
-fplot = FieldPlotter(np.zeros(2), extent=5, npoints=600)
-from pytential.target import PointsTarget
+    fld_in_vol = bind(
+            (slow_qbx, PointsTarget(fplot.points)),
+            op)(queue, sigma=sigma, k=k).get()
 
-fld_in_vol = bind(
-        (slow_qbx, PointsTarget(fplot.points)),
-        op)(queue, sigma=sigma, k=k).get()
+    fmm_fld_in_vol = bind(
+            (qbx, PointsTarget(fplot.points)),
+            op)(queue, sigma=sigma, k=k).get()
 
-fmm_fld_in_vol = bind(
-        (qbx, PointsTarget(fplot.points)),
-        op)(queue, sigma=sigma, k=k).get()
+    err = fmm_fld_in_vol-fld_in_vol
 
-err = fmm_fld_in_vol-fld_in_vol
-im = fplot.show_scalar_in_matplotlib(np.log10(np.abs(err)))
+    import matplotlib
+    matplotlib.use('Agg')
+    im = fplot.show_scalar_in_matplotlib(np.log10(np.abs(err)))
 
-from matplotlib.colors import Normalize
-im.set_norm(Normalize(vmin=-6, vmax=0))
+    from matplotlib.colors import Normalize
+    im.set_norm(Normalize(vmin=-12, vmax=0))
 
-import matplotlib.pyplot as pt
-from matplotlib.ticker import NullFormatter
-pt.gca().xaxis.set_major_formatter(NullFormatter())
-pt.gca().yaxis.set_major_formatter(NullFormatter())
+    import matplotlib.pyplot as pt
+    from matplotlib.ticker import NullFormatter
+    pt.gca().xaxis.set_major_formatter(NullFormatter())
+    pt.gca().yaxis.set_major_formatter(NullFormatter())
 
-cb = pt.colorbar(shrink=0.9)
-cb.set_label(r"$\log_{10}(\mathdefault{Error})$")
+    cb = pt.colorbar(shrink=0.9)
+    cb.set_label(r"$\log_{10}(\mathdefault{Error})$")
 
-pt.savefig("fmm-error-order-%d.pdf" % qbx_order)
+    pt.savefig("fmm-error-order-%d.pdf" % qbx_order)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/blob-2d.step b/examples/geometries/blob-2d.step
similarity index 100%
rename from examples/blob-2d.step
rename to examples/geometries/blob-2d.step
diff --git a/examples/circle.step b/examples/geometries/circle.step
similarity index 100%
rename from examples/circle.step
rename to examples/geometries/circle.step
diff --git a/examples/circles.step b/examples/geometries/circles.step
similarity index 100%
rename from examples/circles.step
rename to examples/geometries/circles.step
diff --git a/examples/ellipsoid.step b/examples/geometries/ellipsoid.step
similarity index 100%
rename from examples/ellipsoid.step
rename to examples/geometries/ellipsoid.step
diff --git a/examples/molecule.step b/examples/geometries/molecule.step
similarity index 100%
rename from examples/molecule.step
rename to examples/geometries/molecule.step
diff --git a/examples/two-balls.step b/examples/geometries/two-balls.step
similarity index 100%
rename from examples/two-balls.step
rename to examples/geometries/two-balls.step
diff --git a/examples/two-cylinders-smooth.step b/examples/geometries/two-cylinders-smooth.step
similarity index 100%
rename from examples/two-cylinders-smooth.step
rename to examples/geometries/two-cylinders-smooth.step
diff --git a/examples/helmholtz-dirichlet.py b/examples/helmholtz-dirichlet.py
index d6d9baa6c609399028305edb249565d0457d878d..22f8fa8a0083da62f1f17d7c79053056ce3ab33c 100644
--- a/examples/helmholtz-dirichlet.py
+++ b/examples/helmholtz-dirichlet.py
@@ -25,7 +25,7 @@ k = 3
 
 def main():
     import logging
-    logging.basicConfig(level=logging.INFO)
+    logging.basicConfig(level=logging.WARNING)  # INFO for more progress info
 
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
@@ -167,7 +167,7 @@ def main():
 
     #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5)
     fplot.write_vtk_file(
-            "potential.vts",
+            "potential-helm.vts",
             [
                 ("potential", fld_in_vol),
                 ("indicator", indicator),
diff --git a/examples/laplace-dirichlet-3d.py b/examples/laplace-dirichlet-3d.py
new file mode 100644
index 0000000000000000000000000000000000000000..4166dddfa8d8e1606885e167fa65a8881fe64484
--- /dev/null
+++ b/examples/laplace-dirichlet-3d.py
@@ -0,0 +1,170 @@
+import numpy as np
+import numpy.linalg as la
+import pyopencl as cl
+import pyopencl.clmath  # noqa
+
+from meshmode.discretization import Discretization
+from meshmode.discretization.poly_element import \
+        InterpolatoryQuadratureSimplexGroupFactory
+
+from pytential import bind, sym, norm  # noqa
+from pytential.target import PointsTarget
+
+# {{{ set some constants for use below
+
+nelements = 20
+bdry_quad_order = 4
+mesh_order = bdry_quad_order
+qbx_order = bdry_quad_order
+bdry_ovsmp_quad_order = 4*bdry_quad_order
+fmm_order = 3
+
+# }}}
+
+
+def main():
+    import logging
+    logging.basicConfig(level=logging.WARNING)  # INFO for more progress info
+
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    from meshmode.mesh.generation import generate_torus
+
+    rout = 10
+    rin = 1
+    if 1:
+        base_mesh = generate_torus(
+                rout, rin, 40, 4,
+                mesh_order)
+
+        from meshmode.mesh.processing import affine_map, merge_disjoint_meshes
+        # nx = 1
+        # ny = 1
+        nz = 1
+        dz = 0
+        meshes = [
+                affine_map(
+                    base_mesh,
+                    A=np.diag([1, 1, 1]),
+                    b=np.array([0, 0, iz*dz]))
+                for iz in range(nz)]
+
+        mesh = merge_disjoint_meshes(meshes, single_group=True)
+
+        if 0:
+            from meshmode.mesh.visualization import draw_curve
+            draw_curve(mesh)
+            import matplotlib.pyplot as plt
+            plt.show()
+
+    pre_density_discr = Discretization(
+            cl_ctx, mesh,
+            InterpolatoryQuadratureSimplexGroupFactory(bdry_quad_order))
+
+    from pytential.qbx import (
+            QBXLayerPotentialSource, QBXTargetAssociationFailedException)
+    qbx, _ = QBXLayerPotentialSource(
+            pre_density_discr, fine_order=bdry_ovsmp_quad_order, qbx_order=qbx_order,
+            fmm_order=fmm_order
+            ).with_refinement()
+    density_discr = qbx.density_discr
+
+    # {{{ describe bvp
+
+    from sumpy.kernel import LaplaceKernel
+    kernel = LaplaceKernel(3)
+
+    cse = sym.cse
+
+    sigma_sym = sym.var("sigma")
+    #sqrt_w = sym.sqrt_jac_q_weight(3)
+    sqrt_w = 1
+    inv_sqrt_w_sigma = cse(sigma_sym/sqrt_w)
+
+    # -1 for interior Dirichlet
+    # +1 for exterior Dirichlet
+    loc_sign = +1
+
+    bdry_op_sym = (loc_sign*0.5*sigma_sym
+            + sqrt_w*(
+                sym.S(kernel, inv_sqrt_w_sigma)
+                + sym.D(kernel, inv_sqrt_w_sigma)
+                ))
+
+    # }}}
+
+    bound_op = bind(qbx, bdry_op_sym)
+
+    # {{{ fix rhs and solve
+
+    nodes = density_discr.nodes().with_queue(queue)
+    source = np.array([rout, 0, 0])
+
+    def u_incoming_func(x):
+        #        return 1/cl.clmath.sqrt( (x[0] - source[0])**2
+        #                                +(x[1] - source[1])**2
+        #                                +(x[2] - source[2])**2 )
+        return 1.0/la.norm(x.get()-source[:, None], axis=0)
+
+    bc = cl.array.to_device(queue, u_incoming_func(nodes))
+
+    bvp_rhs = bind(qbx, sqrt_w*sym.var("bc"))(queue, bc=bc)
+
+    from pytential.solve import gmres
+    gmres_result = gmres(
+            bound_op.scipy_op(queue, "sigma", dtype=np.float64),
+            bvp_rhs, tol=1e-14, progress=True,
+            stall_iterations=0,
+            hard_failure=True)
+
+    sigma = bind(qbx, sym.var("sigma")/sqrt_w)(queue, sigma=gmres_result.solution)
+
+    # }}}
+
+    from meshmode.discretization.visualization import make_visualizer
+    bdry_vis = make_visualizer(queue, density_discr, 20)
+    bdry_vis.write_vtk_file("laplace.vtu", [
+        ("sigma", sigma),
+        ])
+
+    # {{{ postprocess/visualize
+
+    repr_kwargs = dict(qbx_forced_limit=None)
+    representation_sym = (
+            sym.S(kernel, inv_sqrt_w_sigma, **repr_kwargs)
+            + sym.D(kernel, inv_sqrt_w_sigma, **repr_kwargs))
+
+    from sumpy.visualization import FieldPlotter
+    fplot = FieldPlotter(np.zeros(3), extent=20, npoints=50)
+
+    targets = cl.array.to_device(queue, fplot.points)
+
+    qbx_stick_out = qbx.copy(target_stick_out_factor=0.2)
+
+    try:
+        fld_in_vol = bind(
+                (qbx_stick_out, PointsTarget(targets)),
+                representation_sym)(queue, sigma=sigma).get()
+    except QBXTargetAssociationFailedException as e:
+        fplot.write_vtk_file(
+                "failed-targets.vts",
+                [
+                    ("failed", e.failed_target_flags.get(queue))
+                    ]
+                )
+        raise
+
+    #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5)
+    fplot.write_vtk_file(
+            "potential-laplace-3d.vts",
+            [
+                ("potential", fld_in_vol),
+                ]
+            )
+
+    # }}}
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/layerpot-3d.py b/examples/layerpot-3d.py
index 0a35ebd7fa418b8cf517ada57bc7dee5e39f1471..28f0967e8aec28332902a128d0fb1efafb100d4e 100644
--- a/examples/layerpot-3d.py
+++ b/examples/layerpot-3d.py
@@ -9,10 +9,6 @@ import faulthandler
 from six.moves import range
 faulthandler.enable()
 
-import logging
-logging.basicConfig(level=logging.INFO)
-logger = logging.getLogger(__name__)
-
 cl_ctx = cl.create_some_context()
 queue = cl.CommandQueue(cl_ctx)
 
@@ -21,10 +17,10 @@ qbx_order = 3
 mode_nr = 4
 
 if 1:
-    cad_file_name = "ellipsoid.step"
+    cad_file_name = "geometries/ellipsoid.step"
     h = 0.6
 else:
-    cad_file_name = "two-cylinders-smooth.step"
+    cad_file_name = "geometries/two-cylinders-smooth.step"
     h = 0.4
 
 k = 0
@@ -34,76 +30,86 @@ else:
     kernel = LaplaceKernel(3)
 #kernel = OneKernel()
 
-from meshmode.mesh.io import generate_gmsh, FileSource
-mesh = generate_gmsh(
-        FileSource(cad_file_name), 2, order=2,
-        other_options=["-string", "Mesh.CharacteristicLengthMax = %g;" % h])
 
-from meshmode.mesh.processing import perform_flips
-# Flip elements--gmsh generates inside-out geometry.
-mesh = perform_flips(mesh, np.ones(mesh.nelements))
+def main():
+    import logging
+    logger = logging.getLogger(__name__)
+    logging.basicConfig(level=logging.WARNING)  # INFO for more progress info
+
+    from meshmode.mesh.io import generate_gmsh, FileSource
+    mesh = generate_gmsh(
+            FileSource(cad_file_name), 2, order=2,
+            other_options=["-string", "Mesh.CharacteristicLengthMax = %g;" % h])
+
+    from meshmode.mesh.processing import perform_flips
+    # Flip elements--gmsh generates inside-out geometry.
+    mesh = perform_flips(mesh, np.ones(mesh.nelements))
+
+    from meshmode.mesh.processing import find_bounding_box
+    bbox_min, bbox_max = find_bounding_box(mesh)
+    bbox_center = 0.5*(bbox_min+bbox_max)
+    bbox_size = max(bbox_max-bbox_min) / 2
+
+    logger.info("%d elements" % mesh.nelements)
 
-from meshmode.mesh.processing import find_bounding_box
-bbox_min, bbox_max = find_bounding_box(mesh)
-bbox_center = 0.5*(bbox_min+bbox_max)
-bbox_size = max(bbox_max-bbox_min) / 2
+    from pytential.qbx import QBXLayerPotentialSource
+    from meshmode.discretization import Discretization
+    from meshmode.discretization.poly_element import \
+            InterpolatoryQuadratureSimplexGroupFactory
 
-logger.info("%d elements" % mesh.nelements)
+    density_discr = Discretization(
+            cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order))
 
-from pytential.qbx import QBXLayerPotentialSource
-from meshmode.discretization import Discretization
-from meshmode.discretization.poly_element import \
-        InterpolatoryQuadratureSimplexGroupFactory
+    qbx, _ = QBXLayerPotentialSource(density_discr, 4*target_order, qbx_order,
+            fmm_order=qbx_order + 3,
+            target_association_tolerance=0.15).with_refinement()
 
-density_discr = Discretization(
-        cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order))
+    nodes = density_discr.nodes().with_queue(queue)
 
-qbx, _ = QBXLayerPotentialSource(density_discr, 4*target_order, qbx_order,
-        fmm_order=qbx_order + 3,
-        target_association_tolerance=0.15).with_refinement()
+    angle = cl.clmath.atan2(nodes[1], nodes[0])
 
-nodes = density_discr.nodes().with_queue(queue)
+    from pytential import bind, sym
+    #op = sym.d_dx(sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None))
+    op = sym.D(kernel, sym.var("sigma"), qbx_forced_limit=None)
+    #op = sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None)
 
-angle = cl.clmath.atan2(nodes[1], nodes[0])
+    sigma = cl.clmath.cos(mode_nr*angle)
+    if 0:
+        sigma = 0*angle
+        from random import randrange
+        for i in range(5):
+            sigma[randrange(len(sigma))] = 1
 
-from pytential import bind, sym
-#op = sym.d_dx(sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None))
-op = sym.D(kernel, sym.var("sigma"), qbx_forced_limit=None)
-#op = sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None)
+    if isinstance(kernel, HelmholtzKernel):
+        sigma = sigma.astype(np.complex128)
 
-sigma = cl.clmath.cos(mode_nr*angle)
-if 0:
-    sigma = 0*angle
-    from random import randrange
-    for i in range(5):
-        sigma[randrange(len(sigma))] = 1
+    fplot = FieldPlotter(bbox_center, extent=3.5*bbox_size, npoints=150)
 
-if isinstance(kernel, HelmholtzKernel):
-    sigma = sigma.astype(np.complex128)
+    from pytential.target import PointsTarget
+    fld_in_vol = bind(
+            (qbx, PointsTarget(fplot.points)),
+            op)(queue, sigma=sigma, k=k).get()
 
-fplot = FieldPlotter(bbox_center, extent=3.5*bbox_size, npoints=150)
+    #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5)
+    fplot.write_vtk_file(
+            "potential-3d.vts",
+            [
+                ("potential", fld_in_vol)
+                ]
+            )
 
-from pytential.target import PointsTarget
-fld_in_vol = bind(
-        (qbx, PointsTarget(fplot.points)),
-        op)(queue, sigma=sigma, k=k).get()
+    bdry_normals = bind(
+            density_discr,
+            sym.normal(density_discr.ambient_dim))(queue).as_vector(dtype=object)
 
-#fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5)
-fplot.write_vtk_file(
-        "potential.vts",
-        [
-            ("potential", fld_in_vol)
-            ]
-        )
+    from meshmode.discretization.visualization import make_visualizer
+    bdry_vis = make_visualizer(queue, density_discr, target_order)
 
-bdry_normals = bind(
-        density_discr,
-        sym.normal(density_discr.ambient_dim))(queue).as_vector(dtype=object)
+    bdry_vis.write_vtk_file("source-3d.vtu", [
+        ("sigma", sigma),
+        ("bdry_normals", bdry_normals),
+        ])
 
-from meshmode.discretization.visualization import make_visualizer
-bdry_vis = make_visualizer(queue, density_discr, target_order)
 
-bdry_vis.write_vtk_file("source.vtu", [
-    ("sigma", sigma),
-    ("bdry_normals", bdry_normals),
-    ])
+if __name__ == "__main__":
+    main()
diff --git a/examples/layerpot.py b/examples/layerpot.py
index 0371b1ff557f57a7b2f839be6571bded120728bb..7b4737da00d6d1d1cc76e31f39932fc7c12783e8 100644
--- a/examples/layerpot.py
+++ b/examples/layerpot.py
@@ -15,13 +15,6 @@ import faulthandler
 from six.moves import range
 faulthandler.enable()
 
-import logging
-logging.basicConfig(level=logging.INFO)
-logger = logging.getLogger(__name__)
-
-cl_ctx = cl.create_some_context()
-queue = cl.CommandQueue(cl_ctx)
-
 target_order = 16
 qbx_order = 3
 nelements = 60
@@ -36,96 +29,107 @@ else:
     kernel_kwargs = {}
 #kernel = OneKernel()
 
-from meshmode.mesh.generation import (  # noqa
-        make_curve_mesh, starfish, ellipse, drop)
-mesh = make_curve_mesh(
-        #lambda t: ellipse(1, t),
-        starfish,
-        np.linspace(0, 1, nelements+1),
-        target_order)
 
-from pytential.qbx import QBXLayerPotentialSource
-from meshmode.discretization import Discretization
-from meshmode.discretization.poly_element import \
-        InterpolatoryQuadratureSimplexGroupFactory
+def main():
+    import logging
+    logging.basicConfig(level=logging.WARNING)  # INFO for more progress info
+
+    cl_ctx = cl.create_some_context()
+    queue = cl.CommandQueue(cl_ctx)
+
+    from meshmode.mesh.generation import (  # noqa
+            make_curve_mesh, starfish, ellipse, drop)
+    mesh = make_curve_mesh(
+            #lambda t: ellipse(1, t),
+            starfish,
+            np.linspace(0, 1, nelements+1),
+            target_order)
+
+    from pytential.qbx import QBXLayerPotentialSource
+    from meshmode.discretization import Discretization
+    from meshmode.discretization.poly_element import \
+            InterpolatoryQuadratureSimplexGroupFactory
+
+    pre_density_discr = Discretization(
+            cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order))
 
-pre_density_discr = Discretization(
-        cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order))
+    qbx, _ = QBXLayerPotentialSource(pre_density_discr, 4*target_order, qbx_order,
+            fmm_order=qbx_order+3,
+            target_association_tolerance=0.005).with_refinement()
 
-qbx, _ = QBXLayerPotentialSource(pre_density_discr, 4*target_order, qbx_order,
-        fmm_order=qbx_order+3,
-        target_association_tolerance=0.005).with_refinement()
+    density_discr = qbx.density_discr
 
-density_discr = qbx.density_discr
+    nodes = density_discr.nodes().with_queue(queue)
 
-nodes = density_discr.nodes().with_queue(queue)
+    angle = cl.clmath.atan2(nodes[1], nodes[0])
 
-angle = cl.clmath.atan2(nodes[1], nodes[0])
+    def op(**kwargs):
+        kwargs.update(kernel_kwargs)
 
+        #op = sym.d_dx(sym.S(kernel, sym.var("sigma"), **kwargs))
+        return sym.D(kernel, sym.var("sigma"), **kwargs)
+        #op = sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None, **kwargs)
 
-def op(**kwargs):
-    kwargs.update(kernel_kwargs)
+    sigma = cl.clmath.cos(mode_nr*angle)
+    if 0:
+        sigma = 0*angle
+        from random import randrange
+        for i in range(5):
+            sigma[randrange(len(sigma))] = 1
 
-    #op = sym.d_dx(sym.S(kernel, sym.var("sigma"), **kwargs))
-    return sym.D(kernel, sym.var("sigma"), **kwargs)
-    #op = sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None, **kwargs)
+    if isinstance(kernel, HelmholtzKernel):
+        sigma = sigma.astype(np.complex128)
 
+    bound_bdry_op = bind(qbx, op())
+    #mlab.figure(bgcolor=(1, 1, 1))
+    if 1:
+        fplot = FieldPlotter(np.zeros(2), extent=5, npoints=1000)
+        from pytential.target import PointsTarget
 
-sigma = cl.clmath.cos(mode_nr*angle)
-if 0:
-    sigma = 0*angle
-    from random import randrange
-    for i in range(5):
-        sigma[randrange(len(sigma))] = 1
+        targets_dev = cl.array.to_device(queue, fplot.points)
+        fld_in_vol = bind(
+                (qbx, PointsTarget(targets_dev)),
+                op(qbx_forced_limit=None))(queue, sigma=sigma, k=k).get()
 
-if isinstance(kernel, HelmholtzKernel):
-    sigma = sigma.astype(np.complex128)
+        if enable_mayavi:
+            fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5)
+        else:
+            fplot.write_vtk_file(
+                    "potential-2d.vts",
+                    [
+                        ("potential", fld_in_vol)
+                        ]
+                    )
 
-bound_bdry_op = bind(qbx, op())
-#mlab.figure(bgcolor=(1, 1, 1))
-if 1:
-    fplot = FieldPlotter(np.zeros(2), extent=5, npoints=1000)
-    from pytential.target import PointsTarget
+    if 0:
+        def apply_op(density):
+            return bound_bdry_op(
+                    queue, sigma=cl.array.to_device(queue, density), k=k).get()
 
-    targets_dev = cl.array.to_device(queue, fplot.points)
-    fld_in_vol = bind(
-            (qbx, PointsTarget(targets_dev)),
-            op(qbx_forced_limit=None))(queue, sigma=sigma, k=k).get()
+        from sumpy.tools import build_matrix
+        n = len(sigma)
+        mat = build_matrix(apply_op, dtype=np.float64, shape=(n, n))
+
+        import matplotlib.pyplot as pt
+        pt.imshow(mat)
+        pt.colorbar()
+        pt.show()
 
     if enable_mayavi:
-        fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5)
-    else:
-        fplot.write_vtk_file(
-                "potential.vts",
-                [
-                    ("potential", fld_in_vol)
-                    ]
-                )
-
-if 0:
-    def apply_op(density):
-        return bound_bdry_op(
-                queue, sigma=cl.array.to_device(queue, density), k=k).get()
-
-    from sumpy.tools import build_matrix
-    n = len(sigma)
-    mat = build_matrix(apply_op, dtype=np.float64, shape=(n, n))
-
-    import matplotlib.pyplot as pt
-    pt.imshow(mat)
-    pt.colorbar()
-    pt.show()
+        # {{{ plot boundary field
 
-if enable_mayavi:
-    # {{{ plot boundary field
+        fld_on_bdry = bound_bdry_op(queue, sigma=sigma, k=k).get()
 
-    fld_on_bdry = bound_bdry_op(queue, sigma=sigma, k=k).get()
+        nodes_host = density_discr.nodes().get(queue=queue)
+        mlab.points3d(nodes_host[0], nodes_host[1],
+                fld_on_bdry.real, scale_factor=0.03)
 
-    nodes_host = density_discr.nodes().get(queue=queue)
-    mlab.points3d(nodes_host[0], nodes_host[1], fld_on_bdry.real, scale_factor=0.03)
+        # }}}
 
-    # }}}
+    if enable_mayavi:
+        mlab.colorbar()
+        mlab.show()
 
-if enable_mayavi:
-    mlab.colorbar()
-    mlab.show()
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/perf-model.py b/examples/perf-model.py
deleted file mode 100644
index 3a87d63113ecedda8378485a4869e2771f0aaa4e..0000000000000000000000000000000000000000
--- a/examples/perf-model.py
+++ /dev/null
@@ -1,18 +0,0 @@
-# ------------------------------
-nlevels = 6
-nboxes = 1365
-nsources = 60
-ntargets = 12040
-form_mp = 60*p_fmm
-prop_upward = 1365*p_fmm**2
-part_direct = 196560
-m2l = 31920*p_fmm**2
-mp_eval = 0
-form_local = 65000*p_fmm
-prop_downward = 1365*p_fmm**2
-eval_part = 12040*p_fmm
-ncenters = 2040
-qbxl_direct = 2370940*p_qbx
-qbx_m2l = 35339*p_fmm*p_qbx
-qbx_l2l = 2040*p_fmm*p_qbx
-qbx_eval = 1902*p_qbx
diff --git a/examples/scaling-study.py b/examples/scaling-study.py
index 183fc915cb286ccb7731a76e9c1ed6cc3869efd5..3327e3c8c6ce71262018551008a203a04d68e70b 100644
--- a/examples/scaling-study.py
+++ b/examples/scaling-study.py
@@ -16,8 +16,8 @@ bdry_quad_order = 4
 mesh_order = bdry_quad_order
 qbx_order = bdry_quad_order
 bdry_ovsmp_quad_order = 4*bdry_quad_order
-fmm_order = 25
-k = 25
+fmm_order = 10
+k = 0
 
 # }}}
 
@@ -54,7 +54,7 @@ def make_mesh(nx, ny):
 
 def timing_run(nx, ny):
     import logging
-    logging.basicConfig(level=logging.INFO)
+    logging.basicConfig(level=logging.WARNING)  # INFO for more progress info
 
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
@@ -154,6 +154,7 @@ def timing_run(nx, ny):
                 sym_op)(
                 queue, sigma=ones_density).get()
 
+        qbx_stick_out = qbx.copy(target_stick_out_factor=0.1)
         try:
             fld_in_vol = bind(
                     (qbx_stick_out, PointsTarget(targets)),
@@ -169,7 +170,7 @@ def timing_run(nx, ny):
 
         #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5)
         fplot.write_vtk_file(
-                "potential.vts",
+                "potential-scaling.vts",
                 [
                     ("potential", fld_in_vol),
                     ("indicator", indicator)
diff --git a/experiments/README.md b/experiments/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..d0f56efd2f7e29d95e945154962d5f00fdd98cb1
--- /dev/null
+++ b/experiments/README.md
@@ -0,0 +1,7 @@
+# Experiments
+
+What you find in this directory are experiments that *may* have done something
+useful at some point (or not). Unlike `examples`, they are not being tested on
+an ongoing basis.
+
+So if what you find here breaks for you, you get to keep both pieces.
diff --git a/examples/cahn-hilliard.py b/experiments/cahn-hilliard.py
similarity index 100%
rename from examples/cahn-hilliard.py
rename to experiments/cahn-hilliard.py
diff --git a/examples/find-photonic-mode-sk.py b/experiments/find-photonic-mode-sk.py
similarity index 100%
rename from examples/find-photonic-mode-sk.py
rename to experiments/find-photonic-mode-sk.py
diff --git a/examples/find-photonic-mode.py b/experiments/find-photonic-mode.py
similarity index 100%
rename from examples/find-photonic-mode.py
rename to experiments/find-photonic-mode.py
diff --git a/examples/helmholtz-expression-tree.py b/experiments/helmholtz-expression-tree.py
similarity index 100%
rename from examples/helmholtz-expression-tree.py
rename to experiments/helmholtz-expression-tree.py
diff --git a/examples/maxwell.py b/experiments/maxwell.py
similarity index 100%
rename from examples/maxwell.py
rename to experiments/maxwell.py
diff --git a/examples/maxwell_sphere.py b/experiments/maxwell_sphere.py
similarity index 100%
rename from examples/maxwell_sphere.py
rename to experiments/maxwell_sphere.py
diff --git a/examples/poisson.py b/experiments/poisson.py
similarity index 100%
rename from examples/poisson.py
rename to experiments/poisson.py
diff --git a/examples/qbx-tangential-deriv-jump.py b/experiments/qbx-tangential-deriv-jump.py
similarity index 100%
rename from examples/qbx-tangential-deriv-jump.py
rename to experiments/qbx-tangential-deriv-jump.py
diff --git a/examples/stokes-2d-interior.py b/experiments/stokes-2d-interior.py
similarity index 100%
rename from examples/stokes-2d-interior.py
rename to experiments/stokes-2d-interior.py
diff --git a/examples/two-domain-helmholtz.py b/experiments/two-domain-helmholtz.py
similarity index 100%
rename from examples/two-domain-helmholtz.py
rename to experiments/two-domain-helmholtz.py