diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py index 669048077..bb05c4608 100644 --- a/pytential/symbolic/mappers.py +++ b/pytential/symbolic/mappers.py @@ -280,7 +280,12 @@ def map_common_subexpression(self, expr): # {{{ FlattenMapper class FlattenMapper(FlattenMapperBase, IdentityMapper): - pass + def map_int_g(self, expr): + densities, kernel_arguments, changed = rec_int_g_arguments(self, expr) + if not changed: + return expr + + return replace(expr, densities=densities, kernel_arguments=kernel_arguments) def flatten(expr): diff --git a/pytential/symbolic/pde/system_utils.py b/pytential/symbolic/pde/system_utils.py index c8ca2d1a8..2e6d56dcf 100644 --- a/pytential/symbolic/pde/system_utils.py +++ b/pytential/symbolic/pde/system_utils.py @@ -21,27 +21,21 @@ """ import logging -import warnings from collections.abc import Mapping, Sequence -from dataclasses import dataclass -from typing import Any +from dataclasses import dataclass, replace +from typing import Any, cast import numpy as np -import pymbolic -import sumpy.symbolic as sym -from pytools import \ - generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam +import sumpy.symbolic as sp +from pytools import generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam from pytools import memoize_on_first_arg from sumpy.kernel import (AxisSourceDerivative, AxisTargetDerivative, DirectionalSourceDerivative, ExpressionKernel, Kernel, KernelWrapper, TargetPointMultiplier) -from pymbolic.typing import ExpressionT, ArithmeticExpressionT +from pymbolic.typing import ArithmeticExpression -import pytential -from pytential.symbolic.mappers import IdentityMapper -from pytential.symbolic.primitives import (DEFAULT_SOURCE, IntG, - NodeCoordinateComponent, - hashable_kernel_args) +from pytential import sym +from pytential.symbolic.mappers import IdentityMapper, flatten from pytential.utils import chop, solve_from_lu logger = logging.getLogger(__name__) @@ -72,8 +66,8 @@ class RewriteFailedError(RuntimeError): def rewrite_using_base_kernel( - exprs: Sequence[ExpressionT], - base_kernel: Kernel = _NO_ARG_SENTINEL) -> list[ExpressionT]: + exprs: Sequence[ArithmeticExpression], + base_kernel: Kernel = _NO_ARG_SENTINEL) -> list[ArithmeticExpression]: """ Rewrites a list of expressions with :class:`~pytential.symbolic.primitives.IntG` objects using *base_kernel*. @@ -101,7 +95,7 @@ def rewrite_using_base_kernel( raise NotImplementedError mapper = RewriteUsingBaseKernelMapper(base_kernel) - return [mapper(expr) for expr in exprs] + return [cast(ArithmeticExpression, mapper(expr)) for expr in exprs] class RewriteUsingBaseKernelMapper(IdentityMapper): @@ -131,8 +125,9 @@ def map_int_g(self, expr): self.base_kernel) for new_int_g in new_int_gs) -def _get_sympy_kernel_expression(expr: ArithmeticExpressionT, - kernel_arguments: Mapping[str, Any]) -> sym.Basic: +def _get_sympy_kernel_expression( + expr: ArithmeticExpression, + kernel_arguments: Mapping[str, Any]) -> sp.Basic: """Convert a :mod:`pymbolic` expression to :mod:`sympy` expression after substituting kernel arguments. @@ -149,14 +144,14 @@ def _get_sympy_kernel_expression(expr: ArithmeticExpressionT, def _monom_to_expr(monom: Sequence[int], - variables: Sequence[sym.Basic | ArithmeticExpressionT] - ) -> sym.Basic | ArithmeticExpressionT: + variables: Sequence[sp.Basic | ArithmeticExpression] + ) -> sp.Basic | ArithmeticExpression: """Convert a monomial to an expression using given variables. For example, ``[3, 2, 1]`` with variables ``[x, y, z]`` is converted to ``x^3 y^2 z``. """ - prod: ArithmeticExpressionT = 1 + prod: ArithmeticExpression = 1 for i, nrepeats in enumerate(monom): for _ in range(nrepeats): prod *= variables[i] @@ -164,7 +159,7 @@ def _monom_to_expr(monom: Sequence[int], return prod -def convert_target_transformation_to_source(int_g: IntG) -> list[IntG]: +def convert_target_transformation_to_source(int_g: sym.IntG) -> list[sym.IntG]: r"""Convert an ``IntG`` with :class:`~sumpy.kernel.AxisTargetDerivative` or :class:`~sumpy.kernel.TargetPointMultiplier` to a list of ``IntG``\ s without them and only source dependent transformations. @@ -182,8 +177,9 @@ def convert_target_transformation_to_source(int_g: IntG) -> list[IntG]: knl = int_g.target_kernel if not knl.is_translation_invariant: - warnings.warn(f"Translation variant kernel ({knl}) found.", - stacklevel=2) + from warnings import warn + + warn(f"Translation variant kernel ({knl}) found.", stacklevel=2) return [int_g] # we use a symbol for d = (x - y) @@ -204,7 +200,9 @@ def convert_target_transformation_to_source(int_g: IntG) -> list[IntG]: expr = expr.diff(ds[knl.axis]) found = True else: - warnings.warn( + from warnings import warn + + warn( f"Unknown target kernel ({knl}) found. " "Returning IntG expression unchanged.", stacklevel=2) return [int_g] @@ -213,9 +211,9 @@ def convert_target_transformation_to_source(int_g: IntG) -> list[IntG]: if not found: return [int_g] - int_g = int_g.copy(target_kernel=knl) + int_g = replace(int_g, target_kernel=knl) - sources_pymbolic = [NodeCoordinateComponent(i) for i in range(knl.dim)] + sources_pymbolic = sym.nodes(knl.dim).as_vector() expr = expr.expand() # Now the expr is an Add and looks like # u_{d[0], d[1]}(d, y)*d[0]*y[1] + u(d, y) * d[1] @@ -255,7 +253,7 @@ def convert_target_transformation_to_source(int_g: IntG) -> list[IntG]: for _ in range(nrepeats): knl = AxisSourceDerivative(axis, knl) new_source_kernels.append(knl) - new_int_g = int_g.copy(source_kernels=new_source_kernels) + new_int_g = replace(int_g, source_kernels=tuple(new_source_kernels)) (monom, coeff,) = remaining_factors.terms()[0] # Now from d[0]*y[1], we separate the two terms @@ -266,33 +264,39 @@ def convert_target_transformation_to_source(int_g: IntG) -> list[IntG]: * conv(coeff) # since d/d(d) = - d/d(y), we multiply by -1 to get source derivatives density_multiplier *= (-1)**int(sum(nrepeats for _, nrepeats in derivatives)) - new_int_gs = _multiply_int_g(new_int_g, sym.sympify(expr_multiplier), + new_int_gs = _multiply_int_g(new_int_g, sp.sympify(expr_multiplier), density_multiplier) result.extend(new_int_gs) return result -def _multiply_int_g(int_g: IntG, expr_multiplier: sym.Basic, - density_multiplier: ArithmeticExpressionT) -> list[IntG]: +def _multiply_int_g( + int_g: sym.IntG, + expr_multiplier: sp.Basic, + density_multiplier: ArithmeticExpression) -> list[sym.IntG]: """Multiply the expression in ``IntG`` with the *expr_multiplier* which is a symbolic (:mod:`sympy` or :mod:`symengine`) expression and multiply the densities with *density_multiplier* which is a :mod:`pymbolic` expression. """ + from pymbolic import substitute + result = [] base_kernel = int_g.target_kernel.get_base_kernel() - sym_d = sym.make_sym_vector("d", base_kernel.dim) + sym_d = sp.make_sym_vector("d", base_kernel.dim) base_kernel_expr = _get_sympy_kernel_expression(base_kernel.expression, int_g.kernel_arguments) - subst = {pymbolic.var(f"d{i}"): pymbolic.var("d")[i] for i in + subst = {sym.var(f"d{i}"): sym.var("d")[i] for i in range(base_kernel.dim)} - conv = sym.SympyToPymbolicMapper() + conv = sp.SympyToPymbolicMapper() if expr_multiplier == 1: # if there's no expr_multiplier, only multiply the densities - return [int_g.copy(densities=tuple(density*density_multiplier - for density in int_g.densities))] + return [replace( + int_g, + densities=tuple(density*density_multiplier for density in int_g.densities)) + ] for knl, density in zip(int_g.source_kernels, int_g.densities, strict=True): if expr_multiplier == 1: @@ -300,55 +304,58 @@ def _multiply_int_g(int_g: IntG, expr_multiplier: sym.Basic, else: new_expr = conv(knl.postprocess_at_source(base_kernel_expr, sym_d) * expr_multiplier) - new_expr = pymbolic.substitute(new_expr, subst) - new_knl = ExpressionKernel(knl.dim, new_expr, + new_expr = substitute(new_expr, subst) + new_knl = ExpressionKernel(knl.dim, flatten(new_expr), knl.get_base_kernel().global_scaling_const, knl.is_complex_valued) - result.append(int_g.copy(target_kernel=new_knl, + result.append(replace( + int_g, + target_kernel=new_knl, densities=(density*density_multiplier,), - source_kernels=(new_knl,))) + source_kernels=(new_knl,) + )) return result def rewrite_int_g_using_base_kernel( - int_g: IntG, base_kernel: ExpressionKernel) -> ArithmeticExpressionT: + int_g: sym.IntG, base_kernel: ExpressionKernel) -> ArithmeticExpression: r"""Rewrite an ``IntG`` to an expression with ``IntG``\ s having the base kernel *base_kernel*. """ - result: ArithmeticExpressionT = 0 + result: ArithmeticExpression = 0 for knl, density in zip(int_g.source_kernels, int_g.densities, strict=True): result += _rewrite_int_g_using_base_kernel( - int_g.copy(source_kernels=(knl,), densities=(density,)), + replace(int_g, source_kernels=(knl,), densities=(density,)), base_kernel) return result def _rewrite_int_g_using_base_kernel( - int_g: IntG, base_kernel: ExpressionKernel) -> ArithmeticExpressionT: + int_g: sym.IntG, base_kernel: ExpressionKernel) -> ArithmeticExpression: r"""Rewrites an ``IntG`` with only one source kernel to an expression with ``IntG``\ s having the base kernel *base_kernel*. """ target_kernel = int_g.target_kernel.replace_base_kernel(base_kernel) dim = target_kernel.dim - result: ArithmeticExpressionT = 0 + result: ArithmeticExpression = 0 density, = int_g.densities source_kernel, = int_g.source_kernels deriv_relation = get_deriv_relation_kernel(source_kernel.get_base_kernel(), base_kernel, hashable_kernel_arguments=( - hashable_kernel_args(int_g.kernel_arguments))) + sym.hashable_kernel_args(int_g.kernel_arguments))) const = deriv_relation.const # NOTE: we set a dofdesc here to force the evaluation of this integral # on the source instead of the target when using automatic tagging # see :meth:`pytential.symbolic.mappers.LocationTagger._default_dofdesc` if int_g.source.geometry is None: - dd = int_g.source.copy(geometry=DEFAULT_SOURCE) + dd = int_g.source.copy(geometry=sym.DEFAULT_SOURCE) else: dd = int_g.source - const *= pytential.sym.integral(dim, dim-1, density, dofdesc=dd) + const *= sym.integral(dim, dim-1, density, dofdesc=dd) if const != 0 and target_kernel != target_kernel.get_base_kernel(): # There might be some TargetPointMultipliers hanging around. @@ -377,8 +384,13 @@ def _rewrite_int_g_using_base_kernel( for _ in range(val): knl = AxisSourceDerivative(d, knl) c *= -1 - result += int_g.copy(source_kernels=(knl,), target_kernel=target_kernel, - densities=(density * c,), kernel_arguments=new_kernel_args) + result += replace( + int_g, + source_kernels=(knl,), + target_kernel=target_kernel, + densities=(density * c,), + kernel_arguments=new_kernel_args) + return result @@ -394,9 +406,9 @@ class DerivRelation: .. autoattribute:: linear_combination """ - const: ArithmeticExpressionT + const: ArithmeticExpression """A constant to add to the combination.""" - linear_combination: Sequence[tuple[tuple[int, ...], ArithmeticExpressionT]] + linear_combination: Sequence[tuple[tuple[int, ...], ArithmeticExpression]] """A list of pairs ``(mi, coeffs)``.""" @@ -432,7 +444,7 @@ def get_deriv_relation( res = [] for knl in kernels: res.append(get_deriv_relation_kernel(knl, base_kernel, - hashable_kernel_arguments=hashable_kernel_args(kernel_arguments), + hashable_kernel_arguments=sym.hashable_kernel_args(kernel_arguments), tol=tol, order=order)) return res @@ -460,14 +472,14 @@ def get_deriv_relation_kernel( order=order, hashable_kernel_arguments=hashable_kernel_arguments) dim = base_kernel.dim - sym_vec = sym.make_sym_vector("d", dim) - sympy_conv = sym.SympyToPymbolicMapper() + sym_vec = sp.make_sym_vector("d", dim) + sympy_conv = sp.SympyToPymbolicMapper() expr = _get_sympy_kernel_expression(kernel.expression, kernel_arguments) vec = [] for i in range(len(mis)): vec.append(evalf(expr.xreplace(dict(zip(sym_vec, rand[:, i], strict=True))))) - vec = sym.Matrix(vec) + vec = sp.Matrix(vec) result = [] const = 0 logger.debug("%s = ", kernel) @@ -493,8 +505,8 @@ def get_deriv_relation_kernel( @dataclass class LUFactorization: - L: sym.Matrix - U: sym.Matrix + L: sp.Matrix + U: sp.Matrix perm: Sequence[tuple[int, int]] @@ -539,8 +551,8 @@ def _get_base_kernel_matrix_lu_factorization( rand: np.ndarray = rng.integers(1, 10**15, size=(dim, len(mis))).astype(object) for i in range(rand.shape[0]): for j in range(rand.shape[1]): - rand[i, j] = sym.sympify(rand[i, j])/10**15 - sym_vec = sym.make_sym_vector("d", dim) + rand[i, j] = sp.sympify(rand[i, j])/10**15 + sym_vec = sp.make_sym_vector("d", dim) base_expr = _get_sympy_kernel_expression(base_kernel.expression, dict(hashable_kernel_arguments)) @@ -560,7 +572,7 @@ def _get_base_kernel_matrix_lu_factorization( row.append(1) mat.append(row) - sym_mat = sym.Matrix(mat) + sym_mat = sp.Matrix(mat) failed = False try: L, U, perm = sym_mat.LUdecomposition() @@ -569,7 +581,7 @@ def _get_base_kernel_matrix_lu_factorization( # and sympy returns U with last row zero failed = True - if not sym.USE_SYMENGINE and all(expr == 0 for expr in U[-1, :]): + if not sp.USE_SYMENGINE and all(expr == 0 for expr in U[-1, :]): failed = True if failed: diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index 17e244300..62decc28d 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -1625,7 +1625,7 @@ class IntG(Expression): derivatives attached. k-th elements represents the k-th source derivative operator above. """ - densities: tuple[Expression, ...] + densities: tuple[ArithmeticExpression, ...] """A tuple of density expressions. Length of this tuple must match the length of the *source_kernels* arguments. """ diff --git a/test/test_pde_system_utils.py b/test/test_pde_system_utils.py index 877c800cd..85df3bc9f 100644 --- a/test/test_pde_system_utils.py +++ b/test/test_pde_system_utils.py @@ -18,160 +18,205 @@ THE SOFTWARE. """ -from pytential.symbolic.pde.system_utils import ( - convert_target_transformation_to_source, rewrite_int_g_using_base_kernel) -from pytential.symbolic.primitives import IntG -from pytential import sym -import pytential +from functools import partial + import numpy as np +import pymbolic.primitives as prim from sumpy.kernel import ( LaplaceKernel, HelmholtzKernel, ExpressionKernel, BiharmonicKernel, StokesletKernel, AxisTargetDerivative, TargetPointMultiplier, AxisSourceDerivative) from sumpy.symbolic import USE_SYMENGINE -from pymbolic.primitives import make_sym_vector -import pymbolic.primitives as prim +from pytential import sym +from pytential.symbolic.mappers import flatten +from pytential.symbolic.pde.system_utils import ( + convert_target_transformation_to_source, rewrite_int_g_using_base_kernel) def test_convert_target_deriv(): knl = LaplaceKernel(2) - int_g = IntG(AxisTargetDerivative(0, knl), [AxisSourceDerivative(1, knl), knl], - [1, 2], qbx_forced_limit=1) - expected_int_g = IntG(knl, - [AxisSourceDerivative(0, AxisSourceDerivative(1, knl)), - AxisSourceDerivative(0, knl)], [-1, -2], qbx_forced_limit=1) + dsource = partial(AxisSourceDerivative, inner_kernel=knl) - assert sum(convert_target_transformation_to_source(int_g)) == expected_int_g + int_g = sym.IntG( + AxisTargetDerivative(0, knl), + (dsource(1), knl), + (1, 2), qbx_forced_limit=1) + expected_int_g = sym.IntG( + knl, + (AxisSourceDerivative(0, dsource(1)), dsource(0)), + (-1, -2), qbx_forced_limit=1) + expr = sum(convert_target_transformation_to_source(int_g)) + assert flatten(expr) == expected_int_g -def test_convert_target_point_multiplier(): - xs = sym.nodes(3).as_vector() - knl = LaplaceKernel(3) - int_g = IntG(TargetPointMultiplier(0, knl), [AxisSourceDerivative(1, knl), knl], - [1, 2], qbx_forced_limit=1) +def test_convert_target_point_multiplier(): + ambient_dim = 3 + knl = LaplaceKernel(ambient_dim) + dsource = partial(AxisSourceDerivative, inner_kernel=knl) - d = make_sym_vector("d", 3) + int_g = sym.IntG( + TargetPointMultiplier(0, knl), + (dsource(1), knl), + (1, 2), qbx_forced_limit=1) + d = sym.make_sym_vector("d", ambient_dim) r2 = d[2]**2 + d[1]**2 + d[0]**2 - eknl0 = ExpressionKernel(3, d[1]*d[0]*r2**prim.Quotient(-3, 2), - knl.global_scaling_const, False) - eknl2 = ExpressionKernel(3, d[0]*r2**prim.Quotient(-1, 2), - knl.global_scaling_const, False) + + eknl0 = ExpressionKernel(ambient_dim, + d[1]*d[0]*r2**prim.Quotient(-3, 2), + knl.global_scaling_const, is_complex_valued=False) + eknl2 = ExpressionKernel(ambient_dim, + d[0]*r2**prim.Quotient(-1, 2), + knl.global_scaling_const, is_complex_valued=False) r2 = d[0]**2 + d[1]**2 + d[2]**2 - eknl1 = ExpressionKernel(3, d[0]*d[1]*r2**prim.Quotient(-3, 2), - knl.global_scaling_const, False) - eknl3 = ExpressionKernel(3, d[0]*r2**prim.Quotient(-1, 2), - knl.global_scaling_const, False) - possible_int_g1 = IntG(eknl0, [eknl0], [1], qbx_forced_limit=1) + \ - IntG(eknl2, [eknl2], [2], qbx_forced_limit=1) + \ - IntG(knl, [AxisSourceDerivative(1, knl), knl], - [xs[0], 2*xs[0]], qbx_forced_limit=1) + eknl1 = ExpressionKernel(ambient_dim, + d[0]*d[1]*r2**prim.Quotient(-3, 2), + knl.global_scaling_const, is_complex_valued=False) + eknl3 = ExpressionKernel(ambient_dim, + d[0]*r2**prim.Quotient(-1, 2), + knl.global_scaling_const, is_complex_valued=False) - possible_int_g2 = IntG(eknl1, [eknl1], [1], qbx_forced_limit=1) + \ - IntG(eknl3, [eknl3], [2], qbx_forced_limit=1) + \ - IntG(knl, [AxisSourceDerivative(1, knl), knl], - [xs[0], 2*xs[0]], qbx_forced_limit=1) + xs = sym.nodes(3).as_vector() - assert sum(convert_target_transformation_to_source(int_g)) in \ - [possible_int_g1, possible_int_g2] + possible_int_g1 = flatten( + sym.IntG(eknl0, (eknl0,), (1,), qbx_forced_limit=1) + + sym.IntG(eknl2, (eknl2,), (2,), qbx_forced_limit=1) + + sym.IntG(knl, (dsource(1), knl), (xs[0], 2*xs[0]), qbx_forced_limit=1)) + possible_int_g2 = flatten( + sym.IntG(eknl1, (eknl1,), (1,), qbx_forced_limit=1) + + sym.IntG(eknl3, (eknl3,), (2,), qbx_forced_limit=1) + + sym.IntG(knl, (dsource(1), knl), (xs[0], 2*xs[0]), qbx_forced_limit=1)) + + expr = flatten(sum(convert_target_transformation_to_source(int_g))) + assert expr in (possible_int_g1, possible_int_g2) def test_product_rule(): - xs = sym.nodes(3).as_vector() + ambient_dim = 3 + knl = LaplaceKernel(ambient_dim) + dsource = partial(AxisSourceDerivative, inner_kernel=knl) - knl = LaplaceKernel(3) - int_g = IntG(AxisTargetDerivative(0, TargetPointMultiplier(0, knl)), [knl], [1], - qbx_forced_limit=1) + int_g = sym.IntG( + AxisTargetDerivative(0, TargetPointMultiplier(0, knl)), + (knl,), (1,), qbx_forced_limit=1) - d = make_sym_vector("d", 3) + d = sym.make_sym_vector("d", 3) r2 = d[2]**2 + d[1]**2 + d[0]**2 - eknl0 = ExpressionKernel(3, d[0]**2*r2**prim.Quotient(-3, 2), - knl.global_scaling_const, False) - r2 = d[0]**2 + d[1]**2 + d[2]**2 - eknl1 = ExpressionKernel(3, d[0]**2*r2**prim.Quotient(-3, 2), - knl.global_scaling_const, False) - possible_int_g1 = IntG(eknl0, [eknl0], [-1], qbx_forced_limit=1) + \ - IntG(knl, [AxisSourceDerivative(0, knl)], [xs[0]*(-1)], qbx_forced_limit=1) - possible_int_g2 = IntG(eknl1, [eknl1], [-1], qbx_forced_limit=1) + \ - IntG(knl, [AxisSourceDerivative(0, knl)], [xs[0]*(-1)], qbx_forced_limit=1) + eknl0 = ExpressionKernel(ambient_dim, + d[0]**2*r2**prim.Quotient(-3, 2), + knl.global_scaling_const, is_complex_valued=False) - assert sum(convert_target_transformation_to_source(int_g)) in \ - [possible_int_g1, possible_int_g2] + r2 = d[0]**2 + d[1]**2 + d[2]**2 + eknl1 = ExpressionKernel(ambient_dim, + d[0]**2*r2**prim.Quotient(-3, 2), + knl.global_scaling_const, is_complex_valued=False) -def test_convert_helmholtz(): xs = sym.nodes(3).as_vector() - knl = HelmholtzKernel(3) - int_g = IntG(TargetPointMultiplier(0, knl), [knl], [1], - qbx_forced_limit=1, k=1) + possible_int_g1 = flatten( + sym.IntG(eknl0, (eknl0,), (-1,), qbx_forced_limit=1) + + sym.IntG(knl, (dsource(0),), (xs[0]*(-1),), qbx_forced_limit=1) + ) + possible_int_g2 = flatten( + sym.IntG(eknl1, (eknl1,), (-1,), qbx_forced_limit=1) + + sym.IntG(knl, (dsource(0),), (xs[0]*(-1),), qbx_forced_limit=1)) + + expr = flatten(sum(convert_target_transformation_to_source(int_g))) + assert expr in [possible_int_g1, possible_int_g2] + - d = make_sym_vector("d", 3) +def test_convert_helmholtz(): + ambient_dim = 3 + knl = HelmholtzKernel(ambient_dim) + int_g = sym.IntG( + TargetPointMultiplier(0, knl), + (knl,), (1,), qbx_forced_limit=1, kernel_arguments={"k": 1}) + + d = sym.make_sym_vector("d", 3) exp = prim.Variable("exp") if USE_SYMENGINE: r2 = d[2]**2 + d[1]**2 + d[0]**2 - eknl = ExpressionKernel(3, exp(1j*r2**prim.Quotient(1, 2))*d[0] - * r2**prim.Quotient(-1, 2), - knl.global_scaling_const, knl.is_complex_valued) + eknl = ExpressionKernel(ambient_dim, + exp(1j*r2**prim.Quotient(1, 2))*d[0] + * r2**prim.Quotient(-1, 2), + knl.global_scaling_const, + is_complex_valued=knl.is_complex_valued) else: r2 = d[0]**2 + d[1]**2 + d[2]**2 - eknl = ExpressionKernel(3, d[0]*r2**prim.Quotient(-1, 2) - * exp(1j*r2**prim.Quotient(1, 2)), - knl.global_scaling_const, knl.is_complex_valued) + eknl = ExpressionKernel(ambient_dim, + d[0]*r2**prim.Quotient(-1, 2) + * exp(1j*r2**prim.Quotient(1, 2)), + knl.global_scaling_const, + is_complex_valued=knl.is_complex_valued) - expected_int_g = IntG(eknl, [eknl], [1], qbx_forced_limit=1, - kernel_arguments={"k": 1}) + \ - IntG(knl, [knl], [xs[0]], qbx_forced_limit=1, k=1) + xs = sym.nodes(3).as_vector() + expected_int_g = flatten( + sym.IntG(eknl, (eknl,), (1,), qbx_forced_limit=1, kernel_arguments={"k": 1}) + + sym.IntG(knl, (knl,), (xs[0],), qbx_forced_limit=1, kernel_arguments={"k": 1}) + ) - assert expected_int_g == sum(convert_target_transformation_to_source(int_g)) + expr = flatten(sum(convert_target_transformation_to_source(int_g))) + assert expr == expected_int_g def test_convert_int_g_base(): - knl = LaplaceKernel(3) - int_g = IntG(knl, [knl], [1], qbx_forced_limit=1) + ambient_dim = 3 + knl = LaplaceKernel(ambient_dim) + int_g = sym.IntG(knl, (knl,), (1,), qbx_forced_limit=1) - base_knl = BiharmonicKernel(3) - expected_int_g = sum( - IntG(base_knl, [AxisSourceDerivative(d, AxisSourceDerivative(d, base_knl))], - [-1], qbx_forced_limit=1) for d in range(3)) + base_knl = BiharmonicKernel(ambient_dim) + expected_int_g = sum(sym.IntG( + base_knl, + (AxisSourceDerivative(d, AxisSourceDerivative(d, base_knl)),), + (-1,), qbx_forced_limit=1) + for d in range(ambient_dim)) - assert expected_int_g == rewrite_int_g_using_base_kernel(int_g, - base_kernel=base_knl) + expr = flatten(rewrite_int_g_using_base_kernel(int_g, base_kernel=base_knl)) + assert expr == flatten(expected_int_g) def test_convert_int_g_base_with_const(): - knl = StokesletKernel(2, 0, 0) - int_g = IntG(knl, [knl], [1], qbx_forced_limit=1, mu=2) + ambient_dim = 2 + knl = StokesletKernel(ambient_dim, 0, 0) + base_knl = BiharmonicKernel(ambient_dim) - base_knl = BiharmonicKernel(2) - dim = 2 - dd = pytential.sym.DOFDescriptor(geometry=pytential.sym.DEFAULT_SOURCE) + dd = sym.DOFDescriptor(sym.DEFAULT_SOURCE) + int_g = sym.IntG(knl, (knl,), (1,), qbx_forced_limit=1, kernel_arguments={"mu": 2}) - expected_int_g = (-0.1875) / np.pi * \ - pytential.sym.integral(dim, dim-1, 1, dofdesc=dd) + \ - IntG(base_knl, - [AxisSourceDerivative(1, AxisSourceDerivative(1, base_knl))], [0.5], - qbx_forced_limit=1) - assert rewrite_int_g_using_base_kernel(int_g, - base_kernel=base_knl) == expected_int_g + expected_int_g = flatten( + (-0.1875) / np.pi * sym.integral(ambient_dim, ambient_dim - 1, 1, dofdesc=dd) + + sym.IntG(base_knl, + (AxisSourceDerivative(1, AxisSourceDerivative(1, base_knl)),), + (0.5,), qbx_forced_limit=1)) + expr = flatten(rewrite_int_g_using_base_kernel(int_g, base_kernel=base_knl)) + assert expr == expected_int_g -def test_convert_int_g_base_with_const_and_deriv(): - knl = StokesletKernel(2, 0, 0) - int_g = IntG(knl, [AxisSourceDerivative(0, knl)], [1], qbx_forced_limit=1, mu=2) - - base_knl = BiharmonicKernel(2) - expected_int_g = IntG(base_knl, - [AxisSourceDerivative(1, AxisSourceDerivative(1, - AxisSourceDerivative(0, base_knl)))], [0.5], - qbx_forced_limit=1) - assert rewrite_int_g_using_base_kernel(int_g, - base_kernel=base_knl) == expected_int_g +def test_convert_int_g_base_with_const_and_deriv(): + ambient_dim = 2 + knl = StokesletKernel(ambient_dim, 0, 0) + base_knl = BiharmonicKernel(ambient_dim) + + int_g = sym.IntG( + knl, + (AxisSourceDerivative(0, knl),), (1,), qbx_forced_limit=1, + kernel_arguments={"mu": 2}) + + expected_int_g = sym.IntG( + base_knl, + (AxisSourceDerivative( + 1, AxisSourceDerivative( + 1, AxisSourceDerivative(0, base_knl))),), + (0.5,), qbx_forced_limit=1) + + expr = flatten(rewrite_int_g_using_base_kernel(int_g, base_kernel=base_knl)) + assert expr == expected_int_g diff --git a/test/test_stokes.py b/test/test_stokes.py index 0e63e6f3a..a2e63cce3 100644 --- a/test/test_stokes.py +++ b/test/test_stokes.py @@ -612,13 +612,13 @@ def run_stokes_pde(actx_factory, case, identity, resolution, visualize=False): identity.apply_pde_using_pytential())(actx)) m = np.max([np.linalg.norm(p, ord=np.inf) for p in potential_host]) - error = [0.0] * places.ambient_dim - error[:places.ambient_dim] = [np.linalg.norm(x, ord=np.inf)/m for x in result] + error = ( + [np.linalg.norm(x, ord=np.inf)/m for x in result] + + [np.linalg.norm(x, ord=np.inf)/m for x in result_pytential]) - error += [np.linalg.norm(x, ord=np.inf)/m for x in result_pytential] logger.info( - "resolution %4d h_min %.5e h_max %.5e error %.5e %.5e %.5e", - resolution, h_min, h_max, *error) + "resolution %4d h_min %.5e h_max %.5e error %s", + resolution, h_min, h_max, *error[:places.ambient_dim]) return h_max, error @@ -628,6 +628,10 @@ def __init__(self, ambient_dim, wrapper): self.ambient_dim = ambient_dim self.wrapper = wrapper + @property + def dim(self) -> int: + return self.ambient_dim + 1 + def apply_operator(self): dim = self.ambient_dim args = { @@ -675,6 +679,10 @@ def __init__(self, ambient_dim, wrapper): self.ambient_dim = ambient_dim self.wrapper = wrapper + @property + def dim(self) -> int: + return self.ambient_dim + def apply_operator(self): dim = self.ambient_dim args = { @@ -795,7 +803,7 @@ def test_elasticity_pde(actx_factory, dim, method, nu, is_double_layer, case.ambient_dim, mu=1, nu=nu, method=Method[method])) from pytools.convergence import EOCRecorder - eocs = [EOCRecorder() for _ in range(2*case.ambient_dim)] + eocs = [EOCRecorder() for _ in range(2*identity.dim)] for resolution in case.resolutions: h_max, errors = run_stokes_pde(