-
Isuru Fernando authoredIsuru Fernando authored
test_kernels.py 28.17 KiB
__copyright__ = "Copyright (C) 2012 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 numpy as np
import numpy.linalg as la
import sys
import pytest
import pyopencl as cl
from pyopencl.tools import ( # noqa
pytest_generate_tests_for_pyopencl as pytest_generate_tests)
from sumpy.expansion.multipole import (
VolumeTaylorMultipoleExpansion, H2DMultipoleExpansion,
VolumeTaylorMultipoleExpansionBase,
LinearPDEConformingVolumeTaylorMultipoleExpansion)
from sumpy.expansion.local import (
VolumeTaylorLocalExpansion, H2DLocalExpansion,
LinearPDEConformingVolumeTaylorLocalExpansion)
from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, AxisTargetDerivative,
DirectionalSourceDerivative, BiharmonicKernel, StokesletKernel)
import sumpy.symbolic as sym
from pytools.convergence import PConvergenceVerifier
import logging
logger = logging.getLogger(__name__)
try:
import faulthandler
except ImportError:
pass
else:
faulthandler.enable()
@pytest.mark.parametrize("exclude_self", (True, False))
def test_p2p(ctx_factory, exclude_self):
ctx = ctx_factory()
queue = cl.CommandQueue(ctx)
dimensions = 3
n = 5000
from sumpy.p2p import P2P
lknl = LaplaceKernel(dimensions)
knl = P2P(ctx,
[lknl, AxisTargetDerivative(0, lknl)],
exclude_self=exclude_self)
targets = np.random.rand(dimensions, n)
sources = targets if exclude_self else np.random.rand(dimensions, n)
strengths = np.ones(n, dtype=np.float64)
extra_kwargs = {}
if exclude_self:
extra_kwargs["target_to_source"] = np.arange(n, dtype=np.int32)
evt, (potential, x_derivative) = knl(
queue, targets, sources, [strengths],
out_host=True, **extra_kwargs)
potential_ref = np.empty_like(potential)
targets = targets.T
sources = sources.T
for itarg in range(n):
with np.errstate(divide="ignore"):
invdists = np.sum((targets[itarg] - sources) ** 2, axis=-1) ** -0.5
if exclude_self:
assert np.isinf(invdists[itarg])
invdists[itarg] = 0
potential_ref[itarg] = np.sum(strengths * invdists)
potential_ref *= 1/(4*np.pi)
rel_err = la.norm(potential - potential_ref)/la.norm(potential_ref)
print(rel_err)
assert rel_err < 1e-3
@pytest.mark.parametrize(("base_knl", "expn_class"), [
(LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion),
(LaplaceKernel(2), LinearPDEConformingVolumeTaylorMultipoleExpansion),
])
def test_p2e_multiple(ctx_factory, base_knl, expn_class):
from sympy.core.cache import clear_cache
clear_cache()
order = 4
ctx = ctx_factory()
queue = cl.CommandQueue(ctx)
np.random.seed(17)
nsources = 100
extra_kwargs = {}
if isinstance(base_knl, HelmholtzKernel):
if base_knl.allow_evanescent:
extra_kwargs["k"] = 0.2 * (0.707 + 0.707j)
else:
extra_kwargs["k"] = 0.2
if isinstance(base_knl, StokesletKernel):
extra_kwargs["mu"] = 0.2
source_kernels = [
DirectionalSourceDerivative(base_knl, "dir_vec"),
base_knl,
]
knl = base_knl
expn = expn_class(knl, order=order)
from sumpy import P2EFromSingleBox
center = np.array([2, 1, 0][:knl.dim], np.float64)
sources = (0.7*(-0.5+np.random.rand(knl.dim, nsources).astype(np.float64))
+ center[:, np.newaxis])
strengths = [
np.ones(nsources, dtype=np.float64) * (1/nsources),
np.ones(nsources, dtype=np.float64) * (2/nsources)
]
source_boxes = np.array([0], dtype=np.int32)
box_source_starts = np.array([0], dtype=np.int32)
box_source_counts_nonchild = np.array([nsources], dtype=np.int32)
alpha = np.linspace(0, 2*np.pi, nsources, np.float64)
dir_vec = np.vstack([np.cos(alpha), np.sin(alpha)])
from sumpy.expansion.local import LocalExpansionBase
if issubclass(expn_class, LocalExpansionBase):
loc_center = np.array([5.5, 0.0, 0.0][:knl.dim]) + center
centers = np.array(loc_center, dtype=np.float64).reshape(knl.dim, 1)
else:
centers = (np.array([0.0, 0.0, 0.0][:knl.dim],
dtype=np.float64).reshape(knl.dim, 1)
+ center[:, np.newaxis])
rscale = 0.5 # pick something non-1
# apply p2e at the same time
p2e = P2EFromSingleBox(ctx, expn, kernels=source_kernels, strength_usage=[0, 1])
evt, (mpoles,) = p2e(queue,
source_boxes=source_boxes,
box_source_starts=box_source_starts,
box_source_counts_nonchild=box_source_counts_nonchild,
centers=centers,
sources=sources,
strengths=strengths,
nboxes=1,
tgt_base_ibox=0,
rscale=rscale,
#flags="print_hl_cl",
out_host=True,
dir_vec=dir_vec,
**extra_kwargs)
actual_result = mpoles
# apply p2e separately
expected_result = np.zeros_like(actual_result)
for i, source_kernel in enumerate(source_kernels):
extra_source_kwargs = extra_kwargs.copy()
if isinstance(source_kernel, DirectionalSourceDerivative):
extra_source_kwargs["dir_vec"] = dir_vec
p2e = P2EFromSingleBox(ctx, expn,
kernels=[source_kernel], strength_usage=[i])
evt, (mpoles,) = p2e(queue,
source_boxes=source_boxes,
box_source_starts=box_source_starts,
box_source_counts_nonchild=box_source_counts_nonchild,
centers=centers,
sources=sources,
strengths=strengths,
nboxes=1,
tgt_base_ibox=0,
rscale=rscale,
#flags="print_hl_cl",
out_host=True, **extra_source_kwargs)
expected_result += mpoles
norm = la.norm(actual_result - expected_result)/la.norm(expected_result)
assert norm < 1e-12
@pytest.mark.parametrize("order", [4])
@pytest.mark.parametrize(("base_knl", "expn_class"), [
(LaplaceKernel(2), VolumeTaylorLocalExpansion),
(LaplaceKernel(2), VolumeTaylorMultipoleExpansion),
(LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion),
(LaplaceKernel(2), LinearPDEConformingVolumeTaylorMultipoleExpansion),
(HelmholtzKernel(2), VolumeTaylorMultipoleExpansion),
(HelmholtzKernel(2), VolumeTaylorLocalExpansion),
(HelmholtzKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion),
(HelmholtzKernel(2), LinearPDEConformingVolumeTaylorMultipoleExpansion),
(HelmholtzKernel(2), H2DLocalExpansion),
(HelmholtzKernel(2), H2DMultipoleExpansion),
(DirectionalSourceDerivative(BiharmonicKernel(2), "dir_vec"),
VolumeTaylorMultipoleExpansion),
(DirectionalSourceDerivative(BiharmonicKernel(2), "dir_vec"),
VolumeTaylorLocalExpansion),
(HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorMultipoleExpansion),
(HelmholtzKernel(2, allow_evanescent=True), VolumeTaylorLocalExpansion),
(HelmholtzKernel(2, allow_evanescent=True),
LinearPDEConformingVolumeTaylorLocalExpansion),
(HelmholtzKernel(2, allow_evanescent=True),
LinearPDEConformingVolumeTaylorMultipoleExpansion),
(HelmholtzKernel(2, allow_evanescent=True), H2DLocalExpansion),
(HelmholtzKernel(2, allow_evanescent=True), H2DMultipoleExpansion),
])
@pytest.mark.parametrize("with_source_derivative", [
False,
True
])
# Sample: test_p2e2p(cl._csc, LaplaceKernel(2), VolumeTaylorLocalExpansion, 4, False)
def test_p2e2p(ctx_factory, base_knl, expn_class, order, with_source_derivative):
#logging.basicConfig(level=logging.INFO)
from sympy.core.cache import clear_cache
clear_cache()
ctx = ctx_factory()
queue = cl.CommandQueue(ctx)
np.random.seed(17)
res = 100
nsources = 100
extra_kwargs = {}
if isinstance(base_knl, HelmholtzKernel):
if base_knl.allow_evanescent:
extra_kwargs["k"] = 0.2 * (0.707 + 0.707j)
else:
extra_kwargs["k"] = 0.2
if isinstance(base_knl, StokesletKernel):
extra_kwargs["mu"] = 0.2
if with_source_derivative:
knl = DirectionalSourceDerivative(base_knl, "dir_vec")
else:
knl = base_knl
target_kernels = [
knl,
AxisTargetDerivative(0, knl),
]
expn = expn_class(knl, order=order)
from sumpy import P2EFromSingleBox, E2PFromSingleBox, P2P
p2e = P2EFromSingleBox(ctx, expn, kernels=[knl])
e2p = E2PFromSingleBox(ctx, expn, kernels=target_kernels)
p2p = P2P(ctx, target_kernels, exclude_self=False)
from pytools.convergence import EOCRecorder
eoc_rec_pot = EOCRecorder()
eoc_rec_grad_x = EOCRecorder()
from sumpy.expansion.local import LocalExpansionBase
if issubclass(expn_class, LocalExpansionBase):
h_values = [1/5, 1/7, 1/20]
else:
h_values = [1/2, 1/3, 1/5]
center = np.array([2, 1, 0][:knl.dim], np.float64)
sources = (0.7*(-0.5+np.random.rand(knl.dim, nsources).astype(np.float64))
+ center[:, np.newaxis])
strengths = np.ones(nsources, dtype=np.float64) * (1/nsources)
source_boxes = np.array([0], dtype=np.int32)
box_source_starts = np.array([0], dtype=np.int32)
box_source_counts_nonchild = np.array([nsources], dtype=np.int32)
extra_source_kwargs = extra_kwargs.copy()
if isinstance(knl, DirectionalSourceDerivative):
alpha = np.linspace(0, 2*np.pi, nsources, np.float64)
dir_vec = np.vstack([np.cos(alpha), np.sin(alpha)])
extra_source_kwargs["dir_vec"] = dir_vec
from sumpy.visualization import FieldPlotter
for h in h_values:
if issubclass(expn_class, LocalExpansionBase):
loc_center = np.array([5.5, 0.0, 0.0][:knl.dim]) + center
centers = np.array(loc_center, dtype=np.float64).reshape(knl.dim, 1)
fp = FieldPlotter(loc_center, extent=h, npoints=res)
else:
eval_center = np.array([1/h, 0.0, 0.0][:knl.dim]) + center
fp = FieldPlotter(eval_center, extent=0.1, npoints=res)
centers = (np.array([0.0, 0.0, 0.0][:knl.dim],
dtype=np.float64).reshape(knl.dim, 1)
+ center[:, np.newaxis])
targets = fp.points
rscale = 0.5 # pick something non-1
# {{{ apply p2e
evt, (mpoles,) = p2e(queue,
source_boxes=source_boxes,
box_source_starts=box_source_starts,
box_source_counts_nonchild=box_source_counts_nonchild,
centers=centers,
sources=sources,
strengths=(strengths,),
nboxes=1,
tgt_base_ibox=0,
rscale=rscale,
#flags="print_hl_cl",
out_host=True, **extra_source_kwargs)
# }}}
# {{{ apply e2p
ntargets = targets.shape[-1]
box_target_starts = np.array([0], dtype=np.int32)
box_target_counts_nonchild = np.array([ntargets], dtype=np.int32)
evt, (pot, grad_x, ) = e2p(
queue,
src_expansions=mpoles,
src_base_ibox=0,
target_boxes=source_boxes,
box_target_starts=box_target_starts,
box_target_counts_nonchild=box_target_counts_nonchild,
centers=centers,
targets=targets,
rscale=rscale,
#flags="print_hl_cl",
out_host=True, **extra_kwargs)
# }}}
# {{{ compute (direct) reference solution
evt, (pot_direct, grad_x_direct, ) = p2p(
queue,
targets, sources, (strengths,),
out_host=True,
**extra_source_kwargs)
err_pot = la.norm((pot - pot_direct)/res**2)
err_grad_x = la.norm((grad_x - grad_x_direct)/res**2)
if 1:
err_pot = err_pot / la.norm((pot_direct)/res**2)
err_grad_x = err_grad_x / la.norm((grad_x_direct)/res**2)
if 0:
import matplotlib.pyplot as pt
from matplotlib.colors import Normalize
pt.subplot(131)
im = fp.show_scalar_in_matplotlib(pot.real)
im.set_norm(Normalize(vmin=-0.1, vmax=0.1))
pt.subplot(132)
im = fp.show_scalar_in_matplotlib(pot_direct.real)
im.set_norm(Normalize(vmin=-0.1, vmax=0.1))
pt.colorbar()
pt.subplot(133)
im = fp.show_scalar_in_matplotlib(np.log10(1e-15+np.abs(pot-pot_direct)))
im.set_norm(Normalize(vmin=-6, vmax=1))
pt.colorbar()
pt.show()
# }}}
eoc_rec_pot.add_data_point(h, err_pot)
eoc_rec_grad_x.add_data_point(h, err_grad_x)
print(expn_class, knl, order)
print("POTENTIAL:")
print(eoc_rec_pot)
print("X TARGET DERIVATIVE:")
print(eoc_rec_grad_x)
tgt_order = order + 1
if issubclass(expn_class, LocalExpansionBase):
tgt_order_grad = tgt_order - 1
slack = 0.7
grad_slack = 0.5
else:
tgt_order_grad = tgt_order + 1
slack = 0.5
grad_slack = 1
if order <= 2:
slack += 1
grad_slack += 1
if isinstance(knl, DirectionalSourceDerivative):
slack += 1
grad_slack += 2
if isinstance(base_knl, DirectionalSourceDerivative):
slack += 1
grad_slack += 2
if isinstance(base_knl, HelmholtzKernel):
if base_knl.allow_evanescent:
slack += 0.5
grad_slack += 0.5
if issubclass(expn_class, VolumeTaylorMultipoleExpansionBase):
slack += 0.3
grad_slack += 0.3
assert eoc_rec_pot.order_estimate() > tgt_order - slack
assert eoc_rec_grad_x.order_estimate() > tgt_order_grad - grad_slack
@pytest.mark.parametrize("knl, local_expn_class, mpole_expn_class", [
(LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion),
(LaplaceKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion,
LinearPDEConformingVolumeTaylorMultipoleExpansion),
(HelmholtzKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion),
(HelmholtzKernel(2), LinearPDEConformingVolumeTaylorLocalExpansion,
LinearPDEConformingVolumeTaylorMultipoleExpansion),
(HelmholtzKernel(2), H2DLocalExpansion, H2DMultipoleExpansion),
(StokesletKernel(2, 0, 0), VolumeTaylorLocalExpansion,
VolumeTaylorMultipoleExpansion),
(StokesletKernel(2, 0, 0), LinearPDEConformingVolumeTaylorLocalExpansion,
LinearPDEConformingVolumeTaylorMultipoleExpansion),
])
def test_translations(ctx_factory, knl, local_expn_class, mpole_expn_class):
logging.basicConfig(level=logging.INFO)
from sympy.core.cache import clear_cache
clear_cache()
ctx = ctx_factory()
queue = cl.CommandQueue(ctx)
np.random.seed(17)
res = 20
nsources = 15
target_kernels = [knl]
extra_kwargs = {}
if isinstance(knl, HelmholtzKernel):
extra_kwargs["k"] = 0.05
if isinstance(knl, StokesletKernel):
extra_kwargs["mu"] = 0.05
# Just to make sure things also work away from the origin
origin = np.array([2, 1, 0][:knl.dim], np.float64)
sources = (0.7*(-0.5+np.random.rand(knl.dim, nsources).astype(np.float64))
+ origin[:, np.newaxis])
strengths = np.ones(nsources, dtype=np.float64) * (1/nsources)
pconv_verifier_p2m2p = PConvergenceVerifier()
pconv_verifier_p2m2m2p = PConvergenceVerifier()
pconv_verifier_p2m2m2l2p = PConvergenceVerifier()
pconv_verifier_full = PConvergenceVerifier()
from sumpy.visualization import FieldPlotter
eval_offset = np.array([5.5, 0.0, 0][:knl.dim])
centers = (np.array(
[
# box 0: particles, first mpole here
[0, 0, 0][:knl.dim],
# box 1: second mpole here
np.array([-0.2, 0.1, 0][:knl.dim], np.float64),
# box 2: first local here
eval_offset + np.array([0.3, -0.2, 0][:knl.dim], np.float64),
# box 3: second local and eval here
eval_offset
],
dtype=np.float64) + origin).T.copy()
del eval_offset
from sumpy.expansion import VolumeTaylorExpansionBase
if isinstance(knl, HelmholtzKernel) and \
issubclass(local_expn_class, VolumeTaylorExpansionBase):
# FIXME: Embarrassing--but we run out of memory for higher orders.
orders = [2, 3]
else:
orders = [2, 3, 4]
nboxes = centers.shape[-1]
def eval_at(e2p, source_box_nr, rscale):
e2p_target_boxes = np.array([source_box_nr], dtype=np.int32)
# These are indexed by global box numbers.
e2p_box_target_starts = np.array([0, 0, 0, 0], dtype=np.int32)
e2p_box_target_counts_nonchild = np.array([0, 0, 0, 0],
dtype=np.int32)
e2p_box_target_counts_nonchild[source_box_nr] = ntargets
evt, (pot,) = e2p(
queue,
src_expansions=mpoles,
src_base_ibox=0,
target_boxes=e2p_target_boxes,
box_target_starts=e2p_box_target_starts,
box_target_counts_nonchild=e2p_box_target_counts_nonchild,
centers=centers,
targets=targets,
rscale=rscale,
out_host=True, **extra_kwargs
)
return pot
for order in orders:
m_expn = mpole_expn_class(knl, order=order)
l_expn = local_expn_class(knl, order=order)
from sumpy import P2EFromSingleBox, E2PFromSingleBox, P2P, E2EFromCSR
p2m = P2EFromSingleBox(ctx, m_expn)
m2m = E2EFromCSR(ctx, m_expn, m_expn)
m2p = E2PFromSingleBox(ctx, m_expn, target_kernels)
m2l = E2EFromCSR(ctx, m_expn, l_expn)
l2l = E2EFromCSR(ctx, l_expn, l_expn)
l2p = E2PFromSingleBox(ctx, l_expn, target_kernels)
p2p = P2P(ctx, target_kernels, exclude_self=False)
fp = FieldPlotter(centers[:, -1], extent=0.3, npoints=res)
targets = fp.points
# {{{ compute (direct) reference solution
evt, (pot_direct,) = p2p(
queue,
targets, sources, (strengths,),
out_host=True, **extra_kwargs)
# }}}
m1_rscale = 0.5
m2_rscale = 0.25
l1_rscale = 0.5
l2_rscale = 0.25
# {{{ apply P2M
p2m_source_boxes = np.array([0], dtype=np.int32)
# These are indexed by global box numbers.
p2m_box_source_starts = np.array([0, 0, 0, 0], dtype=np.int32)
p2m_box_source_counts_nonchild = np.array([nsources, 0, 0, 0],
dtype=np.int32)
evt, (mpoles,) = p2m(queue,
source_boxes=p2m_source_boxes,
box_source_starts=p2m_box_source_starts,
box_source_counts_nonchild=p2m_box_source_counts_nonchild,
centers=centers,
sources=sources,
strengths=(strengths,),
nboxes=nboxes,
rscale=m1_rscale,
tgt_base_ibox=0,
#flags="print_hl_wrapper",
out_host=True, **extra_kwargs)
# }}}
ntargets = targets.shape[-1]
pot = eval_at(m2p, 0, m1_rscale)
err = la.norm((pot - pot_direct)/res**2)
err = err / (la.norm(pot_direct) / res**2)
pconv_verifier_p2m2p.add_data_point(order, err)
# {{{ apply M2M
m2m_target_boxes = np.array([1], dtype=np.int32)
m2m_src_box_starts = np.array([0, 1], dtype=np.int32)
m2m_src_box_lists = np.array([0], dtype=np.int32)
evt, (mpoles,) = m2m(queue,
src_expansions=mpoles,
src_base_ibox=0,
tgt_base_ibox=0,
ntgt_level_boxes=mpoles.shape[0],
target_boxes=m2m_target_boxes,
src_box_starts=m2m_src_box_starts,
src_box_lists=m2m_src_box_lists,
centers=centers,
src_rscale=m1_rscale,
tgt_rscale=m2_rscale,
#flags="print_hl_cl",
out_host=True, **extra_kwargs)
# }}}
pot = eval_at(m2p, 1, m2_rscale)
err = la.norm((pot - pot_direct)/res**2)
err = err / (la.norm(pot_direct) / res**2)
pconv_verifier_p2m2m2p.add_data_point(order, err)
# {{{ apply M2L
m2l_target_boxes = np.array([2], dtype=np.int32)
m2l_src_box_starts = np.array([0, 1], dtype=np.int32)
m2l_src_box_lists = np.array([1], dtype=np.int32)
evt, (mpoles,) = m2l(queue,
src_expansions=mpoles,
src_base_ibox=0,
tgt_base_ibox=0,
ntgt_level_boxes=mpoles.shape[0],
target_boxes=m2l_target_boxes,
src_box_starts=m2l_src_box_starts,
src_box_lists=m2l_src_box_lists,
centers=centers,
src_rscale=m2_rscale,
tgt_rscale=l1_rscale,
#flags="print_hl_cl",
out_host=True, **extra_kwargs)
# }}}
pot = eval_at(l2p, 2, l1_rscale)
err = la.norm((pot - pot_direct)/res**2)
err = err / (la.norm(pot_direct) / res**2)
pconv_verifier_p2m2m2l2p.add_data_point(order, err)
# {{{ apply L2L
l2l_target_boxes = np.array([3], dtype=np.int32)
l2l_src_box_starts = np.array([0, 1], dtype=np.int32)
l2l_src_box_lists = np.array([2], dtype=np.int32)
evt, (mpoles,) = l2l(queue,
src_expansions=mpoles,
src_base_ibox=0,
tgt_base_ibox=0,
ntgt_level_boxes=mpoles.shape[0],
target_boxes=l2l_target_boxes,
src_box_starts=l2l_src_box_starts,
src_box_lists=l2l_src_box_lists,
centers=centers,
src_rscale=l1_rscale,
tgt_rscale=l2_rscale,
#flags="print_hl_wrapper",
out_host=True, **extra_kwargs)
# }}}
pot = eval_at(l2p, 3, l2_rscale)
err = la.norm((pot - pot_direct)/res**2)
err = err / (la.norm(pot_direct) / res**2)
pconv_verifier_full.add_data_point(order, err)
for name, verifier in [
("p2m2p", pconv_verifier_p2m2p),
("p2m2m2p", pconv_verifier_p2m2m2p),
("p2m2m2l2p", pconv_verifier_p2m2m2l2p),
("full", pconv_verifier_full),
]:
print(30*"-")
print(name)
print(30*"-")
print(verifier)
print(30*"-")
verifier()
@pytest.mark.parametrize("order", [4])
@pytest.mark.parametrize(("base_knl", "local_expn_class", "mpole_expn_class"), [
(LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion),
])
@pytest.mark.parametrize("with_source_derivative", [
False,
True
])
def test_m2m_and_l2l_exprs_simpler(base_knl, local_expn_class, mpole_expn_class,
order, with_source_derivative):
from sympy.core.cache import clear_cache
clear_cache()
np.random.seed(17)
extra_kwargs = {}
if isinstance(base_knl, HelmholtzKernel):
if base_knl.allow_evanescent:
extra_kwargs["k"] = 0.2 * (0.707 + 0.707j)
else:
extra_kwargs["k"] = 0.2
if isinstance(base_knl, StokesletKernel):
extra_kwargs["mu"] = 0.2
if with_source_derivative:
knl = DirectionalSourceDerivative(base_knl, "dir_vec")
else:
knl = base_knl
mpole_expn = mpole_expn_class(knl, order=order)
local_expn = local_expn_class(knl, order=order)
from sumpy.symbolic import make_sym_vector, Symbol, USE_SYMENGINE
dvec = make_sym_vector("d", knl.dim)
src_coeff_exprs = [Symbol("src_coeff%d" % i) for i in range(len(mpole_expn))]
src_rscale = 3
tgt_rscale = 2
faster_m2m = mpole_expn.translate_from(mpole_expn, src_coeff_exprs, src_rscale,
dvec, tgt_rscale)
slower_m2m = mpole_expn.translate_from(mpole_expn, src_coeff_exprs, src_rscale,
dvec, tgt_rscale, _fast_version=False)
def _check_equal(expr1, expr2):
if USE_SYMENGINE:
return float((expr1 - expr2).expand()) == 0.0
else:
# with sympy we are using UnevaluatedExpr and expand doesn't expand it
# Running doit replaces UnevaluatedExpr with evaluated exprs
return float((expr1 - expr2).doit().expand()) == 0.0
for expr1, expr2 in zip(faster_m2m, slower_m2m):
assert _check_equal(expr1, expr2)
faster_l2l = local_expn.translate_from(local_expn, src_coeff_exprs, src_rscale,
dvec, tgt_rscale)
slower_l2l = local_expn.translate_from(local_expn, src_coeff_exprs, src_rscale,
dvec, tgt_rscale, _fast_version=False)
for expr1, expr2 in zip(faster_l2l, slower_l2l):
assert _check_equal(expr1, expr2)
# {{{ test toeplitz
def _m2l_translate_simple(tgt_expansion, src_expansion, src_coeff_exprs, src_rscale,
dvec, tgt_rscale):
if not tgt_expansion.use_rscale:
src_rscale = 1
tgt_rscale = 1
from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansionBase
if not isinstance(src_expansion, VolumeTaylorMultipoleExpansionBase):
return 1
# We know the general form of the multipole expansion is:
#
# coeff0 * diff(kernel, mi0) + coeff1 * diff(kernel, mi1) + ...
#
# To get the local expansion coefficients, we take derivatives of
# the multipole expansion.
taker = src_expansion.kernel.get_derivative_taker(dvec, src_rscale, sac=None)
from sumpy.tools import add_mi
result = []
for deriv in tgt_expansion.get_coefficient_identifiers():
local_result = []
for coeff, term in zip(
src_coeff_exprs,
src_expansion.get_coefficient_identifiers()):
kernel_deriv = taker.diff(add_mi(deriv, term)) / src_rscale**sum(deriv)
local_result.append(
coeff * kernel_deriv * tgt_rscale**sum(deriv))
result.append(sym.Add(*local_result))
return result
def test_m2l_toeplitz():
dim = 3
knl = LaplaceKernel(dim)
local_expn_class = LinearPDEConformingVolumeTaylorLocalExpansion
mpole_expn_class = LinearPDEConformingVolumeTaylorMultipoleExpansion
local_expn = local_expn_class(knl, order=5)
mpole_expn = mpole_expn_class(knl, order=5)
dvec = sym.make_sym_vector("d", dim)
src_coeff_exprs = list(1 + np.random.randn(len(mpole_expn)))
src_rscale = 2.0
tgt_rscale = 1.0
expected_output = _m2l_translate_simple(local_expn, mpole_expn, src_coeff_exprs,
src_rscale, dvec, tgt_rscale)
actual_output = local_expn.translate_from(mpole_expn, src_coeff_exprs,
src_rscale, dvec, tgt_rscale, sac=None)
replace_dict = dict((d, np.random.rand(1)[0]) for d in dvec)
for sym_a, sym_b in zip(expected_output, actual_output):
num_a = sym_a.xreplace(replace_dict)
num_b = sym_b.xreplace(replace_dict)
assert abs(num_a - num_b)/abs(num_a) < 1e-10
# }}}
# You can test individual routines by typing
# $ python test_kernels.py 'test_p2p(cl.create_some_context)'
if __name__ == "__main__":
if len(sys.argv) > 1:
exec(sys.argv[1])
else:
from pytest import main
main([__file__])
# vim: fdm=marker