Welcome to sumpy’s documentation!¶
Sumpy is mainly a ‘scaffolding’ package for Fast Multipole and quadrature methods. If you’re building one of those and need code generation for the required Multipole and local expansions, come right on in. Together with boxtree, there is a full, symbolically kernel-independent FMM implementation here.
Contents¶
Indices and tables¶
Example¶
import numpy as np
import numpy.linalg as la
import pyopencl as cl
try:
import matplotlib.pyplot as plt
USE_MATPLOTLIB = True
except ImportError:
USE_MATPLOTLIB = False
try:
from mayavi import mlab
USE_MAYAVI = True
except ImportError:
USE_MAYAVI = False
import logging
logging.basicConfig(level=logging.INFO)
def process_kernel(knl, what_operator):
target_knl = knl
source_knl = knl
if what_operator == "S":
pass
elif what_operator == "S0":
from sumpy.kernel import AxisTargetDerivative
target_knl = AxisTargetDerivative(0, knl)
elif what_operator == "S1":
from sumpy.kernel import AxisTargetDerivative
target_knl = AxisTargetDerivative(1, knl)
elif what_operator == "D":
from sumpy.kernel import DirectionalSourceDerivative
source_knl = DirectionalSourceDerivative(knl)
# DirectionalTargetDerivative (temporarily?) removed
# elif what_operator == "S'":
# from sumpy.kernel import DirectionalTargetDerivative
# knl = DirectionalTargetDerivative(knl)
else:
raise RuntimeError(f"unrecognized operator '{what_operator}'")
return source_knl, target_knl
def draw_pot_figure(aspect_ratio,
nsrc=100, novsmp=None, helmholtz_k=0,
what_operator="S",
what_operator_lpot=None,
order=4,
ovsmp_center_exp=0.66,
force_center_side=None):
if novsmp is None:
novsmp = 4*nsrc
if what_operator_lpot is None:
what_operator_lpot = what_operator
from sumpy.array_context import PyOpenCLArrayContext
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
actx = PyOpenCLArrayContext(queue)
# {{{ make plot targets
center = np.asarray([0, 0], dtype=np.float64)
from sumpy.visualization import FieldPlotter
fp = FieldPlotter(center, npoints=1000, extent=6)
# }}}
# {{{ make p2p kernel calculator
from sumpy.expansion.local import H2DLocalExpansion, LineTaylorLocalExpansion
from sumpy.kernel import HelmholtzKernel, LaplaceKernel
from sumpy.p2p import P2P
if helmholtz_k:
if isinstance(helmholtz_k, complex):
knl = HelmholtzKernel(2, allow_evanescent=True)
expn_class = H2DLocalExpansion
knl_kwargs = {"k": helmholtz_k}
else:
knl = HelmholtzKernel(2)
expn_class = H2DLocalExpansion
knl_kwargs = {"k": helmholtz_k}
else:
knl = LaplaceKernel(2)
expn_class = LineTaylorLocalExpansion
knl_kwargs = {}
vol_source_knl, vol_target_knl = process_kernel(knl, what_operator)
p2p = P2P(actx.context,
source_kernels=(vol_source_knl,),
target_kernels=(vol_target_knl,),
exclude_self=False,
value_dtypes=np.complex128)
lpot_source_knl, lpot_target_knl = process_kernel(knl, what_operator_lpot)
from sumpy.qbx import LayerPotential
lpot = LayerPotential(actx.context,
expansion=expn_class(knl, order=order),
source_kernels=(lpot_source_knl,),
target_kernels=(lpot_target_knl,),
value_dtypes=np.complex128)
# }}}
# {{{ set up geometry
# r,a,b match the corresponding letters from G. J. Rodin and O. Steinbach,
# Boundary Element Preconditioners for problems defined on slender domains.
# https://dx.doi.org/10.1137/S1064827500372067
a = 1
b = 1/aspect_ratio
def map_to_curve(t):
t = t*(2*np.pi)
x = a*np.cos(t)
y = b*np.sin(t)
w = (np.zeros_like(t)+1)/len(t)
return x, y, w
from curve import CurveGrid
native_t = np.linspace(0, 1, nsrc, endpoint=False)
native_x, native_y, native_weights = map_to_curve(native_t)
native_curve = CurveGrid(native_x, native_y)
ovsmp_t = np.linspace(0, 1, novsmp, endpoint=False)
ovsmp_x, ovsmp_y, ovsmp_weights = map_to_curve(ovsmp_t)
ovsmp_curve = CurveGrid(ovsmp_x, ovsmp_y)
curve_len = np.sum(ovsmp_weights * ovsmp_curve.speed)
hovsmp = curve_len/novsmp
center_dist = 5*hovsmp
if force_center_side is not None:
center_side = force_center_side*np.ones(len(native_curve))
else:
center_side = -np.sign(native_curve.mean_curvature)
centers = (native_curve.pos
+ center_side[:, np.newaxis]
* center_dist*native_curve.normal)
if 0:
native_curve.plot()
plt.show()
volpot_kwargs = knl_kwargs.copy()
lpot_kwargs = knl_kwargs.copy()
if what_operator == "D":
volpot_kwargs["src_derivative_dir"] = actx.from_numpy(native_curve.normal)
if what_operator_lpot == "D":
lpot_kwargs["src_derivative_dir"] = actx.from_numpy(ovsmp_curve.normal)
if what_operator_lpot == "S'":
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 scipy.sparse.linalg import LinearOperator
from sumpy.tools import build_matrix
def apply_lpot(x):
xovsmp = np.dot(fim, x)
_evt, (y,) = lpot(actx.queue,
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 actx.to_numpy(y)
op = LinearOperator((nsrc, nsrc), apply_lpot)
mat = build_matrix(op, dtype=np.complex128)
w, _v = la.eig(mat)
plt.plot(w.real, "o-")
# import sys; sys.exit(0)
return
# }}}
# {{{ compute potentials
mode_nr = 0
density = np.cos(mode_nr*2*np.pi*native_t).astype(np.complex128)
strength = actx.from_numpy(native_curve.speed * native_weights * density)
_evt, (vol_pot,) = p2p(actx.queue,
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,
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)
# }}}
if USE_MATPLOTLIB:
# {{{ plot on-surface potential in 2D
plt.plot(curve_pot, label="pot")
plt.plot(density, label="dens")
plt.legend()
plt.show()
# }}}
fp.write_vtk_file("potential.vts", [
("potential", vol_pot.real)
])
if USE_MATPLOTLIB:
# {{{ 2D false-color plot
plt.clf()
plotval = np.log10(1e-20+np.abs(vol_pot))
im = fp.show_scalar_in_matplotlib(plotval.real)
from matplotlib.colors import Normalize
im.set_norm(Normalize(vmin=-2, vmax=1))
src = native_curve.pos
plt.plot(src[:, 0], src[:, 1], "o-k")
# close the curve
plt.plot(src[-1::-len(src)+1, 0], src[-1::-len(src)+1, 1], "o-k")
cb = plt.colorbar(shrink=0.9)
cb.set_label(r"$\log_{10}(\mathdefault{Error})$")
fp.set_matplotlib_limits()
# }}}
else:
# {{{ 3D plots
plotval_vol = vol_pot.real
plotval_c = curve_pot.real
scale = 1
if 0:
# crop singularities--doesn't work very well
neighbors = [
np.roll(plotval_vol, 3, 0),
np.roll(plotval_vol, -3, 0),
np.roll(plotval_vol, 6, 0),
np.roll(plotval_vol, -6, 0),
]
avg = np.average(np.abs(plotval_vol))
outlier_flag = sum(
np.abs(plotval_vol-nb) for nb in neighbors) > avg
plotval_vol[outlier_flag] = sum(
nb[outlier_flag] for nb in neighbors)/len(neighbors)
if USE_MAYAVI:
fp.show_scalar_in_mayavi(scale*plotval_vol, max_val=1)
mlab.colorbar()
if 1:
mlab.points3d(
native_curve.pos[0],
native_curve.pos[1],
scale*plotval_c, scale_factor=0.02)
mlab.show()
# }}}
if __name__ == "__main__":
draw_pot_figure(
aspect_ratio=1, nsrc=100, novsmp=100, helmholtz_k=(35+4j)*0.3,
what_operator="D", what_operator_lpot="D", force_center_side=1)
if USE_MATPLOTLIB:
plt.savefig("eigvals-ext-nsrc100-novsmp100.pdf")
plt.clf()
# draw_pot_figure(
# aspect_ratio=1, nsrc=100, novsmp=100, helmholtz_k=0,
# what_operator="D", what_operator_lpot="D", force_center_side=-1)
# plt.savefig("eigvals-int-nsrc100-novsmp100.pdf")
# plt.clf()
# draw_pot_figure(
# aspect_ratio=1, nsrc=100, novsmp=200, helmholtz_k=0,
# what_operator="D", what_operator_lpot="D", force_center_side=-1)
# plt.savefig("eigvals-int-nsrc100-novsmp200.pdf")
# plt.clf()
# vim: fdm=marker