Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

no_numpy in sumpy kernels #184

Merged
merged 4 commits into from
Aug 4, 2023
Merged

no_numpy in sumpy kernels #184

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
632b386
examples: only send cl arrays to kernels
alexfikl Aug 3, 2023
0de148a
toys: only send cl arrays to kernels
alexfikl Aug 3, 2023
8b41573
tests: only send cl arrays to kernels
alexfikl Aug 4, 2023
312771e
do not use out_host=True
alexfikl Aug 4, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
51 changes: 32 additions & 19 deletions examples/curve-pot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,8 @@ def draw_pot_figure(aspect_ratio,
knl_kwargs = {}

vol_source_knl, vol_target_knl = process_kernel(knl, what_operator)
p2p = P2P(actx.context, source_kernels=(vol_source_knl,),
p2p = P2P(actx.context,
source_kernels=(vol_source_knl,),
target_kernels=(vol_target_knl,),
exclude_self=False,
value_dtypes=np.complex128)
Expand Down Expand Up @@ -157,35 +158,40 @@ def map_to_curve(t):
lpot_kwargs = knl_kwargs.copy()

if what_operator == "D":
volpot_kwargs["src_derivative_dir"] = native_curve.normal
volpot_kwargs["src_derivative_dir"] = actx.from_numpy(native_curve.normal)

if what_operator_lpot == "D":
lpot_kwargs["src_derivative_dir"] = ovsmp_curve.normal
lpot_kwargs["src_derivative_dir"] = actx.from_numpy(ovsmp_curve.normal)

if what_operator_lpot == "S'":
lpot_kwargs["tgt_derivative_dir"] = native_curve.normal
lpot_kwargs["tgt_derivative_dir"] = actx.from_numpy(native_curve.normal)

# }}}

targets = actx.from_numpy(fp.points)
sources = actx.from_numpy(native_curve.pos)
ovsmp_sources = actx.from_numpy(ovsmp_curve.pos)

if 0:
# {{{ build matrix

from fourier import make_fourier_interp_matrix
fim = make_fourier_interp_matrix(novsmp, nsrc)

from sumpy.tools import build_matrix
from scipy.sparse.linalg import LinearOperator

def apply_lpot(x):
xovsmp = np.dot(fim, x)
evt, (y,) = lpot(actx.queue,
native_curve.pos,
ovsmp_curve.pos,
centers,
[xovsmp * ovsmp_curve.speed * ovsmp_weights],
expansion_radii=np.ones(centers.shape[1]),
sources,
ovsmp_sources,
actx.from_numpy(centers),
[actx.from_numpy(xovsmp * ovsmp_curve.speed * ovsmp_weights)],
expansion_radii=actx.from_numpy(np.ones(centers.shape[1])),
**lpot_kwargs)

return y
return actx.to_numpy(y)

op = LinearOperator((nsrc, nsrc), apply_lpot)
mat = build_matrix(op, dtype=np.complex128)
Expand All @@ -200,19 +206,26 @@ def apply_lpot(x):

mode_nr = 0
density = np.cos(mode_nr*2*np.pi*native_t).astype(np.complex128)
ovsmp_density = np.cos(mode_nr*2*np.pi*ovsmp_t).astype(np.complex128)
strength = actx.from_numpy(native_curve.speed * native_weights * density)

evt, (vol_pot,) = p2p(actx.queue,
fp.points,
native_curve.pos,
[native_curve.speed*native_weights*density], **volpot_kwargs)
targets,
sources,
[strength], **volpot_kwargs)
vol_pot = actx.to_numpy(vol_pot)

ovsmp_density = np.cos(mode_nr*2*np.pi*ovsmp_t).astype(np.complex128)
ovsmp_strength = actx.from_numpy(
ovsmp_curve.speed * ovsmp_weights * ovsmp_density)

evt, (curve_pot,) = lpot(actx.queue,
native_curve.pos,
ovsmp_curve.pos,
centers,
[ovsmp_density * ovsmp_curve.speed * ovsmp_weights],
expansion_radii=np.ones(centers.shape[1]),
sources,
ovsmp_sources,
actx.from_numpy(centers),
[ovsmp_strength],
expansion_radii=actx.from_numpy(np.ones(centers.shape[1])),
**lpot_kwargs)
curve_pot = actx.to_numpy(curve_pot)

# }}}

Expand Down
100 changes: 58 additions & 42 deletions sumpy/toys.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -201,48 +201,57 @@ def get_l2l(self, from_order, to_order):
# {{{ helpers

def _p2e(psource, center, rscale, order, p2e, expn_class, expn_kwargs):
source_boxes = np.array([0], dtype=np.int32)
box_source_starts = np.array([0], dtype=np.int32)
box_source_counts_nonchild = np.array(
[psource.points.shape[-1]], dtype=np.int32)

toy_ctx = psource.toy_ctx
queue = toy_ctx.queue

source_boxes = cl.array.to_device(
queue, np.array([0], dtype=np.int32))
box_source_starts = cl.array.to_device(
queue, np.array([0], dtype=np.int32))
box_source_counts_nonchild = cl.array.to_device(
queue, np.array([psource.points.shape[-1]], dtype=np.int32))

center = np.asarray(center)
centers = np.array(center, dtype=np.float64).reshape(
toy_ctx.kernel.dim, 1)
centers = cl.array.to_device(
queue,
np.array(center, dtype=np.float64).reshape(toy_ctx.kernel.dim, 1))

evt, (coeffs,) = p2e(toy_ctx.queue,
source_boxes=source_boxes,
box_source_starts=box_source_starts,
box_source_counts_nonchild=box_source_counts_nonchild,
centers=centers,
sources=psource.points,
strengths=(psource.weights,),
sources=cl.array.to_device(queue, psource.points),
strengths=(cl.array.to_device(queue, psource.weights),),
rscale=rscale,
nboxes=1,
tgt_base_ibox=0,

#flags="print_hl_cl",
out_host=True,
**toy_ctx.extra_source_and_kernel_kwargs)

return expn_class(toy_ctx, center, rscale, order, coeffs[0],
return expn_class(toy_ctx, center, rscale, order, coeffs[0].get(queue),
derived_from=psource, **expn_kwargs)


def _e2p(psource, targets, e2p):
ntargets = targets.shape[-1]
toy_ctx = psource.toy_ctx
queue = toy_ctx.queue

boxes = np.array([0], dtype=np.int32)
ntargets = targets.shape[-1]
boxes = cl.array.to_device(
queue, np.array([0], dtype=np.int32))
box_target_starts = cl.array.to_device(
queue, np.array([0], dtype=np.int32))
box_target_counts_nonchild = cl.array.to_device(
queue, np.array([ntargets], dtype=np.int32))

box_target_starts = np.array([0], dtype=np.int32)
box_target_counts_nonchild = np.array([ntargets], dtype=np.int32)
centers = cl.array.to_device(
queue,
np.array(psource.center, dtype=np.float64).reshape(toy_ctx.kernel.dim, 1))

toy_ctx = psource.toy_ctx
centers = np.array(psource.center, dtype=np.float64).reshape(
toy_ctx.kernel.dim, 1)
from pytools.obj_array import make_obj_array
from sumpy.tools import vector_to_device

coeffs = np.array([psource.coeffs])
coeffs = cl.array.to_device(queue, np.array([psource.coeffs]))
evt, (pot,) = e2p(
toy_ctx.queue,
src_expansions=coeffs,
Expand All @@ -252,31 +261,36 @@ def _e2p(psource, targets, e2p):
box_target_counts_nonchild=box_target_counts_nonchild,
centers=centers,
rscale=psource.rscale,
targets=targets,
#flags="print_hl_cl",
out_host=True, **toy_ctx.extra_kernel_kwargs)
targets=vector_to_device(queue, make_obj_array(targets)),
**toy_ctx.extra_kernel_kwargs)

return pot
return pot.get(queue)


def _e2e(psource, to_center, to_rscale, to_order, e2e, expn_class, expn_kwargs):
toy_ctx = psource.toy_ctx

target_boxes = np.array([1], dtype=np.int32)
src_box_starts = np.array([0, 1], dtype=np.int32)
src_box_lists = np.array([0], dtype=np.int32)

centers = (np.array(
queue = toy_ctx.queue

target_boxes = cl.array.to_device(
queue, np.array([1], dtype=np.int32))
src_box_starts = cl.array.to_device(
queue, np.array([0, 1], dtype=np.int32))
src_box_lists = cl.array.to_device(
queue, np.array([0], dtype=np.int32))

centers = cl.array.to_device(
queue,
np.array(
[
# box 0: source
psource.center,

# box 1: target
to_center,
],
dtype=np.float64)).T.copy()

coeffs = np.array([psource.coeffs])
],
dtype=np.float64).T.copy()
)
coeffs = cl.array.to_device(queue, np.array([psource.coeffs]))

evt, (to_coeffs,) = e2e(
toy_ctx.queue,
Expand All @@ -293,11 +307,10 @@ def _e2e(psource, to_center, to_rscale, to_order, e2e, expn_class, expn_kwargs):

src_rscale=psource.rscale,
tgt_rscale=to_rscale,
**toy_ctx.extra_kernel_kwargs)

#flags="print_hl_cl",
out_host=True, **toy_ctx.extra_kernel_kwargs)

return expn_class(toy_ctx, to_center, to_rscale, to_order, to_coeffs[1],
return expn_class(
toy_ctx, to_center, to_rscale, to_order, to_coeffs[1].get(queue),
derived_from=psource, **expn_kwargs)

# }}}
Expand Down Expand Up @@ -443,12 +456,15 @@ def __init__(self,
self._center = center

def eval(self, targets: np.ndarray) -> np.ndarray:
queue = self.toy_ctx.queue
evt, (potential,) = self.toy_ctx.get_p2p()(
self.toy_ctx.queue, targets, self.points, [self.weights],
out_host=True,
queue,
cl.array.to_device(queue, targets),
cl.array.to_device(queue, self.points),
[cl.array.to_device(queue, self.weights)],
**self.toy_ctx.extra_source_and_kernel_kwargs)

return potential
return potential.get(queue)

@property
def center(self):
Expand Down
21 changes: 14 additions & 7 deletions test/test_fmm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,9 +62,13 @@

# {{{ test_sumpy_fmm

@pytest.mark.parametrize("use_translation_classes, use_fft, fft_backend",
[(False, False, None), (True, False, None), (True, True, "loopy"),
(True, True, "pyvkfft")])
@pytest.mark.parametrize(
("use_translation_classes", "use_fft", "fft_backend"), [
(False, False, None),
(True, False, None),
(True, True, "loopy"),
(True, True, "pyvkfft"),
])
@pytest.mark.parametrize(
("knl", "local_expn_class", "mpole_expn_class",
"order_varies_with_level"), [
Expand Down Expand Up @@ -92,6 +96,9 @@
def test_sumpy_fmm(actx_factory, knl, local_expn_class, mpole_expn_class,
order_varies_with_level, use_translation_classes, use_fft,
fft_backend, visualize=False):
if fft_backend == "pyvkfft":
pytest.importorskip("pyvkfft")

if visualize:
logging.basicConfig(level=logging.INFO)

Expand Down Expand Up @@ -399,7 +406,7 @@ def test_unified_single_and_double(actx_factory, visualize=False):
strength_usages = [[0], [1], [0, 1]]

alpha = np.linspace(0, 2*np.pi, nsources, np.float64)
dir_vec = np.vstack([np.cos(alpha), np.sin(alpha)])
dir_vec = actx.from_numpy(np.vstack([np.cos(alpha), np.sin(alpha)]))

results = []
for source_kernels, strength_usage in zip(source_kernel_vecs, strength_usages):
Expand Down Expand Up @@ -528,7 +535,7 @@ def test_sumpy_fmm_exclude_self(actx_factory, visualize=False):
rng = np.random.default_rng(44)
weights = actx.from_numpy(rng.random(nsources, dtype=np.float64))

target_to_source = np.arange(tree.ntargets, dtype=np.int32)
target_to_source = actx.from_numpy(np.arange(tree.ntargets, dtype=np.int32))
self_extra_kwargs = {"target_to_source": target_to_source}

target_kernels = [knl]
Expand Down Expand Up @@ -597,7 +604,7 @@ def test_sumpy_axis_source_derivative(actx_factory, visualize=False):
rng = np.random.default_rng(12)
weights = actx.from_numpy(rng.random(nsources, dtype=np.float64))

target_to_source = np.arange(tree.ntargets, dtype=np.int32)
target_to_source = actx.from_numpy(np.arange(tree.ntargets, dtype=np.int32))
self_extra_kwargs = {"target_to_source": target_to_source}

from sumpy.kernel import AxisTargetDerivative, AxisSourceDerivative
Expand Down Expand Up @@ -665,7 +672,7 @@ def test_sumpy_target_point_multiplier(actx_factory, deriv_axes, visualize=False
rng = np.random.default_rng(12)
weights = actx.from_numpy(rng.random(nsources, dtype=np.float64))

target_to_source = np.arange(tree.ntargets, dtype=np.int32)
target_to_source = actx.from_numpy(np.arange(tree.ntargets, dtype=np.int32))
self_extra_kwargs = {"target_to_source": target_to_source}

from sumpy.kernel import TargetPointMultiplier, AxisTargetDerivative
Expand Down
Loading