From 0bc2d916a416bce55c1347088321d77d3a74643f Mon Sep 17 00:00:00 2001
From: Andreas Kloeckner <inform@tiker.net>
Date: Tue, 19 May 2015 11:19:46 -0500
Subject: [PATCH] Enable evanescent Helmholtz
---
examples/curve-pot.py | 44 ++++++++++++++++++++++++++-----------------
sumpy/codegen.py | 16 +++++++++++++++-
sumpy/qbx.py | 4 ----
3 files changed, 42 insertions(+), 22 deletions(-)
diff --git a/examples/curve-pot.py b/examples/curve-pot.py
index 15d14552..06ce7aeb 100644
--- a/examples/curve-pot.py
+++ b/examples/curve-pot.py
@@ -31,10 +31,13 @@ def draw_pot_figure(aspect_ratio,
nsrc=100, novsmp=None, helmholtz_k=0,
what_operator="S",
what_operator_lpot=None,
- order=5,
+ order=4,
ovsmp_center_exp=0.66,
force_center_side=None):
+ import logging
+ logging.basicConfig(level=logging.INFO)
+
if novsmp is None:
novsmp = 4*nsrc
@@ -48,7 +51,7 @@ def draw_pot_figure(aspect_ratio,
center = np.asarray([0, 0], dtype=np.float64)
from sumpy.visualization import FieldPlotter
- fp = FieldPlotter(center, npoints=1000, extent=3)
+ fp = FieldPlotter(center, npoints=1000, extent=6)
# }}}
@@ -58,9 +61,15 @@ def draw_pot_figure(aspect_ratio,
from sumpy.kernel import LaplaceKernel, HelmholtzKernel
from sumpy.expansion.local import H2DLocalExpansion, LineTaylorLocalExpansion
if helmholtz_k:
- knl = HelmholtzKernel(2)
- expn_class = H2DLocalExpansion
- knl_kwargs = {"k": 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
@@ -166,8 +175,9 @@ def draw_pot_figure(aspect_ratio,
# {{{ compute potentials
- density = np.cos(3*2*np.pi*native_t).astype(np.complex128)
- ovsmp_density = np.cos(3*2*np.pi*ovsmp_t).astype(np.complex128)
+ 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)
evt, (vol_pot,) = p2p(queue, fp.points, native_curve.pos,
[native_curve.speed*native_weights*density], **volpot_kwargs)
@@ -178,7 +188,7 @@ def draw_pot_figure(aspect_ratio,
# }}}
- if 1:
+ if 0:
# {{{ plot on-surface potential in 2D
pt.plot(curve_pot, label="pot")
@@ -247,16 +257,16 @@ def draw_pot_figure(aspect_ratio,
if __name__ == "__main__":
- draw_pot_figure(aspect_ratio=1, nsrc=100, novsmp=100, helmholtz_k=0,
+ draw_pot_figure(aspect_ratio=1, nsrc=100, novsmp=100, helmholtz_k=(15+4j)*0.3,
what_operator="D", what_operator_lpot="D", force_center_side=1)
pt.savefig("eigvals-ext-nsrc100-novsmp100.pdf")
- pt.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)
- pt.savefig("eigvals-int-nsrc100-novsmp100.pdf")
- pt.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)
- pt.savefig("eigvals-int-nsrc100-novsmp200.pdf")
+ #pt.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)
+ #pt.savefig("eigvals-int-nsrc100-novsmp100.pdf")
+ #pt.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)
+ #pt.savefig("eigvals-int-nsrc100-novsmp200.pdf")
# vim: fdm=marker
diff --git a/sumpy/codegen.py b/sumpy/codegen.py
index 29137877..f2382282 100644
--- a/sumpy/codegen.py
+++ b/sumpy/codegen.py
@@ -48,7 +48,7 @@ class SympyToPymbolicMapper(SympyToPymbolicMapperBase):
self.assignments = assignments
self.derivative_cse_names = set()
- def map_Derivative(self, expr):
+ def map_Derivative(self, expr): # noqa
# Sympy has picked up the habit of picking arguments out of derivatives
# and pronounce them common subexpressions. Me no like. Undo it, so
# that the bessel substitutor further down can do its job.
@@ -70,9 +70,12 @@ class SympyToPymbolicMapper(SympyToPymbolicMapperBase):
BESSEL_PREAMBLE = """//CL//
#include <pyopencl-bessel-j.cl>
#include <pyopencl-bessel-y.cl>
+#include <pyopencl-bessel-j-complex.cl>
"""
HANKEL_PREAMBLE = """//CL//
+#include <pyopencl-hankel-complex.cl>
+
typedef struct hank1_01_result_str
{
cdouble_t order0, order1;
@@ -85,6 +88,13 @@ hank1_01_result hank1_01(cdouble_t z)
result.order1 = cdouble_new(bessel_j1(z.x), bessel_y1(z.x));
return result;
}
+
+hank1_01_result hank1_01_complex(cdouble_t z)
+{
+ hank1_01_result result;
+ hankel_01_complex(z, &result.order0, &result.order1, 1);
+ return result;
+}
"""
@@ -120,6 +130,9 @@ def bessel_mangler(target, identifier, arg_dtypes):
raise NotImplementedError("Only the PyOpenCLTarget is supported as of now")
if identifier == "hank1_01":
+ if arg_dtypes[0].kind == "c":
+ identifier = "hank1_01_complex"
+
return (np.dtype(hank1_01_result_dtype),
identifier, (np.dtype(np.complex128),))
elif identifier == "bessel_jv":
@@ -227,6 +240,7 @@ class BesselDerivativeReplacer(IdentityMapper):
import sympy as sp
# AS (9.1.31)
+ # http://dlmf.nist.gov/10.6.7
if order >= 0:
order_str = str(order)
else:
diff --git a/sumpy/qbx.py b/sumpy/qbx.py
index 36b268dc..ef6b3d15 100644
--- a/sumpy/qbx.py
+++ b/sumpy/qbx.py
@@ -158,10 +158,6 @@ class LayerPotentialBase(KernelComputation):
loopy_knl = expn.prepare_loopy_kernel(loopy_knl)
loopy_knl = lp.tag_data_axes(loopy_knl, "center", "sep,C")
- # loopy_knl = lp.set_options(loopy_knl,
- # #write_cl=True, highlight_cl=True,
- # #trace_assignment_values=True
- # )
return loopy_knl
--
GitLab