diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py
index 1431f75b1d01ea0e5b6055c9352d1f2d48aec25e..b3e18650cd86f9aa0a6ab3b52db60446ac872ec8 100644
--- a/pytential/symbolic/pde/maxwell/__init__.py
+++ b/pytential/symbolic/pde/maxwell/__init__.py
@@ -152,7 +152,8 @@ class PECChargeCurrentMFIEOperator:
         return xyz_to_tangential(sym.n_cross(Hinc_xyz))
 
     def rho_operator(self, loc, rho):
-        return loc*0.5*rho+sym.Sp(self.kernel, rho, k=self.k)
+        return loc*0.5*rho+sym.Sp(
+                self.kernel, rho, k=self.k, qbx_forced_limit="avg")
 
     def rho_rhs(self, Jt, Einc_xyz):
         Jxyz = cse(tangential_to_xyz(Jt), "Jxyz")
diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py
index 7f700fa104b2bc10e99416d594aecd090eca69eb..af340805f7af66c381eae1ce64b8fe6fb6f3ca0d 100644
--- a/test/test_global_qbx.py
+++ b/test/test_global_qbx.py
@@ -79,8 +79,8 @@ def iter_elements(discr):
             discr_nodes_idx += discr_group.nunit_nodes
 
 
-def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None):
-    cl_ctx = ctx_getter()
+def run_source_refinement_test(ctx_factory, mesh, order, helmholtz_k=None):
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     from meshmode.discretization import Discretization
@@ -201,30 +201,30 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None):
     ("20-to-1 ellipse", partial(ellipse, 20), 100),
     ("horseshoe", horseshoe, 64),
     ])
-def test_source_refinement_2d(ctx_getter, curve_name, curve_f, nelements):
+def test_source_refinement_2d(ctx_factory, curve_name, curve_f, nelements):
     helmholtz_k = 10
     order = 8
 
     mesh = make_curve_mesh(curve_f, np.linspace(0, 1, nelements+1), order)
-    run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k)
+    run_source_refinement_test(ctx_factory, mesh, order, helmholtz_k)
 
 
 @pytest.mark.parametrize(("surface_name", "surface_f", "order"), [
     ("sphere", partial(generate_icosphere, 1), 4),
     ("torus", partial(generate_torus, 3, 1, n_inner=10, n_outer=7), 6),
     ])
-def test_source_refinement_3d(ctx_getter, surface_name, surface_f, order):
+def test_source_refinement_3d(ctx_factory, surface_name, surface_f, order):
     mesh = surface_f(order=order)
-    run_source_refinement_test(ctx_getter, mesh, order)
+    run_source_refinement_test(ctx_factory, mesh, order)
 
 
 @pytest.mark.parametrize(("curve_name", "curve_f", "nelements"), [
     ("20-to-1 ellipse", partial(ellipse, 20), 100),
     ("horseshoe", horseshoe, 64),
     ])
-def test_target_association(ctx_getter, curve_name, curve_f, nelements,
+def test_target_association(ctx_factory, curve_name, curve_f, nelements,
         visualize=False):
-    cl_ctx = ctx_getter()
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     # {{{ generate lpot source
@@ -403,8 +403,8 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements,
     # }}}
 
 
-def test_target_association_failure(ctx_getter):
-    cl_ctx = ctx_getter()
+def test_target_association_failure(ctx_factory):
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     # {{{ generate circle
diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py
index 4b5d70b675871031023d3f5cf40c049c008d0f5b..c99f3073c8a9be4eb535153411eca05ad470e6dc 100644
--- a/test/test_layer_pot.py
+++ b/test/test_layer_pot.py
@@ -46,8 +46,8 @@ circle = partial(ellipse, 1)
 
 # {{{ geometry test
 
-def test_geometry(ctx_getter):
-    cl_ctx = ctx_getter()
+def test_geometry(ctx_factory):
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     nelements = 30
@@ -79,10 +79,10 @@ def test_geometry(ctx_getter):
 # {{{ test off-surface eval
 
 @pytest.mark.parametrize("use_fmm", [True, False])
-def test_off_surface_eval(ctx_getter, use_fmm, do_plot=False):
+def test_off_surface_eval(ctx_factory, use_fmm, do_plot=False):
     logging.basicConfig(level=logging.INFO)
 
-    cl_ctx = ctx_getter()
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     # prevent cache 'splosion
@@ -146,10 +146,10 @@ def test_off_surface_eval(ctx_getter, use_fmm, do_plot=False):
 
 # {{{ test off-surface eval vs direct
 
-def test_off_surface_eval_vs_direct(ctx_getter,  do_plot=False):
+def test_off_surface_eval_vs_direct(ctx_factory,  do_plot=False):
     logging.basicConfig(level=logging.INFO)
 
-    cl_ctx = ctx_getter()
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     # prevent cache 'splosion
@@ -227,8 +227,8 @@ def test_off_surface_eval_vs_direct(ctx_getter,  do_plot=False):
 # {{{ unregularized tests
 
 
-def test_unregularized_with_ones_kernel(ctx_getter):
-    cl_ctx = ctx_getter()
+def test_unregularized_with_ones_kernel(ctx_factory):
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     nelements = 10
@@ -268,8 +268,8 @@ def test_unregularized_with_ones_kernel(ctx_getter):
     assert np.allclose(result_nonself.get(), 2 * np.pi)
 
 
-def test_unregularized_off_surface_fmm_vs_direct(ctx_getter):
-    cl_ctx = ctx_getter()
+def test_unregularized_off_surface_fmm_vs_direct(ctx_factory):
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     nelements = 300
@@ -317,8 +317,8 @@ def test_unregularized_off_surface_fmm_vs_direct(ctx_getter):
 
 # {{{ test performance data gathering
 
-def test_perf_data_gathering(ctx_getter, n_arms=5):
-    cl_ctx = ctx_getter()
+def test_perf_data_gathering(ctx_factory, n_arms=5):
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     # prevent cache 'splosion
diff --git a/test/test_layer_pot_eigenvalues.py b/test/test_layer_pot_eigenvalues.py
index b4e986d4af6e16cc149effac6a08403d11668b4c..fdfdc93738d823dae9bcf24828362b54b4bc84ff 100644
--- a/test/test_layer_pot_eigenvalues.py
+++ b/test/test_layer_pot_eigenvalues.py
@@ -61,14 +61,14 @@ except ImportError:
 
             (2, 7, 5, True),
             ])
-def test_ellipse_eigenvalues(ctx_getter, ellipse_aspect, mode_nr, qbx_order,
+def test_ellipse_eigenvalues(ctx_factory, ellipse_aspect, mode_nr, qbx_order,
         force_direct):
     logging.basicConfig(level=logging.INFO)
 
     print("ellipse_aspect: %s, mode_nr: %d, qbx_order: %d" % (
             ellipse_aspect, mode_nr, qbx_order))
 
-    cl_ctx = ctx_getter()
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     target_order = 8
@@ -148,7 +148,7 @@ def test_ellipse_eigenvalues(ctx_getter, ellipse_aspect, mode_nr, qbx_order,
 
         sigma = cl.clmath.cos(mode_nr*angle)/J
 
-        s_sigma_op = bind(qbx, sym.S(lap_knl, sym.var("sigma")))
+        s_sigma_op = bind(qbx, sym.S(lap_knl, sym.var("sigma"), qbx_forced_limit=+1))
         s_sigma = s_sigma_op(queue=queue, sigma=sigma)
 
         # SIGN BINGO! :)
@@ -249,13 +249,13 @@ def test_ellipse_eigenvalues(ctx_getter, ellipse_aspect, mode_nr, qbx_order,
     "sumpy",
     "fmmlib",
     ])
-def test_sphere_eigenvalues(ctx_getter, mode_m, mode_n, qbx_order,
+def test_sphere_eigenvalues(ctx_factory, mode_m, mode_n, qbx_order,
         fmm_backend):
     logging.basicConfig(level=logging.INFO)
 
     special = pytest.importorskip("scipy.special")
 
-    cl_ctx = ctx_getter()
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     target_order = 8
@@ -310,7 +310,7 @@ def test_sphere_eigenvalues(ctx_getter, mode_m, mode_n, qbx_order,
 
         # {{{ single layer
 
-        s_sigma_op = bind(qbx, sym.S(lap_knl, sym.var("sigma")))
+        s_sigma_op = bind(qbx, sym.S(lap_knl, sym.var("sigma"), qbx_forced_limit=+1))
         s_sigma = s_sigma_op(queue=queue, sigma=ymn)
         s_eigval = 1/(2*mode_n + 1)
         s_eoc_rec.add_data_point(qbx.h_max, rel_err(s_sigma, s_eigval*ymn))
@@ -319,7 +319,8 @@ def test_sphere_eigenvalues(ctx_getter, mode_m, mode_n, qbx_order,
 
         # {{{ double layer
 
-        d_sigma_op = bind(qbx, sym.D(lap_knl, sym.var("sigma")))
+        d_sigma_op = bind(qbx,
+                sym.D(lap_knl, sym.var("sigma"), qbx_forced_limit="avg"))
         d_sigma = d_sigma_op(queue=queue, sigma=ymn)
         d_eigval = -1/(2*(2*mode_n + 1))
         d_eoc_rec.add_data_point(qbx.h_max, rel_err(d_sigma, d_eigval*ymn))
@@ -328,7 +329,8 @@ def test_sphere_eigenvalues(ctx_getter, mode_m, mode_n, qbx_order,
 
         # {{{ S'
 
-        sp_sigma_op = bind(qbx, sym.Sp(lap_knl, sym.var("sigma")))
+        sp_sigma_op = bind(qbx,
+                 sym.Sp(lap_knl, sym.var("sigma"), qbx_forced_limit="avg"))
         sp_sigma = sp_sigma_op(queue=queue, sigma=ymn)
         sp_eigval = -1/(2*(2*mode_n + 1))
         sp_eoc_rec.add_data_point(qbx.h_max, rel_err(sp_sigma, sp_eigval*ymn))
@@ -337,7 +339,8 @@ def test_sphere_eigenvalues(ctx_getter, mode_m, mode_n, qbx_order,
 
         # {{{ D'
 
-        dp_sigma_op = bind(qbx, sym.Dp(lap_knl, sym.var("sigma")))
+        dp_sigma_op = bind(qbx,
+                sym.Dp(lap_knl, sym.var("sigma"), qbx_forced_limit="avg"))
         dp_sigma = dp_sigma_op(queue=queue, sigma=ymn)
         dp_eigval = -(mode_n*(mode_n+1))/(2*mode_n + 1)
         dp_eoc_rec.add_data_point(qbx.h_max, rel_err(dp_sigma, dp_eigval*ymn))
diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py
index c376fdf3e5072aba8604147c70a4f5d918c15f33..6168ae249f712211fd6b7c5e7984edc91165dd62 100644
--- a/test/test_layer_pot_identity.py
+++ b/test/test_layer_pot_identity.py
@@ -159,9 +159,14 @@ class ZeroCalderonExpr(object):
 
         u_sym = sym.var("u")
 
+        from functools import partial
+        S = partial(sym.S, qbx_forced_limit=+1)
+        Dp = partial(sym.Dp, qbx_forced_limit="avg")
+        Sp = partial(sym.Sp, qbx_forced_limit="avg")
+
         return (
-                    -sym.Dp(kernel, sym.S(kernel, u_sym))
-                    - 0.25*u_sym + sym.Sp(kernel, sym.Sp(kernel, u_sym))
+                    -Dp(kernel, S(kernel, u_sym))
+                    - 0.25*u_sym + Sp(kernel, sym.Sp(kernel, u_sym))
                     )
 
     order_drop = 1
@@ -259,12 +264,12 @@ class DynamicTestCase(object):
         DynamicTestCase(geom, GreenExpr(), 0, fmm_backend="fmmlib"),
         DynamicTestCase(geom, GreenExpr(), 1.2, fmm_backend="fmmlib"),
         ]])
-def test_identity_convergence(ctx_getter,  case, visualize=False):
+def test_identity_convergence(ctx_factory,  case, visualize=False):
     logging.basicConfig(level=logging.INFO)
 
     case.check()
 
-    cl_ctx = ctx_getter()
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     # prevent cache 'splosion
diff --git a/test/test_maxwell.py b/test/test_maxwell.py
index fa9383f4821a179c6ad088608e31d0ed88486c0b..69bdc81e19193d3e62a932c20e7e362db96fea93 100644
--- a/test/test_maxwell.py
+++ b/test/test_maxwell.py
@@ -218,14 +218,14 @@ class EHField(object):
     #tc_int,
     tc_ext,
     ])
-def test_pec_mfie_extinction(ctx_getter, case, visualize=False):
+def test_pec_mfie_extinction(ctx_factory, case, visualize=False):
     """For (say) is_interior=False (the 'exterior' MFIE), this test verifies
     extinction of the combined (incoming + scattered) field on the interior
     of the scatterer.
     """
     logging.basicConfig(level=logging.INFO)
 
-    cl_ctx = ctx_getter()
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     np.random.seed(12)
diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py
index b65b5263f20173c5911c1077c2ae46e6d8bd01cc..811a912de0ebdf5cc1813ff6e063b1bf1a506d51 100644
--- a/test/test_scalar_int_eq.py
+++ b/test/test_scalar_int_eq.py
@@ -885,10 +885,10 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize):
     ])
 # Sample test run:
 # 'test_integral_equation(cl._csc, EllipseIntEqTestCase(0, "dirichlet", +1), visualize=True)'  # noqa: E501
-def test_integral_equation(ctx_getter, case, visualize=False):
+def test_integral_equation(ctx_factory, case, visualize=False):
     logging.basicConfig(level=logging.INFO)
 
-    cl_ctx = ctx_getter()
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     if USE_SYMENGINE and case.fmm_backend is None:
diff --git a/test/too_slow_test_helmholtz.py b/test/too_slow_test_helmholtz.py
index 25eeb075ee2091110a6de0199938835d68bf2b2f..52c304ebd6dd61e4c0392138b3413c67a5c64ff8 100644
--- a/test/too_slow_test_helmholtz.py
+++ b/test/too_slow_test_helmholtz.py
@@ -386,8 +386,8 @@ def run_dielectric_test(cl_ctx, queue, nelements, qbx_order,
     "tm",
     "tem",
     ])
-def test_dielectric(ctx_getter, qbx_order, op_class, mode, visualize=False):
-    cl_ctx = ctx_getter()
+def test_dielectric(ctx_factory, qbx_order, op_class, mode, visualize=False):
+    cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
     import logging