From dd22feb5ddfda9897f31ce1d74335ad3bbf63e49 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 16 Jan 2018 19:26:44 -0600 Subject: [PATCH 1/6] Clean up and test examples dir (Closes #86) --- .gitlab-ci.yml | 16 +- examples/fmm-error.py | 129 ++++++------- examples/{ => geometries}/blob-2d.step | 0 examples/{ => geometries}/circle.step | 0 examples/{ => geometries}/circles.step | 0 examples/{ => geometries}/ellipsoid.step | 0 examples/{ => geometries}/molecule.step | 0 examples/{ => geometries}/two-balls.step | 0 .../two-cylinders-smooth.step | 0 examples/helmholtz-dirichlet.py | 2 +- examples/laplace-dirichlet-3d.py | 170 +++++++++++++++++ examples/layerpot-3d.py | 125 +++++++------ examples/layerpot.py | 173 +++++++++--------- examples/perf-model.py | 18 -- examples/scaling-study.py | 4 +- experiments/README.md | 7 + {examples => experiments}/cahn-hilliard.py | 0 .../find-photonic-mode-sk.py | 0 .../find-photonic-mode.py | 0 .../helmholtz-expression-tree.py | 0 {examples => experiments}/maxwell.py | 0 {examples => experiments}/maxwell_sphere.py | 0 {examples => experiments}/poisson.py | 0 .../qbx-tangential-deriv-jump.py | 0 .../stokes-2d-interior.py | 0 .../two-domain-helmholtz.py | 0 26 files changed, 416 insertions(+), 228 deletions(-) rename examples/{ => geometries}/blob-2d.step (100%) rename examples/{ => geometries}/circle.step (100%) rename examples/{ => geometries}/circles.step (100%) rename examples/{ => geometries}/ellipsoid.step (100%) rename examples/{ => geometries}/molecule.step (100%) rename examples/{ => geometries}/two-balls.step (100%) rename examples/{ => geometries}/two-cylinders-smooth.step (100%) create mode 100644 examples/laplace-dirichlet-3d.py delete mode 100644 examples/perf-model.py create mode 100644 experiments/README.md rename {examples => experiments}/cahn-hilliard.py (100%) rename {examples => experiments}/find-photonic-mode-sk.py (100%) rename {examples => experiments}/find-photonic-mode.py (100%) rename {examples => experiments}/helmholtz-expression-tree.py (100%) rename {examples => experiments}/maxwell.py (100%) rename {examples => experiments}/maxwell_sphere.py (100%) rename {examples => experiments}/poisson.py (100%) rename {examples => experiments}/qbx-tangential-deriv-jump.py (100%) rename {examples => experiments}/stokes-2d-interior.py (100%) rename {examples => experiments}/two-domain-helmholtz.py (100%) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 334a41d0..a05b560d 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -26,6 +26,20 @@ Python 3.6 POCL: except: - tags +Python 3.6 POCL Examples: + script: + - export PY_EXE=python3.6 + - export PYOPENCL_TEST=portable + - export EXTRA_INSTALL="numpy mako" + - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-py-project-and-run-examples.sh + - ". ./build-py-project-and-run-examples.sh" + tags: + - python3.6 + - pocl + - large-node + except: + - tags + Python 3.5 Conda: script: - export SUMPY_FORCE_SYMBOLIC_BACKEND=symengine @@ -66,7 +80,7 @@ Documentation: Flake8: script: - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/prepare-and-run-flake8.sh - - ". ./prepare-and-run-flake8.sh pytential test" + - ". ./prepare-and-run-flake8.sh pytential test examples" tags: - python3.5 except: diff --git a/examples/fmm-error.py b/examples/fmm-error.py index fea97c99..4a15f18f 100644 --- a/examples/fmm-error.py +++ b/examples/fmm-error.py @@ -6,89 +6,90 @@ from meshmode.mesh.generation import ( # noqa from sumpy.visualization import FieldPlotter from sumpy.kernel import LaplaceKernel, HelmholtzKernel -import faulthandler -faulthandler.enable() -import logging -logging.basicConfig(level=logging.INFO) +def main(): + import logging + logging.basicConfig(level=logging.WARNING) # INFO for more progress info -cl_ctx = cl.create_some_context() -queue = cl.CommandQueue(cl_ctx) + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) -target_order = 16 -qbx_order = 3 -nelements = 60 -mode_nr = 0 + target_order = 16 + qbx_order = 3 + nelements = 60 + mode_nr = 0 -k = 0 -if k: - kernel = HelmholtzKernel("k") -else: - kernel = LaplaceKernel() -#kernel = OneKernel() + k = 0 + if k: + kernel = HelmholtzKernel(2) + else: + kernel = LaplaceKernel(2) + #kernel = OneKernel() -mesh = make_curve_mesh( - #lambda t: ellipse(1, t), - starfish, - np.linspace(0, 1, nelements+1), - target_order) + mesh = make_curve_mesh( + #lambda t: ellipse(1, t), + starfish, + np.linspace(0, 1, nelements+1), + target_order) -from pytential.qbx import QBXLayerPotentialSource -from meshmode.discretization import Discretization -from meshmode.discretization.poly_element import \ - InterpolatoryQuadratureSimplexGroupFactory + from pytential.qbx import QBXLayerPotentialSource + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory -density_discr = Discretization( - cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(target_order)) + pre_density_discr = Discretization( + cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) -qbx = QBXLayerPotentialSource( - density_discr, fine_order=2*target_order, - qbx_order=qbx_order, fmm_order=qbx_order) -slow_qbx = QBXLayerPotentialSource( - density_discr, fine_order=2*target_order, - qbx_order=qbx_order, fmm_order=False) + slow_qbx, _ = QBXLayerPotentialSource( + pre_density_discr, fine_order=2*target_order, + qbx_order=qbx_order, fmm_order=False, + target_association_tolerance=.05 + ).with_refinement() + qbx = slow_qbx.copy(fmm_order=10) + density_discr = slow_qbx.density_discr -nodes = density_discr.nodes().with_queue(queue) + nodes = density_discr.nodes().with_queue(queue) -angle = cl.clmath.atan2(nodes[1], nodes[0]) + angle = cl.clmath.atan2(nodes[1], nodes[0]) -from pytential import bind, sym -d = sym.Derivative() -#op = d.nabla[0] * d(sym.S(kernel, sym.var("sigma"))) -#op = sym.D(kernel, sym.var("sigma")) -op = sym.S(kernel, sym.var("sigma")) + from pytential import bind, sym + #op = sym.d_dx(sym.S(kernel, sym.var("sigma")), qbx_forced_limit=None) + #op = sym.D(kernel, sym.var("sigma"), qbx_forced_limit=None) + op = sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None) -sigma = cl.clmath.cos(mode_nr*angle) + sigma = cl.clmath.cos(mode_nr*angle) -if isinstance(kernel, HelmholtzKernel): - sigma = sigma.astype(np.complex128) + if isinstance(kernel, HelmholtzKernel): + sigma = sigma.astype(np.complex128) -bound_bdry_op = bind(qbx, op) + fplot = FieldPlotter(np.zeros(2), extent=5, npoints=600) + from pytential.target import PointsTarget -fplot = FieldPlotter(np.zeros(2), extent=5, npoints=600) -from pytential.target import PointsTarget + fld_in_vol = bind( + (slow_qbx, PointsTarget(fplot.points)), + op)(queue, sigma=sigma, k=k).get() -fld_in_vol = bind( - (slow_qbx, PointsTarget(fplot.points)), - op)(queue, sigma=sigma, k=k).get() + fmm_fld_in_vol = bind( + (qbx, PointsTarget(fplot.points)), + op)(queue, sigma=sigma, k=k).get() -fmm_fld_in_vol = bind( - (qbx, PointsTarget(fplot.points)), - op)(queue, sigma=sigma, k=k).get() + err = fmm_fld_in_vol-fld_in_vol + im = fplot.show_scalar_in_matplotlib(np.log10(np.abs(err))) -err = fmm_fld_in_vol-fld_in_vol -im = fplot.show_scalar_in_matplotlib(np.log10(np.abs(err))) + from matplotlib.colors import Normalize + im.set_norm(Normalize(vmin=-12, vmax=0)) -from matplotlib.colors import Normalize -im.set_norm(Normalize(vmin=-6, vmax=0)) + import matplotlib.pyplot as pt + from matplotlib.ticker import NullFormatter + pt.gca().xaxis.set_major_formatter(NullFormatter()) + pt.gca().yaxis.set_major_formatter(NullFormatter()) -import matplotlib.pyplot as pt -from matplotlib.ticker import NullFormatter -pt.gca().xaxis.set_major_formatter(NullFormatter()) -pt.gca().yaxis.set_major_formatter(NullFormatter()) + cb = pt.colorbar(shrink=0.9) + cb.set_label(r"$\log_{10}(\mathdefault{Error})$") -cb = pt.colorbar(shrink=0.9) -cb.set_label(r"$\log_{10}(\mathdefault{Error})$") + pt.savefig("fmm-error-order-%d.pdf" % qbx_order) -pt.savefig("fmm-error-order-%d.pdf" % qbx_order) + +if __name__ == "__main__": + main() diff --git a/examples/blob-2d.step b/examples/geometries/blob-2d.step similarity index 100% rename from examples/blob-2d.step rename to examples/geometries/blob-2d.step diff --git a/examples/circle.step b/examples/geometries/circle.step similarity index 100% rename from examples/circle.step rename to examples/geometries/circle.step diff --git a/examples/circles.step b/examples/geometries/circles.step similarity index 100% rename from examples/circles.step rename to examples/geometries/circles.step diff --git a/examples/ellipsoid.step b/examples/geometries/ellipsoid.step similarity index 100% rename from examples/ellipsoid.step rename to examples/geometries/ellipsoid.step diff --git a/examples/molecule.step b/examples/geometries/molecule.step similarity index 100% rename from examples/molecule.step rename to examples/geometries/molecule.step diff --git a/examples/two-balls.step b/examples/geometries/two-balls.step similarity index 100% rename from examples/two-balls.step rename to examples/geometries/two-balls.step diff --git a/examples/two-cylinders-smooth.step b/examples/geometries/two-cylinders-smooth.step similarity index 100% rename from examples/two-cylinders-smooth.step rename to examples/geometries/two-cylinders-smooth.step diff --git a/examples/helmholtz-dirichlet.py b/examples/helmholtz-dirichlet.py index d6d9baa6..d5498fe3 100644 --- a/examples/helmholtz-dirichlet.py +++ b/examples/helmholtz-dirichlet.py @@ -25,7 +25,7 @@ k = 3 def main(): import logging - logging.basicConfig(level=logging.INFO) + logging.basicConfig(level=logging.WARNING) # INFO for more progress info cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) diff --git a/examples/laplace-dirichlet-3d.py b/examples/laplace-dirichlet-3d.py new file mode 100644 index 00000000..b2e895ca --- /dev/null +++ b/examples/laplace-dirichlet-3d.py @@ -0,0 +1,170 @@ +import numpy as np +import numpy.linalg as la +import pyopencl as cl +import pyopencl.clmath # noqa + +from meshmode.discretization import Discretization +from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + +from pytential import bind, sym, norm # noqa +from pytential.target import PointsTarget + +# {{{ set some constants for use below + +nelements = 20 +bdry_quad_order = 4 +mesh_order = bdry_quad_order +qbx_order = bdry_quad_order +bdry_ovsmp_quad_order = 4*bdry_quad_order +fmm_order = 3 + +# }}} + + +def main(): + import logging + logging.basicConfig(level=logging.WARNING) # INFO for more progress info + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + from meshmode.mesh.generation import generate_torus + + rout = 10 + rin = 1 + if 1: + base_mesh = generate_torus( + rout, rin, 40, 4, + mesh_order) + + from meshmode.mesh.processing import affine_map, merge_disjoint_meshes + # nx = 1 + # ny = 1 + nz = 1 + dz = 0 + meshes = [ + affine_map( + base_mesh, + A=np.diag([1, 1, 1]), + b=np.array([0, 0, iz*dz])) + for iz in range(nz)] + + mesh = merge_disjoint_meshes(meshes, single_group=True) + + if 0: + from meshmode.mesh.visualization import draw_curve + draw_curve(mesh) + import matplotlib.pyplot as plt + plt.show() + + pre_density_discr = Discretization( + cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(bdry_quad_order)) + + from pytential.qbx import ( + QBXLayerPotentialSource, QBXTargetAssociationFailedException) + qbx, _ = QBXLayerPotentialSource( + pre_density_discr, fine_order=bdry_ovsmp_quad_order, qbx_order=qbx_order, + fmm_order=fmm_order + ).with_refinement() + density_discr = qbx.density_discr + + # {{{ describe bvp + + from sumpy.kernel import LaplaceKernel + kernel = LaplaceKernel(3) + + cse = sym.cse + + sigma_sym = sym.var("sigma") + #sqrt_w = sym.sqrt_jac_q_weight(3) + sqrt_w = 1 + inv_sqrt_w_sigma = cse(sigma_sym/sqrt_w) + + # -1 for interior Dirichlet + # +1 for exterior Dirichlet + loc_sign = +1 + + bdry_op_sym = (loc_sign*0.5*sigma_sym + + sqrt_w*( + sym.S(kernel, inv_sqrt_w_sigma) + + sym.D(kernel, inv_sqrt_w_sigma) + )) + + # }}} + + bound_op = bind(qbx, bdry_op_sym) + + # {{{ fix rhs and solve + + nodes = density_discr.nodes().with_queue(queue) + source = np.array([rout, 0, 0]) + + def u_incoming_func(x): + # return 1/cl.clmath.sqrt( (x[0] - source[0])**2 + # +(x[1] - source[1])**2 + # +(x[2] - source[2])**2 ) + return 1.0/la.norm(x.get()-source[:, None], axis=0) + + bc = cl.array.to_device(queue, u_incoming_func(nodes)) + + bvp_rhs = bind(qbx, sqrt_w*sym.var("bc"))(queue, bc=bc) + + from pytential.solve import gmres + gmres_result = gmres( + bound_op.scipy_op(queue, "sigma", dtype=np.float64), + bvp_rhs, tol=1e-14, progress=True, + stall_iterations=0, + hard_failure=True) + + sigma = bind(qbx, sym.var("sigma")/sqrt_w)(queue, sigma=gmres_result.solution) + + # }}} + + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, density_discr, 20) + bdry_vis.write_vtk_file("laplace.vtu", [ + ("sigma", sigma), + ]) + + # {{{ postprocess/visualize + + repr_kwargs = dict(qbx_forced_limit=None) + representation_sym = ( + sym.S(kernel, inv_sqrt_w_sigma, **repr_kwargs) + + sym.D(kernel, inv_sqrt_w_sigma, **repr_kwargs)) + + from sumpy.visualization import FieldPlotter + fplot = FieldPlotter(np.zeros(3), extent=20, npoints=50) + + targets = cl.array.to_device(queue, fplot.points) + + qbx_stick_out = qbx.copy(target_stick_out_factor=0.2) + + try: + fld_in_vol = bind( + (qbx_stick_out, PointsTarget(targets)), + representation_sym)(queue, sigma=sigma).get() + except QBXTargetAssociationFailedException as e: + fplot.write_vtk_file( + "failed-targets.vts", + [ + ("failed", e.failed_target_flags.get(queue)) + ] + ) + raise + + #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) + fplot.write_vtk_file( + "potential.vts", + [ + ("potential", fld_in_vol), + ] + ) + + # }}} + + +if __name__ == "__main__": + main() diff --git a/examples/layerpot-3d.py b/examples/layerpot-3d.py index 0a35ebd7..cc620997 100644 --- a/examples/layerpot-3d.py +++ b/examples/layerpot-3d.py @@ -21,10 +21,10 @@ qbx_order = 3 mode_nr = 4 if 1: - cad_file_name = "ellipsoid.step" + cad_file_name = "geometries/ellipsoid.step" h = 0.6 else: - cad_file_name = "two-cylinders-smooth.step" + cad_file_name = "geometries/two-cylinders-smooth.step" h = 0.4 k = 0 @@ -34,76 +34,85 @@ else: kernel = LaplaceKernel(3) #kernel = OneKernel() -from meshmode.mesh.io import generate_gmsh, FileSource -mesh = generate_gmsh( - FileSource(cad_file_name), 2, order=2, - other_options=["-string", "Mesh.CharacteristicLengthMax = %g;" % h]) -from meshmode.mesh.processing import perform_flips -# Flip elements--gmsh generates inside-out geometry. -mesh = perform_flips(mesh, np.ones(mesh.nelements)) +def main(): + import logging + logging.basicConfig(level=logging.WARNING) # INFO for more progress info -from meshmode.mesh.processing import find_bounding_box -bbox_min, bbox_max = find_bounding_box(mesh) -bbox_center = 0.5*(bbox_min+bbox_max) -bbox_size = max(bbox_max-bbox_min) / 2 + from meshmode.mesh.io import generate_gmsh, FileSource + mesh = generate_gmsh( + FileSource(cad_file_name), 2, order=2, + other_options=["-string", "Mesh.CharacteristicLengthMax = %g;" % h]) -logger.info("%d elements" % mesh.nelements) + from meshmode.mesh.processing import perform_flips + # Flip elements--gmsh generates inside-out geometry. + mesh = perform_flips(mesh, np.ones(mesh.nelements)) -from pytential.qbx import QBXLayerPotentialSource -from meshmode.discretization import Discretization -from meshmode.discretization.poly_element import \ - InterpolatoryQuadratureSimplexGroupFactory + from meshmode.mesh.processing import find_bounding_box + bbox_min, bbox_max = find_bounding_box(mesh) + bbox_center = 0.5*(bbox_min+bbox_max) + bbox_size = max(bbox_max-bbox_min) / 2 -density_discr = Discretization( - cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) + logger.info("%d elements" % mesh.nelements) -qbx, _ = QBXLayerPotentialSource(density_discr, 4*target_order, qbx_order, - fmm_order=qbx_order + 3, - target_association_tolerance=0.15).with_refinement() + from pytential.qbx import QBXLayerPotentialSource + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory -nodes = density_discr.nodes().with_queue(queue) + density_discr = Discretization( + cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) -angle = cl.clmath.atan2(nodes[1], nodes[0]) + qbx, _ = QBXLayerPotentialSource(density_discr, 4*target_order, qbx_order, + fmm_order=qbx_order + 3, + target_association_tolerance=0.15).with_refinement() -from pytential import bind, sym -#op = sym.d_dx(sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None)) -op = sym.D(kernel, sym.var("sigma"), qbx_forced_limit=None) -#op = sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None) + nodes = density_discr.nodes().with_queue(queue) -sigma = cl.clmath.cos(mode_nr*angle) -if 0: - sigma = 0*angle - from random import randrange - for i in range(5): - sigma[randrange(len(sigma))] = 1 + angle = cl.clmath.atan2(nodes[1], nodes[0]) -if isinstance(kernel, HelmholtzKernel): - sigma = sigma.astype(np.complex128) + from pytential import bind, sym + #op = sym.d_dx(sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None)) + op = sym.D(kernel, sym.var("sigma"), qbx_forced_limit=None) + #op = sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None) -fplot = FieldPlotter(bbox_center, extent=3.5*bbox_size, npoints=150) + sigma = cl.clmath.cos(mode_nr*angle) + if 0: + sigma = 0*angle + from random import randrange + for i in range(5): + sigma[randrange(len(sigma))] = 1 -from pytential.target import PointsTarget -fld_in_vol = bind( - (qbx, PointsTarget(fplot.points)), - op)(queue, sigma=sigma, k=k).get() + if isinstance(kernel, HelmholtzKernel): + sigma = sigma.astype(np.complex128) -#fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) -fplot.write_vtk_file( - "potential.vts", - [ - ("potential", fld_in_vol) - ] - ) + fplot = FieldPlotter(bbox_center, extent=3.5*bbox_size, npoints=150) -bdry_normals = bind( - density_discr, - sym.normal(density_discr.ambient_dim))(queue).as_vector(dtype=object) + from pytential.target import PointsTarget + fld_in_vol = bind( + (qbx, PointsTarget(fplot.points)), + op)(queue, sigma=sigma, k=k).get() -from meshmode.discretization.visualization import make_visualizer -bdry_vis = make_visualizer(queue, density_discr, target_order) + #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) + fplot.write_vtk_file( + "potential.vts", + [ + ("potential", fld_in_vol) + ] + ) -bdry_vis.write_vtk_file("source.vtu", [ - ("sigma", sigma), - ("bdry_normals", bdry_normals), - ]) + bdry_normals = bind( + density_discr, + sym.normal(density_discr.ambient_dim))(queue).as_vector(dtype=object) + + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, density_discr, target_order) + + bdry_vis.write_vtk_file("source.vtu", [ + ("sigma", sigma), + ("bdry_normals", bdry_normals), + ]) + + +if __name__ == "__main__": + main() diff --git a/examples/layerpot.py b/examples/layerpot.py index 0371b1ff..f5e56de3 100644 --- a/examples/layerpot.py +++ b/examples/layerpot.py @@ -36,96 +36,101 @@ else: kernel_kwargs = {} #kernel = OneKernel() -from meshmode.mesh.generation import ( # noqa - make_curve_mesh, starfish, ellipse, drop) -mesh = make_curve_mesh( - #lambda t: ellipse(1, t), - starfish, - np.linspace(0, 1, nelements+1), - target_order) -from pytential.qbx import QBXLayerPotentialSource -from meshmode.discretization import Discretization -from meshmode.discretization.poly_element import \ - InterpolatoryQuadratureSimplexGroupFactory +def main(): + from meshmode.mesh.generation import ( # noqa + make_curve_mesh, starfish, ellipse, drop) + mesh = make_curve_mesh( + #lambda t: ellipse(1, t), + starfish, + np.linspace(0, 1, nelements+1), + target_order) + + from pytential.qbx import QBXLayerPotentialSource + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + pre_density_discr = Discretization( + cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + qbx, _ = QBXLayerPotentialSource(pre_density_discr, 4*target_order, qbx_order, + fmm_order=qbx_order+3, + target_association_tolerance=0.005).with_refinement() + + density_discr = qbx.density_discr + + nodes = density_discr.nodes().with_queue(queue) + + angle = cl.clmath.atan2(nodes[1], nodes[0]) + + def op(**kwargs): + kwargs.update(kernel_kwargs) + + #op = sym.d_dx(sym.S(kernel, sym.var("sigma"), **kwargs)) + return sym.D(kernel, sym.var("sigma"), **kwargs) + #op = sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None, **kwargs) + + sigma = cl.clmath.cos(mode_nr*angle) + if 0: + sigma = 0*angle + from random import randrange + for i in range(5): + sigma[randrange(len(sigma))] = 1 + + if isinstance(kernel, HelmholtzKernel): + sigma = sigma.astype(np.complex128) + + bound_bdry_op = bind(qbx, op()) + #mlab.figure(bgcolor=(1, 1, 1)) + if 1: + fplot = FieldPlotter(np.zeros(2), extent=5, npoints=1000) + from pytential.target import PointsTarget + + targets_dev = cl.array.to_device(queue, fplot.points) + fld_in_vol = bind( + (qbx, PointsTarget(targets_dev)), + op(qbx_forced_limit=None))(queue, sigma=sigma, k=k).get() + + if enable_mayavi: + fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) + else: + fplot.write_vtk_file( + "potential.vts", + [ + ("potential", fld_in_vol) + ] + ) + + if 0: + def apply_op(density): + return bound_bdry_op( + queue, sigma=cl.array.to_device(queue, density), k=k).get() + + from sumpy.tools import build_matrix + n = len(sigma) + mat = build_matrix(apply_op, dtype=np.float64, shape=(n, n)) + + import matplotlib.pyplot as pt + pt.imshow(mat) + pt.colorbar() + pt.show() -pre_density_discr = Discretization( - cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) - -qbx, _ = QBXLayerPotentialSource(pre_density_discr, 4*target_order, qbx_order, - fmm_order=qbx_order+3, - target_association_tolerance=0.005).with_refinement() - -density_discr = qbx.density_discr - -nodes = density_discr.nodes().with_queue(queue) - -angle = cl.clmath.atan2(nodes[1], nodes[0]) - - -def op(**kwargs): - kwargs.update(kernel_kwargs) - - #op = sym.d_dx(sym.S(kernel, sym.var("sigma"), **kwargs)) - return sym.D(kernel, sym.var("sigma"), **kwargs) - #op = sym.S(kernel, sym.var("sigma"), qbx_forced_limit=None, **kwargs) - - -sigma = cl.clmath.cos(mode_nr*angle) -if 0: - sigma = 0*angle - from random import randrange - for i in range(5): - sigma[randrange(len(sigma))] = 1 + if enable_mayavi: + # {{{ plot boundary field -if isinstance(kernel, HelmholtzKernel): - sigma = sigma.astype(np.complex128) + fld_on_bdry = bound_bdry_op(queue, sigma=sigma, k=k).get() -bound_bdry_op = bind(qbx, op()) -#mlab.figure(bgcolor=(1, 1, 1)) -if 1: - fplot = FieldPlotter(np.zeros(2), extent=5, npoints=1000) - from pytential.target import PointsTarget + nodes_host = density_discr.nodes().get(queue=queue) + mlab.points3d(nodes_host[0], nodes_host[1], + fld_on_bdry.real, scale_factor=0.03) - targets_dev = cl.array.to_device(queue, fplot.points) - fld_in_vol = bind( - (qbx, PointsTarget(targets_dev)), - op(qbx_forced_limit=None))(queue, sigma=sigma, k=k).get() + # }}} if enable_mayavi: - fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) - else: - fplot.write_vtk_file( - "potential.vts", - [ - ("potential", fld_in_vol) - ] - ) - -if 0: - def apply_op(density): - return bound_bdry_op( - queue, sigma=cl.array.to_device(queue, density), k=k).get() - - from sumpy.tools import build_matrix - n = len(sigma) - mat = build_matrix(apply_op, dtype=np.float64, shape=(n, n)) - - import matplotlib.pyplot as pt - pt.imshow(mat) - pt.colorbar() - pt.show() - -if enable_mayavi: - # {{{ plot boundary field + mlab.colorbar() + mlab.show() - fld_on_bdry = bound_bdry_op(queue, sigma=sigma, k=k).get() - nodes_host = density_discr.nodes().get(queue=queue) - mlab.points3d(nodes_host[0], nodes_host[1], fld_on_bdry.real, scale_factor=0.03) - - # }}} - -if enable_mayavi: - mlab.colorbar() - mlab.show() +if __name__ == "__main__": + main() diff --git a/examples/perf-model.py b/examples/perf-model.py deleted file mode 100644 index 3a87d631..00000000 --- a/examples/perf-model.py +++ /dev/null @@ -1,18 +0,0 @@ -# ------------------------------ -nlevels = 6 -nboxes = 1365 -nsources = 60 -ntargets = 12040 -form_mp = 60*p_fmm -prop_upward = 1365*p_fmm**2 -part_direct = 196560 -m2l = 31920*p_fmm**2 -mp_eval = 0 -form_local = 65000*p_fmm -prop_downward = 1365*p_fmm**2 -eval_part = 12040*p_fmm -ncenters = 2040 -qbxl_direct = 2370940*p_qbx -qbx_m2l = 35339*p_fmm*p_qbx -qbx_l2l = 2040*p_fmm*p_qbx -qbx_eval = 1902*p_qbx diff --git a/examples/scaling-study.py b/examples/scaling-study.py index 183fc915..4d20bcc1 100644 --- a/examples/scaling-study.py +++ b/examples/scaling-study.py @@ -16,8 +16,8 @@ bdry_quad_order = 4 mesh_order = bdry_quad_order qbx_order = bdry_quad_order bdry_ovsmp_quad_order = 4*bdry_quad_order -fmm_order = 25 -k = 25 +fmm_order = 10 +k = 0 # }}} diff --git a/experiments/README.md b/experiments/README.md new file mode 100644 index 00000000..d0f56efd --- /dev/null +++ b/experiments/README.md @@ -0,0 +1,7 @@ +# Experiments + +What you find in this directory are experiments that *may* have done something +useful at some point (or not). Unlike `examples`, they are not being tested on +an ongoing basis. + +So if what you find here breaks for you, you get to keep both pieces. diff --git a/examples/cahn-hilliard.py b/experiments/cahn-hilliard.py similarity index 100% rename from examples/cahn-hilliard.py rename to experiments/cahn-hilliard.py diff --git a/examples/find-photonic-mode-sk.py b/experiments/find-photonic-mode-sk.py similarity index 100% rename from examples/find-photonic-mode-sk.py rename to experiments/find-photonic-mode-sk.py diff --git a/examples/find-photonic-mode.py b/experiments/find-photonic-mode.py similarity index 100% rename from examples/find-photonic-mode.py rename to experiments/find-photonic-mode.py diff --git a/examples/helmholtz-expression-tree.py b/experiments/helmholtz-expression-tree.py similarity index 100% rename from examples/helmholtz-expression-tree.py rename to experiments/helmholtz-expression-tree.py diff --git a/examples/maxwell.py b/experiments/maxwell.py similarity index 100% rename from examples/maxwell.py rename to experiments/maxwell.py diff --git a/examples/maxwell_sphere.py b/experiments/maxwell_sphere.py similarity index 100% rename from examples/maxwell_sphere.py rename to experiments/maxwell_sphere.py diff --git a/examples/poisson.py b/experiments/poisson.py similarity index 100% rename from examples/poisson.py rename to experiments/poisson.py diff --git a/examples/qbx-tangential-deriv-jump.py b/experiments/qbx-tangential-deriv-jump.py similarity index 100% rename from examples/qbx-tangential-deriv-jump.py rename to experiments/qbx-tangential-deriv-jump.py diff --git a/examples/stokes-2d-interior.py b/experiments/stokes-2d-interior.py similarity index 100% rename from examples/stokes-2d-interior.py rename to experiments/stokes-2d-interior.py diff --git a/examples/two-domain-helmholtz.py b/experiments/two-domain-helmholtz.py similarity index 100% rename from examples/two-domain-helmholtz.py rename to experiments/two-domain-helmholtz.py -- GitLab From fc1afe493f5346d6676562e7643b12bb88c69cbf Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 17 Jan 2018 00:45:45 -0600 Subject: [PATCH 2/6] More examples cleanup --- .gitlab-ci.yml | 2 +- examples/layerpot.py | 13 ++++++------- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index a05b560d..8d25d57a 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -30,7 +30,7 @@ Python 3.6 POCL Examples: script: - export PY_EXE=python3.6 - export PYOPENCL_TEST=portable - - export EXTRA_INSTALL="numpy mako" + - export EXTRA_INSTALL="numpy mako pyvisfile" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-py-project-and-run-examples.sh - ". ./build-py-project-and-run-examples.sh" tags: diff --git a/examples/layerpot.py b/examples/layerpot.py index f5e56de3..222606d2 100644 --- a/examples/layerpot.py +++ b/examples/layerpot.py @@ -15,13 +15,6 @@ import faulthandler from six.moves import range faulthandler.enable() -import logging -logging.basicConfig(level=logging.INFO) -logger = logging.getLogger(__name__) - -cl_ctx = cl.create_some_context() -queue = cl.CommandQueue(cl_ctx) - target_order = 16 qbx_order = 3 nelements = 60 @@ -38,6 +31,12 @@ else: def main(): + import logging + logging.basicConfig(level=logging.WARNING) # INFO for more progress info + + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + from meshmode.mesh.generation import ( # noqa make_curve_mesh, starfish, ellipse, drop) mesh = make_curve_mesh( -- GitLab From 05276326baecf79bcf15e48d6d278bb7cea090f7 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 17 Jan 2018 11:33:34 -0600 Subject: [PATCH 3/6] Tweak examples --- examples/helmholtz-dirichlet.py | 2 +- examples/laplace-dirichlet-3d.py | 2 +- examples/layerpot-3d.py | 9 +++------ examples/layerpot.py | 2 +- examples/scaling-study.py | 3 ++- 5 files changed, 8 insertions(+), 10 deletions(-) diff --git a/examples/helmholtz-dirichlet.py b/examples/helmholtz-dirichlet.py index d5498fe3..22f8fa8a 100644 --- a/examples/helmholtz-dirichlet.py +++ b/examples/helmholtz-dirichlet.py @@ -167,7 +167,7 @@ def main(): #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) fplot.write_vtk_file( - "potential.vts", + "potential-helm.vts", [ ("potential", fld_in_vol), ("indicator", indicator), diff --git a/examples/laplace-dirichlet-3d.py b/examples/laplace-dirichlet-3d.py index b2e895ca..4166dddf 100644 --- a/examples/laplace-dirichlet-3d.py +++ b/examples/laplace-dirichlet-3d.py @@ -157,7 +157,7 @@ def main(): #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) fplot.write_vtk_file( - "potential.vts", + "potential-laplace-3d.vts", [ ("potential", fld_in_vol), ] diff --git a/examples/layerpot-3d.py b/examples/layerpot-3d.py index cc620997..28f0967e 100644 --- a/examples/layerpot-3d.py +++ b/examples/layerpot-3d.py @@ -9,10 +9,6 @@ import faulthandler from six.moves import range faulthandler.enable() -import logging -logging.basicConfig(level=logging.INFO) -logger = logging.getLogger(__name__) - cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) @@ -37,6 +33,7 @@ else: def main(): import logging + logger = logging.getLogger(__name__) logging.basicConfig(level=logging.WARNING) # INFO for more progress info from meshmode.mesh.io import generate_gmsh, FileSource @@ -95,7 +92,7 @@ def main(): #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) fplot.write_vtk_file( - "potential.vts", + "potential-3d.vts", [ ("potential", fld_in_vol) ] @@ -108,7 +105,7 @@ def main(): from meshmode.discretization.visualization import make_visualizer bdry_vis = make_visualizer(queue, density_discr, target_order) - bdry_vis.write_vtk_file("source.vtu", [ + bdry_vis.write_vtk_file("source-3d.vtu", [ ("sigma", sigma), ("bdry_normals", bdry_normals), ]) diff --git a/examples/layerpot.py b/examples/layerpot.py index 222606d2..7b4737da 100644 --- a/examples/layerpot.py +++ b/examples/layerpot.py @@ -95,7 +95,7 @@ def main(): fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) else: fplot.write_vtk_file( - "potential.vts", + "potential-2d.vts", [ ("potential", fld_in_vol) ] diff --git a/examples/scaling-study.py b/examples/scaling-study.py index 4d20bcc1..ce6fc88f 100644 --- a/examples/scaling-study.py +++ b/examples/scaling-study.py @@ -154,6 +154,7 @@ def timing_run(nx, ny): sym_op)( queue, sigma=ones_density).get() + qbx_stick_out = qbx.copy(target_stick_out_factor=0.1) try: fld_in_vol = bind( (qbx_stick_out, PointsTarget(targets)), @@ -169,7 +170,7 @@ def timing_run(nx, ny): #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) fplot.write_vtk_file( - "potential.vts", + "potential-scaling.vts", [ ("potential", fld_in_vol), ("indicator", indicator) -- GitLab From 94291e4c2d75b25f5f58cbf1dbf1264aaf5452d3 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 17 Jan 2018 12:09:30 -0600 Subject: [PATCH 4/6] Install matplotlib for example execution --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 8d25d57a..ad62d954 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -30,7 +30,7 @@ Python 3.6 POCL Examples: script: - export PY_EXE=python3.6 - export PYOPENCL_TEST=portable - - export EXTRA_INSTALL="numpy mako pyvisfile" + - export EXTRA_INSTALL="numpy mako pyvisfile matplotlib" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-py-project-and-run-examples.sh - ". ./build-py-project-and-run-examples.sh" tags: -- GitLab From f0ddfea87f36c5f2a83e13b8b40099e014caca41 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 18 Jan 2018 17:31:57 -0600 Subject: [PATCH 5/6] Make matplotlib go without a display --- examples/fmm-error.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/fmm-error.py b/examples/fmm-error.py index 4a15f18f..110ec66b 100644 --- a/examples/fmm-error.py +++ b/examples/fmm-error.py @@ -75,6 +75,9 @@ def main(): op)(queue, sigma=sigma, k=k).get() err = fmm_fld_in_vol-fld_in_vol + + import matplotlib + matplotlib.use('Agg') im = fplot.show_scalar_in_matplotlib(np.log10(np.abs(err))) from matplotlib.colors import Normalize -- GitLab From 344ad5e1603222163a04672ebe1601da1d5077fd Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 18 Jan 2018 17:59:53 -0600 Subject: [PATCH 6/6] Switch scaling study to use less verbose logging --- examples/scaling-study.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/scaling-study.py b/examples/scaling-study.py index ce6fc88f..3327e3c8 100644 --- a/examples/scaling-study.py +++ b/examples/scaling-study.py @@ -54,7 +54,7 @@ def make_mesh(nx, ny): def timing_run(nx, ny): import logging - logging.basicConfig(level=logging.INFO) + logging.basicConfig(level=logging.WARNING) # INFO for more progress info cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) -- GitLab