1+ r"""
2+ With the functionality in this module, we compute layer potentials
3+ using a recurrence for one-dimensional derivatives of the corresponding
4+ Green's function. See recurrence.py.
5+
6+ .. autofunction:: recurrence_qbx_lp
7+ """
8+ from __future__ import annotations
9+
10+
11+ __copyright__ = """
12+ Copyright (C) 2024 Hirish Chandrasekaran
13+ Copyright (C) 2024 Andreas Kloeckner
14+ """
15+
16+ __license__ = """
17+ Permission is hereby granted, free of charge, to any person obtaining a copy
18+ of this software and associated documentation files (the "Software"), to deal
19+ in the Software without restriction, including without limitation the rights
20+ to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
21+ copies of the Software, and to permit persons to whom the Software is
22+ furnished to do so, subject to the following conditions:
23+
24+ The above copyright notice and this permission notice shall be included in
25+ all copies or substantial portions of the Software.
26+
27+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
28+ IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
29+ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
30+ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
31+ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
32+ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
33+ THE SOFTWARE.
34+ """
35+
36+ import math
37+ from typing import Sequence
38+
39+ import numpy as np
40+ import sympy as sp
41+
42+ from sumpy .recurrence import _make_sympy_vec , get_reindexed_and_center_origin_on_axis_recurrence
43+
44+
45+ # ================ Transform/Rotate =================
46+ def _produce_orthogonal_basis (normals : np .ndarray ) -> Sequence [np .ndarray ]:
47+ ndim , ncenters = normals .shape
48+ orth_coordsys = [normals ]
49+ for i in range (1 , ndim ):
50+ v = np .random .rand (ndim , ncenters ) # noqa: NPY002
51+ v = v / np .linalg .norm (v , 2 , axis = 0 )
52+ for j in range (i ):
53+ v = v - np .einsum ("dc,dc->c" , v ,
54+ orth_coordsys [j ]) * orth_coordsys [j ]
55+ v = v / np .linalg .norm (v , 2 , axis = 0 )
56+ orth_coordsys .append (v )
57+
58+ return orth_coordsys
59+
60+
61+ def _compute_rotated_shifted_coordinates (
62+ sources : np .ndarray ,
63+ centers : np .ndarray ,
64+ normals : np .ndarray
65+ ) -> np .ndarray :
66+ cts = sources [:, None ] - centers [:, :, None ]
67+ orth_coordsys = _produce_orthogonal_basis (normals )
68+ cts_rotated_shifted = np .einsum ("idc,dcs->ics" , orth_coordsys , cts )
69+
70+ return cts_rotated_shifted
71+
72+
73+ # ================ Recurrence LP Eval =================
74+ def recurrence_qbx_lp (sources , centers , normals , strengths , radius , pde , g_x_y ,
75+ ndim , p ) -> np .ndarray :
76+ r"""
77+ A function that computes a single-layer potential using a recurrence.
78+
79+ :arg sources: a (ndim, nsources) array of source locations
80+ :arg centers: a (ndim, ncenters) array of center locations
81+ :arg normals: a (ndim, ncenters) array of normals
82+ :arg strengths: array corresponding to quadrature weight multiplied by
83+ density
84+ :arg radius: expansion radius
85+ :arg pde: pde that we are computing layer potential for
86+ :arg g_x_y: a green's function in (x0, x1, ...) source and
87+ (t0, t1, ...) target
88+ :arg ndim: number of spatial variables
89+ :arg p: order of expansion computed
90+ """
91+
92+ # ------------- 2. Compute rotated/shifted coordinates
93+ cts_r_s = _compute_rotated_shifted_coordinates (sources , centers , normals )
94+
95+ # ------------- 4. Define input variables for green's function expression
96+ var = _make_sympy_vec ("x" , ndim )
97+ var_t = _make_sympy_vec ("t" , ndim )
98+
99+ # ------------ 5. Compute recurrence
100+ n_initial , order , recurrence = get_reindexed_and_center_origin_on_axis_recurrence (pde )
101+
102+ # ------------ 6. Set order p = 5
103+ n_p = sources .shape [1 ]
104+ storage = [np .zeros ((n_p , n_p ))] * order
105+
106+ s = sp .Function ("s" )
107+ n = sp .symbols ("n" )
108+
109+ def generate_lamb_expr (i , n_initial ):
110+ arg_list = []
111+ for j in range (order , 0 , - 1 ):
112+ # pylint: disable-next=not-callable
113+ arg_list .append (s (i - j ))
114+ for j in range (ndim ):
115+ arg_list .append (var [j ])
116+
117+ lamb_expr_symb_deriv = sp .diff (g_x_y , var_t [0 ], i )
118+ for j in range (ndim ):
119+ lamb_expr_symb_deriv = lamb_expr_symb_deriv .subs (var_t [j ], 0 )
120+
121+ if i < n_initial :
122+ lamb_expr_symb = lamb_expr_symb_deriv
123+ else :
124+ lamb_expr_symb = recurrence .subs (n , i )
125+ #print("=============== ORDER = " + str(i))
126+ #print(lamb_expr_symb)
127+ return sp .lambdify (arg_list , lamb_expr_symb ), sp .lambdify (arg_list , lamb_expr_symb_deriv )
128+
129+ interactions = 0
130+ coord = [cts_r_s [j ] for j in range (ndim )]
131+ for i in range (p + 1 ):
132+ lamb_expr , true_lamb_expr = generate_lamb_expr (i , n_initial )
133+ a = [* storage , * coord ]
134+ s_new = lamb_expr (* a )
135+ s_new_true = true_lamb_expr (* a )
136+ arg_max = np .argmax (abs (s_new - s_new_true )/ abs (s_new_true ))
137+ print ((s_new - s_new_true ).reshape (- 1 )[arg_max ]/ s_new_true .reshape (- 1 )[arg_max ])
138+ print ("x:" , coord [0 ].reshape (- 1 )[arg_max ], "y:" , coord [1 ].reshape (- 1 )[arg_max ],
139+ "s_recur:" , s_new .reshape (- 1 )[arg_max ], "s_true:" , s_new_true .reshape (- 1 )[arg_max ], "order: " , i )
140+ interactions += s_new * radius ** i / math .factorial (i )
141+
142+ storage .pop (0 )
143+ storage .append (s_new )
144+
145+ exp_res = (interactions * strengths [None , :]).sum (axis = 1 )
146+
147+ return exp_res
0 commit comments