diff --git a/sumpy/kernel_rewrite.py b/sumpy/kernel_rewrite.py new file mode 100644 index 00000000..d7e37e51 --- /dev/null +++ b/sumpy/kernel_rewrite.py @@ -0,0 +1,622 @@ +from __future__ import annotations + + +__copyright__ = """ +Copyright (C) 2012 Andreas Kloeckner +Copyright (C) 2020 Isuru Fernando +Copyright (C) 2026 Alexandru Fikl +""" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import logging +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any, NamedTuple + +import numpy as np + +from pytools import ( + generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, +) + +import sumpy.symbolic as sym + + +if TYPE_CHECKING: + from collections.abc import Sequence + + import optype.numpy as onp + + from pymbolic.typing import ArithmeticExpression + + from sumpy.expansions.diff_op import MultiIndex + from sumpy.kernel import ScalarKernel + +logger = logging.getLogger(__name__) + + +# {{{ rewrite_using_base_kernel + + +class RewriteFailedError(Exception): + pass + + +@dataclass(frozen=True) +class LinearOperatorRepresentation: + r"""Expresses a target kernel as a linear operator acting on a base kernel. + + .. math:: + + G(\boldsymbol{r}) = C + \sum_{|\boldsymbol{\alpha}| < p} + c_{\boldsymbol{\alpha}} + \frac{\partial^{|\boldsymbol{\alpha}|} G_0} + {\partial \boldsymbol{r}^{\boldsymbol{\alpha}}} + + .. autoattribute:: target_kernel + .. autoattribute:: base_kernel + .. autoattribute:: mis + .. autoattribute:: coeffs + """ + + target_kernel: ScalarKernel + """The target kernel expressed in terms of :attr:`base_kernel`.""" + base_kernel: ScalarKernel + """A base kernel used to express the target kernel.""" + + mis: Sequence[MultiIndex] + """Multi-indices for each non-zero term in the linear combination.""" + coeffs: Sequence[ArithmeticExpression] + """Constant coefficients in the linear combination. Note that the coefficients + have length ``len(mis) + 1``, where the first element is always the constant term. + """ + + def pretty(self) -> str: + from sumpy.kernel import AxisTargetDerivative + + terms = [] + if self.coeffs[0] != 0: + terms.append(str(self.coeffs[0])) + + for mi, c in zip(self.mis, self.coeffs[1:], strict=True): + expr = self.base_kernel + for d, n in enumerate(mi): + for _ in range(n): + expr = AxisTargetDerivative(d, expr) + + terms.append(f"{c} * {expr}") + + return f"{self.target_kernel} = " + " + ".join(terms) + + +def rewrite_using_base_kernel( + target_kernel: ScalarKernel, + base_kernel: ScalarKernel, + *, + order: int | None = None, + atol: float = 1.0e-10, + ) -> LinearOperatorRepresentation: + pde = base_kernel.get_pde_as_diff_op() + if order is None: + order = pde.order + + # TODO: pick the best algorithm here + return rewrite_using_base_kernel_lu( + target_kernel, base_kernel, + min_order=order, + ) + + +# }}} + + +# {{{ rewrite_using_base_kernel_fourier + + +def rewrite_using_base_kernel_fourier( + target_kernel: ScalarKernel, + base_kernel: ScalarKernel, +) -> LinearOperatorRepresentation: + r"""Find a relation between *target_kernel* and *base_kernel* using the + Fourier symbol of their respective PDEs. + + The algorithm works by computing the Fourier symbols :math:`P_{\text{base}}` + and :math:`P_{\text{target}}` of the scalar PDEs satisfied by the two + kernels. The target kernel can be expressed as a differential operator + applied to the base kernel if and only if the ratio + + .. math:: + + D(i \boldsymbol{k}) = \frac{P_{\text{base}}(\boldsymbol{k})} + {P_{\text{target}}(\boldsymbol{k})} + + is a polynomial in :math:`\boldsymbol{k}`. When it is, each monomial + :math:`\prod_j (i k_j)^{\alpha_j}` corresponds to the derivative + :math:`\partial^{|\alpha|} / \partial \boldsymbol{r}^{\alpha}`. + """ + import sympy as sp + + # FIXME: want to add a general check for this? + try: + _ = base_kernel.get_pde_system_kernel() + except TypeError: + pass + else: + raise ValueError( + f"'base_kernel' cannot be part of a system: {type(base_kernel)}" + ) + + dim = base_kernel.dim + if target_kernel.dim != dim: + raise ValueError( + f"kernel dimension mismatch: {target_kernel.dim} (target_kernel) and " + f"{dim} (base_kernel)" + ) + + from sumpy.expansion.diff_op import to_fourier_matrix + + ks = sp.Matrix([sp.Symbol(f"_k{j}") for j in range(dim)]) + + pde_base = base_kernel.get_pde_as_diff_op() + fourier_base = to_fourier_matrix(pde_base, ks).inv() + + try: + target_system_kernel, idx = target_kernel.get_pde_system_kernel() + except TypeError: + target_system_kernel, idx = None, None + + if target_system_kernel is None: + pde_target = target_kernel.get_pde_as_diff_op() + fourier_target = to_fourier_matrix(pde_target, ks).inv() + else: + assert idx is not None + + pde_target_system = target_system_kernel.get_pde_as_diff_op() + fourier_target_system = to_fourier_matrix(pde_target_system, ks) + fourier_target_system_inv = sp.simplify(fourier_target_system.inv()) + + fourier_target = sp.Matrix([fourier_target_system_inv[idx]]) + + p_base: sp.Expr = sp.Integer(1) + for entry in fourier_base: + p_base = p_base * entry + + p_target: sp.Expr = sp.Integer(1) + for entry in fourier_target: + p_target = p_target * entry + + quotient = sp.simplify(p_target / p_base) + if not quotient.is_polynomial(*ks): + raise RewriteFailedError( + f"cannot rewrite {target_kernel} in terms of {base_kernel}" + ) + + mis: list[MultiIndex] = [] + coeffs: list[sp.Expr] = [] + + for exponent, coeff in sp.Poly(quotient, *ks).as_dict().items(): + coeff = sp.simplify(coeff) + if coeff == sp.Integer(0): + continue + + mis.append(exponent) + coeffs.append(coeff) + + if not mis: + raise RewriteFailedError( + f"cannot rewrite {target_kernel} in terms of {base_kernel}" + ) + + # compute constant term, if any + dvec = sym.make_sym_vector("d", dim) + base_scaled = ( + base_kernel.get_global_scaling_const() + * base_kernel.get_expression(dvec) + ) + + # build the differential operator part in real space + const: sp.Expr = ( + target_kernel.get_global_scaling_const() + * target_kernel.get_expression(dvec) + ) + + for mi, c in zip(mis, coeffs, strict=True): + term = c * base_scaled + for i, n in enumerate(mi): + term = term.diff(dvec[i], n) + const = const - sp.simplify(term) + const = sp.simplify(const) + + return LinearOperatorRepresentation( + target_kernel, + base_kernel, + mis, + [sym.to_pymbolic(expr) for expr in [const, *coeffs]], + ) + + +# }}} + + +# {{{ rewrite_using_base_kernel_lu + + +class FactorizationFailedError(Exception): + pass + + +class _LUDecomposition(NamedTuple): + L: sym.Matrix + U: sym.Matrix + permutation: Sequence[tuple[int, int]] + + mis: Sequence[MultiIndex] + """The multi-indices for which the derivatives were computed. These correspond + to rows in the matrix and should be used to recover the expansion of the + target kernel in terms of the base kernel. + """ + points: onp.Array2D[Any] + """An array of shape ``(dim, npoints)`` of points where the base kernel was + evaluated to compute the current LU factorization. + """ + + +def evalf(expr: sym.Expr, prec: int = 100) -> sym.Expr: + """Evaluate an expression numerically using ``prec`` number of bits.""" + from sumpy.symbolic import USE_SYMENGINE + + if USE_SYMENGINE: + return expr.n(prec=prec) + else: + import sympy + dps = int(sympy.log(2**prec, 10)) + return expr.n(n=dps) + + +def round_expr( + expr: sym.Basic, atol: float = 1.0e-8, rtol: float = 1.0e-5 + ) -> sym.Basic: + """Round all numeric values in *expr* to the nearest integer or fraction. + + This function clips all numbers close to zero (effectively removing them + from the expression due to SymPy's automatic simplifications) and rounds all + numbers to the nearest integer using the given *atol* and *rtol*. + """ + from fractions import Fraction + + nums = expr.atoms(sym.Float) + replace_dict: dict[Any, float] = {} + + for num in nums: + value = float(num) + nearest_int = round(value) + + if abs(value - nearest_int) < atol + rtol * abs(value): + replace_dict[num] = sym.Integer(nearest_int) + else: + frac = Fraction(value).limit_denominator(1000) + if abs(float(frac) - value) < atol + rtol * abs(value): + replace_dict[num] = sym.Rational(frac.numerator, frac.denominator) + else: + replace_dict[num] = value + + return expr.xreplace(replace_dict) + + +def simplify(expr: sym.Basic) -> sym.Basic: + if sym.USE_SYMENGINE: + return expr.simplify() + else: + import sympy as sp + + return sp.simplify(expr) + + +def _generate_points_shells( + dim: int, + npoints: int, *, + rmin: float = 0.25, + rmax: float = 2.0, + rng: np.random.Generator | None = None + ) -> onp.Array2D[Any]: + if dim < 1: + raise ValueError(f"'dim' must be >= 1: {dim!r}") + + if npoints < 1: + raise ValueError(f"'npoints' must be >= 1: {npoints!r}") + + if rmin >= rmax: + raise ValueError(f"'rmin' must be smaller than 'rmax': {rmin} >= {rmax}") + + if rng is None: + rng = np.random.default_rng() + + # make log spaced shell radii + radii = np.exp(rng.uniform(np.log(rmin), np.log(rmax), npoints)) + + # generate points on the unit sphere + p = rng.standard_normal((dim, npoints)) + p /= np.linalg.norm(p, axis=0) + + return radii * p + + +def _make_expr_derivatives( + base_expr: sym.Expr, dvec: sym.Matrix, mis: Sequence[MultiIndex] + ) -> dict[MultiIndex, sym.Basic]: + mi_to_derivative: dict[MultiIndex, sym.Basic] = {} + for mi in mis: + expr = base_expr + for i, nderivs in enumerate(mi): + if nderivs == 0: + continue + expr = expr.diff(dvec[i], nderivs) + + mi_to_derivative[mi] = expr + + return mi_to_derivative + + +def _make_derivative_matrix( + points: onp.Array2D[Any], + dvec: sym.Matrix, + mis: Sequence[MultiIndex], + mi_to_derivative: dict[MultiIndex, sym.Basic] + ) -> sym.Matrix: + # evaluate derivatives at points and construct matrix + entries: list[list[sym.Basic]] = [] + for i in range(points.shape[1]): + row: list[sym.Basic] = [sym.Integer(1)] + + for mi in mis: + expr = mi_to_derivative[mi].xreplace( + dict(zip(dvec, points[:, i], strict=True)) + ) + row.append(evalf(expr)) + + entries.append(row) + + return sym.Matrix(entries) + + +def _check_linear_combination( + target_expr: sym.Expr, + base_expr: sym.Expr, + dvec: sym.Matrix, + mis: Sequence[MultiIndex], + coeffs: Sequence[sym.Basic], + *, + rng: np.random.Generator, + rtol: float = 1.0e-8, + ) -> bool: + dim = len(dvec) + mi_to_derivative = _make_expr_derivatives(base_expr, dvec, mis) + points = _generate_points_shells(dim, len(mis) + 1, rng=rng) + + max_lhs = 0.0 + max_err = 0.0 + for i in range(points.shape[1]): + subst = dict(zip(dvec, points[:, i], strict=True)) + lhs = evalf(target_expr.xreplace(subst)) + + rhs = evalf(coeffs[0]) + for c, mi in zip(coeffs[1:], mis, strict=True): + rhs += evalf(c * mi_to_derivative[mi].xreplace(subst)) + + if lhs.is_number and rhs.is_number: + max_lhs = max(max_lhs, abs(lhs)) + max_err = max(max_err, abs(lhs - rhs)) + else: + max_err = round_expr(simplify(lhs - rhs), rtol=rtol) + if not max_err.is_zero: + return False + + return max_err <= rtol * max(max_lhs, 1.0) + + +def _get_base_kernel_matrix_lu_factorization( + base_kernel: ScalarKernel, + order: int, + *, + rng: np.random.Generator, + retries: int, + ) -> _LUDecomposition: + pde = base_kernel.get_pde_as_diff_op() + if order > pde.order: + raise NotImplementedError( + "Rewriting when the base kernel's derivatives are linearly dependent " + "is not implemented") + + dim = base_kernel.dim + + mis = list(gnitstam(order, dim)) + if order == pde.order: + pde_mis = [ident.mi for eq in pde.eqs for ident in eq] + pde_mis = [mi for mi in pde_mis if sum(mi) == order] + mis.remove(pde_mis[-1]) + + logger.debug("Removing %s to avoid linear dependent mis", pde_mis[-1]) + + # get sympy expression for the base kernel + dvec = sym.make_sym_vector("d", dim) + base_expr = base_kernel.get_expression(dvec) + + # evaluate all the needed derivatives + mi_to_derivative = _make_expr_derivatives(base_expr, dvec, mis) + + # try to LU factorize on random points + for _ in range(retries): + points = _generate_points_shells(dim, len(mis) + 1, rng=rng) + mat = _make_derivative_matrix(points, dvec, mis, mi_to_derivative) + + # NOTE: this can only happen if the points are somehow degenerate, e.g. + # points[:, 0] == points[:, 1], so it shouldn't happen under normal + # operating conditions? + try: + L, U, perm = mat.LUdecomposition() # ruff: ignore[non-lowercase-variable-in-function] + except RuntimeError: + continue + else: + # NOTE: and sympy seems to set the last row of U to 0 + if not sym.USE_SYMENGINE and all(expr == 0 for expr in U[-1, :]): + continue + + return _LUDecomposition(L, U, perm, mis, points) + + raise FactorizationFailedError( + f"failed to compute LU factorization to order {order} for {base_kernel} " + f"after {retries} retries" + ) + + +def rewrite_using_base_kernel_lu( + target_kernel: ScalarKernel, + base_kernel: ScalarKernel, + *, + min_order: int | None = None, + retries: int = 5, + rng: np.random.Generator | None = None, + ) -> LinearOperatorRepresentation: + """Find a relation between the *target_kernel* and the *base_kernel* using + a numerical LU-based algorithm. + + The algorithm samples the *base_kernel* and its derivatives at random + points to get a matrix ``A``. It also samples the target kernel at the same + points to get a vector ``b`` and solving for the system ``A c = b`` using + an LU factorization of ``A``. The solution ``c`` is the vector of coefficients + in the linear combination :class:`LinearOperatorRepresentation`. + + :arg min_order: starting minimum derivative order to use when attempting the + decomposition. By default, this will be the order of the PDE solved by + *base_kernel*. + :arg retries: maximum number of retries for each order. If the LU decomposition + fails due to a poor choice of random points, it is retried several times. + """ + try: + _ = base_kernel.get_pde_system_kernel() + except TypeError: + pass + else: + raise ValueError( + f"'base_kernel' cannot be part of a system: {type(base_kernel)}" + ) + + pde = base_kernel.get_pde_as_diff_op() + if min_order is None: + min_order = pde.order + + if min_order > pde.order: + raise NotImplementedError( + "Rewriting when the base kernel's derivatives are linearly dependent " + "is not implemented") + + if rng is None: + rng = np.random.default_rng() + + coeffs: list[sym.Basic] = [] + mis: list[MultiIndex] = [] + + dim = base_kernel.dim + dvec = sym.make_sym_vector("d", dim) + target_expr = target_kernel.get_expression(dvec) + base_expr = base_kernel.get_expression(dvec) + + target_scaling = target_kernel.get_global_scaling_const() + base_scaling = base_kernel.get_global_scaling_const() + + order = min_order + while order <= pde.order: + try: + lu = _get_base_kernel_matrix_lu_factorization( + base_kernel, order, rng=rng, retries=retries + ) + except FactorizationFailedError as exc: + if order == pde.order: + raise RewriteFailedError( + f"failed to compute LU factorization for orders in " + f"[{min_order}, {pde.order}] for base kernel {base_kernel} " + f"after {retries} retries" + ) from exc + + order += 1 + continue + + # evaluate right-hand side + b = sym.Matrix([ + target_expr.xreplace(dict(zip(dvec, lu.points[:, i], strict=True))) + for i in range(lu.points.shape[1]) + ]) + + # solve + all_coeffs = sym.solve_lu(lu.L, lu.U, lu.permutation, b, + postprocess=lambda x: x.expand()) + + # gather all non-zero coefficients from the result + const = sym.Integer(0) + coeffs = [] + mis = [] + for i, coeff in enumerate(all_coeffs): + coeff = round_expr(evalf(coeff)) + if coeff == 0: + continue + + if i == 0: + const = coeff + logger.debug(" %s", const) + else: + mis.append(lu.mis[i - 1]) + coeffs.append(coeff) + logger.debug(" + %s*%s.diff%s", coeff, base_kernel, lu.mis[-1]) + + if coeffs: + coeffs.insert(0, const) + + success = _check_linear_combination( + target_expr, + base_expr, + dvec, + mis, + coeffs, + rng=rng, + ) + if coeffs and success: + break + + coeffs = [] + order += 1 + + if not coeffs: + raise RewriteFailedError( + f"could not express {target_kernel} in terms of {base_kernel}" + ) + + return LinearOperatorRepresentation( + target_kernel, base_kernel, mis, + tuple( + sym.to_pymbolic(simplify( + c * (target_scaling if i == 0 else (target_scaling / base_scaling)) + ), symbols=True) + for i, c in enumerate(coeffs) + ), + ) + + +# }}} diff --git a/sumpy/test/test_kernel_rewrite.py b/sumpy/test/test_kernel_rewrite.py new file mode 100644 index 00000000..5a9f800b --- /dev/null +++ b/sumpy/test/test_kernel_rewrite.py @@ -0,0 +1,297 @@ +from __future__ import annotations + + +__copyright__ = """ +Copyright (C) 2020 Isuru Fernando +Copyright (C) 2026 Alexandru Fikl +""" + +import logging + +import numpy as np +import pytest + +import sumpy.symbolic as sym +from sumpy.kernel import ( + BiharmonicKernel, + BrinkmanletComponentKernel, + ElasticityComponentKernel, + LaplaceKernel, + StokesletComponentKernel, + StressletComponentKernel, + YukawaKernel, +) +from sumpy.kernel_rewrite import ( + LinearOperatorRepresentation, + rewrite_using_base_kernel_fourier, + rewrite_using_base_kernel_lu, +) + + +logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO) + + +def mi_derivative(expr: sym.Expr, x: sym.Matrix, mi: tuple[int, ...]) -> sym.Expr: + result = expr + for i, n in enumerate(mi): + result = result.diff(x[i], n) + + return result + + +def check_kernel_rewrite(op: LinearOperatorRepresentation) -> None: + from sumpy.kernel_rewrite import evalf, simplify + + dim = op.target_kernel.dim + dvec = sym.make_sym_vector("d", dim) + + target_expr = ( + op.target_kernel.get_global_scaling_const() + * op.target_kernel.get_expression(dvec)) + base_expr = ( + op.base_kernel.get_global_scaling_const() + * op.base_kernel.get_expression(dvec)) + + expr = sym.to_symbolic(op.coeffs[0], symbols=True) + sum(( + sym.to_symbolic(c, symbols=True) * mi_derivative(base_expr, dvec, mi) + for c, mi in zip(op.coeffs[1:], op.mis, strict=True) + ), sym.Integer(0)) + + result = evalf(simplify(target_expr - expr)) + assert abs(result) < 3.0 * 1.0e-16 + + +# {{{ test_rewrite_using_base_kernel_lu_conditioning + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_conditioning(dim: int) -> None: + from pytools import ( + generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, + ) + + from sumpy.kernel_rewrite import ( + _generate_points_shells, + _make_derivative_matrix, + _make_expr_derivatives, + ) + + rng = np.random.default_rng(42) + base_kernel = BiharmonicKernel(dim) + + dvec = sym.make_sym_vector("d", dim) + base_expr = base_kernel.get_expression(dvec) + + pde = base_kernel.get_pde_as_diff_op() + mis = list(gnitstam(pde.order, dim)) + pde_mis = [ident.mi for eq in pde.eqs for ident in eq] + pde_mis = [mi for mi in pde_mis if sum(mi) == pde.order] + mis.remove(pde_mis[-1]) + + mi_to_derivative = _make_expr_derivatives(base_expr, dvec, mis) + + nruns = 16 + kappa = np.empty(nruns) + + for i in range(nruns): + points = _generate_points_shells(dim, len(mis) + 1, rng=rng) + mat = _make_derivative_matrix(points, dvec, mis, mi_to_derivative) + + mat = np.array([ + [float(mat[i, j]) for j in range(mat.shape[1])] + for i in range(mat.shape[0]) + ]) + + kappa[i] = np.linalg.cond(mat) + logger.info("kappa = %.8e", kappa[i]) + + logger.info("median: %.8e max %.8e", np.median(kappa), np.max(kappa)) + assert np.max(kappa) < 3.0e+5 + +# }}} + + +# {{{ test_rewrite_using_base_kernel_lu_laplace_biharmonic + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_laplace_biharmonic(dim: int) -> None: + rng = np.random.default_rng(seed=42) + + base_kernel = BiharmonicKernel(dim) + target_kernel = LaplaceKernel(dim) + result = rewrite_using_base_kernel_lu(target_kernel, base_kernel, rng=rng) + + print(result.pretty()) + check_kernel_rewrite(result) + +# }}} + + +# {{{ test_rewrite_using_base_kernel_lu_stokeslet_biharmonic + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_stokeslet_biharmonic(dim: int) -> None: + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + rng = np.random.default_rng(seed=42) + + base_kernel = BiharmonicKernel(dim) + for i, j in gnitb(dim, 2): + target_kernel = StokesletComponentKernel(dim, i, j, viscosity_mu_name="mu") + result = rewrite_using_base_kernel_lu(target_kernel, base_kernel, rng=rng) + print(result.pretty()) + check_kernel_rewrite(result) + + +# }}} + + +# {{{ test_rewrite_using_base_kernel_lu_brinkman_yukawa + +@pytest.mark.skip(reason="needs two bases (Laplace and Yukawa)") +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_brinkman_yukawa(dim: int) -> None: + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + rng = np.random.default_rng(seed=42) + + base_kernel = YukawaKernel(dim, yukawa_lambda_name="k") + for i, j in gnitb(dim, 2): + target_kernel = BrinkmanletComponentKernel( + dim, i, j, viscosity_mu_name="mu", darcy_impermeability_name="k") + result = rewrite_using_base_kernel_lu(target_kernel, base_kernel, rng=rng) + print(result.pretty()) + check_kernel_rewrite(result) + + +# }}} + + +# {{{ test_rewrite_using_base_kernel_lu_stresslet_biharmonic + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_stresslet_biharmonic(dim: int) -> None: + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + rng = np.random.default_rng(seed=42) + + base_kernel = BiharmonicKernel(dim) + for i, j, k in gnitb(dim, 3): + target_kernel = StressletComponentKernel(dim, i, j, k, viscosity_mu_name="mu") + result = rewrite_using_base_kernel_lu(target_kernel, base_kernel, rng=rng) + print(result.pretty()) + check_kernel_rewrite(result) + + +# }}} + + +# {{{ test_rewrite_using_base_kernel_lu_elasticity_biharmonic + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_elasticity_biharmonic(dim: int) -> None: + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + rng = np.random.default_rng(seed=42) + + base_kernel = BiharmonicKernel(dim) + for i, j in gnitb(dim, 2): + target_kernel = ElasticityComponentKernel( + dim, i, j, viscosity_mu_name="mu", poisson_ratio_name="nu" + ) + result = rewrite_using_base_kernel_lu(target_kernel, base_kernel, rng=rng) + print(result.pretty()) + check_kernel_rewrite(result) + +# }}} + + +# {{{ test_rewrite_using_base_kernel_fourier_laplace_biharmonic + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_fourier_laplace_biharmonic(dim: int) -> None: + """Test that the Fourier-based algorithm recovers Laplace from biharmonic.""" + base_kernel = BiharmonicKernel(dim) + target_kernel = LaplaceKernel(dim) + result = rewrite_using_base_kernel_fourier(target_kernel, base_kernel) + + logger.info(result.pretty()) + check_kernel_rewrite(result) + + +# }}} + + +# {{{ test_rewrite_using_base_kernel_fourier_stokeslet_biharmonic + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_fourier_stokeslet_biharmonic(dim: int) -> None: + """Test that the Fourier-based algorithm recovers the Stokeslet from biharmonic.""" + from itertools import product + + base_kernel = BiharmonicKernel(dim) + + for i, j in product(range(dim), repeat=2): + target_kernel = StokesletComponentKernel(dim, i, j, viscosity_mu_name="mu") + result = rewrite_using_base_kernel_fourier(target_kernel, base_kernel) + + logger.info(result.pretty()) + check_kernel_rewrite(result) + +# }}} + + +# {{{ test_rewrite_using_base_kernel_fourier_brinkman_yukawa + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_fourier_brinkman_yukawa(dim: int) -> None: + from itertools import product + + base_kernel = YukawaKernel(dim, yukawa_lambda_name="k") + + for i, j in product(range(dim), repeat=2): + target_kernel = BrinkmanletComponentKernel( + dim, i, j, viscosity_mu_name="mu", darcy_impermeability_name="k") + result = rewrite_using_base_kernel_fourier(target_kernel, base_kernel) + + logger.info(result.pretty()) + check_kernel_rewrite(result) + +# }}} + + +# {{{ test_rewrite_using_base_kernel_fourier_indivisible + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_fourier_biharmonic_laplace(dim: int) -> None: + from sumpy.kernel_rewrite import ( + RewriteFailedError, + rewrite_using_base_kernel_fourier, + ) + + # Laplace Fourier symbol (-|k|^2) is not divisible by biharmonic (|k|^4) + base_kernel = LaplaceKernel(dim) + target_kernel = BiharmonicKernel(dim) + + with pytest.raises(RewriteFailedError, match="cannot rewrite"): + rewrite_using_base_kernel_fourier(target_kernel, base_kernel) + + +# }}} + + +if __name__ == "__main__": + import sys + + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + pytest.main([__file__]) + +# vim: fdm=marker