diff --git a/doc/conf.py b/doc/conf.py index ee91574f1..31bc6933f 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -30,4 +30,5 @@ "pytools": ("https://documen.tician.de/pytools", None), "scipy": ("https://docs.scipy.org/doc/scipy", None), "sumpy": ("https://documen.tician.de/sumpy", None), + "sympy": ("https://docs.sympy.org/latest/", None), } diff --git a/pytential/linalg/direct_solver_symbolic.py b/pytential/linalg/direct_solver_symbolic.py index c5f523930..fc6322f5f 100644 --- a/pytential/linalg/direct_solver_symbolic.py +++ b/pytential/linalg/direct_solver_symbolic.py @@ -102,10 +102,12 @@ def map_int_g(self, expr): if name not in source_args } - return expr.copy(target_kernel=target_kernel, - source_kernels=source_kernels, - densities=self.rec(expr.densities), - kernel_arguments=kernel_arguments) + from dataclasses import replace + return replace(expr, + target_kernel=target_kernel, + source_kernels=source_kernels, + densities=self.rec(expr.densities), + kernel_arguments=kernel_arguments) # }}} diff --git a/pytential/linalg/skeletonization.py b/pytential/linalg/skeletonization.py index bfe252dab..e9bb73782 100644 --- a/pytential/linalg/skeletonization.py +++ b/pytential/linalg/skeletonization.py @@ -29,6 +29,7 @@ from arraycontext import PyOpenCLArrayContext, Array from pytential import GeometryCollection, sym +from pytential.symbolic.matrix import ClusterMatrixBuilderBase from pytential.linalg.utils import IndexList, TargetAndSourceClusterList from pytential.linalg.proxy import ProxyGeneratorBase, ProxyClusterGeometryData from pytential.linalg.direct_solver_symbolic import ( @@ -136,7 +137,7 @@ def prg(): lang_version=lp.MOST_RECENT_LANGUAGE_VERSION, ) - return knl + return knl.executor(actx.context) waa = bind(places, sym.weights_and_area_elements( places.ambient_dim, dofdesc=domain))(actx) @@ -253,17 +254,17 @@ class SkeletonizationWrangler: domains: tuple[sym.DOFDescriptor, ...] context: dict[str, Any] - neighbor_cluster_builder: Callable[..., np.ndarray] + neighbor_cluster_builder: type[ClusterMatrixBuilderBase] # target skeletonization weighted_targets: bool target_proxy_exprs: np.ndarray - proxy_target_cluster_builder: Callable[..., np.ndarray] + proxy_target_cluster_builder: type[ClusterMatrixBuilderBase] # source skeletonization weighted_sources: bool source_proxy_exprs: np.ndarray - proxy_source_cluster_builder: Callable[..., np.ndarray] + proxy_source_cluster_builder: type[ClusterMatrixBuilderBase] @property def nrows(self) -> int: @@ -386,9 +387,9 @@ def make_skeletonization_wrangler( # internal _weighted_proxy: bool | tuple[bool, bool] | None = None, - _proxy_source_cluster_builder: Callable[..., np.ndarray] | None = None, - _proxy_target_cluster_builder: Callable[..., np.ndarray] | None = None, - _neighbor_cluster_builder: Callable[..., np.ndarray] | None = None, + _proxy_source_cluster_builder: type[ClusterMatrixBuilderBase] | None = None, + _proxy_target_cluster_builder: type[ClusterMatrixBuilderBase] | None = None, + _neighbor_cluster_builder: type[ClusterMatrixBuilderBase] | None = None, ) -> SkeletonizationWrangler: if context is None: context = {} @@ -396,13 +397,14 @@ def make_skeletonization_wrangler( # {{{ setup expressions try: - exprs = list(exprs) + lpot_exprs = list(exprs) except TypeError: - exprs = [exprs] + lpot_exprs = [exprs] try: input_exprs = list(input_exprs) except TypeError: + assert not isinstance(input_exprs, Sequence) input_exprs = [input_exprs] from pytential.symbolic.execution import _prepare_auto_where, _prepare_domains @@ -410,11 +412,11 @@ def make_skeletonization_wrangler( auto_where = _prepare_auto_where(auto_where, places) domains = _prepare_domains(len(input_exprs), places, domains, auto_where[0]) - exprs = prepare_expr(places, exprs, auto_where) + prepared_lpot_exprs = prepare_expr(places, lpot_exprs, auto_where) source_proxy_exprs = prepare_proxy_expr( - places, exprs, (auto_where[0], PROXY_SKELETONIZATION_TARGET)) + places, prepared_lpot_exprs, (auto_where[0], PROXY_SKELETONIZATION_TARGET)) target_proxy_exprs = prepare_proxy_expr( - places, exprs, (PROXY_SKELETONIZATION_SOURCE, auto_where[1])) + places, prepared_lpot_exprs, (PROXY_SKELETONIZATION_SOURCE, auto_where[1])) # }}} @@ -449,7 +451,7 @@ def make_skeletonization_wrangler( return SkeletonizationWrangler( # operator - exprs=exprs, + exprs=prepared_lpot_exprs, input_exprs=tuple(input_exprs), domains=tuple(domains), context=context, diff --git a/pytential/linalg/utils.py b/pytential/linalg/utils.py index 17895b125..db6da7c35 100644 --- a/pytential/linalg/utils.py +++ b/pytential/linalg/utils.py @@ -282,7 +282,7 @@ def prg(): lang_version=MOST_RECENT_LANGUAGE_VERSION) knl = lp.split_iname(knl, "icluster", 128, outer_tag="g.0") - return knl + return knl.executor(actx.context) @memoize_in(mindex, (make_index_cluster_cartesian_product, "index_product")) def _product(): diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 092b785bd..705c33302 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -414,8 +414,8 @@ def preprocess_optemplate(self, name, discretizations, expr): def op_group_features(self, expr): from pytential.utils import sort_arrays_together result = ( - expr.source, *sort_arrays_together(expr.source_kernels, - expr.densities, key=str) + expr.source, + *sort_arrays_together(expr.source_kernels, expr.densities, key=str) ) return result diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py index ea4bef3b9..7491aee1d 100644 --- a/pytential/qbx/refinement.py +++ b/pytential/qbx/refinement.py @@ -268,7 +268,7 @@ def element_prop_threshold_checker(self): lang_version=MOST_RECENT_LANGUAGE_VERSION) knl = lp.split_iname(knl, "ielement", 128, inner_tag="l.0", outer_tag="g.0") - return knl + return knl.executor(self.array_context.context) def get_wrangler(self): return RefinerWrangler(self.array_context, self) @@ -388,8 +388,8 @@ def check_sufficient_source_quadrature_resolution(self, sym.ElementwiseMax( sym._source_danger_zone_radii( stage2_density_discr.ambient_dim, - dofdesc=sym.QBX_SOURCE_STAGE2), - dofdesc=sym.GRANULARITY_ELEMENT) + dofdesc=sym.as_dofdesc(sym.QBX_SOURCE_STAGE2)), + dofdesc=sym.as_dofdesc(sym.GRANULARITY_ELEMENT)) )(self.array_context), self.array_context) unwrap_args = AreaQueryElementwiseTemplate.unwrap_args @@ -633,7 +633,8 @@ def _refine_qbx_stage1(lpot_source, density_discr, quad_resolution_by_element = bind(stage1_density_discr, sym.ElementwiseMax( sym._quad_resolution(stage1_density_discr.ambient_dim), - dofdesc=sym.GRANULARITY_ELEMENT))(actx) + dofdesc=sym.as_dofdesc(sym.GRANULARITY_ELEMENT) + ))(actx) violates_kernel_length_scale = \ wrangler.check_element_prop_threshold( @@ -653,7 +654,8 @@ def _refine_qbx_stage1(lpot_source, density_discr, scaled_max_curvature_by_element = bind(stage1_density_discr, sym.ElementwiseMax( sym._scaled_max_curvature(stage1_density_discr.ambient_dim), - dofdesc=sym.GRANULARITY_ELEMENT))(actx) + dofdesc=sym.as_dofdesc(sym.GRANULARITY_ELEMENT) + ))(actx) violates_scaled_max_curv = \ wrangler.check_element_prop_threshold( diff --git a/pytential/qbx/target_assoc.py b/pytential/qbx/target_assoc.py index d3345e83b..56d3bd3d5 100644 --- a/pytential/qbx/target_assoc.py +++ b/pytential/qbx/target_assoc.py @@ -810,7 +810,7 @@ def make_target_flags(self, target_discrs_and_qbx_sides): return target_flags def make_default_target_association(self, ntargets): - target_to_center = self.array_context.zeros(ntargets, dtype=np.int32) + target_to_center = self.array_context.np.zeros(ntargets, dtype=np.int32) target_to_center.fill(-1) target_to_center.finish() diff --git a/pytential/symbolic/compiler.py b/pytential/symbolic/compiler.py index f303b5741..0f2bba39f 100644 --- a/pytential/symbolic/compiler.py +++ b/pytential/symbolic/compiler.py @@ -27,7 +27,7 @@ import numpy as np -from pymbolic.primitives import cse_scope, Expression, Variable +from pymbolic.primitives import cse_scope, Expression, Variable, Subscript from sumpy.kernel import Kernel from pytential.symbolic.primitives import ( @@ -44,6 +44,7 @@ class Statement: .. attribute:: exprs .. attribute:: priority """ + names: list[str] exprs: list[Expression] priority: int @@ -52,7 +53,7 @@ def get_assignees(self) -> set[str]: raise NotImplementedError( f"get_assignees for '{self.__class__.__name__}'") - def get_dependencies(self, dep_mapper: DependencyMapper) -> set[Expression]: + def get_dependencies(self, dep_mapper: DependencyMapper) -> set[Variable]: raise NotImplementedError( f"get_dependencies for '{self.__class__.__name__}'") @@ -80,14 +81,19 @@ def __post_init__(self): def get_assignees(self): return set(self.names) - def get_dependencies(self, dep_mapper: DependencyMapper) -> set[Expression]: + def get_dependencies(self, dep_mapper: DependencyMapper) -> set[Variable]: from operator import or_ - deps = reduce(or_, (dep_mapper(expr) for expr in self.exprs)) + all_deps = reduce(or_, (dep_mapper(expr) for expr in self.exprs)) + + deps: set[Variable] = set() + for dep in all_deps: + if isinstance(dep, Variable): + if dep.name not in self.names: + deps.add(dep) + else: + raise TypeError(f"Unsupported dependency type: {type(dep)}") - return { - dep - for dep in deps - if dep.name not in self.names} + return deps def __str__(self): comment = self.comment @@ -189,13 +195,16 @@ class ComputePotential(Statement): def get_assignees(self): return {o.name for o in self.outputs} - def get_dependencies(self, dep_mapper: DependencyMapper) -> set[Expression]: - result = set(dep_mapper(self.densities[0])) - for density in self.densities[1:]: - result.update(dep_mapper(density)) + def get_dependencies(self, dep_mapper: DependencyMapper) -> set[Variable]: + from itertools import chain - for arg_expr in self.kernel_arguments.values(): - result.update(dep_mapper(arg_expr)) + result: set[Variable] = set() + for expr in chain(self.densities, self.kernel_arguments.values()): + for dep in dep_mapper(expr): + if isinstance(dep, Variable): + result.add(dep) + else: + raise TypeError(f"Unsupported dependency type: {type(dep)}") return result @@ -546,9 +555,7 @@ def make_assign( def assign_to_new_var( self, expr: Expression, priority: int = 0, prefix: str | None = None, - ) -> Variable: - from pymbolic.primitives import Subscript - + ) -> Variable | Subscript: # Observe that the only things that can be legally subscripted # are variables. All other expressions are broken down into # their scalar components. diff --git a/pytential/symbolic/dof_desc.py b/pytential/symbolic/dof_desc.py index b7190576f..044eeb685 100644 --- a/pytential/symbolic/dof_desc.py +++ b/pytential/symbolic/dof_desc.py @@ -275,6 +275,8 @@ def as_dofdesc(desc: DOFDescriptorLike) -> DOFDescriptor: # {{{ type annotations +DEFAULT_DOFDESC = DOFDescriptor() + DiscretizationStages = ( type[QBX_SOURCE_STAGE1] | type[QBX_SOURCE_STAGE2] diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index e533c1f0c..13f3e78d1 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -20,6 +20,7 @@ THE SOFTWARE. """ +from dataclasses import replace from functools import reduce from pymbolic.mapper.stringifier import ( @@ -28,7 +29,7 @@ from pymbolic.mapper import ( Mapper, CachedMapper, - CSECachingMapperMixin + CSECachingMapperMixin, ) from pymbolic.mapper.dependency import ( DependencyMapper as DependencyMapperBase) @@ -50,6 +51,7 @@ as DerivativeSourceFinderBase, GraphvizMapper as GraphvizMapperBase) +from pymbolic.typing import ExpressionT import pytential.symbolic.primitives as prim @@ -72,7 +74,7 @@ def rec_int_g_arguments(mapper, expr): # {{{ IdentityMapper -class IdentityMapper(IdentityMapperBase): +class IdentityMapper(IdentityMapperBase[[]]): def map_node_sum(self, expr): operand = self.rec(expr.operand) if operand is expr.operand: @@ -130,9 +132,7 @@ def map_int_g(self, expr): if not changed: return expr - return expr.copy( - densities=densities, - kernel_arguments=kernel_arguments) + return replace(expr, densities=densities, kernel_arguments=kernel_arguments) def map_interpolation(self, expr): operand = self.rec(expr.operand) @@ -262,10 +262,7 @@ def map_int_g(self, expr): if not changed: return expr - return expr.copy( - densities=densities, - kernel_arguments=kernel_arguments, - ) + return replace(expr, densities=densities, kernel_arguments=kernel_arguments) def map_common_subexpression(self, expr): child = self.rec(expr.child) @@ -294,15 +291,20 @@ def flatten(expr): # {{{ LocationTagger -class LocationTagger(CSECachingMapperMixin, IdentityMapper): +class LocationTagger(CSECachingMapperMixin[ExpressionT, []], + IdentityMapper): """Used internally by :class:`ToTargetTagger`.""" def __init__(self, default_target, default_source): self.default_source = default_source self.default_target = default_target - map_common_subexpression_uncached = \ - IdentityMapper.map_common_subexpression + def map_common_subexpression_uncached(self, expr) -> ExpressionT: + # Mypy 1.13 complains about this: + # error: Too few arguments for "map_common_subexpression" of "IdentityMapper" [call-arg] # noqa: E501 + # error: Argument 1 to "map_common_subexpression" of "IdentityMapper" has incompatible type "LocationTagger"; expected "IdentityMapper[P]" [arg-type] # noqa: E501 + # This seems spurious? + return IdentityMapper.map_common_subexpression(self, expr) # type: ignore[arg-type, call-arg] def _default_dofdesc(self, dofdesc): if dofdesc.geometry is None: @@ -539,9 +541,9 @@ def map_product(self, expr): def map_int_g(self, expr): from sumpy.kernel import AxisTargetDerivative - return expr.copy( - target_kernel=AxisTargetDerivative( - self.ambient_axis, expr.target_kernel)) + + target_kernel = AxisTargetDerivative(self.ambient_axis, expr.target_kernel) + return replace(expr, target_kernel=target_kernel) class DerivativeSourceAndNablaComponentCollector( @@ -587,15 +589,15 @@ def map_int_g(self, expr): raise ValueError( "Unregularized evaluation does not support one-sided limits") - expr = expr.copy( - qbx_forced_limit=None, - densities=self.rec(expr.densities), - kernel_arguments={ - name: self.rec(arg_expr) - for name, arg_expr in expr.kernel_arguments.items() - }) - - return expr + return replace( + expr, + qbx_forced_limit=None, + densities=self.rec(expr.densities), + kernel_arguments={ + name: self.rec(arg_expr) + for name, arg_expr in expr.kernel_arguments.items() + } + ) # }}} @@ -643,7 +645,7 @@ def map_num_reference_derivative(self, expr): def map_int_g(self, expr): if expr.target.discr_stage is None: - expr = expr.copy(target=expr.target.to_stage1()) + expr = replace(expr, target=expr.target.to_stage1()) if expr.source.discr_stage is not None: return expr @@ -655,8 +657,9 @@ def map_int_g(self, expr): from_dd = expr.source.to_stage1() to_dd = from_dd.to_quad_stage2() - densities = [prim.interp(from_dd, to_dd, self.rec(density)) for - density in expr.densities] + densities = tuple( + prim.interp(from_dd, to_dd, self.rec(density)) for + density in expr.densities) from_dd = from_dd.copy(discr_stage=self.from_discr_stage) kernel_arguments = { @@ -664,7 +667,8 @@ def map_int_g(self, expr): self.rec(self.tagger(arg_expr))) for name, arg_expr in expr.kernel_arguments.items()} - return expr.copy( + return replace( + expr, densities=densities, kernel_arguments=kernel_arguments, source=to_dd) @@ -695,7 +699,8 @@ def map_int_g(self, expr): is_self = source_discr is target_discr - expr = expr.copy( + expr = replace( + expr, densities=self.rec(expr.densities), kernel_arguments={ name: self.rec(arg_expr) @@ -724,8 +729,8 @@ def map_int_g(self, expr): if expr.qbx_forced_limit == "avg": return 0.5*( - expr.copy(qbx_forced_limit=+1) - + expr.copy(qbx_forced_limit=-1)) + replace(expr, qbx_forced_limit=+1) + + replace(expr, qbx_forced_limit=-1)) else: return expr diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 6531e207c..969e1696b 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -20,27 +20,29 @@ THE SOFTWARE. """ -from sys import intern +from dataclasses import field from warnings import warn from functools import partial -from typing import ClassVar +from typing import Any, Union, Literal import numpy as np +import modepy as mp from pymbolic.primitives import ( # noqa: N813 Expression as ExpressionBase, Variable, Variable as var, cse_scope as cse_scope_base, - make_common_subexpression as cse) + make_common_subexpression as cse, + expr_dataclass) from pymbolic.geometric_algebra import MultiVector, componentwise from pymbolic.geometric_algebra.primitives import ( NablaComponent, Derivative as DerivativeBase) from pymbolic.primitives import make_sym_vector from pytools.obj_array import make_obj_array, flat_obj_array -from sumpy.kernel import SpatialConstant +from sumpy.kernel import Kernel, SpatialConstant from pytential.symbolic.dof_desc import ( - DEFAULT_SOURCE, DEFAULT_TARGET, + DEFAULT_SOURCE, DEFAULT_TARGET, DEFAULT_DOFDESC, QBX_SOURCE_STAGE1, QBX_SOURCE_STAGE2, QBX_SOURCE_QUAD_STAGE2, GRANULARITY_NODE, GRANULARITY_CENTER, GRANULARITY_ELEMENT, DOFDescriptor, DOFDescriptorLike, @@ -87,10 +89,18 @@ :class:`~meshmode.dof_array.DOFArray` is used and otherwise :class:`~pyopencl.array.Array` is used. +.. autoclass:: Expression + :show-inheritance: + :undoc-members: + :members: mapper_method + Diagnostics ^^^^^^^^^^^ .. autoclass:: ErrorExpression + :show-inheritance: + :undoc-members: + :members: mapper_method .. _placeholders: @@ -98,7 +108,12 @@ ^^^^^^^^^^^^ .. autoclass:: var + .. autoclass:: SpatialConstant + :show-inheritance: + :undoc-members: + :members: mapper_method + .. autofunction:: make_sym_mv .. autofunction:: make_sym_surface_mv @@ -135,8 +150,21 @@ Discretization properties ^^^^^^^^^^^^^^^^^^^^^^^^^ +.. autoclass:: DiscretizationProperty + :show-inheritance: + :undoc-members: + :members: mapper_method + .. autoclass:: IsShapeClass + :show-inheritance: + :undoc-members: + :members: mapper_method + .. autoclass:: QWeight + :show-inheritance: + :undoc-members: + :members: mapper_method + .. autofunction:: nodes .. autofunction:: parametrization_derivative .. autofunction:: parametrization_derivative_matrix @@ -158,28 +186,66 @@ ^^^^^^^^^^^^^^^^^^^ .. autoclass:: NumReferenceDerivative + :show-inheritance: + :undoc-members: + :members: mapper_method + .. autoclass:: NodeSum + :undoc-members: + :members: mapper_method + .. autoclass:: NodeMax + :undoc-members: + :members: mapper_method + .. autoclass:: NodeMin + :undoc-members: + :members: mapper_method + .. autoclass:: ElementwiseSum + :undoc-members: + :members: mapper_method + +.. autoclass:: ElementwiseMin + :undoc-members: + :members: mapper_method + .. autoclass:: ElementwiseMax + :undoc-members: + :members: mapper_method + .. autofunction:: integral + .. autoclass:: Ones + :show-inheritance: + :undoc-members: + :members: mapper_method + .. autofunction:: ones_vec .. autofunction:: area .. autofunction:: mean + .. autoclass:: IterativeInverse + :show-inheritance: + :undoc-members: + :members: mapper_method + Operators ^^^^^^^^^ .. autoclass:: Interpolation + :show-inheritance: + :undoc-members: + :members: mapper_method + .. autofunction:: interp Geometric Calculus (based on Geometric/Clifford Algebra) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. autoclass:: Derivative + :undoc-members: Conventional Calculus ^^^^^^^^^^^^^^^^^^^^^ @@ -196,6 +262,10 @@ ^^^^^^^^^^^^^^^^ .. autoclass:: IntG + :show-inheritance: + :undoc-members: + :members: mapper_method + .. autofunction:: int_g_dsource .. autofunction:: int_g_vec @@ -226,6 +296,8 @@ """ __all__ = ( + "Expression", + "ErrorExpression", "var", "SpatialConstant", "make_sym_mv", "make_sym_surface_mv", @@ -275,6 +347,10 @@ ) +Operand = Union["Expression", np.ndarray, MultiVector] +QBXForcedLimit = int | Literal["avg"] | None + + class _NoArgSentinel: pass @@ -298,35 +374,35 @@ def array_to_tuple(ary): # }}} +@expr_dataclass() class Expression(ExpressionBase): + """A subclass of :class:`pymbolic.primitives.Expression` for use with + :mod:`pytential` mappers. + """ + def make_stringifier(self, originating_stringifier=None): from pytential.symbolic.mappers import StringifyMapper return StringifyMapper() +@expr_dataclass() class NamedIntermediateResult(Variable): - # These are inserted by the pytential 'compiler'. - pass + """Internal variables used by ``pytential.compiler``.""" +@expr_dataclass() class ErrorExpression(Expression): """An expression that, if evaluated, causes an error with the supplied *message*. - .. automethod:: __init__ + .. autoattribute:: message """ - init_arg_names = ("message",) - def __init__(self, message): - self.message = message + message: str + """The error message to raise when this expression is encountered.""" - def __getinitargs__(self): - return (self.message,) - mapper_method = intern("map_error_expression") - - -def make_sym_mv(name, num_components): +def make_sym_mv(name: str, num_components: int) -> MultiVector[Expression]: return MultiVector(make_sym_vector(name, num_components)) @@ -334,12 +410,13 @@ def make_sym_surface_mv(name, ambient_dim, dim, dofdesc=None): par_grad = parametrization_derivative_matrix(ambient_dim, dim, dofdesc) return sum( - var("%s%d" % (name, i)) - * cse(MultiVector(vec), "tangent%d" % i, cse_scope.DISCRETIZATION) + var(f"{name}{i}") + * cse(MultiVector(vec), f"tangent{i}", cse_scope.DISCRETIZATION) for i, vec in enumerate(par_grad.T)) -class Function(var): +@expr_dataclass() +class Function(Variable): def __call__(self, operand, *args, **kwargs): # If the call is handed an object array full of operands, # return an object array of the operator applied to each of the @@ -355,10 +432,10 @@ def make_op(operand_i): return var.__call__(self, operand, *args, **kwargs) +@expr_dataclass() class NumpyMathFunction(Function): """A math function named within the numpy naming convention and with numpy-like semantics.""" - pass real = NumpyMathFunction("real") @@ -391,146 +468,151 @@ class NumpyMathFunction(Function): # {{{ discretization properties +@expr_dataclass() class DiscretizationProperty(Expression): """A quantity that depends exclusively on the discretization. - .. attribute:: dofdesc + .. autoattribute:: dofdesc """ - init_arg_names: ClassVar[tuple[str, ...]] = ("dofdesc",) - - def __init__(self, dofdesc=None): - """ - :arg dofdesc: |dofdesc-blurb| - """ + dofdesc: DOFDescriptor + """The descriptor for this quantity that selects its geometry on evaluation.""" - self.dofdesc = as_dofdesc(dofdesc) + def __post_init__(self) -> None: + if not isinstance(self.dofdesc, DOFDescriptor): + warn("Passing a 'dofdesc' that is not a 'DOFDescriptor' to " + f"{type(self).__name__!r} is deprecated and will stop working " + "in 2025. Use 'as_dofdesc' to convert the descriptor.", + DeprecationWarning, stacklevel=2) - def __getinitargs__(self): - return (self.dofdesc,) + object.__setattr__(self, "dofdesc", as_dofdesc(self.dofdesc)) +@expr_dataclass(init=False) class IsShapeClass(DiscretizationProperty): """A predicate that is *True* if the elements of the discretization have a unique type that matches :attr:`shape`. - .. attribute:: shape - - A :class:`modepy.Shape` subclass. + .. autoattribute:: shape """ - init_arg_names = ("shape", "dofdesc") - - def __init__(self, shape, dofdesc) -> None: - super().__init__(dofdesc) - self.shape = shape - def __getinitargs__(self): - return (self.shape, self.dofdesc) + shape: mp.Shape + """A :class:`modepy.Shape` subclass.""" - mapper_method = intern("map_is_shape_class") + # FIXME: this is added for backwards compatibility with pre-dataclass expressions + def __init__(self, shape: mp.Shape, dofdesc: DOFDescriptorLike) -> None: + object.__setattr__(self, "shape", shape) + super().__init__(dofdesc) # type: ignore[arg-type] +@expr_dataclass() class QWeight(DiscretizationProperty): """Bare quadrature weights (without Jacobians).""" - mapper_method = intern("map_q_weight") - +@expr_dataclass(init=False) class NodeCoordinateComponent(DiscretizationProperty): + """ + .. autoattribute:: ambient_axis + """ - init_arg_names = ("ambient_axis", "dofdesc") - - def __init__(self, ambient_axis, dofdesc=None): - """ - :arg dofdesc: |dofdesc-blurb| - """ - self.ambient_axis = ambient_axis - DiscretizationProperty.__init__(self, dofdesc) - - def __getinitargs__(self): - return (self.ambient_axis, self.dofdesc) + ambient_axis: int + """The axis index this node coordinate represents, i.e. 0 for $x$, etc.""" - mapper_method = intern("map_node_coordinate_component") + # FIXME: this is added for backwards compatibility with pre-dataclass expressions + def __init__(self, ambient_axis: int, dofdesc: DOFDescriptorLike) -> None: + object.__setattr__(self, "ambient_axis", ambient_axis) + super().__init__(dofdesc) # type: ignore[arg-type] def nodes(ambient_dim, dofdesc=None): """Return a :class:`pymbolic.geometric_algebra.MultiVector` of node locations. """ - + dofdesc = as_dofdesc(dofdesc) return MultiVector( make_obj_array([ NodeCoordinateComponent(i, dofdesc) for i in range(ambient_dim)])) +@expr_dataclass() class NumReferenceDerivative(DiscretizationProperty): - """An operator that takes a derivative - of *operand* with respect to the the element - reference coordinates. + """An operator that takes a derivative of *operand* with respect to the the + element reference coordinates. + + .. autoattribute:: ref_axes + .. autoattribute:: operand """ - init_arg_names = ("ref_axes", "operand", "dofdesc") + ref_axes: tuple[tuple[int, int], ...] + """A tuple of pairs ``(axis, derivative_order)`` that define the reference + derivatives taken on the given *operand*. The tuple must be sorted with + respect to the axis index. For example, ``((0, 2), (1, 1))`` is a correct + input as it is sorted and each axis is unique. It denotes a second derivative + with respect to $x$ (0) and a first derivative with respect to $y$ (1). + """ + operand: Operand + """An operand to differentiate.""" - def __new__(cls, ref_axes=None, operand=None, dofdesc=None): + def __new__(cls, + ref_axes: int | tuple[tuple[int, int], ...] | None = None, + operand: Operand | None = None, + dofdesc: DOFDescriptor | None = None) -> "NumReferenceDerivative": # If the constructor is handed a multivector object, return an # object array of the operator applied to each of the # coefficients in the multivector. if isinstance(operand, np.ndarray): def make_op(operand_i): - return cls(ref_axes, operand_i, dofdesc=dofdesc) + return cls(ref_axes, operand_i, as_dofdesc(dofdesc)) return componentwise(make_op, operand) else: return DiscretizationProperty.__new__(cls) - def __init__(self, ref_axes, operand, dofdesc=None): - """ - :arg ref_axes: a :class:`tuple` of tuples indicating indices of - coordinate axes of the reference element to the number of derivatives - which will be taken. For example, the value ``((0, 2), (1, 1))`` - indicates that Each axis must occur at most once. The tuple must be - sorted by the axis index. - - May also be a singile integer *i*, which is viewed as equivalent - to ``((i, 1),)``. - :arg dofdesc: |dofdesc-blurb| - """ - - if isinstance(ref_axes, int): - ref_axes = ((ref_axes, 1),) + # FIXME: this is added for backwards compatibility with pre-dataclass expressions + def __init__(self, + ref_axes: tuple[tuple[int, int], ...], + operand: Expression, + dofdesc: DOFDescriptorLike) -> None: + object.__setattr__(self, "ref_axes", ref_axes) + object.__setattr__(self, "operand", operand) + super().__init__(dofdesc) # type: ignore[arg-type] - if not isinstance(ref_axes, tuple): - raise ValueError("ref_axes must be a tuple") + if isinstance(self.ref_axes, int): + warn(f"Passing an 'int' as 'ref_axes' to {type(self).__name__!r} " + "is deprecated and will be removed in 2025. Pass the " + "well-formatted tuple '((ref_axes, 1),)' instead.", + DeprecationWarning, stacklevel=2) - if tuple(sorted(ref_axes)) != ref_axes: - raise ValueError("ref_axes must be sorted") + object.__setattr__(self, "ref_axes", ((self.ref_axes, 1),)) - if len(dict(ref_axes)) != len(ref_axes): - raise ValueError("ref_axes must not contain an axis more than once") + if not isinstance(self.ref_axes, tuple): + raise ValueError(f"'ref_axes' must be a tuple: {type(self)}") - self.ref_axes = ref_axes + if tuple(sorted(self.ref_axes)) != self.ref_axes: + raise ValueError( + f"'ref_axes' must be sorted by axis index: {self.ref_axes}" + ) - self.operand = operand - DiscretizationProperty.__init__(self, dofdesc) - - def __getinitargs__(self): - return (self.ref_axes, self.operand, self.dofdesc) - - mapper_method = intern("map_num_reference_derivative") + if len(dict(self.ref_axes)) != len(self.ref_axes): + raise ValueError( + f"'ref_axes' must not contain an axis more than once: {self.ref_axes}" + ) def reference_jacobian(func, output_dim, dim, dofdesc=None): """Return a :class:`numpy.ndarray` representing the Jacobian of a vector function with respect to the reference coordinates. """ + dofdesc = as_dofdesc(dofdesc) jac = np.zeros((output_dim, dim), object) for i in range(output_dim): func_component = func[i] for j in range(dim): - jac[i, j] = NumReferenceDerivative(j, func_component, dofdesc) + jac[i, j] = NumReferenceDerivative(((j, 1),), func_component, dofdesc) return jac @@ -540,6 +622,7 @@ def parametrization_derivative_matrix(ambient_dim, dim, dofdesc=None): reference-to-global parametrization. """ + dofdesc = as_dofdesc(dofdesc) return cse( reference_jacobian( [NodeCoordinateComponent(i, dofdesc) for i in range(ambient_dim)], @@ -578,6 +661,7 @@ def area_element(ambient_dim, dim=None, dofdesc=None): def sqrt_jac_q_weight(ambient_dim, dim=None, dofdesc=None): + dofdesc = as_dofdesc(dofdesc) return cse( sqrt( area_element(ambient_dim, dim, dofdesc) @@ -591,6 +675,7 @@ def normal(ambient_dim, dim=None, dofdesc=None): # Don't be tempted to add a sign here. As it is, it produces # exterior normals for positively oriented curves and surfaces. + dofdesc = as_dofdesc(dofdesc) pder = ( pseudoscalar(ambient_dim, dim, dofdesc) / area_element(ambient_dim, dim, dofdesc)) @@ -650,6 +735,7 @@ def second_fundamental_form(ambient_dim, dim=None, dofdesc=None): if not (ambient_dim == 3 and dim == 2): raise NotImplementedError("only available for surfaces in 3D") + dofdesc = as_dofdesc(dofdesc) r = nodes(ambient_dim, dofdesc=dofdesc).as_vector() # https://en.wikipedia.org/w/index.php?title=Second_fundamental_form&oldid=821047433#Classical_notation @@ -1020,18 +1106,35 @@ def weights_and_area_elements(ambient_dim, dim=None, dofdesc=None): # {{{ operators +@expr_dataclass() class Interpolation(Expression): """Interpolate quantity from a DOF described by *from_dd* to a DOF - described by *to_dd*.""" + described by *to_dd*." + + .. autoattribute:: from_dd + .. autoattribute:: to_dd + .. autoattribute:: operand + """ - init_arg_names = ("from_dd", "to_dd", "operand") + from_dd: DOFDescriptor + """A descriptor for the geometry on which *operand* is defined.""" + to_dd: DOFDescriptor + """A descriptor for the geometry to which to interpolate *operand* to.""" + operand: Operand + """An expression or array of expressions to interpolate. Arrays are + interpolated componentwise. + """ - def __new__(cls, from_dd, to_dd, operand): + def __new__(cls, + from_dd: DOFDescriptorLike, + to_dd: DOFDescriptorLike, + operand: Operand) -> "Interpolation": from_dd = as_dofdesc(from_dd) to_dd = as_dofdesc(to_dd) if from_dd == to_dd: - return operand + # FIXME: __new__ should return a class instance + return operand # type: ignore[return-value] if isinstance(operand, np.ndarray): def make_op(operand_i): @@ -1041,26 +1144,39 @@ def make_op(operand_i): else: return Expression.__new__(cls) - def __init__(self, from_dd, to_dd, operand): - self.from_dd = as_dofdesc(from_dd) - self.to_dd = as_dofdesc(to_dd) - self.operand = operand + def __post_init__(self) -> None: + if not isinstance(self.from_dd, DOFDescriptor): + warn("Passing a 'from_dd' that is not a 'DOFDescriptor' to " + f"{type(self).__name__!r} is deprecated and will stop working " + "in 2025. Use 'as_dofdesc' to convert the descriptor.", + DeprecationWarning, stacklevel=2) + + object.__setattr__(self, "from_dd", as_dofdesc(self.from_dd)) - def __getinitargs__(self): - return (self.from_dd, self.to_dd, self.operand) + if not isinstance(self.to_dd, DOFDescriptor): + warn("Passing a 'to_dd' that is not a 'DOFDescriptor' to " + f"{type(self).__name__!r} is deprecated and will stop working " + "in 2025. Use 'as_dofdesc' to convert the descriptor.", + DeprecationWarning, stacklevel=2) - mapper_method = intern("map_interpolation") + object.__setattr__(self, "to_dd", as_dofdesc(self.to_dd)) def interp(from_dd, to_dd, operand): - return Interpolation(from_dd, to_dd, operand) + return Interpolation(as_dofdesc(from_dd), as_dofdesc(to_dd), operand) +@expr_dataclass() class SingleScalarOperandExpression(Expression): + """ + .. autoattribute:: operand + """ - init_arg_names = ("operand",) + operand: Operand + """An expression or an array on which to apply the operation.""" - def __new__(cls, operand=None): + def __new__(cls, + operand: Operand | None = None) -> "SingleScalarOperandExpression": # If the constructor is handed a multivector object, return an # object array of the operator applied to each of the # coefficients in the multivector. @@ -1073,112 +1189,134 @@ def make_op(operand_i): else: return Expression.__new__(cls) - def __init__(self, operand): - self.operand = operand - - def __getinitargs__(self): - return (self.operand,) - +@expr_dataclass() class NodeSum(SingleScalarOperandExpression): - """Implements a global sum over all discretization nodes.""" + """Bases: :class:`~pytential.symbolic.primitives.Expression`. - mapper_method = "map_node_sum" + Implements a global sum over all discretization nodes. + """ +@expr_dataclass() class NodeMax(SingleScalarOperandExpression): - """Implements a global maximum over all discretization nodes.""" + """Bases: :class:`~pytential.symbolic.primitives.Expression`. - mapper_method = "map_node_max" + Implements a global maximum over all discretization nodes. + """ +@expr_dataclass() class NodeMin(SingleScalarOperandExpression): - """Implements a global minimum over all discretization nodes.""" + """Bases: :class:`~pytential.symbolic.primitives.Expression`. - mapper_method = "map_node_min" + Implements a global minimum over all discretization nodes. + """ def integral(ambient_dim, dim, operand, dofdesc=None): """A volume integral of *operand*.""" + dofdesc = as_dofdesc(dofdesc) return NodeSum( area_element(ambient_dim, dim, dofdesc) * QWeight(dofdesc) * operand) +@expr_dataclass() class SingleScalarOperandExpressionWithWhere(Expression): + """ + .. autoattribute:: operand + .. autoattribute:: dofdesc + """ - init_arg_names = ("operand", "dofdesc") + operand: Operand + """An expression or an array on which to apply the operation.""" + dofdesc: DOFDescriptor + """The descriptor for the geometry where the *operand* is defined.""" - def __new__(cls, operand=None, dofdesc=None): + def __new__(cls, + operand: Operand | None = None, + dofdesc: DOFDescriptorLike | None = None, + ) -> "SingleScalarOperandExpressionWithWhere": # If the constructor is handed a multivector object, return an # object array of the operator applied to each of the # coefficients in the multivector. if isinstance(operand, np.ndarray | MultiVector): def make_op(operand_i): - return cls(operand_i, dofdesc) + return cls(operand_i, as_dofdesc(dofdesc)) return componentwise(make_op, operand) else: return Expression.__new__(cls) - def __init__(self, operand, dofdesc=None): - self.operand = operand - self.dofdesc = as_dofdesc(dofdesc) + def __post_init__(self) -> None: + if not isinstance(self.dofdesc, DOFDescriptor): + warn("Passing a 'dofdesc' that is not a 'DOFDescriptor' to " + f"{type(self).__name__!r} is deprecated and will stop working " + "in 2025. Use 'as_dofdesc' to convert the descriptor.", + DeprecationWarning, stacklevel=2) - def __getinitargs__(self): - return (self.operand, self.dofdesc) + object.__setattr__(self, "dofdesc", as_dofdesc(self.dofdesc)) +@expr_dataclass() class ElementwiseSum(SingleScalarOperandExpressionWithWhere): - """Returns a vector of DOFs with all entries on each element set + """Bases: :class:`~pytential.symbolic.primitives.Expression`. + + Returns a vector of DOFs with all entries on each element set to the sum of DOFs on that element. """ - mapper_method = "map_elementwise_sum" - +@expr_dataclass() class ElementwiseMin(SingleScalarOperandExpressionWithWhere): - """Returns a vector of DOFs with all entries on each element set + """Bases: :class:`~pytential.symbolic.primitives.Expression`. + + Returns a vector of DOFs with all entries on each element set to the minimum of DOFs on that element. """ - mapper_method = "map_elementwise_min" - +@expr_dataclass() class ElementwiseMax(SingleScalarOperandExpressionWithWhere): - """Returns a vector of DOFs with all entries on each element set + """Bases: :class:`~pytential.symbolic.primitives.Expression`. + + Returns a vector of DOFs with all entries on each element set to the maximum of DOFs on that element. """ - mapper_method = "map_elementwise_max" - +@expr_dataclass() class Ones(Expression): - """A DOF-vector that is constant *one* on the whole - discretization. + """A DOF-vector that is constant *one* on the whole discretization. """ - init_arg_names = ("dofdesc",) - - def __init__(self, dofdesc=None): - self.dofdesc = as_dofdesc(dofdesc) + # pylint: disable-next=invalid-field-call + dofdesc: DOFDescriptor = field(default_factory=lambda: DEFAULT_DOFDESC) + """A descriptor for the discretization where the array is defined.""" - def __getinitargs__(self): - return (self.dofdesc,) + def __post_init__(self) -> None: + if not isinstance(self.dofdesc, DOFDescriptor): + warn("Passing a 'dofdesc' that is not a 'DOFDescriptor' to " + f"{type(self).__name__!r} is deprecated and will stop working " + "in 2025. Use 'as_dofdesc' to convert the descriptor.", + DeprecationWarning, stacklevel=2) - mapper_method = intern("map_ones") + object.__setattr__(self, "dofdesc", as_dofdesc(self.dofdesc)) def ones_vec(dim, dofdesc=None): from pytools.obj_array import make_obj_array - return MultiVector( - make_obj_array(dim*[Ones(dofdesc)])) + + dofdesc = as_dofdesc(dofdesc) + return MultiVector(make_obj_array(dim*[Ones(dofdesc)])) def area(ambient_dim, dim, dofdesc=None): + dofdesc = as_dofdesc(dofdesc) return cse(integral(ambient_dim, dim, Ones(dofdesc), dofdesc), "area", cse_scope.DISCRETIZATION) @@ -1189,36 +1327,48 @@ def mean(ambient_dim, dim, operand, dofdesc=None): / area(ambient_dim, dim, dofdesc)) +@expr_dataclass() class IterativeInverse(Expression): + """A symbolic :math:`A x = b` linear solve expression. - init_arg_names = ("expression", "rhs", "variable_name", "extra_vars", "dofdesc") - - def __init__(self, expression, rhs, variable_name, extra_vars=None, - dofdesc=None): - if extra_vars is None: - extra_vars = {} - self.expression = expression - self.rhs = rhs - self.variable_name = variable_name - self.extra_vars = extra_vars - self.dofdesc = as_dofdesc(dofdesc) + .. autoattribute:: expression + .. autoattribute:: rhs + .. autoattribute:: variable_name + .. autoattribute:: extra_vars + .. autoattribute:: dofdesc + """ - def __getinitargs__(self): - return (self.expression, self.rhs, self.variable_name, - self.extra_vars, self.dofdesc) + expression: Expression + """The operator *A* used in the linear solve.""" + rhs: Expression + """The right-hand side variable used in the linear solve.""" + variable_name: str + """The name of the variable to solve for.""" + extra_vars: dict[str, Variable] + """A dictionary of additional variables required to define the operator.""" + dofdesc: DOFDescriptor + """A descriptor for the geometry on which the solution is defined.""" - def get_hash(self): - return hash((self.__class__, - self.expression, - self.rhs, - self.variable_name, - frozenset(self.extra_vars.items()), - self.dofdesc)) + def __post_init__(self) -> None: + if not isinstance(self.dofdesc, DOFDescriptor): + warn("Passing a 'dofdesc' that is not a 'DOFDescriptor' to " + f"{type(self).__name__!r} is deprecated and will stop working " + "in 2025. Use 'as_dofdesc' to convert the descriptor.", + DeprecationWarning, stacklevel=2) - mapper_method = intern("map_inverse") + object.__setattr__(self, "dofdesc", as_dofdesc(self.dofdesc)) class Derivative(DerivativeBase): + """A symbolic derivative. + + This mechanism cannot be used to take more than one derivative at a time. + + .. automethod:: __call__ + .. automethod:: dnabla + .. automethod:: resolve + """ + @property def nabla(self): raise ValueError("Derivative.nabla should not be used" @@ -1297,171 +1447,189 @@ def hashable_kernel_args(kernel_arguments): return tuple(hashable_args) +@expr_dataclass(hash=False) class IntG(Expression): r""" .. math:: - \int_\Gamma T (\sum S_k G(x-y) \sigma_k(y)) dS_y + \int_\Gamma T \left[\sum S_k[G](x-y) \sigma_k(y)\right] \,\mathrm{d} S_y where :math:`\sigma_k` is the k-th *density*, :math:`G` is a Green's function, :math:`S_k` are source derivative operators and :math:`T` is a target derivative operator. - .. attribute:: target_kernel - .. attribute:: source_kernels - .. attribute:: densities - .. attribute:: qbx_forced_limit - .. attribute:: kernel_arguments + .. autoattribute:: target_kernel + .. autoattribute:: source_kernels + .. autoattribute:: densities + .. autoattribute:: qbx_forced_limit + .. autoattribute:: source + .. autoattribute:: target + .. autoattribute:: kernel_arguments """ - init_arg_names = ("target_kernel", "source_kernels", "densities", - "qbx_forced_limit", "source", "target", "kernel_arguments") + target_kernel: Kernel + """An instance of :class:`~sumpy.kernel.Kernel` with only target dervatives + attached. This represents the target derivative operator :math:`T` above. - def __init__(self, target_kernel, source_kernels, densities, - qbx_forced_limit, source=None, target=None, - kernel_arguments=None, - **kwargs): - """*target_derivatives* and later arguments should be considered - keyword-only. - - :arg source_kernels: a tuple of instances of :class:`sumpy.kernel.Kernel` - with only source derivatives attached. k-th elements represents the - k-th source derivative operator above. - :arg target_kernel: an instance of :class:`sumpy.kernel.Kernel` with only - target dervatives attached. This represents the target derivative - operator :math:`T` above. Note that the term ``target_kernel`` is - bad as it's not a kernel and merely represents a target derivative - operator. This name will change once :mod:`sumpy` properly - supports derivative operators. This also means that the user has to - make sure that base kernels of all the kernels passed are the same. - :arg densities: a tuple of density expressions. Length of this tuple - must match the length of the source_kernels arguments. - :arg qbx_forced_limit: +1 if the output is required to originate from a - QBX center on the "+" side of the boundary. -1 for the other side. - Evaluation at a target with a value of +/- 1 in *qbx_forced_limit* - will fail if no QBX center is found. - - +2 may be used to *allow* evaluation QBX center on the "+" side of the - (but disallow evaluation using a center on the "-" side). Potential - evaluation at the target still succeeds if no applicable QBX center - is found. (-2 for the analogous behavior on the "-" side.) - - *None* may be used to avoid expressing a side preference for close - evaluation. - - ``'avg'`` may be used as a shorthand to evaluate this potential - as an average of the ``+1`` and the ``-1`` value. - - :arg kernel_arguments: A dictionary mapping named - :class:`sumpy.kernel.Kernel` arguments - (see :meth:`sumpy.kernel.Kernel.get_args` - and :meth:`sumpy.kernel.Kernel.get_source_args`) - to expressions that determine them - - :arg source: The symbolic name of the source discretization. This name - is bound to a concrete :class:`pytential.source.LayerPotentialSourceBase` - by :func:`pytential.bind`. - - :arg target: The symbolic name of the set of targets. This name gets - assigned to a concrete target set by :func:`pytential.bind`. - - *kwargs* has the same meaning as *kernel_arguments* can be used as a - more user-friendly interface. - """ - - if kernel_arguments is None: - kernel_arguments = {} - - if isinstance(kernel_arguments, tuple): - kernel_arguments = dict(kernel_arguments) - - if qbx_forced_limit not in [-1, +1, -2, +2, "avg", None]: - raise ValueError("invalid value (%s) of qbx_forced_limit" - % qbx_forced_limit) - - source_kernels = tuple(source_kernels) - densities = tuple(densities) - kernel_arg_names = set() + Note that the term ``target_kernel`` is bad as it's not a kernel and merely + represents a target derivative operator. This name will change once :mod:`sumpy` + properly supports derivative operators. This also means that the user has to + make sure that base kernels of all the kernels passed are the same. + """ + source_kernels: tuple[Kernel, ...] + """A tuple of instances of :class:`~sumpy.kernel.Kernel` with only source + derivatives attached. k-th elements represents the k-th source derivative + operator above. + """ + densities: tuple[Expression, ...] + """A tuple of density expressions. Length of this tuple must match the length + of the *source_kernels* arguments. + """ + qbx_forced_limit: QBXForcedLimit + """Limit used for the QBX expansions. Can take one of the values + + * *None*: may be used to avoid expressing a side preference for close + evaluation. + * ``+1``: if the output is required to originate from a QBX center on the "+" + side of the boundary. + * ``-1``: for the "-" side. + * ``'avg'``: may be used as a shorthand to evaluate this potential as an + average of the ``+1`` and the ``-1`` value. + * ``+2`` may be used to *allow* evaluation QBX center on the "+" side of the + (but disallow evaluation using a center on the "-" side). + * ``-2``: for the "-" side. + + Evaluation at a target with a value of ``±1`` in *qbx_forced_limit* will + fail if no QBX center is found. To allow potential evaluation at the target + to succeeds even if no applicable QBX center is found use ``±2``. + """ - for kernel in (*source_kernels, target_kernel): - for karg in (kernel.get_args() + kernel.get_source_args()): - kernel_arg_names.add(karg.loopy_arg.name) + # pylint: disable-next=invalid-field-call + source: DOFDescriptor = field(default_factory=lambda: DEFAULT_DOFDESC) + """The symbolic name of the source discretization. This name is bound to a + concrete :class:`~pytential.source.LayerPotentialSourceBase` + by :func:`pytential.bind`. + """ - from pytools import single_valued + # pylint: disable-next=invalid-field-call + target: DOFDescriptor = field(default_factory=lambda: DEFAULT_DOFDESC) + """The symbolic name of the set of targets. This name gets assigned to a + concrete target set by :func:`pytential.bind`. + """ - single_valued(kernel.get_base_kernel() for - kernel in (*source_kernels, target_kernel)) + # pylint: disable-next=invalid-field-call + kernel_arguments: dict[str, Operand] = field(default_factory=dict) + """A dictionary mapping named :class:`~sumpy.kernel.Kernel` arguments + (see :meth:`~sumpy.kernel.Kernel.get_args` and + :meth:`~sumpy.kernel.Kernel.get_source_args`) to expressions that determine + them. + """ - kernel_arguments = kernel_arguments.copy() - if kwargs: - for name, val in kwargs.items(): - if name in kernel_arguments: - raise ValueError("'%s' already set in kernel_arguments" - % name) + def __post_init__(self) -> None: + if self.qbx_forced_limit not in {-1, +1, -2, +2, "avg", None}: + raise ValueError( + f"Invalid value for 'qbx_forced_limit': {self.qbx_forced_limit}" + ) - if name not in kernel_arg_names: - raise TypeError("'%s' not recognized as kernel argument" - % name) + if not isinstance(self.source_kernels, tuple): + warn(f"'source_kernels' is not tuple ({type(self.source_kernels)}). " + "Passing a different type is deprecated and will stop working in " + "2025.", DeprecationWarning, stacklevel=2) - kernel_arguments[name] = val + object.__setattr__(self, "source_kernels", tuple(self.source_kernels)) - provided_arg_names = set(kernel_arguments.keys()) - missing_args = kernel_arg_names - provided_arg_names - if missing_args: - raise TypeError("kernel argument(s) '%s' not supplied" - % ", ".join(missing_args)) + if not isinstance(self.densities, tuple): + warn(f"'densities' is not tuple ({type(self.densities)}). " + "Passing a different type is deprecated and will stop working in " + "2025.", DeprecationWarning, stacklevel=2) - extraneous_args = provided_arg_names - kernel_arg_names - if missing_args: - raise TypeError("kernel arguments '%s' not recognized" - % ", ".join(extraneous_args)) + object.__setattr__(self, "densities", tuple(self.densities)) - self.target_kernel = target_kernel - self.source_kernels = source_kernels - self.densities = tuple(densities) - self.qbx_forced_limit = qbx_forced_limit - self.source = as_dofdesc(source) - self.target = as_dofdesc(target) - self.kernel_arguments = kernel_arguments + if not isinstance(self.source, DOFDescriptor): + warn("Passing a 'source' descriptor that is not a 'DOFDescriptor' to " + f"{type(self).__name__!r} is deprecated and will stop working " + "in 2025. Use 'as_dofdesc' to convert the descriptor.", + DeprecationWarning, stacklevel=2) - def copy(self, target_kernel=None, source_kernels=None, densities=None, - qbx_forced_limit=_NoArgSentinel, - source=_NoArgSentinel, target=_NoArgSentinel, - kernel_arguments=None): - if target_kernel is None: - target_kernel = self.target_kernel + object.__setattr__(self, "source", as_dofdesc(self.source)) - if source_kernels is None: - source_kernels = self.source_kernels + if not isinstance(self.target, DOFDescriptor): + warn("Passing a 'target' descriptor that is not a 'DOFDescriptor' to " + f"{type(self).__name__!r} is deprecated and will stop working " + "in 2025. Use 'as_dofdesc' to convert the descriptor.", + DeprecationWarning, stacklevel=2) - if densities is None: - densities = self.densities + object.__setattr__(self, "target", as_dofdesc(self.target)) - if kernel_arguments is None: - kernel_arguments = self.kernel_arguments + if not isinstance(self.kernel_arguments, dict): + warn(f"'kernel_arguments' is not a dict ({type(self.kernel_arguments)}). " + "Passing a different type is deprecated and will stop being " + "supported in 2025.", DeprecationWarning, stacklevel=2) - if qbx_forced_limit is _NoArgSentinel: - qbx_forced_limit = self.qbx_forced_limit + kernel_arguments = self.kernel_arguments if self.kernel_arguments else {} + object.__setattr__(self, "kernel_arguments", dict(kernel_arguments)) - source = self.source if source is _NoArgSentinel else as_dofdesc(source) - target = self.target if target is _NoArgSentinel else as_dofdesc(target) - return type(self)(target_kernel, source_kernels, densities, - qbx_forced_limit=qbx_forced_limit, - source=source, - target=target, - kernel_arguments=kernel_arguments) + from pytools import single_valued - def __getinitargs__(self): - return (self.target_kernel, self.source_kernels, self.densities, - self.qbx_forced_limit, self.source, self.target, - hashable_kernel_args(self.kernel_arguments)) + kernels = (*self.source_kernels, self.target_kernel) + single_valued(kernel.get_base_kernel() for kernel in kernels) - def __setstate__(self, state): - # Overwrite pymbolic.Expression.__setstate__ - assert len(self.init_arg_names) == len(state), type(self) - self.__init__(*state) + kernel_arg_names = set() + for kernel in kernels: + for karg in (kernel.get_args() + kernel.get_source_args()): + kernel_arg_names.add(karg.loopy_arg.name) + + provided_arg_names = set(self.kernel_arguments.keys()) # pylint: disable=no-member + missing_args = kernel_arg_names - provided_arg_names + if missing_args: + raise ValueError( + "Kernel argument(s) not supplied: '{}'".format(", ".join(missing_args)) + ) - mapper_method = intern("map_int_g") + # FIXME: this check is clearly wrong :( + extra_args = provided_arg_names - kernel_arg_names + if missing_args: + raise ValueError( + "Kernel argument(s) not recognized: '{}'".format(", ".join(extra_args)) + ) + + def copy(self, **kwargs) -> "IntG": + warn(f"'{type(self).__name__}.copy' is deprecated and will be removed in " + f"2025. {type(self)} is a dataclass now and can use " + "'dataclasses.replace'.", DeprecationWarning, stacklevel=2) + + from dataclasses import replace + return replace(self, **kwargs) + + def __eq__(self, other: Any) -> bool: + if self is other: + return True + if self.__class__ is not other.__class__: + return False + if hash(self) != hash(other): + return False + + return ( + self.__class__ == other.__class__ + and self.target_kernel == other.target_kernel + and self.source_kernels == other.source_kernels + and self.densities == other.densities + and self.source == other.source + and self.target == other.target + and ( + hashable_kernel_args(self.kernel_arguments) + == hashable_kernel_args(other.kernel_arguments)) + ) + + def __hash__(self) -> int: + return hash(( + self.target_kernel, + self.source_kernels, + self.densities, + self.qbx_forced_limit, + self.source, self.target, + hashable_kernel_args(self.kernel_arguments) + )) _DIR_VEC_NAME = "dsource_vec" @@ -1515,17 +1683,10 @@ def int_g_dsource(ambient_dim, dsource, kernel, density, r""" .. math:: - \int_\Gamma \operatorname{dsource} \dot \nabla_y - \dot g(x-y) \sigma(y) dS_y + \int_\Gamma \operatorname{dsource} \dot \nabla_y G(x-y) \sigma(y) dS_y - where :math:`\sigma` is *density*, and - *dsource*, a multivector. - Note that the first product in the integrand - is a geometric product. - - .. attribute:: dsource - - A :class:`pymbolic.geometric_algebra.MultiVector`. + where :math:`\sigma` is *density* and *dsource* is a multivector. + Note that the first product in the integrand is a geometric product. """ if kernel_arguments is None: @@ -1587,13 +1748,23 @@ def int_g_vec(kernel, density, qbx_forced_limit, source=None, target=None, txr = TargetTransformationRemover() target_kernel = sxr(kernel) - source_kernels = [txr(kernel)] + source_kernels = (txr(kernel),) + + if kernel_arguments is None: + kernel_arguments = {} + + if kwargs is not None: + kernel_arguments = {**kernel_arguments, **kwargs} def make_op(operand_i): - return IntG(target_kernel=target_kernel, densities=[operand_i], + return IntG( + target_kernel=target_kernel, source_kernels=source_kernels, - qbx_forced_limit=qbx_forced_limit, source=source, target=target, - kernel_arguments=kernel_arguments, **kwargs) + densities=(operand_i,), + qbx_forced_limit=qbx_forced_limit, + source=as_dofdesc(source), + target=as_dofdesc(target), + kernel_arguments=kernel_arguments) if isinstance(density, np.ndarray | MultiVector): return componentwise(make_op, density) @@ -1650,7 +1821,6 @@ def Sp(kernel, *args, **kwargs): kwargs["qbx_forced_limit"] = "avg" ambient_dim = kwargs.get("ambient_dim") - from sumpy.kernel import Kernel if ambient_dim is None and isinstance(kernel, Kernel): ambient_dim = kernel.dim if ambient_dim is None: @@ -1666,7 +1836,6 @@ def Sp(kernel, *args, **kwargs): def Spp(kernel, *args, **kwargs): ambient_dim = kwargs.get("ambient_dim") - from sumpy.kernel import Kernel if ambient_dim is None and isinstance(kernel, Kernel): ambient_dim = kernel.dim if ambient_dim is None: @@ -1683,7 +1852,6 @@ def Spp(kernel, *args, **kwargs): def D(kernel, *args, **kwargs): ambient_dim = kwargs.get("ambient_dim") - from sumpy.kernel import Kernel if ambient_dim is None and isinstance(kernel, Kernel): ambient_dim = kernel.dim if ambient_dim is None: @@ -1706,7 +1874,6 @@ def D(kernel, *args, **kwargs): def Dp(kernel, *args, **kwargs): ambient_dim = kwargs.get("ambient_dim") - from sumpy.kernel import Kernel if ambient_dim is None and isinstance(kernel, Kernel): ambient_dim = kernel.dim if ambient_dim is None: @@ -1750,9 +1917,9 @@ def tangential_onb(ambient_dim, dim=None, dofdesc=None): q = avec for j in range(k): q = q - np.dot(avec, orth_pd_mat[:, j])*orth_pd_mat[:, j] - q = cse(q, "q%d" % k) + q = cse(q, f"q{k}") - orth_pd_mat[:, k] = cse(q/sqrt(np.sum(q**2)), "orth_pd_vec%d_" % k) + orth_pd_mat[:, k] = cse(q/sqrt(np.sum(q**2)), f"orth_pd_vec{k}_") # }}} @@ -1778,7 +1945,8 @@ def tangential_to_xyz(tangential_vec, dofdesc=None): def project_to_tangential(xyz_vec, dofdesc=None): return tangential_to_xyz( - cse(xyz_to_tangential(xyz_vec, dofdesc), dofdesc)) + cse(xyz_to_tangential(xyz_vec, dofdesc)), + dofdesc) def n_dot(vec, dofdesc=None): diff --git a/test/test_cost_model.py b/test/test_cost_model.py index 9ca455d54..7e7ebeba2 100644 --- a/test/test_cost_model.py +++ b/test/test_cost_model.py @@ -336,7 +336,7 @@ def test_timing_data_gathering(ctx_factory): cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx, properties=cl.command_queue_properties.PROFILING_ENABLE) - actx = PyOpenCLArrayContext(queue, force_device_scalars=True) + actx = PyOpenCLArrayContext(queue) lpot_source = get_lpot_source(actx, 2) places = GeometryCollection(lpot_source) diff --git a/test/test_symbolic.py b/test/test_symbolic.py index 9abd50c13..cb24fd10b 100644 --- a/test/test_symbolic.py +++ b/test/test_symbolic.py @@ -349,7 +349,7 @@ def randrange_like(xi, offset): for xi in arys )) r = bind( - discr, func(sym.var("x"), dofdesc=sym.GRANULARITY_NODE) + discr, func(sym.var("x"), dofdesc=sym.as_dofdesc(sym.GRANULARITY_NODE)) )(actx, x=x) assert actx.to_numpy(flat_norm(r - expected)) < 1.0e-15 @@ -358,7 +358,7 @@ def randrange_like(xi, offset): for xi in arys )) r = bind( - discr, func(sym.var("x"), dofdesc=sym.GRANULARITY_ELEMENT) + discr, func(sym.var("x"), dofdesc=sym.as_dofdesc(sym.GRANULARITY_ELEMENT)) )(actx, x=x) assert actx.to_numpy(flat_norm(r - expected)) < 1.0e-15 @@ -391,8 +391,8 @@ def principal_directions(ambient_dim, dim=None, dofdesc=None): k1, k2 = principal_curvatures(ambient_dim, dim=dim, dofdesc=dofdesc) from pytools.obj_array import make_obj_array - d1 = sym.cse(make_obj_array([s12, -(s11 - k1)])) - d2 = sym.cse(make_obj_array([-(s22 - k2), s12])) + d1 = sym.cse(make_obj_array([s12, -(s11 - k1)]), scope=sym.cse_scope.EVALUATION) + d2 = sym.cse(make_obj_array([-(s22 - k2), s12]), scope=sym.cse_scope.EVALUATION) form1 = sym.first_fundamental_form(ambient_dim, dim=dim, dofdesc=dofdesc) return make_obj_array([ @@ -489,9 +489,10 @@ def test_mapper_int_g_term_collector(op_name, k=0): expr_only_intgs = IntGTermCollector()(expr) # FIXME: how to check this did something? - sigma = sym.cse(op.get_density_var("sigma") / op.get_sqrt_weight()) + sigma = sym.cse(op.get_density_var("sigma") / op.get_sqrt_weight(), + scope=sym.cse_scope.EVALUATION) if op_name == "dirichlet": - expected_expr = -1 * sym.D(op.kernel, sigma) + expected_expr = -1 * sym.D(op.kernel, sigma, qbx_forced_limit="avg") elif op_name == "neumann": int_g = sym.S(op.kernel, sigma, qbx_forced_limit="avg") expected_expr = sym.div([int_g] * ambient_dim)