diff --git a/doc/expansion.rst b/doc/expansion.rst index 5d72d735a..ea2680340 100644 --- a/doc/expansion.rst +++ b/doc/expansion.rst @@ -27,3 +27,8 @@ Estimating Expansion Orders --------------------------- .. automodule:: sumpy.expansion.level_to_order + +Recurrences +----------- + +.. automodule:: sumpy.recurrence diff --git a/sumpy/recurrence.py b/sumpy/recurrence.py new file mode 100644 index 000000000..916b9d756 --- /dev/null +++ b/sumpy/recurrence.py @@ -0,0 +1,401 @@ +r""" +With the functionality in this module, we aim to compute a recurrence for +one-dimensional derivatives of functions :math:`f:\mathbb R^n \to \mathbb R` +for functions satisfying two assumptions: + +- :math:`f` satisfies a PDE that is linear and has coefficients polynomial + in the coordinates. +- :math:`f` only depends on the radius :math:`r`, + i.e. :math:`f(\boldsymbol x)=f(|\boldsymbol x|_2)`. + +This process proceeds in multiple steps: + +- Convert from the PDE to an ODE in :math:`r`, using :func:`pde_to_ode_in_r`. +- Convert from an ODE in :math:`r` to one in :math:`x`, +using :func:`ode_in_r_to_x`. +- Sort general-form ODE in :math:`x` into a coefficient array, using + :func:`ode_in_x_to_coeff_array`. +- Finally, get an expression for the recurrence, using + :func:`recurrence_from_coeff_array`. + +The whole process can be automated using :func:`recurrence_from_pde`. + +.. autofunction:: pde_to_ode_in_r +.. autofunction:: ode_in_r_to_x +.. autofunction:: ode_in_x_to_coeff_array +.. autofunction:: recurrence_from_coeff_array +.. autofunction:: recurrence_from_pde +.. autofunction:: process_recurrence_relation +.. autofunction:: shift_recurrence +.. autofunction:: get_processed_and_shifted_recurrence +""" + +from __future__ import annotations + + +__copyright__ = """ +Copyright (C) 2024 Hirish Chandrasekaran +Copyright (C) 2024 Andreas Kloeckner +""" + +__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 math +from typing import TypeVar + +import numpy as np +import sympy as sp + +from pytools.obj_array import make_obj_array + +from sumpy.expansion.diff_op import ( + DerivativeIdentifier, + LinearPDESystemOperator, +) + + +# similar to make_sym_vector in sumpy.symbolic, but returns an object array +# instead of a sympy.Matrix. +def _make_sympy_vec(name, n): + return make_obj_array([sp.Symbol(f"{name}{i}") for i in range(n)]) + + +def pde_to_ode_in_r(pde: LinearPDESystemOperator) -> tuple[ + sp.Expr, np.ndarray, int +]: + r""" + Returns an ODE satisfied by the radial derivatives of a function + :math:`f:\mathbb R^n \to \mathbb R` satisfying + :math:`f(\boldsymbol x)=f(|\boldsymbol x|_2)` and *pde*. + + :arg pde: must satisfy ``pde.eqs == 1`` and have polynomial coefficients. + + :returns: a tuple ``(ode_in_r, var, ode_order)``, where + - *ode_in_r* with derivatives given as :class:`sympy.Derivative` + - *var* is an object array of :class:`sympy.Symbol`, with successive + variables + representing the Cartesian coordinate directions. + - *ode_order* the order of ODE that is returned + """ + if len(pde.eqs) != 1: + raise ValueError("PDE must be scalar") + + dim = pde.dim + ode_order = pde.order + pde_eqn, = pde.eqs + + var = _make_sympy_vec("x", dim) + r = sp.sqrt(sum(var**2)) + eps = sp.symbols("epsilon") + rval = r + eps + f = sp.Function("f") + + def apply_deriv_id(expr: sp.Expr, + deriv_id: DerivativeIdentifier) -> sp.Expr: + for i, nderivs in enumerate(deriv_id.mi): + expr = expr.diff(var[i], nderivs) + return expr + + ode_in_r = sum( + # pylint: disable-next=not-callable + coeff * apply_deriv_id(f(rval), deriv_id) + for deriv_id, coeff in pde_eqn.items() + ) + + f_r_derivs = _make_sympy_vec("f_r", ode_order+1) + # pylint: disable-next=not-callable + f_derivs = [sp.diff(f(rval), eps, i) for i in range(ode_order+1)] + + # PDE ORDER = ODE ORDER + for i in range(ode_order+1): + ode_in_r = ode_in_r.subs(f_derivs[i], f_r_derivs[i]) + + return ode_in_r, var, ode_order + + +def _generate_nd_derivative_relations(var: np.ndarray, ode_order: int) -> dict: + r""" + Using the chain rule outputs a vector that gives in each component + respectively + :math:`[f(r), f'(r), \dots, f^{(ode_order)}(r)]` as a linear combination of + :math:`[f(x), f'(x), \dots, f^{(ode_order)}(x)]` + + :arg var: array of sympy variables math:`[x_0, x_1, \dots]` + :arg ode_order: the order of the ODE that we will be translating + """ + f_r_derivs = _make_sympy_vec("f_r", ode_order+1) + f_x_derivs = _make_sympy_vec("f_x", ode_order+1) + f = sp.Function("f") + eps = sp.symbols("epsilon") + rval = sp.sqrt(sum(var**2)) + eps + # pylint: disable=not-callable + f_derivs_x = [sp.diff(f(rval), var[0], i) for i in range(ode_order+1)] + f_derivs = [sp.diff(f(rval), eps, i) for i in range(ode_order+1)] + # pylint: disable=not-callable + for i in range(len(f_derivs_x)): + for j in range(len(f_derivs)): + f_derivs_x[i] = f_derivs_x[i].subs(f_derivs[j], f_r_derivs[j]) + system = [f_x_derivs[i] - f_derivs_x[i] for i in range(ode_order+1)] + return sp.solve(system, *f_r_derivs, dict=True)[0] + + +def ode_in_r_to_x(ode_in_r: sp.Expr, var: np.ndarray, + ode_order: int) -> sp.Expr: + r""" + Translates an ode in the variable r into an ode in the variable x + by replacing the terms :math:`f, f_r, f_{rr}, \dots` as a linear + combinations of + :math:`f, f_x, f_{xx}, \dots` using the chain rule. + + :arg ode_in_r: a linear combination of :math:`f, f_r, f_{rr}, \dots` + represented by the sympy variables :math:`f_{r0}, f_{r1}, f_{r2}, \dots` + :arg var: array of sympy variables :math:`[x_0, x_1, \dots]` + :arg ode_order: the order of the input ODE + + :returns: *ode_in_x* a linear combination of :math:`f, f_x, f_{xx}, \dots` + represented by the sympy variables :math:`f_{x0}, f_{x1}, f_{x2}, + \dots` with coefficients as rational functions in + :math:`x_0, x_1, \dots` + """ + subme = _generate_nd_derivative_relations(var, ode_order+1) + ode_in_x = ode_in_r + f_r_derivs = _make_sympy_vec("f_r", ode_order+1) + for i in range(ode_order+1): + ode_in_x = ode_in_x.subs(f_r_derivs[i], subme[f_r_derivs[i]]) + return ode_in_x + + +ODECoefficients = list[list[sp.Expr]] + + +def ode_in_x_to_coeff_array(poly: sp.Poly, ode_order: int, var: + np.ndarray) -> ODECoefficients: + r""" + Organizes the coefficients of an ODE in the :math:`x_0` variable into a + 2D array. + + :arg poly: a sympy polynomial in + :math:`\partial_{x_0}^0 f, \partial_{x_0}^1 f,\cdots` of the form + :math:`(b_{00} x_0^0 + b_{01} x_0^1 + \cdots) \partial_{x_0}^0 f + + (b_{10} x_0^0 + b_{11} x_0^1 +\cdots) \partial_x^1 f` + + :arg var: array of sympy variables :math:`[x_0, x_1, \dots]` + :arg ode_order: the order of the input ODE we return a sequence + + :returns: *coeffs* a sequence of of sequences, with the outer sequence + iterating over derivative orders, and each inner sequence iterating + over powers of :math:`x_0`, so that, in terms of the above form, + coeffs is :math:`[[b_{00}, b_{01}, ...], [b_{10}, b_{11}, ...], ...]` + """ + return [ + # recast ODE coefficient obtained below as polynomial in x0 + sp.Poly( + # get coefficient of deriv_ind'th derivative + poly.coeff_monomial(poly.gens[deriv_ind]), + + var[0]) + # get poly coefficients in /ascending/ order + .all_coeffs()[::-1] + for deriv_ind in range(ode_order+1)] + + +NumberT = TypeVar("NumberT", int, float, complex) + + +def _falling_factorial(arg: NumberT, num_terms: int) -> NumberT: + result = 1 + for i in range(num_terms): + result = result * (arg - i) + return result + + +def _auto_product_rule_single_term(p: int, m: int, var: np.ndarray) -> sp.Expr: + r""" + We assume that we are given the expression :math:`x_0^p f^(m)(x_0)`. We + then output the nth order derivative of the expression where :math:`n` is + a symbolic variable. + We let :math:`s(i)` represent the ith order derivative of f when + we output the final result. + :arg var: array of sympy variables :math:`[x_0, x_1, \dots]` + """ + n = sp.symbols("n") + s = sp.Function("s") + + return sum( + # pylint: disable=not-callable + _falling_factorial(n, i) + * math.comb(p, i) * s(n-i+m) * var[0]**(p-i) + for i in range(p+1) + ) + + +def recurrence_from_coeff_array(coeffs: list, var: np.ndarray) -> sp.Expr: + r""" + A function that takes in as input an organized 2D coefficient array (see + above) and outputs a recurrence relation. + + :arg coeffs: a sequence of of sequences, described in + :func:`ode_in_x_to_coeff_array` + :arg var: array of sympy variables :math:`[x_0, x_1, \dots]` + """ + final_recurrence = 0 + # Outer loop is derivative direction + # Inner is polynomial order of x_0 + for m, _ in enumerate(coeffs): + for p, _ in enumerate(coeffs[m]): + final_recurrence += coeffs[m][p] * _auto_product_rule_single_term( + p, m, var) + return final_recurrence + + +def recurrence_from_pde(pde: LinearPDESystemOperator) -> sp.Expr: + r""" + A function that takes in as input a sympy PDE and outputs a recurrence + relation. + + :arg pde: a :class:`sumpy.expansion.diff_op.LinearSystemPDEOperator` + that must satisfy ``pde.eqs == 1`` and have polynomial coefficients + in the coordinates. + :arg var: array of sympy variables :math:`[x_0, x_1, \dots]` + """ + ode_in_r, var, ode_order = pde_to_ode_in_r(pde) + ode_in_x = ode_in_r_to_x(ode_in_r, var, ode_order).simplify() + ode_in_x_cleared = (ode_in_x * var[0]**(ode_order+1)).simplify() + # ode_in_x_cleared shouldn't have rational function coefficients + assert sp.together(ode_in_x_cleared) == ode_in_x_cleared + f_x_derivs = _make_sympy_vec("f_x", ode_order+1) + poly = sp.Poly(ode_in_x_cleared, *f_x_derivs) + coeffs = ode_in_x_to_coeff_array(poly, ode_order, var) + return recurrence_from_coeff_array(coeffs, var) + + +def process_recurrence_relation(r: sp.Expr) -> tuple[int, sp.Expr]: + r""" + A function that takes in as input a recurrence and outputs a recurrence + relation that has the nth term in terms of the n-1th, n-2th etc. + Also returns the order of the recurrence relation. + + :arg recurrence: a recurrence relation in :math:`s(n)` + """ + idx_l, terms = _extract_idx_terms_from_recurrence(r) + # Order is the max difference between highest/lowest in idx_l + order = max(idx_l) - min(idx_l) + 1 + + # How much do we need to shift the recurrence relation + shift_idx = max(idx_l) + + # Get the respective coefficients in the recurrence relation from r + n = sp.symbols("n") + coeffs = sp.poly(r, list(terms)).coeffs() + + # Re-arrange the recurrence relation so we get s(n) = ____ + # in terms of s(n-1), ... + true_recurrence = sum(coeffs[i]/coeffs[-1] * terms[i] + for i in range(0, len(terms)-1)) + true_recurrence1 = true_recurrence.subs(n, n-shift_idx) + + return order, true_recurrence1 + + +def _extract_idx_terms_from_recurrence(r: sp.Expr) -> tuple[np.ndarray, + np.ndarray]: + r""" + Given a recurrence extracts the variables in the recurrence + as well as the indexes, both in sorted order. + + :arg r: recurrence to extract terms from + """ + # We're assuming here that s(...) are the only function calls. + terms = list(r.atoms(sp.Function)) + terms = np.array(terms) + + idx_l = [] + for i in range(len(terms)): + tms = list(terms[i].atoms(sp.Number)) + if len(tms) == 1: + idx_l.append(tms[0]) + else: + idx_l.append(0) + idx_l = np.array(idx_l, dtype="int") + idx_sort = idx_l.argsort() + idx_l = idx_l[idx_sort] + terms = terms[idx_sort] + + return idx_l, terms + + +def _check_neg_ind(r_n): + r""" + Simply checks if a negative index exists in a recurrence relation. + """ + + idx_l, _ = _extract_idx_terms_from_recurrence(r_n) + + return np.any(idx_l < 0) + + +def _get_initial_c(recurrence): + r""" + For a given recurrence checks how many initial conditions by + checking for non-negative indexed terms. + """ + n = sp.symbols("n") + + i = 0 + r_c = recurrence.subs(n, i) + while _check_neg_ind(r_c): + i += 1 + r_c = recurrence.subs(n, i) + return i + + +def shift_recurrence(r: sp.Expr) -> sp.Expr: + r""" + A function that "shifts" the recurrence so it's center is placed + at the origin and source is the input for the recurrence generated. + + :arg recurrence: a recurrence relation in :math:`s(n)` + """ + idx_l, terms = _extract_idx_terms_from_recurrence(r) + + r_ret = r + + n = sp.symbols("n") + for i in range(len(idx_l)): + r_ret = r_ret.subs(terms[i], (-1)**(n+idx_l[i])*terms[i]) + + return r_ret*((-1)**(n+1)) + + +def get_processed_and_shifted_recurrence(pde) -> tuple[int, int, + sp.Expr]: + r""" + A function that "shifts" the recurrence so the expansion center is placed + at the origin and source is the input for the recurrence generated. + + Also processes the recurrence so s(n) is in terms of s(n-1), etc. + + :arg recurrence: a recurrence relation in :math:`s(n)` + """ + r = recurrence_from_pde(pde) + order, r_p = process_recurrence_relation(r) + n_initial = _get_initial_c(r_p) + r_s = shift_recurrence(r_p) + return n_initial, order, r_s diff --git a/sumpy/recurrenceqbx.py b/sumpy/recurrenceqbx.py new file mode 100644 index 000000000..0407ba02b --- /dev/null +++ b/sumpy/recurrenceqbx.py @@ -0,0 +1,115 @@ +r""" +With the functionality in this module, we aim to compute layer potentials +using a recurrence for one-dimensional derivatives of the corresponding +Green's function. See recurrence.py. + +.. autofunction:: recurrence_qbx_lp +""" +from __future__ import annotations # noqa: I001 + +import math +from typing import Sequence + +import numpy as np +import sympy as sp + +from sumpy.recurrence import ( + _make_sympy_vec, + get_processed_and_shifted_recurrence +) + + +# ================ Transform/Rotate ================= +def _produce_orthogonal_basis(normals: np.ndarray) -> Sequence[np.ndarray]: + ndim, ncenters = normals.shape + orth_coordsys = [normals] + for i in range(1, ndim): + v = np.random.rand(ndim, ncenters) # noqa: NPY002 + v = v/np.linalg.norm(v, 2, axis=0) + for j in range(i): + v = v - np.einsum("dc,dc->c", v, + orth_coordsys[j]) * orth_coordsys[j] + v = v/np.linalg.norm(v, 2, axis=0) + orth_coordsys.append(v) + + return orth_coordsys + + +def _compute_rotated_shifted_coordinates( + sources: np.ndarray, + centers: np.ndarray, + normals: np.ndarray + ) -> np.ndarray: + cts = sources[:, None] - centers[:, :, None] + orth_coordsys = _produce_orthogonal_basis(normals) + cts_rotated_shifted = np.einsum("idc,dcs->ics", orth_coordsys, cts) + + return cts_rotated_shifted + + +# ================ Recurrence LP Eval ================= +def recurrence_qbx_lp(sources, centers, normals, strengths, radius, pde, g_x_y, + ndim, p) -> np.ndarray: + r""" + A function that computes a single-layer potential using a recurrence. + + :arg sources: a (ndim, nsources) array of source locations + :arg centers: a (ndim, ncenters) array of center locations + :arg normals: a (ndim, ncenters) array of normals + :arg strengths: array corresponding to quadrature weight multiplied by + density + :arg radius: expansion radius + :arg pde: pde that we are computing layer potential for + :arg g_x_y: a green's function in (x0, x1, ...) source and + (t0, t1, ...) target + :arg ndim: number of spatial variables + :arg p: order of expansion computed + """ + + # ------------- 2. Compute rotated/shifted coordinates + cts_r_s = _compute_rotated_shifted_coordinates(sources, centers, normals) + + # ------------- 4. Define input variables for green's function expression + var = _make_sympy_vec("x", ndim) + var_t = _make_sympy_vec("t", ndim) + + # ------------ 5. Compute recurrence + n_initial, order, recurrence = get_processed_and_shifted_recurrence(pde) + + # ------------ 6. Set order p = 5 + n_p = sources.shape[1] + storage = [np.zeros((n_p, n_p))] * order + + s = sp.Function("s") + r, n = sp.symbols("r,n") + + def generate_lamb_expr(i, n_initial): + arg_list = [] + for j in range(order, 0, -1): + # pylint: disable-next=not-callable + arg_list.append(s(i-j)) + arg_list.append(var[0]) + arg_list.append(var[1]) + arg_list.append(r) + + if i < n_initial: + lamb_expr = sp.diff(g_x_y, var_t[0], i) + for j in range(ndim): + lamb_expr = lamb_expr.subs(var_t[j], 0) + else: + lamb_expr = recurrence.subs(n, i) + return sp.lambdify(arg_list, lamb_expr) + + interactions = 0 + for i in range(p+1): + lamb_expr = generate_lamb_expr(i, n_initial) + a = [*storage, cts_r_s[0], cts_r_s[1], radius] + s_new = lamb_expr(*a) + interactions += s_new * radius**i/math.factorial(i) + + storage.pop(0) + storage.append(s_new) + + exp_res = (interactions * strengths[None, :]).sum(axis=1) + + return exp_res diff --git a/test/test_recurrence.py b/test/test_recurrence.py new file mode 100644 index 000000000..8ad80f3db --- /dev/null +++ b/test/test_recurrence.py @@ -0,0 +1,195 @@ +r""" +With the functionality in this module, we aim to test recurrence +code. + +.. autofunction:: test_laplace3d +.. autofunction:: test_helmholtz3d +.. autofunction:: test_laplace2d +""" +from __future__ import annotations + +import numpy as np +import sympy as sp +from sympy import hankel1 + +from sumpy.expansion.diff_op import ( + laplacian, + make_identity_diff_op, +) +from sumpy.recurrence import _make_sympy_vec, get_processed_and_shifted_recurrence + + +def test_laplace3d(): + r""" + Tests recurrence code for orders up to 6 laplace3d. + """ + w = make_identity_diff_op(3) + laplace3d = laplacian(w) + _, _, r = get_processed_and_shifted_recurrence(laplace3d) + n = sp.symbols("n") + s = sp.Function("s") + + var = _make_sympy_vec("x", 3) + var_t = _make_sympy_vec("t", 3) + abs_dist = sp.sqrt((var[0]-var_t[0])**2 + + (var[1]-var_t[1])**2 + (var[2]-var_t[2])**2) + g_x_y = 1/abs_dist + derivs = [sp.diff(g_x_y, + var_t[0], i).subs(var_t[0], 0).subs(var_t[1], 0).subs(var_t[2], 0) + for i in range(6)] + + # pylint: disable-next=not-callable + subs_dict = {s(0): derivs[0], s(1): derivs[1]} + check_2_s = r.subs(n, 2).subs(subs_dict) - derivs[2] + # pylint: disable-next=not-callable + subs_dict[s(2)] = derivs[2] + check_3_s = r.subs(n, 3).subs(subs_dict) - derivs[3] + # pylint: disable-next=not-callable + subs_dict[s(3)] = derivs[3] + check_4_s = r.subs(n, 4).subs(subs_dict) - derivs[4] + # pylint: disable-next=not-callable + subs_dict[s(4)] = derivs[4] + check_5_s = r.subs(n, 5).subs(subs_dict) - derivs[5] + + x_coord = np.random.rand() # noqa: NPY002 + y_coord = np.random.rand() # noqa: NPY002 + z_coord = np.random.rand() # noqa: NPY002 + coord_dict = {var[0]: x_coord, var[1]: y_coord, var[2]: z_coord} + + assert abs(check_2_s.subs(coord_dict)) <= 1e-15 + assert abs(check_3_s.subs(coord_dict)) <= 1e-14 + assert abs(check_4_s.subs(coord_dict)) <= 1e-12 + assert abs(check_5_s.subs(coord_dict)) <= 1e-12 + + +def test_helmholtz3d(): + r""" + Tests recurrence code for orders up to 6 helmholtz3d. + """ + w = make_identity_diff_op(3) + helmholtz3d = laplacian(w) + w + _, _, r = get_processed_and_shifted_recurrence(helmholtz3d) + + n = sp.symbols("n") + s = sp.Function("s") + + var = _make_sympy_vec("x", 3) + var_t = _make_sympy_vec("t", 3) + abs_dist = sp.sqrt((var[0]-var_t[0])**2 + + (var[1]-var_t[1])**2 + (var[2]-var_t[2])**2) + g_x_y = sp.exp(1j * abs_dist) / abs_dist + derivs = [sp.diff(g_x_y, + var_t[0], i).subs(var_t[0], 0).subs(var_t[1], 0).subs(var_t[2], 0) + for i in range(6)] + + # pylint: disable-next=not-callable + subs_dict = {s(0): derivs[0], s(1): derivs[1]} + check_2_s = r.subs(n, 2).subs(subs_dict) - derivs[2] + # pylint: disable-next=not-callable + subs_dict[s(2)] = derivs[2] + check_3_s = r.subs(n, 3).subs(subs_dict) - derivs[3] + # pylint: disable-next=not-callable + subs_dict[s(3)] = derivs[3] + check_4_s = r.subs(n, 4).subs(subs_dict) - derivs[4] + # pylint: disable-next=not-callable + subs_dict[s(4)] = derivs[4] + check_5_s = r.subs(n, 5).subs(subs_dict) - derivs[5] + + x_coord = np.random.rand() # noqa: NPY002 + y_coord = np.random.rand() # noqa: NPY002 + z_coord = np.random.rand() # noqa: NPY002 + coord_dict = {var[0]: x_coord, var[1]: y_coord, var[2]: z_coord} + + assert abs(abs(check_2_s.subs(coord_dict))) <= 1e-15 + assert abs(abs(check_3_s.subs(coord_dict))) <= 1e-14 + assert abs(abs(check_4_s.subs(coord_dict))) <= 1e-12 + assert abs(abs(check_5_s.subs(coord_dict))) <= 1e-12 + + +def test_helmholtz2d(): + r""" + Tests recurrence code for orders up to 6 helmholtz2d. + """ + w = make_identity_diff_op(2) + helmholtz2d = laplacian(w) + w + _, _, r = get_processed_and_shifted_recurrence(helmholtz2d) + + n = sp.symbols("n") + s = sp.Function("s") + + var = _make_sympy_vec("x", 2) + var_t = _make_sympy_vec("t", 2) + abs_dist = sp.sqrt((var[0]-var_t[0])**2 + + (var[1]-var_t[1])**2) + k = 1 + g_x_y = (1j/4) * hankel1(0, k * abs_dist) + derivs = [sp.diff(g_x_y, + var_t[0], i).subs(var_t[0], 0).subs(var_t[1], 0) + for i in range(6)] + x_coord = np.random.rand() # noqa: NPY002 + y_coord = np.random.rand() # noqa: NPY002 + coord_dict = {var[0]: x_coord, var[1]: y_coord} + derivs = [derivs[i].subs(coord_dict) for i in range(6)] + + # pylint: disable-next=not-callable + subs_dict = {s(0): derivs[0], s(1): derivs[1]} + check_2_s = r.subs(n, 2).subs(subs_dict) - derivs[2] + # pylint: disable-next=not-callable + subs_dict[s(2)] = derivs[2] + check_3_s = r.subs(n, 3).subs(subs_dict) - derivs[3] + # pylint: disable-next=not-callable + subs_dict[s(3)] = derivs[3] + check_4_s = r.subs(n, 4).subs(subs_dict) - derivs[4] + # pylint: disable-next=not-callable + subs_dict[s(4)] = derivs[4] + check_5_s = r.subs(n, 5).subs(subs_dict) - derivs[5] + + f2 = sp.lambdify([var[0], var[1]], check_2_s) + assert abs(f2(x_coord, y_coord)) <= 1e-13 + f3 = sp.lambdify([var[0], var[1]], check_3_s) + assert abs(f3(x_coord, y_coord)) <= 1e-13 + f4 = sp.lambdify([var[0], var[1]], check_4_s) + assert abs(f4(x_coord, y_coord)) <= 1e-13 + f5 = sp.lambdify([var[0], var[1]], check_5_s) + assert abs(f5(x_coord, y_coord)) <= 1e-12 + + +def test_laplace2d(): + r""" + Tests recurrence code for orders up to 6 laplace2d. + """ + w = make_identity_diff_op(2) + laplace2d = laplacian(w) + _, _, r = get_processed_and_shifted_recurrence(laplace2d) + + n = sp.symbols("n") + s = sp.Function("s") + + var = _make_sympy_vec("x", 2) + var_t = _make_sympy_vec("t", 2) + g_x_y = sp.log(sp.sqrt((var[0]-var_t[0])**2 + (var[1]-var_t[1])**2)) + derivs = [sp.diff(g_x_y, + var_t[0], i).subs(var_t[0], 0).subs(var_t[1], 0) + for i in range(6)] + + # pylint: disable-next=not-callable + subs_dict = {s(0): derivs[0], s(1): derivs[1]} + check_2_s = r.subs(n, 2).subs(subs_dict) - derivs[2] + # pylint: disable-next=not-callable + subs_dict[s(2)] = derivs[2] + check_3_s = r.subs(n, 3).subs(subs_dict) - derivs[3] + # pylint: disable-next=not-callable + subs_dict[s(3)] = derivs[3] + check_4_s = r.subs(n, 4).subs(subs_dict) - derivs[4] + # pylint: disable-next=not-callable + subs_dict[s(4)] = derivs[4] + check_5_s = r.subs(n, 5).subs(subs_dict) - derivs[5] + + x_coord = np.random.rand() # noqa: NPY002 + y_coord = np.random.rand() # noqa: NPY002 + coord_dict = {var[0]: x_coord, var[1]: y_coord} + + assert abs(abs(check_2_s.subs(coord_dict))) <= 1e-15 + assert abs(abs(check_3_s.subs(coord_dict))) <= 1e-14 + assert abs(abs(check_4_s.subs(coord_dict))) <= 1e-12 + assert abs(abs(check_5_s.subs(coord_dict))) <= 1e-12 diff --git a/test/test_recurrenceqbx.py b/test/test_recurrenceqbx.py new file mode 100644 index 000000000..788f41494 --- /dev/null +++ b/test/test_recurrenceqbx.py @@ -0,0 +1,151 @@ +r""" +With the functionality in this module, we aim to test recurrence ++ qbx code. +""" +from __future__ import annotations + +import numpy as np +import sympy as sp + +from sympy import hankel1 +from sumpy.array_context import _acf +from sumpy.expansion.diff_op import ( + laplacian, + make_identity_diff_op, +) +from sumpy.expansion.local import LineTaylorLocalExpansion +from sumpy.kernel import HelmholtzKernel, LaplaceKernel +from sumpy.qbx import LayerPotential +from sumpy.recurrenceqbx import _make_sympy_vec, recurrence_qbx_lp + + +actx_factory = _acf +ExpnClass = LineTaylorLocalExpansion + +actx = actx_factory() +lknl = LaplaceKernel(2) +hlknl = HelmholtzKernel(2, "k") + + +def _qbx_lp_helmholtz_general(sources, targets, centers, radius, strengths, order): + lpot = LayerPotential(actx.context, + expansion=ExpnClass(hlknl, order), + target_kernels=(hlknl,), + source_kernels=(hlknl,)) + + # print(lpot.get_kernel()) + expansion_radii = actx.from_numpy(radius * np.ones(sources.shape[1])) + sources = actx.from_numpy(sources) + targets = actx.from_numpy(targets) + centers = actx.from_numpy(centers) + + strengths = (strengths,) + extra_kernel_kwargs = {"k": 1} + _evt, (result_qbx,) = lpot( + actx.queue, + targets, sources, centers, strengths, + expansion_radii=expansion_radii, + k=1) + result_qbx = actx.to_numpy(result_qbx) + + return result_qbx + + +def _qbx_lp_laplace_general(sources, targets, centers, radius, strengths, order): + lpot = LayerPotential(actx.context, + expansion=ExpnClass(lknl, order), + target_kernels=(lknl,), + source_kernels=(lknl,)) + + # print(lpot.get_kernel()) + expansion_radii = actx.from_numpy(radius * np.ones(sources.shape[1])) + sources = actx.from_numpy(sources) + targets = actx.from_numpy(targets) + centers = actx.from_numpy(centers) + + strengths = (strengths,) + + _evt, (result_qbx,) = lpot( + actx.queue, + targets, sources, centers, strengths, + expansion_radii=expansion_radii) + result_qbx = actx.to_numpy(result_qbx) + + return result_qbx + + +def _create_ellipse(n_p): + h = 9.688 / n_p + radius = 7*h + t = np.linspace(0, 2 * np.pi, n_p, endpoint=False) + + unit_circle_param = np.exp(1j * t) + unit_circle = np.array([2 * unit_circle_param.real, unit_circle_param.imag]) + + sources = unit_circle + normals = np.array([unit_circle_param.real, 2*unit_circle_param.imag]) + normals = normals / np.linalg.norm(normals, axis=0) + centers = sources - normals * radius + + mode_nr = 25 + density = np.cos(mode_nr * t) + + return sources, centers, normals, density, h, radius + + +def test_recurrence_laplace_2d_ellipse(): + r""" + Tests recurrence + qbx code. + """ + + # ------------- 1. Define PDE, Green's Function + w = make_identity_diff_op(2) + laplace2d = laplacian(w) + + var = _make_sympy_vec("x", 2) + var_t = _make_sympy_vec("t", 2) + g_x_y = (-1/(2*np.pi)) * sp.log(sp.sqrt((var[0]-var_t[0])**2 + + (var[1]-var_t[1])**2)) + + p = 4 + err = [] + for n_p in range(200, 1001, 200): + sources, centers, normals, density, h, radius = _create_ellipse(n_p) + strengths = h * density + exp_res = recurrence_qbx_lp(sources, centers, normals, + strengths, radius, laplace2d, + g_x_y, 2, p) + qbx_res = _qbx_lp_laplace_general(sources, sources, centers, + radius, strengths, p) + # qbx_res,_ = lpot_eval_circle(sources.shape[1], p) + err.append(np.max(np.abs(exp_res - qbx_res))) + assert np.max(err) <= 1e-13 + + +def test_recurrence_helmholtz_2d_ellipse(): + r""" + Tests recurrence + qbx code. + """ + # ------------- 1. Define PDE, Green's Function + w = make_identity_diff_op(2) + helmholtz2d = laplacian(w) + w + + var = _make_sympy_vec("x", 2) + var_t = _make_sympy_vec("t", 2) + k = 1 + abs_dist = sp.sqrt((var[0]-var_t[0])**2 + (var[1]-var_t[1])**2) + g_x_y = (1j/4) * hankel1(0, k * abs_dist) + + p = 4 + err = [] + for n_p in range(200, 1001, 200): + sources, centers, normals, density, h, radius = _create_ellipse(n_p) + strengths = h * density + exp_res = recurrence_qbx_lp(sources, centers, normals, strengths, + radius, helmholtz2d, g_x_y, 2, p) + qbx_res = _qbx_lp_helmholtz_general(sources, sources, centers, radius, strengths, p) + #qbx_res,_ = lpot_eval_circle(sources.shape[1], p) + err.append(np.max(np.abs(exp_res - qbx_res))) + assert np.max(err) <= 1e-13 + +# test_recurrence_helmholtz_2d_ellipse()