Fix some more deprecations from primitives

examples/cost.py

@@ -115,7 +115,7 @@ def get_test_density(actx, density_discr):
 
 def calibrate_cost_model(ctx):
     queue = cl.CommandQueue(ctx)
-    actx = PyOpenCLArrayContext(queue, force_device_scalars=True)
+    actx = PyOpenCLArrayContext(queue)
     cost_model = QBXCostModel()
 
     model_results = []
@@ -153,7 +153,7 @@ def calibrate_cost_model(ctx):
 
 def test_cost_model(ctx, calibration_params):
     queue = cl.CommandQueue(ctx)
-    actx = PyOpenCLArrayContext(queue, force_device_scalars=True)
+    actx = PyOpenCLArrayContext(queue)
     cost_model = QBXCostModel()
 
     for lpot_source in test_geometries(actx):

examples/fmm-error.py

@@ -15,7 +15,7 @@ def main():
     import pyopencl as cl
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue, force_device_scalars=True)
+    actx = PyOpenCLArrayContext(queue)
 
     target_order = 16
     qbx_order = 3

examples/helmholtz-dirichlet.py

@@ -29,7 +29,7 @@ def main(mesh_name="ellipse", visualize=False):
     import pyopencl as cl
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue, force_device_scalars=True)
+    actx = PyOpenCLArrayContext(queue)
 
     from meshmode.mesh.generation import ellipse, make_curve_mesh
     from functools import partial

examples/laplace-dirichlet-3d.py

@@ -27,7 +27,7 @@ def main(mesh_name="torus", visualize=False):
     import pyopencl as cl
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue, force_device_scalars=True)
+    actx = PyOpenCLArrayContext(queue)
 
     if mesh_name == "torus":
         rout = 10

examples/laplace-dirichlet-simple.py

@@ -28,7 +28,7 @@ def main(mesh_name="starfish", visualize=False):
     import pyopencl as cl
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue, force_device_scalars=True)
+    actx = PyOpenCLArrayContext(queue)
 
     from meshmode.mesh.generation import ellipse, make_curve_mesh, starfish
     from functools import partial

examples/layerpot-3d.py

@@ -21,7 +21,7 @@ def main(mesh_name="ellipsoid"):
     import pyopencl as cl
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue, force_device_scalars=True)
+    actx = PyOpenCLArrayContext(queue)
 
     if mesh_name == "ellipsoid":
         cad_file_name = "geometries/ellipsoid.step"

examples/layerpot.py

@@ -23,7 +23,7 @@ def main(curve_fn=starfish, visualize=True):
     import pyopencl as cl
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue, force_device_scalars=True)
+    actx = PyOpenCLArrayContext(queue)
 
     from meshmode.mesh.generation import make_curve_mesh
     mesh = make_curve_mesh(

examples/scaling-study.py

@@ -58,7 +58,7 @@ def timing_run(nx, ny, visualize=False):
     import pyopencl as cl
     cl_ctx = cl.create_some_context()
     queue = cl.CommandQueue(cl_ctx)
-    actx = PyOpenCLArrayContext(queue, force_device_scalars=True)
+    actx = PyOpenCLArrayContext(queue)
 
     mesh = make_mesh(nx=nx, ny=ny, visualize=visualize)
 

pytential/qbx/target_specific/impl.pyx

@@ -129,7 +129,7 @@ def h3dall_wrapper(nterms, z, scale, hs, hders):
 
 # {{{ helpers
 
-cdef void legvals(double x, int n, double[] vals, double[] derivs) nogil:
+cdef void legvals(double x, int n, double[] vals, double[] derivs) noexcept nogil:
     """Compute the values of the Legendre polynomial up to order n at x.
     Optionally, if derivs is non-NULL, compute the values of the derivative too.
 
@@ -189,7 +189,7 @@ cdef void ts_helmholtz_precompute(
         double complex k,
         double complex[] jvals,
         double complex[] jderivs,
-        double *jscale) nogil:
+        double *jscale) noexcept nogil:
     """Evaluate the source-invariant Bessel terms of the Helmholtz target-specific
     expansion."""
 
@@ -224,7 +224,7 @@ cdef double complex ts_laplace_s(
         double[3] center,
         double[3] target,
         double complex charge,
-        int order) nogil:
+        int order) noexcept nogil:
     """Evaluate the target-specific expansion of the Laplace single-layer kernel."""
 
     cdef:
@@ -278,7 +278,7 @@ cdef void ts_laplace_sp(
         double[3] center,
         double[3] target,
         double complex charge,
-        int order) nogil:
+        int order) noexcept nogil:
     """Evaluate the target-specific expansion of the gradient of the Laplace
     single-layer kernel."""
 
@@ -324,7 +324,7 @@ cdef double complex ts_laplace_d(
         double[3] target,
         double[3] dipole,
         double complex dipstr,
-        int order) nogil:
+        int order) noexcept nogil:
     """Evaluate the target-specific expansion of the Laplace double-layer kernel."""
 
     cdef:
@@ -370,7 +370,7 @@ cdef double complex ts_helmholtz_s(
         int order,
         double complex k,
         double complex[] jvals,
-        double jscale) nogil:
+        double jscale) noexcept nogil:
     """Evaluate the target-specific expansion of the Helmholtz single-layer
     kernel."""
 
@@ -429,7 +429,7 @@ cdef void ts_helmholtz_sp(
         double complex k,
         double complex[] jvals,
         double complex[] jderivs,
-        double jscale) nogil:
+        double jscale) noexcept nogil:
     """Evaluate the target-specific expansion of the gradient of the Helmholtz
     single-layer kernel."""
 
@@ -501,7 +501,7 @@ cdef double complex ts_helmholtz_d(
         int order,
         double complex k,
         double complex[] jvals,
-        double jscale) nogil:
+        double jscale) noexcept nogil:
     """Evaluate the target-specific expansion of the Helmholtz double-layer
     kernel."""
 

pytential/symbolic/execution.py

@@ -249,7 +249,7 @@ def map_interpolation(self, expr):
 
         if isinstance(operand,
                 (*self.array_context.array_types,
-                 list, np.ndarray, DOFArray)):
+                 list, tuple, np.ndarray, DOFArray)):
             conn = self.places.get_connection(expr.from_dd, expr.to_dd)
             return conn(operand)
         elif isinstance(operand, int | float | complex | np.number):

pytential/symbolic/mappers.py

@@ -296,8 +296,8 @@ class LocationTagger(CSECachingMapperMixin[Expression, []],
     """Used internally by :class:`ToTargetTagger`."""
 
     def __init__(self, default_target, default_source):
-        self.default_source = default_source
-        self.default_target = default_target
+        self.default_source = prim.as_dofdesc(default_source)
+        self.default_target = prim.as_dofdesc(default_target)
 
     def map_common_subexpression_uncached(self, expr) -> Expression:
         # Mypy 1.13 complains about this:

pytential/symbolic/primitives.py

@@ -104,6 +104,8 @@
 
 .. autoclass:: ArithmeticExpressionT
 
+.. autoclass:: QBXForcedLimit
+
 .. class:: P
 
     See :class:`pytools.P`
@@ -329,6 +331,10 @@
 
 __all__ = (
     "Expression",
+    "Operand",
+    "ArithmeticExpressionT",
+    "QBXForcedLimit",
+    "for_each_expression",
 
     "ErrorExpression",
 
@@ -395,7 +401,7 @@ def make_stringifier(self, originating_stringifier=None):
 
 Operand: TypeAlias = (
     ArithmeticExpression | np.ndarray[Any, np.dtype[Any]] | MultiVector)
-QBXForcedLimit = int | Literal["avg"] | None
+QBXForcedLimit: TypeAlias = int | Literal["avg"] | None
 
 # NOTE: this will likely live in pymbolic at some point, but for now we take it!
 ArithmeticExpressionT = TypeVar("ArithmeticExpressionT", bound=ArithmeticExpression)
@@ -611,14 +617,6 @@ def __new__(cls,
                 operand: Operand | None = None,
                 dofdesc: DOFDescriptor | None = None,
                 ) -> "NumReferenceDerivative":
-        if isinstance(ref_axes, int):
-            warn(f"Passing an 'int' as 'ref_axes' to {cls.__name__!r} "
-                 "is deprecated and will result in an error in 2025. Pass the "
-                 "well-formatted tuple '((ref_axes, 1),)' instead.",
-                 DeprecationWarning, stacklevel=2)
-
-            ref_axes = ((ref_axes, 1),)
-
         if isinstance(operand, np.ndarray | MultiVector):
             warn(f"Passing {type(operand)} directly to {cls.__name__!r} "
                  "is deprecated and will result in an error from 2025. Use "
@@ -641,6 +639,14 @@ def __init__(self,
                  ref_axes: tuple[tuple[int, int], ...],
                  operand: ArithmeticExpression,
                  dofdesc: DOFDescriptorLike) -> None:
+        if isinstance(ref_axes, int):
+            warn(f"Passing an 'int' as 'ref_axes' to {type(self).__name__!r} "
+                 "is deprecated and will result in an error in 2025. Pass the "
+                 "well-formatted tuple '((ref_axes, 1),)' instead.",
+                 DeprecationWarning, stacklevel=2)
+
+            ref_axes = ((ref_axes, 1),)
+
         if not isinstance(ref_axes, tuple):
             raise ValueError(f"'ref_axes' must be a tuple: {type(ref_axes)}")
 
@@ -856,12 +862,13 @@ def shape_operator(ambient_dim, dim=None, dofdesc=None):
 
 def _element_size(ambient_dim, dim=None, dofdesc=None):
     # A broken quasi-1D approximation of 1D element size. Do not use.
-
     if dim is None:
         dim = ambient_dim - 1
 
+    dofdesc = as_dofdesc(dofdesc)
     return elementwise_sum(
-            area_element(ambient_dim=ambient_dim, dim=dim) * QWeight(),
+            area_element(ambient_dim=ambient_dim, dim=dim, dofdesc=dofdesc)
+            * QWeight(dofdesc),
             dofdesc)**(1/dim)
 
 
@@ -1137,10 +1144,10 @@ def expansion_centers(ambient_dim, side, dim=None, dofdesc=None):
 
 
 def interleaved_expansion_centers(ambient_dim, dim=None, dofdesc=None):
-    centers = [
+    centers = (
             expansion_centers(ambient_dim, -1, dim=dim, dofdesc=dofdesc),
             expansion_centers(ambient_dim, +1, dim=dim, dofdesc=dofdesc)
-            ]
+            )
 
     source = as_dofdesc(dofdesc)
     target = source.copy(granularity=GRANULARITY_CENTER)
@@ -1167,8 +1174,9 @@ def h_min(ambient_dim, dim=None, dofdesc=None):
 
 def weights_and_area_elements(ambient_dim, dim=None, dofdesc=None):
     """Combines :func:`area_element` and :class:`QWeight`."""
-
+    dofdesc = as_dofdesc(dofdesc)
     area = area_element(ambient_dim, dim=dim, dofdesc=dofdesc)
+
     return cse(area * QWeight(dofdesc=dofdesc),
             "weights_area_elements",
             cse_scope.DISCRETIZATION)