diff --git a/examples/curve-pot.py b/examples/curve-pot.py index 7720ef5bb11e35e21dc6eae168279405ef0736cd..4a09b403f755d9bdf35a250aff75453a7b60047c 100644 --- a/examples/curve-pot.py +++ b/examples/curve-pot.py @@ -6,7 +6,6 @@ import matplotlib.pyplot as pt from pytools import Record - def process_kernel(knl, what_operator): if what_operator == "S": pass @@ -19,14 +18,16 @@ def process_kernel(knl, what_operator): elif what_operator == "D": from sumpy.kernel import SourceDerivative knl = SourceDerivative(knl) - elif what_operator == "S'": - from sumpy.kernel import DirectionalTargetDerivative - knl = DirectionalTargetDerivative(knl) + # DirectionalTargetDerivative (temporarily?) removed + # elif what_operator == "S'": + # from sumpy.kernel import DirectionalTargetDerivative + # knl = DirectionalTargetDerivative(knl) else: raise RuntimeError("unrecognized operator '%s'" % what_operator) return knl + def draw_pot_figure(aspect_ratio, nsrc=100, novsmp=None, helmholtz_k=0, what_operator="S", @@ -46,7 +47,7 @@ def draw_pot_figure(aspect_ratio, # {{{ make plot targets - center = np.asarray([0,0], dtype=np.float64) + center = np.asarray([0, 0], dtype=np.float64) from sumpy.visualization import FieldPlotter fp = FieldPlotter(center, points=1000, extent=3) @@ -113,7 +114,6 @@ def draw_pot_figure(aspect_ratio, curve_len = np.sum(ovsmp_weights * ovsmp_curve.speed) hovsmp = curve_len/novsmp - hnative = curve_len/nsrc center_dist = 5*hovsmp if force_center_side is not None: @@ -123,7 +123,7 @@ def draw_pot_figure(aspect_ratio, centers = (native_curve.pos + center_side[:, np.newaxis] - * center_dist*native_curve.normal) + * center_dist*native_curve.normal) #native_curve.plot() #pt.show() @@ -142,13 +142,14 @@ def draw_pot_figure(aspect_ratio, # }}} - if 1: + 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(queue, native_curve.pos, ovsmp_curve.pos, centers, @@ -188,7 +189,7 @@ def draw_pot_figure(aspect_ratio, # }}} - if 1: + if 0: # {{{ apply jump terms class JumpTermArgs(Record): @@ -255,23 +256,22 @@ def draw_pot_figure(aspect_ratio, 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) + plotval_vol[outlier_flag] = sum( + nb[outlier_flag] for nb in neighbors)/len(neighbors) fp.show_scalar_in_mayavi(scale*plotval_vol, maxval=1) from mayavi import mlab mlab.colorbar() if 1: mlab.points3d( - native_curve.pos[:,0], - native_curve.pos[:,1], + 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=0, what_operator="D", what_operator_lpot="D", force_center_side=1) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index c2d0a42551aba8cd21399600ccc7c9e9797409b5..1513595af0365d8209dd9dc7b9868436f1399a6e 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -19,8 +19,8 @@ class ExpansionBase(object): return self.kernel.dimensions @property - def is_complex(self): - return self.kernel.is_complex + def is_complex_valued(self): + return self.kernel.is_complex_valued def prepare_loopy_kernel(self, loopy_knl): return self.kernel.prepare_loopy_kernel(loopy_knl) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 2baaefe7bbbb60be4921d190b6d207d29a356063..3d42c4520d30cbc3dd108bc4417c365adb37a254 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -26,6 +26,7 @@ THE SOFTWARE. import loopy as lp import sympy as sp import numpy as np +from pymbolic.primitives import Expression from pymbolic.mapper import IdentityMapper @@ -34,7 +35,7 @@ from pymbolic.mapper import IdentityMapper class Kernel(object): """Basic kernel interface. - .. attribute:: is_complex + .. attribute:: is_complex_valued .. attribute:: dimensions *dimensions* is allowed to be *None* if the dimensionality is not yet @@ -112,7 +113,7 @@ class Kernel(object): # {{{ PDE kernels class LaplaceKernel(Kernel): - is_complex = False + is_complex_valued = False def __repr__(self): if self.dimensions is not None: @@ -170,7 +171,7 @@ class HelmholtzKernel(Kernel): return HelmholtzKernel(dimensions, self.helmholtz_k_name, self.allow_evanescent) - is_complex = True + is_complex_valued = True def prepare_loopy_kernel(self, loopy_knl): # does loopy_knl already know about hank1_01? @@ -222,13 +223,18 @@ class HelmholtzKernel(Kernel): class DifferenceKernel(Kernel): def __init__(self, kernel_plus, kernel_minus): - self.kernel_plus = kernel_plus - self.kernel_minus = kernel_minus + if (kernel_plus.get_base_kernel() is not kernel_plus + or kernel_minus.get_base_kernel() is not kernel_minus): + raise ValueError("kernels in difference kernel " + "must be basic, unwrapped PDE kernels") - if self.kernel_plus.dimensions != self.kernel_minus.dimensions: + if kernel_plus.dimensions != kernel_minus.dimensions: raise ValueError( "kernels in difference kernel have different dimensions") + self.kernel_plus = kernel_plus + self.kernel_minus = kernel_minus + Kernel.__init__(self, self.kernel_plus.dimensions) def fix_dimensions(self, dimensions): @@ -258,6 +264,25 @@ def normalize_kernel(kernel): return kernel +def count_derivatives(kernel): + if isinstance(kernel, DerivativeBase): + return 1 + count_derivatives(kernel.kernel) + elif isinstance(kernel, KernelWrapper): + return count_derivatives(kernel.kernel) + else: + return 0 + + +def remove_axis_target_derivatives(kernel): + if isinstance(kernel, AxisTargetDerivative): + return remove_axis_target_derivatives(kernel.kernel) + elif isinstance(kernel, KernelWrapper): + kernel.replace_inner_kernel( + remove_axis_target_derivatives(kernel.kernel)) + else: + return kernel + + # {{{ a kernel defined as wrapping another one--e.g., derivatives class KernelWrapper(Kernel): @@ -272,8 +297,8 @@ class KernelWrapper(Kernel): return self.kernel.prepare_loopy_kernel(loopy_knl) @property - def is_complex(self): - return self.kernel.is_complex + def is_complex_valued(self): + return self.kernel.is_complex_valued def get_expression(self, dist_vec): return self.kernel.get_expression(dist_vec) @@ -301,15 +326,33 @@ class KernelWrapper(Kernel): # {{{ derivatives +class KernelPartialDerivative(Expression): + """A placeholder for a partial derivative of a kernel, to be used + in the *derivative_expression* arguments of :class:`DirectionalSourceDerivative` + and :class:`DirectionalTargetDerivative`. + """ + + def __init__(self, axis): + self.axis = axis + + def __getinitargs__(self): + return (self.axis,) + + mapper_method = "map_kernel_partial_derivative" + + class DerivativeBase(KernelWrapper): pass -class TargetDerivative(DerivativeBase): +class AxisTargetDerivative(DerivativeBase): def __init__(self, axis, kernel): KernelWrapper.__init__(self, kernel) self.axis = axis + def replace_inner_kernel(self, new_inner_kernel): + return AxisTargetDerivative(new_inner_kernel, self.derivative_expression) + def postprocess_at_target(self, expr, bvec): expr = self.kernel.postprocess_at_target(expr, bvec) return expr.diff(bvec[self.axis]) @@ -331,11 +374,18 @@ class _VectorIndexPrefixer(IdentityMapper): class DirectionalTargetDerivative(DerivativeBase): - def __init__(self, kernel, dir_vec_name="tgt_derivative_dir", - dir_vec_dtype=np.float64): + def __init__(self, kernel, derivative_expression): + """ + :arg derivative_expression: A scalar expression for the directional + derivative to be computed. :class:`KernelPartialDerivative` is + used as a placeholder for partial derivatives of the kernel. + """ KernelWrapper.__init__(self, kernel) - self.dir_vec_name = dir_vec_name - self.dir_vec_dtype = dir_vec_dtype + self.derivative_expression = derivative_expression + + def replace_inner_kernel(self, new_inner_kernel): + return DirectionalTargetDerivative( + new_inner_kernel, self.derivative_expression) def transform_to_code(self, expr): from sumpy.codegen import VectorComponentRewriter @@ -357,18 +407,21 @@ class DirectionalTargetDerivative(DerivativeBase): return sum(dir_vec[axis]*expr.diff(bvec[axis]) for axis in range(dimensions)) - def get_args(self): - return self.kernel.get_args() + [ - lp.GlobalArg(self.dir_vec_name, self.dir_vec_dtype, - shape=("ntgt", self.dimensions), order="C")] +class DirectionalSourceDerivative(DerivativeBase): + def __init__(self, kernel, derivative_expression): + """ + :arg derivative_expression: A scalar expression for the directional + derivative to be computed. :class:`KernelPartialDerivative` is + used as a placeholder for partial derivatives of the kernel. + """ -class SourceDerivative(DerivativeBase): - def __init__(self, kernel, dir_vec_name="src_derivative_dir", - dir_vec_dtype=np.float64): KernelWrapper.__init__(self, kernel) - self.dir_vec_name = dir_vec_name - self.dir_vec_dtype = dir_vec_dtype + self.derivative_expression = derivative_expression + + def replace_inner_kernel(self, new_inner_kernel): + return DirectionalSourceDerivative( + new_inner_kernel, self.derivative_expression) def transform_to_code(self, expr): from sumpy.codegen import VectorComponentRewriter @@ -390,11 +443,6 @@ class SourceDerivative(DerivativeBase): return sum(-dir_vec[axis]*expr.diff(avec[axis]) for axis in range(dimensions)) - def get_args(self): - return self.kernel.get_args() + [ - lp.GlobalArg(self.dir_vec_name, self.dir_vec_dtype, - shape=("nsrc", self.dimensions), order="C")] - # }}} # vim: fdm=marker diff --git a/sumpy/layerpot.py b/sumpy/layerpot.py index 37e4470a0f09cb6c4ce2cadd6c36cd0cee05c840..d7209ff1c3ca9d3f92709ebf8506235419f304cb 100644 --- a/sumpy/layerpot.py +++ b/sumpy/layerpot.py @@ -257,7 +257,8 @@ class LayerPotentialMatrixGenerator(LayerPotentialBase): def find_jump_term(kernel, arg_provider): from sumpy.kernel import ( - TargetDerivative, SourceDerivative, + AxisTargetDerivative, + DirectionalSourceDerivative, DirectionalTargetDerivative, DerivativeBase) @@ -265,13 +266,13 @@ def find_jump_term(kernel, arg_provider): src_derivatives = [] while isinstance(kernel, DerivativeBase): - if isinstance(kernel, TargetDerivative): + if isinstance(kernel, AxisTargetDerivative): tgt_derivatives.append(kernel.axis) kernel = kernel.kernel elif isinstance(kernel, DirectionalTargetDerivative): tgt_derivatives.append(kernel.dir_vec_name) kernel = kernel.kernel - elif isinstance(kernel, SourceDerivative): + elif isinstance(kernel, DirectionalSourceDerivative): src_derivatives.append(kernel.dir_vec_name) kernel = kernel.kernel else: @@ -436,83 +437,6 @@ class _JumpTermSymbolicArgumentProvider(object): # }}} - -# {{{ jump term applier - -class JumpTermApplier(KernelComputation): - def __init__(self, ctx, kernels, density_usage=None, - value_dtypes=None, density_dtypes=None, - geometry_dtype=None, - options=[], name="jump_term", device=None): - KernelComputation.__init__(self, ctx, kernels, density_usage, - value_dtypes, density_dtypes, geometry_dtype, - name, options, device) - - from pytools import single_valued - self.dimensions = single_valued(knl.dimensions for knl in self.kernels) - - def get_kernel(self): - data_args = {} - - exprs = [find_jump_term( - k, - _JumpTermSymbolicArgumentProvider(data_args, - self.dimensions, - "density_%d" % self.strength_usage[i], - self.strength_dtypes[self.strength_usage[i]], - self.geometry_dtype)) - for i, k in enumerate(self.kernels)] - - data_args = list(data_args.itervalues()) - - operand_args = [ - lp.GlobalArg("result_%d" % i, dtype, shape="ntgt", order="C") - for i, dtype in enumerate(self.value_dtypes) - ]+[ - lp.GlobalArg("limit_%d" % i, dtype, shape="ntgt", order="C") - for i, dtype in enumerate(self.value_dtypes)] - - # FIXME (fast) special case for jump term == 0 - loopy_knl = lp.make_kernel(self.device, - "[ntgt] -> {[itgt]: 0<=itgt