From 990ee00faa3c06d4b11c5e4d756a2a4251504730 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 27 Jun 2012 00:07:52 -0400 Subject: [PATCH] Make layer potential drawing toy more flexible. --- examples/curve-pot.py | 130 +++++++++++++++++++++++++++--------------- 1 file changed, 85 insertions(+), 45 deletions(-) diff --git a/examples/curve-pot.py b/examples/curve-pot.py index e1e57378..901e42f8 100644 --- a/examples/curve-pot.py +++ b/examples/curve-pot.py @@ -6,26 +6,48 @@ from pytools import Record +def process_kernel(knl, what_operator): + if what_operator == "S": + pass + elif what_operator == "S0": + from sumpy.kernel import TargetDerivative + knl = TargetDerivative(0, knl) + elif what_operator == "S1": + from sumpy.kernel import TargetDerivative + knl = TargetDerivative(1, knl) + elif what_operator == "D": + from sumpy.kernel import SourceDerivative + knl = SourceDerivative(knl) + 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", order=5, + nsrc=100, novsmp=None, helmholtz_k=0, + what_operator="S", + what_operator_lpot=None, + order=5, ovsmp_center_exp=0.66): if novsmp is None: novsmp = 4*nsrc + if what_operator_lpot is None: + what_operator_lpot = what_operator + ctx = cl.create_some_context() queue = cl.CommandQueue(ctx) - # {{{ make plot targets center = np.asarray([0,0], dtype=np.float64) from sumpy.visualization import FieldPlotter fp = FieldPlotter(center, points=1000, extent=3) - tgt = fp.points.copy() - # }}} # {{{ make p2p kernel calculator @@ -42,23 +64,18 @@ def draw_pot_figure(aspect_ratio, expn_class = LineTaylorLocalExpansion knl_kwargs = {} - if what_operator == "S": - pass - elif what_operator == "D": - from sumpy.kernel import SourceDerivative - knl = SourceDerivative(knl) - else: - raise RuntimeError("unrecognized operator '%s'" % what_operator) - - p2p = P2P(ctx, [knl], exclude_self=False, + vol_knl = process_kernel(knl, what_operator) + p2p = P2P(ctx, [vol_knl], exclude_self=False, value_dtypes=np.complex128) + lpot_knl = process_kernel(knl, what_operator_lpot) + from sumpy.layerpot import LayerPotential, JumpTermApplier lpot = LayerPotential(ctx, [ - expn_class(knl, order=order)], + expn_class(lpot_knl, order=order)], value_dtypes=np.complex128) - jt = JumpTermApplier(ctx, [knl], + jt = JumpTermApplier(ctx, [lpot_knl], value_dtypes=np.complex128) # }}} @@ -105,35 +122,45 @@ def draw_pot_figure(aspect_ratio, #native_curve.plot() #pt.show() - kwargs = knl_kwargs.copy() - if what_operator == "D": - kwargs["src_derivative_dir"] = native_curve.normal - ovsmp_kwargs = knl_kwargs.copy() + volpot_kwargs = knl_kwargs.copy() + lpot_kwargs = knl_kwargs.copy() + if what_operator == "D": - ovsmp_kwargs["src_derivative_dir"] = ovsmp_curve.normal + lpot_kwargs["src_derivative_dir"] = ovsmp_curve.normal + + if what_operator_lpot == "D": + lpot_kwargs["src_derivative_dir"] = ovsmp_curve.normal + + if what_operator_lpot == "S'": + lpot_kwargs["tgt_derivative_dir"] = native_curve.normal # }}} len(native_curve) 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) - evt, (vol_pot,) = p2p(queue, tgt, native_curve.pos, [density], **kwargs) + evt, (vol_pot,) = p2p(queue, fp.points, native_curve.pos, + [native_curve.speed*native_weights*density], **volpot_kwargs) evt, (curve_pot,) = lpot(queue, native_curve.pos, ovsmp_curve.pos, centers, [ovsmp_density], ovsmp_curve.speed, ovsmp_weights, - **ovsmp_kwargs) - class JumpTermArgs(Record): - pass - evt, (curve_pot,) = jt(queue, [curve_pot], JumpTermArgs( - src_derivative_dir=native_curve.normal, - density_0=density, - side=center_side, - normal=native_curve.normal, - )) - pt.plot(curve_pot, label="pot") - pt.plot(density, label="dens") - pt.legend() - pt.show() + **lpot_kwargs) + + if 1: + class JumpTermArgs(Record): + pass + evt, (curve_pot,) = jt(queue, [curve_pot], JumpTermArgs( + src_derivative_dir=native_curve.normal, + tgt_derivative_dir=native_curve.normal, + density_0=density, + side=center_side, + normal=native_curve.normal, + )) + if 0: + pt.plot(curve_pot, label="pot") + pt.plot(density, label="dens") + pt.legend() + pt.show() if 0: pt.clf() @@ -155,26 +182,39 @@ def draw_pot_figure(aspect_ratio, #pt.gca().yaxis.set_major_formatter(NullFormatter()) fp.set_matplotlib_limits() else: - if 0: - plotval_vol = np.log10(1e-20+np.abs(vol_pot)) - else: - plotval_vol = vol_pot.real - plotval_c = curve_pot.real + plotval_vol = vol_pot.real + plotval_c = curve_pot.real + + scale = 1 - scale = 3e-2 + 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) fp.show_scalar_in_mayavi(scale*plotval_vol, maxval=1) from mayavi import mlab mlab.colorbar() - mlab.points3d( - native_curve.pos[:,0], - native_curve.pos[:,1], - scale*plotval_c) + 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=300, helmholtz_k=0, what_operator="D") + draw_pot_figure(aspect_ratio=1, nsrc=1000, helmholtz_k=0, + what_operator="S", + what_operator_lpot="S'") pt.show() -- GitLab