Skip to content

Commit 2a7adcd

Browse files
ShawnL00inducer
ShawnL00
authored and
inducer
committed
Add tests for compressed coefficient difference formula for L2L, M2M
1 parent 118ae83 commit 2a7adcd

3 files changed

Lines changed: 516 additions & 0 deletions

File tree

‎sumpy/test/coeff_test_tools.py‎

Lines changed: 68 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,68 @@
1+
from __future__ import annotations
2+
3+
4+
__copyright__ = """
5+
Copyright (C) 2026 Shawn/Chaoqi Lin
6+
"""
7+
8+
__license__ = """
9+
Permission is hereby granted, free of charge, to any person obtaining a copy
10+
of this software and associated documentation files (the "Software"), to deal
11+
in the Software without restriction, including without limitation the rights
12+
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13+
copies of the Software, and to permit persons to whom the Software is
14+
furnished to do so, subject to the following conditions:
15+
16+
The above copyright notice and this permission notice shall be included in
17+
all copies or substantial portions of the Software.
18+
19+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20+
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21+
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22+
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23+
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24+
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
25+
THE SOFTWARE.
26+
"""
27+
28+
import numpy as np
29+
import sympy as sp
30+
31+
32+
def to_scalar(val):
33+
"""Convert symbolic or array value to scalar."""
34+
if hasattr(val, "evalf"):
35+
val = val.evalf()
36+
if hasattr(val, "item"):
37+
val = val.item()
38+
return complex(val)
39+
40+
41+
class NumericMatVecOperator:
42+
"""Wrapper for symbolic matrix-vector operator with numeric
43+
substitution."""
44+
45+
def __init__(self, symbolic_op, repl_dict):
46+
self.symbolic_op = symbolic_op
47+
self.repl_dict = repl_dict
48+
self.shape = symbolic_op.shape
49+
50+
def matvec(self, vec):
51+
result = self.symbolic_op.matvec(vec)
52+
out = []
53+
for expr in result:
54+
if hasattr(expr, "xreplace"):
55+
out.append(complex(expr.xreplace(self.repl_dict).evalf()))
56+
else:
57+
out.append(complex(expr))
58+
return np.array(out)
59+
60+
61+
def get_repl_dict(kernel, extra_kwargs):
62+
"""Numeric substitution for symbolic kernel parameters."""
63+
repl_dict = {}
64+
if "lam" in extra_kwargs:
65+
repl_dict[sp.Symbol("lam")] = extra_kwargs["lam"]
66+
if "k" in extra_kwargs:
67+
repl_dict[sp.Symbol("k")] = extra_kwargs["k"]
68+
return repl_dict

‎sumpy/test/test_l2l_coeffs.py‎

Lines changed: 211 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,211 @@
1+
from __future__ import annotations
2+
3+
4+
__copyright__ = """
5+
Copyright (C) 2026 Shawn/Chaoqi Lin
6+
"""
7+
8+
__license__ = """
9+
Permission is hereby granted, free of charge, to any person obtaining a copy
10+
of this software and associated documentation files (the "Software"), to deal
11+
in the Software without restriction, including without limitation the rights
12+
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13+
copies of the Software, and to permit persons to whom the Software is
14+
furnished to do so, subject to the following conditions:
15+
16+
The above copyright notice and this permission notice shall be included in
17+
all copies or substantial portions of the Software.
18+
19+
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20+
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21+
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22+
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23+
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24+
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
25+
THE SOFTWARE.
26+
"""
27+
28+
29+
import math
30+
import sys
31+
from typing import TYPE_CHECKING
32+
33+
import numpy as np
34+
import pytest
35+
import scipy.special as spsp
36+
37+
from arraycontext import (
38+
pytest_generate_tests_for_array_contexts,
39+
)
40+
41+
import sumpy.toys as t
42+
from .coeff_test_tools import NumericMatVecOperator, get_repl_dict, to_scalar
43+
from sumpy.array_context import PytestPyOpenCLArrayContextFactory, _acf # noqa: F401
44+
from sumpy.expansion.local import (
45+
LinearPDEConformingVolumeTaylorLocalExpansion,
46+
VolumeTaylorLocalExpansion,
47+
)
48+
from sumpy.kernel import (
49+
BiharmonicKernel,
50+
HelmholtzKernel,
51+
Kernel,
52+
LaplaceKernel,
53+
YukawaKernel,
54+
)
55+
from sumpy.tools import build_matrix
56+
57+
58+
if TYPE_CHECKING:
59+
from collections.abc import Mapping
60+
61+
from arraycontext import ArrayContextFactory
62+
63+
64+
pytest_generate_tests = pytest_generate_tests_for_array_contexts([
65+
PytestPyOpenCLArrayContextFactory,
66+
])
67+
68+
69+
@pytest.mark.parametrize("knl,extra_kwargs", [
70+
(LaplaceKernel(2), {}),
71+
(YukawaKernel(2), {"lam": 0.1}),
72+
(HelmholtzKernel(2), {"k": 0.5}),
73+
(BiharmonicKernel(2), {}),
74+
])
75+
def test_l2l_coefficient_differences(
76+
actx_factory: ArrayContextFactory,
77+
knl: Kernel,
78+
extra_kwargs: Mapping[str, float],
79+
verbose: bool = True,
80+
):
81+
"""
82+
Tests that the expression for the difference between compressed and uncompressed
83+
translation in the compressed-expansions paper matches implemented reality.
84+
"""
85+
order = 7
86+
dim = 2
87+
repl_dict = get_repl_dict(knl, extra_kwargs)
88+
89+
# Setup sources and centers
90+
source = np.array([[5.0], [5.0]])
91+
c1 = np.array([0.0, 0.0])
92+
c2 = c1 + np.array([-0.5, 1.0])
93+
strength = np.array([1.0])
94+
95+
actx = actx_factory()
96+
toy_ctx = t.ToyContext(
97+
knl,
98+
local_expn_class=LinearPDEConformingVolumeTaylorLocalExpansion,
99+
extra_kernel_kwargs=extra_kwargs
100+
)
101+
toy_ctx_full = t.ToyContext(
102+
knl,
103+
local_expn_class=VolumeTaylorLocalExpansion,
104+
extra_kernel_kwargs=extra_kwargs
105+
)
106+
107+
# Compute expansions
108+
p = t.PointSources(toy_ctx, source, weights=strength)
109+
p_full = t.PointSources(toy_ctx_full, source, weights=strength)
110+
111+
p2l = t.local_expand(actx, p, c1, order=order, rscale=1.0)
112+
p2l2l = t.local_expand(actx, p2l, c2, order=order, rscale=1.0)
113+
p2l_full = t.local_expand(actx, p_full, c1, order=order, rscale=1.0)
114+
p2l2l_full = t.local_expand(actx, p2l_full, c2, order=order, rscale=1.0)
115+
116+
# Build matrix M
117+
p2l2l_expn = LinearPDEConformingVolumeTaylorLocalExpansion(knl, order)
118+
wrangler = p2l2l_expn.expansion_terms_wrangler
119+
M_symbolic = wrangler.get_projection_matrix(rscale=1.0) # noqa: N806
120+
numeric_op = NumericMatVecOperator(M_symbolic, repl_dict)
121+
M = build_matrix(numeric_op, dtype=np.complex128) # noqa: N806
122+
123+
# Get compressed coefficients
124+
mu_c_symbolic = wrangler.get_full_kernel_derivatives_from_stored(
125+
p2l2l.coeffs, rscale=1.0
126+
)
127+
mu_c = []
128+
for coeff in mu_c_symbolic:
129+
if hasattr(coeff, "xreplace"):
130+
mu_c.append(to_scalar(coeff.xreplace(repl_dict)))
131+
else:
132+
mu_c.append(to_scalar(coeff))
133+
134+
# Get identifiers
135+
stored_identifiers = p2l2l_expn.get_coefficient_identifiers()
136+
full_identifiers = p2l2l_expn.get_full_coefficient_identifiers()
137+
lexpn = VolumeTaylorLocalExpansion(knl, order)
138+
lexpn_idx = lexpn.get_full_coefficient_identifiers()
139+
140+
h = c2 - c1
141+
global_const = to_scalar(knl.get_global_scaling_const())
142+
143+
if verbose:
144+
print(f'\n{"="*104}')
145+
print(f"L2L Verification: {type(knl).__name__} (order={order})")
146+
print(f'{"="*104}')
147+
print(f"c1 = {c1}, c2 = {c2}, h = {h}")
148+
print()
149+
print(f"{'i':>3s} | {'ν(i)':>15s} | {'|ν|':4s} | " # noqa: RUF001
150+
f"{'formula':>31s} | {'direct':>31s} | {'abs err':>10s}")
151+
print("-" * 104)
152+
153+
max_abs_error = 0.0
154+
155+
for i, nu_i in enumerate(full_identifiers):
156+
i_card = sum(np.array(nu_i))
157+
158+
# Compute error by formula
159+
error = 0.0 + 0.0j
160+
for k, nu_jk in enumerate(stored_identifiers):
161+
jk_card = sum(np.array(nu_jk))
162+
if jk_card >= i_card:
163+
continue
164+
165+
start_idx = math.comb(order - i_card + dim, dim)
166+
end_idx = math.comb(order - jk_card + dim, dim)
167+
168+
for q_idx in range(start_idx, end_idx):
169+
nu_q = full_identifiers[q_idx]
170+
nu_sum = tuple(map(sum, zip(nu_q, nu_jk, strict=True)))
171+
if nu_sum not in full_identifiers:
172+
continue
173+
174+
deriv_idx = full_identifiers.index(nu_sum)
175+
gamma_deriv = to_scalar(p2l_full.coeffs[deriv_idx])
176+
h_pow = np.prod(h ** np.array(nu_q))
177+
fact_nu_q = np.prod(spsp.factorial(nu_q))
178+
179+
error += -M[i, k] * gamma_deriv * h_pow / fact_nu_q
180+
181+
error /= np.prod(spsp.factorial(nu_i))
182+
error *= global_const
183+
184+
# Compute direct difference
185+
true_i_idx = lexpn_idx.index(nu_i)
186+
mu_full = to_scalar(p2l2l_full.coeffs[true_i_idx])
187+
direct_diff = (mu_full - mu_c[i]) / np.prod(spsp.factorial(nu_i))
188+
direct_diff *= global_const
189+
190+
# Compute errors
191+
abs_err = abs(error - direct_diff)
192+
max_abs_error = max(max_abs_error, abs_err)
193+
194+
if verbose:
195+
print(f"{i:3d} | {nu_i!s:>15s} | {i_card:4d} | "
196+
f"{error.real: .8e}{error.imag:+.8e}j | "
197+
f"{direct_diff.real: .8e}{direct_diff.imag:+.8e}j | "
198+
f"{abs_err:9.2e}")
199+
200+
if verbose:
201+
print(f"\nMaximum absolute error: {max_abs_error:.2e}")
202+
203+
assert max_abs_error < 1e-10, \
204+
f"{type(knl).__name__}: error {max_abs_error:.2e}"
205+
206+
207+
if __name__ == "__main__":
208+
if len(sys.argv) > 1:
209+
exec(sys.argv[1])
210+
else:
211+
pytest.main([__file__])

0 commit comments

Comments
 (0)