diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 7deb4fae9ab9e7df8ec06142658d8cd6fca53356..8373a93d2b77d803ba2d7d54c8d84cfe91ebd078 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -11,6 +11,18 @@ Python 3.5 POCL:
except:
- tags
+Python 3.5 Conda:
+ script:
+ - export SUMPY_NO_CACHE=1
+ - CONDA_ENVIRONMENT=.test-conda-env-py3.yml
+ - REQUIREMENTS_TXT=.test-conda-env-py3-requirements.txt
+ - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project-within-miniconda.sh
+ - ". ./build-and-test-py-project-within-miniconda.sh"
+ tags:
+ - linux
+ except:
+ - tags
+
Python 2.7 POCL:
script:
- export PY_EXE=python2.7
diff --git a/.test-conda-env-py3-requirements.txt b/.test-conda-env-py3-requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7cefac6831754070d40eb8427b09f70a91b498cd
--- /dev/null
+++ b/.test-conda-env-py3-requirements.txt
@@ -0,0 +1,5 @@
+git+https://github.com/inducer/boxtree
+git+https://github.com/inducer/pymbolic
+git+https://github.com/inducer/loopy
+git+https://github.com/inducer/sumpy
+git+https://github.com/inducer/meshmode
diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml
new file mode 100644
index 0000000000000000000000000000000000000000..c098bd19824933d4891d473594390398f659dee2
--- /dev/null
+++ b/.test-conda-env-py3.yml
@@ -0,0 +1,16 @@
+name: test-conda-env-py3
+channels:
+- inducer
+- symengine/label/dev
+- conda-forge
+- defaults
+dependencies:
+- git
+- conda-forge::numpy
+- conda-forge::sympy
+- pocl=0.13
+- islpy
+- pyopencl
+- python=3.5
+- python-symengine=0.2.0.53.g83912b7=py35_1
+# things not in here: loopy boxtree pymbolic pyfmmlib meshmode sumpy
diff --git a/doc/qbx.rst b/doc/qbx.rst
index d297db5eeb540327dc159c0b02e0316f86f5939b..46244da1086be469757fed4251ebdd946ed83af3 100644
--- a/doc/qbx.rst
+++ b/doc/qbx.rst
@@ -3,6 +3,11 @@
QBX internals
=============
+Refinement
+----------
+
+.. automodule:: pytential.qbx.refinement
+
Data structures describing geometry
-----------------------------------
diff --git a/examples/cahn-hilliard.py b/examples/cahn-hilliard.py
new file mode 100644
index 0000000000000000000000000000000000000000..f7cd6559f06878bbfc7811fec1f6c5d814afc193
--- /dev/null
+++ b/examples/cahn-hilliard.py
@@ -0,0 +1,132 @@
+import numpy as np
+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 = 8
+
+# }}}
+
+
+def main():
+ import logging
+ logging.basicConfig(level=logging.INFO)
+
+ cl_ctx = cl.create_some_context()
+ queue = cl.CommandQueue(cl_ctx)
+
+ from meshmode.mesh.generation import ellipse, make_curve_mesh
+ from functools import partial
+
+ mesh = make_curve_mesh(
+ partial(ellipse, 2),
+ np.linspace(0, 1, nelements+1),
+ mesh_order)
+
+ 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,
+ expansion_disks_in_tree_have_extent=True,
+ ).with_refinement()
+ density_discr = qbx.density_discr
+
+ from pytential.symbolic.pde.cahn_hilliard import CahnHilliardOperator
+ chop = CahnHilliardOperator(
+ # FIXME: Constants?
+ lambda1=1.5,
+ lambda2=1.25,
+ c=1)
+
+ unk = chop.make_unknown("sigma")
+ bound_op = bind(qbx, chop.operator(unk))
+
+ # {{{ fix rhs and solve
+
+ nodes = density_discr.nodes().with_queue(queue)
+
+ def g(xvec):
+ x, y = xvec
+ return cl.clmath.atan2(y, x)
+
+ bc = sym.make_obj_array([
+ # FIXME: Realistic BC
+ g(nodes),
+ -g(nodes),
+ ])
+
+ from pytential.solve import gmres
+ gmres_result = gmres(
+ bound_op.scipy_op(queue, "sigma", dtype=np.complex128),
+ bc, tol=1e-8, progress=True,
+ stall_iterations=0,
+ hard_failure=True)
+
+ # }}}
+
+ # {{{ postprocess/visualize
+
+ sigma = gmres_result.solution
+
+ from sumpy.visualization import FieldPlotter
+ fplot = FieldPlotter(np.zeros(2), extent=5, npoints=500)
+
+ targets = cl.array.to_device(queue, fplot.points)
+
+ qbx_stick_out = qbx.copy(target_stick_out_factor=0.05)
+
+ indicator_qbx = qbx_stick_out.copy(qbx_order=2)
+
+ from sumpy.kernel import LaplaceKernel
+ ones_density = density_discr.zeros(queue)
+ ones_density.fill(1)
+ indicator = bind(
+ (indicator_qbx, PointsTarget(targets)),
+ sym.D(LaplaceKernel(2), sym.var("sigma")))(
+ queue, sigma=ones_density).get()
+
+ try:
+ fld_in_vol = bind(
+ (qbx_stick_out, PointsTarget(targets)),
+ chop.representation(unk))(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),
+ ("indicator", indicator),
+ ]
+ )
+
+ # }}}
+
+
+if __name__ == "__main__":
+ main()
diff --git a/examples/qbx-tangential-deriv-jump.py b/examples/qbx-tangential-deriv-jump.py
index a2d9aee155d050e3c19bd253b051d8767a21e77a..6af8bd6805884ea97df69fd703fdeedc6dc1b843 100644
--- a/examples/qbx-tangential-deriv-jump.py
+++ b/examples/qbx-tangential-deriv-jump.py
@@ -1,8 +1,8 @@
import pyopencl as cl
import numpy as np
-import numpy.linalg as la
+import numpy.linalg as la # noqa: F401
-from pytential import bind, sym, norm
+from pytential import bind, sym, norm # noqa: F401
def main():
@@ -38,15 +38,15 @@ def main():
sig_sym = sym.var("sig")
knl = LaplaceKernel(2)
op = join_fields(
- sym.tangential_derivative(
- sym.D(knl, sig_sym, qbx_forced_limit=+1)).a.as_scalar(),
- sym.tangential_derivative(
- sym.D(knl, sig_sym, qbx_forced_limit=-1)).a.as_scalar(),
+ sym.tangential_derivative(mesh.ambient_dim,
+ sym.D(knl, sig_sym, qbx_forced_limit=+1)).as_scalar(),
+ sym.tangential_derivative(mesh.ambient_dim,
+ sym.D(knl, sig_sym, qbx_forced_limit=-1)).as_scalar(),
)
nodes = density_discr.nodes().with_queue(queue)
angle = cl.clmath.atan2(nodes[1], nodes[0])
- n = 30
+ n = 10
sig = cl.clmath.sin(n*angle)
dt_sig = n*cl.clmath.cos(n*angle)
diff --git a/pytential/__init__.py b/pytential/__init__.py
index 98e62c7164105a4c8cb0d1225706ad666ce233d1..7a4e1f4ec91fed3ae83e72396a9e703271a5908b 100644
--- a/pytential/__init__.py
+++ b/pytential/__init__.py
@@ -48,7 +48,14 @@ def _set_up_logging_from_environment():
set_up_logging(pytential_log_var.split(":"), level=level)
+def _set_up_errors():
+ import warnings
+ from pytential.qbx.refinement import RefinerNotConvergedWarning
+ warnings.filterwarnings("error", category=RefinerNotConvergedWarning)
+
+
_set_up_logging_from_environment()
+_set_up_errors()
@memoize_on_first_arg
diff --git a/pytential/muller.py b/pytential/muller.py
index 22fe5458e895236e3ec0fcc6ab7634290e0a84ed..5fc6a82378af0e8c86ac5858187014399cf8c42d 100644
--- a/pytential/muller.py
+++ b/pytential/muller.py
@@ -25,7 +25,7 @@ THE SOFTWARE.
import numpy as np
-def muller_deflate(f, n, maxiter=100, eps=1e-14):
+def muller_deflate(f, n, maxiter=100, eps=1e-14, z_start=None):
"""
:arg n: number of zeros sought
:returns: (roots, niter) - roots of the given function; number of
@@ -52,7 +52,7 @@ def muller_deflate(f, n, maxiter=100, eps=1e-14):
niter.append(niter0)
while (np.isnan(roots[i]) or niter[i] == maxiter) and miter < 50:
- roots0, niter0 = muller(f_deflated, maxiter, eps)
+ roots0, niter0 = muller(f_deflated, maxiter, eps, z_start=z_start)
roots[i] = roots0
niter[i] = niter0
miter = miter+1
diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py
index a18196f401ce355cd5c42f2df0416cf5bc08eb88..4ecf5253e6c0f62f3d0717d7952932f146fc899f 100644
--- a/pytential/qbx/__init__.py
+++ b/pytential/qbx/__init__.py
@@ -25,11 +25,9 @@ THE SOFTWARE.
import six
-import loopy as lp
import numpy as np
from pytools import memoize_method
from meshmode.discretization import Discretization
-from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
from pytential.qbx.target_assoc import QBXTargetAssociationFailedException
from pytential.source import PotentialSource
@@ -46,45 +44,6 @@ __doc__ = """
"""
-# {{{ jump term interface helper
-
-class _JumpTermArgumentProvider(object):
- def __init__(self, discr, density, ds_direction, side=None):
- self.discr = discr
- self.density = density
- self.ds_direction = ds_direction
- self.side = side
-
- @property
- def normal(self):
- return self.discr.curve.normals.reshape(2, -1).T
-
- @property
- def tangent(self):
- return self.discr.curve.tangents.reshape(2, -1).T
-
- @property
- def src_derivative_dir(self):
- return self.ds_direction
-
- @property
- def mean_curvature(self):
- return self.discr.curve.curvature.reshape(-1)
-
- @property
- def density_0(self):
- return self.density.reshape(-1)
-
- @property
- @memoize_method
- def density_0_prime(self):
- diff_mat = self.discr.curve.expansion.get_differentiation_matrix()
- return (2 * np.dot(diff_mat, self.density.T).T.reshape(-1)
- / self.discr.curve.speed.reshape(-1))
-
-# }}}
-
-
class LayerPotentialSource(PotentialSource):
pass
@@ -129,29 +88,37 @@ class QBXLayerPotentialSource(LayerPotentialSource):
.. attribute :: qbx_order
.. attribute :: fmm_order
.. attribute :: cl_context
- .. automethod :: centers
- .. automethod :: panel_sizes
.. automethod :: weights_and_area_elements
.. automethod :: with_refinement
See :ref:`qbxguts` for some information on the inner workings of this.
"""
+
+ # {{{ constructor / copy
+
def __init__(self, density_discr, fine_order,
qbx_order=None,
fmm_order=None,
fmm_level_to_order=None,
target_stick_out_factor=1e-10,
+ base_resampler=None,
# begin undocumented arguments
# FIXME default debug=False once everything works
debug=True,
refined_for_global_qbx=False,
+ expansion_disks_in_tree_have_extent=False,
performance_data_file=None):
"""
:arg fine_order: The total degree to which the (upsampled)
- underlying quadrature is exact.
+ underlying quadrature is exact.
+ :arg base_resampler: A connection used for resampling from
+ *density_discr* the fine density discretization. It is assumed
+ that the fine density discretization given by
+ *base_resampler.to_discr* is *not* already upsampled. May
+ be *None*.
:arg fmm_order: `False` for direct calculation. ``None`` will set
- a reasonable(-ish?) default.
+ a reasonable(-ish?) default.
"""
if fmm_level_to_order is None:
@@ -174,8 +141,13 @@ class QBXLayerPotentialSource(LayerPotentialSource):
self.fmm_level_to_order = fmm_level_to_order
self.target_stick_out_factor = target_stick_out_factor
+ # Default values are lazily provided if these are None
+ self._base_resampler = base_resampler
+
self.debug = debug
self.refined_for_global_qbx = refined_for_global_qbx
+ self.expansion_disks_in_tree_have_extent = \
+ expansion_disks_in_tree_have_extent
self.performance_data_file = performance_data_file
def copy(
@@ -185,6 +157,7 @@ class QBXLayerPotentialSource(LayerPotentialSource):
qbx_order=None,
fmm_level_to_order=None,
target_stick_out_factor=None,
+ base_resampler=None,
debug=None,
refined_for_global_qbx=None,
@@ -199,69 +172,98 @@ class QBXLayerPotentialSource(LayerPotentialSource):
fmm_level_to_order=(
fmm_level_to_order or self.fmm_level_to_order),
target_stick_out_factor=(
- target_stick_out_factor or self.target_stick_out_factor),
+ target_stick_out_factor
+ if target_stick_out_factor is not None
+ else self.target_stick_out_factor),
+ base_resampler=base_resampler or self._base_resampler,
debug=(
debug if debug is not None else self.debug),
refined_for_global_qbx=(
refined_for_global_qbx if refined_for_global_qbx is not None
else self.refined_for_global_qbx),
+ expansion_disks_in_tree_have_extent=(
+ self.expansion_disks_in_tree_have_extent),
performance_data_file=self.performance_data_file)
+ # }}}
+
+ @property
+ def base_fine_density_discr(self):
+ """The refined, interpolation-focused density discretization (no oversampling).
+ """
+ # FIXME: Maybe rename refined_interp_density_discr
+ return (self._base_resampler.to_discr
+ if self._base_resampler is not None
+ else self.density_discr)
+
@property
@memoize_method
def fine_density_discr(self):
+ """The refined, quadrature-focused density discretization (with upsampling).
+ """
+ # FIXME: Maybe rename refined_quad_density_discr
+ from meshmode.discretization.poly_element import (
+ QuadratureSimplexGroupFactory)
+
return Discretization(
- self.density_discr.cl_context, self.density_discr.mesh,
- QuadratureSimplexGroupFactory(self.fine_order), self.real_dtype)
+ self.density_discr.cl_context, self.base_fine_density_discr.mesh,
+ QuadratureSimplexGroupFactory(self.fine_order),
+ self.real_dtype)
@property
@memoize_method
def resampler(self):
- from meshmode.discretization.connection import make_same_mesh_connection
- return make_same_mesh_connection(
- self.fine_density_discr, self.density_discr)
-
- def el_view(self, discr, group_nr, global_array):
- """Return a view of *global_array* of shape
- ``(..., discr.groups[group_nr].nelements)``
- where *global_array* is of shape ``(..., nelements)``,
- where *nelements* is the global (per-discretization) node count.
- """
+ from meshmode.discretization.connection import (
+ make_same_mesh_connection, ChainedDiscretizationConnection)
- group = discr.groups[group_nr]
- el_nr_base = sum(group.nelements for group in discr.groups[:group_nr])
+ conn = make_same_mesh_connection(
+ self.fine_density_discr, self.base_fine_density_discr)
- return global_array[
- ..., el_nr_base:el_nr_base + group.nelements] \
- .reshape(
- global_array.shape[:-1]
- + (group.nelements,))
+ if self._base_resampler is not None:
+ return ChainedDiscretizationConnection([self._base_resampler, conn])
+
+ return conn
+
+ @property
+ @memoize_method
+ def refiner_code_container(self):
+ from pytential.qbx.refinement import RefinerCodeContainer
+ return RefinerCodeContainer(self.cl_context)
@memoize_method
- def with_refinement(self, target_order=None, maxiter=3):
+ def with_refinement(self, target_order=None, kernel_length_scale=None,
+ maxiter=10):
"""
:returns: a tuple ``(lpot_src, cnx)``, where ``lpot_src`` is a
:class:`QBXLayerPotentialSource` and ``cnx`` is a
:class:`meshmode.discretization.connection.DiscretizationConnection`
from the originally given to the refined geometry.
"""
- from pytential.qbx.refinement import QBXLayerPotentialSourceRefiner
- refiner = QBXLayerPotentialSourceRefiner(self.cl_context)
+ from pytential.qbx.refinement import refine_for_global_qbx
+
from meshmode.discretization.poly_element import (
- InterpolatoryQuadratureSimplexGroupFactory)
+ InterpolatoryQuadratureSimplexGroupFactory,
+ QuadratureSimplexGroupFactory)
+
if target_order is None:
target_order = self.density_discr.groups[0].order
- lpot, connection = refiner(self,
+
+ lpot, connection = refine_for_global_qbx(
+ self,
+ self.refiner_code_container,
InterpolatoryQuadratureSimplexGroupFactory(target_order),
+ QuadratureSimplexGroupFactory(self.fine_order),
+ kernel_length_scale=kernel_length_scale,
maxiter=maxiter)
+
return lpot, connection
@property
@memoize_method
def h_max(self):
with cl.CommandQueue(self.cl_context) as queue:
- panel_sizes = self.panel_sizes("npanels").with_queue(queue)
+ panel_sizes = self._panel_sizes("npanels").with_queue(queue)
return np.asscalar(cl.array.max(panel_sizes).get())
@property
@@ -284,106 +286,72 @@ class QBXLayerPotentialSource(LayerPotentialSource):
def complex_dtype(self):
return self.density_discr.complex_dtype
- @memoize_method
- def panel_centers_of_mass(self):
- knl = lp.make_kernel(
- """{[dim,k,i]:
- 0<=dim<ndims and
- 0<=k<nelements and
- 0<=i<nunit_nodes}""",
- """
- panels[dim, k] = sum(i, nodes[dim, k, i])/nunit_nodes
- """,
- default_offset=lp.auto, name="find_panel_centers_of_mass")
-
- knl = lp.fix_parameters(knl, ndims=self.ambient_dim)
-
- knl = lp.split_iname(knl, "k", 128, inner_tag="l.0", outer_tag="g.0")
- knl = lp.tag_inames(knl, dict(dim="ilp"))
+ # {{{ internal API
- with cl.CommandQueue(self.cl_context) as queue:
- mesh = self.density_discr.mesh
- panels = cl.array.empty(queue, (mesh.ambient_dim, mesh.nelements),
- dtype=self.density_discr.real_dtype)
- for group_nr, group in enumerate(self.density_discr.groups):
- _, (result,) = knl(queue,
- nelements=group.nelements,
- nunit_nodes=group.nunit_nodes,
- nodes=group.view(self.density_discr.nodes()),
- panels=self.el_view(self.density_discr, group_nr, panels))
- panels.finish()
- panels = panels.with_queue(None)
- return tuple(panels[d, :] for d in range(mesh.ambient_dim))
+ @memoize_method
+ def _panel_centers_of_mass(self):
+ import pytential.qbx.utils as utils
+ return utils.element_centers_of_mass(self.density_discr)
@memoize_method
- def panel_sizes(self, last_dim_length="nsources"):
- assert last_dim_length in ("nsources", "ncenters", "npanels")
- # To get the panel size this does the equivalent of (∫ 1 ds)**(1/dim).
- # FIXME: Kernel optimizations
+ def _fine_panel_centers_of_mass(self):
+ import pytential.qbx.utils as utils
+ return utils.element_centers_of_mass(self.base_fine_density_discr)
- discr = self.density_discr
+ @memoize_method
+ def _expansion_radii(self, last_dim_length):
+ if last_dim_length == "npanels":
+ # FIXME: Make this an error
- if last_dim_length == "nsources" or last_dim_length == "ncenters":
- knl = lp.make_kernel(
- "{[i,j,k]: 0<=i<nelements and 0<=j,k<nunit_nodes}",
- "panel_sizes[i,j] = sum(k, ds[i,k])**(1/dim)",
- name="compute_size")
+ from warnings import warn
+ warn("Passing 'npanels' as last_dim_length to _expansion_radii is "
+ "deprecated. Expansion radii should be allowed to vary "
+ "within a panel.", stacklevel=3)
- def panel_size_view(discr, group_nr):
- return discr.groups[group_nr].view
+ with cl.CommandQueue(self.cl_context) as queue:
+ return (self._panel_sizes(last_dim_length).with_queue(queue) * 0.5
+ ).with_queue(None)
- elif last_dim_length == "npanels":
- knl = lp.make_kernel(
- "{[i,j]: 0<=i<nelements and 0<=j<nunit_nodes}",
- "panel_sizes[i] = sum(j, ds[i,j])**(1/dim)",
- name="compute_size")
- from functools import partial
+ # _expansion_radii should not be needed for the fine discretization
- def panel_size_view(discr, group_nr):
- return partial(self.el_view, discr, group_nr)
+ @memoize_method
+ def _close_target_tunnel_radius(self, last_dim_length):
+ with cl.CommandQueue(self.cl_context) as queue:
+ return (self._panel_sizes(last_dim_length).with_queue(queue) * 0.5
+ ).with_queue(None)
- else:
- raise ValueError("unknown dim length specified")
+ @memoize_method
+ def _panel_sizes(self, last_dim_length="npanels"):
+ import pytential.qbx.utils as utils
+ return utils.panel_sizes(self.density_discr, last_dim_length)
- knl = lp.fix_parameters(knl, dim=self.dim)
+ @memoize_method
+ def _fine_panel_sizes(self, last_dim_length="npanels"):
+ if last_dim_length != "npanels":
+ raise NotImplementedError()
- with cl.CommandQueue(self.cl_context) as queue:
- from pytential import bind, sym
- ds = bind(
- discr,
- sym.area_element(ambient_dim=discr.ambient_dim, dim=discr.dim)
- * sym.QWeight()
- )(queue)
- panel_sizes = cl.array.empty(
- queue, discr.nnodes
- if last_dim_length in ("nsources", "ncenters")
- else discr.mesh.nelements, discr.real_dtype)
- for group_nr, group in enumerate(discr.groups):
- _, (result,) = knl(queue,
- nelements=group.nelements,
- nunit_nodes=group.nunit_nodes,
- ds=group.view(ds),
- panel_sizes=panel_size_view(
- discr, group_nr)(panel_sizes))
- panel_sizes.finish()
- if last_dim_length == "ncenters":
- from pytential.qbx.utils import get_interleaver_kernel
- knl = get_interleaver_kernel(discr.real_dtype)
- _, (panel_sizes,) = knl(queue, dstlen=2*discr.nnodes,
- src1=panel_sizes, src2=panel_sizes)
- return panel_sizes.with_queue(None)
+ import pytential.qbx.utils as utils
+ return utils.panel_sizes(self.base_fine_density_discr, last_dim_length)
@memoize_method
- def centers(self, sign):
- adim = self.density_discr.ambient_dim
- dim = self.density_discr.dim
+ def qbx_fmm_geometry_data(self, target_discrs_and_qbx_sides):
+ """
+ :arg target_discrs_and_qbx_sides:
+ a tuple of *(discr, qbx_forced_limit)*
+ tuples, where *discr* is a
+ :class:`meshmode.discretization.Discretization`
+ or
+ :class:`pytential.target.TargetBase`
+ instance
+ """
+ from pytential.qbx.geometry import QBXFMMGeometryData
- from pytential import sym, bind
- with cl.CommandQueue(self.cl_context) as queue:
- nodes = bind(self.density_discr, sym.nodes(adim))(queue)
- normals = bind(self.density_discr, sym.normal(adim, dim=dim))(queue)
- panel_sizes = self.panel_sizes().with_queue(queue)
- return (nodes + normals * sign * panel_sizes / 2).as_vector(np.object)
+ return QBXFMMGeometryData(self.qbx_fmm_code_getter,
+ self, target_discrs_and_qbx_sides,
+ target_stick_out_factor=self.target_stick_out_factor,
+ debug=self.debug)
+
+ # }}}
@memoize_method
def weights_and_area_elements(self):
@@ -404,6 +372,8 @@ class QBXLayerPotentialSource(LayerPotentialSource):
return (area_element.with_queue(queue)*qweight).with_queue(None)
+ # {{{ helpers for symbolic operator processing
+
def preprocess_optemplate(self, name, discretizations, expr):
"""
:arg name: The symbolic name for *self*, which the preprocessor
@@ -421,6 +391,10 @@ class QBXLayerPotentialSource(LayerPotentialSource):
return result
+ # }}}
+
+ # {{{ internal functionality for execution
+
def exec_compute_potential_insn(self, queue, insn, bound_expr, evaluate):
from pytools.obj_array import with_object_array_or_scalar
from functools import partial
@@ -450,24 +424,6 @@ class QBXLayerPotentialSource(LayerPotentialSource):
return QBXFMMGeometryCodeGetter(self.cl_context, self.ambient_dim,
debug=self.debug)
- @memoize_method
- def qbx_fmm_geometry_data(self, target_discrs_and_qbx_sides):
- """
- :arg target_discrs_and_qbx_sides:
- a tuple of *(discr, qbx_forced_limit)*
- tuples, where *discr* is a
- :class:`meshmode.discretization.Discretization`
- or
- :class:`pytential.target.TargetBase`
- instance
- """
- from pytential.qbx.geometry import QBXFMMGeometryData
-
- return QBXFMMGeometryData(self.qbx_fmm_code_getter,
- self, target_discrs_and_qbx_sides,
- target_stick_out_factor=self.target_stick_out_factor,
- debug=self.debug)
-
# {{{ fmm-based execution
@memoize_method
@@ -591,7 +547,7 @@ class QBXLayerPotentialSource(LayerPotentialSource):
# }}}
- if len(geo_data.global_qbx_centers()) != geo_data.center_info().ncenters:
+ if len(geo_data.global_qbx_centers()) != geo_data.ncenters:
raise NotImplementedError("geometry has centers requiring local QBX")
from pytential.qbx.geometry import target_state
@@ -709,6 +665,8 @@ class QBXLayerPotentialSource(LayerPotentialSource):
strengths = (evaluate(insn.density).with_queue(queue)
* self.weights_and_area_elements())
+ import pytential.qbx.utils as utils
+
# FIXME: Do this all at once
result = []
for o in insn.outputs:
@@ -724,7 +682,7 @@ class QBXLayerPotentialSource(LayerPotentialSource):
evt, output_for_each_kernel = lpot_applier(
queue, target_discr.nodes(),
self.fine_density_discr.nodes(),
- self.centers(o.qbx_forced_limit),
+ utils.get_centers_on_side(self, o.qbx_forced_limit),
[strengths], **kernel_args)
result.append((o.name, output_for_each_kernel[o.kernel_index]))
else:
@@ -748,12 +706,11 @@ class QBXLayerPotentialSource(LayerPotentialSource):
(target_discr, qbx_forced_limit),
))
- # center_info is independent of targets
- center_info = geo_data.center_info()
+ # center-related info is independent of targets
# First ncenters targets are the centers
tgt_to_qbx_center = (
- geo_data.user_target_to_center()[center_info.ncenters:]
+ geo_data.user_target_to_center()[geo_data.ncenters:]
.copy(queue=queue))
qbx_tgt_numberer = self.get_qbx_target_numberer(
@@ -786,7 +743,7 @@ class QBXLayerPotentialSource(LayerPotentialSource):
queue,
targets=target_discr.nodes(),
sources=self.fine_density_discr.nodes(),
- centers=center_info.centers,
+ centers=geo_data.centers(),
strengths=[strengths],
qbx_tgt_numbers=qbx_tgt_numbers,
qbx_center_numbers=qbx_center_numbers,
@@ -798,6 +755,8 @@ class QBXLayerPotentialSource(LayerPotentialSource):
# }}}
+ # }}}
+
# }}}
diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py
index e3508b7089d2875c12a77079e763201300b5a30a..426de45d3789a02fd392fe5fadc92ab286e72b33 100644
--- a/pytential/qbx/fmm.py
+++ b/pytential/qbx/fmm.py
@@ -153,7 +153,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
return cl.array.zeros(
self.queue,
- (self.geo_data.center_info().ncenters,
+ (self.geo_data.ncenters,
len(qbx_l_expn)),
dtype=self.dtype)
@@ -212,7 +212,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
self.queue,
global_qbx_centers=geo_data.global_qbx_centers(),
qbx_center_to_target_box=geo_data.qbx_center_to_target_box(),
- qbx_centers=geo_data.center_info().centers,
+ qbx_centers=geo_data.centers(),
source_box_starts=starts,
source_box_lists=lists,
@@ -230,7 +230,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
qbx_expansions = self.qbx_local_expansion_zeros()
geo_data = self.geo_data
- if geo_data.center_info().ncenters == 0:
+ if geo_data.ncenters == 0:
return qbx_expansions
traversal = geo_data.traversal()
@@ -249,7 +249,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
qbx_center_to_target_box=geo_data.qbx_center_to_target_box(),
centers=self.tree.box_centers,
- qbx_centers=geo_data.center_info().centers,
+ qbx_centers=geo_data.centers(),
src_expansions=source_mpoles_view,
src_base_ibox=source_level_start_ibox,
@@ -273,7 +273,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
qbx_expansions = self.qbx_local_expansion_zeros()
geo_data = self.geo_data
- if geo_data.center_info().ncenters == 0:
+ if geo_data.ncenters == 0:
return qbx_expansions
trav = geo_data.traversal()
@@ -294,7 +294,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
target_base_ibox=target_level_start_ibox,
centers=self.tree.box_centers,
- qbx_centers=geo_data.center_info().centers,
+ qbx_centers=geo_data.centers(),
expansions=target_locals_view,
qbx_expansions=qbx_expansions,
@@ -322,7 +322,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
qbxl2p = self.code.qbxl2p(self.qbx_order)
evt, pot_res = qbxl2p(self.queue,
- qbx_centers=geo_data.center_info().centers,
+ qbx_centers=geo_data.centers(),
global_qbx_centers=geo_data.global_qbx_centers(),
center_to_targets_starts=ctt.starts,
@@ -369,7 +369,7 @@ def drive_fmm(expansion_wrangler, src_weights):
# Interface guidelines: Attributes of the tree are assumed to be known
# to the expansion wrangler and should not be passed.
- logger.debug("start qbx fmm")
+ logger.info("start qbx fmm")
logger.debug("reorder source weights")
src_weights = wrangler.reorder_sources(src_weights)
@@ -522,7 +522,7 @@ def drive_fmm(expansion_wrangler, src_weights):
# }}}
- logger.debug("qbx fmm complete")
+ logger.info("qbx fmm complete")
return result
@@ -652,7 +652,7 @@ def write_performance_model(outf, geo_data):
qbx_center_to_target_box = geo_data.qbx_center_to_target_box()
center_to_targets_starts = geo_data.center_to_tree_targets().starts
- ncenters = geo_data.center_info().ncenters
+ ncenters = geo_data.ncenters
outf.write("ncenters = {cost}\n"
.format(cost=ncenters))
diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py
index ce93e851d61741af3951096d3f8e131c76b221fc..f2669274bcd018cbc33f5bda5dbd338b173c0446 100644
--- a/pytential/qbx/geometry.py
+++ b/pytential/qbx/geometry.py
@@ -118,72 +118,6 @@ class QBXFMMGeometryCodeGetter(object):
self.ambient_dim = ambient_dim
self.debug = debug
- @property
- @memoize_method
- def pick_expansion_centers(self):
- knl = lp.make_kernel(
- """{[dim,k,i]:
- 0<=dim<ndims and
- 0<=k<nelements and
- 0<=i<nout_nodes}""",
- """
- centers[dim, k, i] = all_centers[dim, k, kept_center_indices[i]]
- radii[k, i] = all_radii[k, kept_center_indices[i]]
- """,
- [
- lp.GlobalArg("all_centers", None,
- shape="ndims,nelements,nunit_nodes"),
- lp.GlobalArg("all_radii", None, shape="nelements,nunit_nodes"),
- lp.ValueArg("nunit_nodes", np.int32),
- "..."
- ],
- default_offset=lp.auto,
- name="center_pick")
-
- knl = lp.fix_parameters(knl, ndims=self.ambient_dim)
-
- knl = lp.tag_array_axes(knl, "centers,all_centers", "sep, C, C")
-
- knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
- return lp.tag_inames(knl, dict(k="g.0", dim="ilp"))
-
- @property
- @memoize_method
- def find_element_centers(self):
- knl = lp.make_kernel(
- """{[dim,k,i]:
- 0<=dim<ndims and
- 0<=k<nelements and
- 0<=i<nunit_nodes}""",
- """
- el_centers[dim, k] = sum(i, nodes[dim, k, i])/nunit_nodes
- """,
- default_offset=lp.auto, name="find_element_centers")
-
- knl = lp.fix_parameters(knl, ndims=self.ambient_dim)
-
- knl = lp.split_iname(knl, "k", 128, inner_tag="l.0", outer_tag="g.0")
- return lp.tag_inames(knl, dict(dim="ilp"))
-
- @property
- @memoize_method
- def find_element_radii(self):
- knl = lp.make_kernel(
- """{[dim,k,i]:
- 0<=dim<ndims and
- 0<=k<nelements and
- 0<=i<nunit_nodes}""",
- """
- el_radii[k] = max(dim, max(i, \
- fabs(nodes[dim, k, i] - el_centers[dim, k])))
- """,
- default_offset=lp.auto, name="find_element_radii")
-
- knl = lp.fix_parameters(knl, ndims=self.ambient_dim)
-
- knl = lp.split_iname(knl, "k", 128, inner_tag="l.0", outer_tag="g.0")
- return lp.tag_inames(knl, dict(dim="unr"))
-
@memoize_method
def copy_targets_kernel(self):
knl = lp.make_kernel(
@@ -236,8 +170,10 @@ class QBXFMMGeometryCodeGetter(object):
# This write is race-free because each center only belongs
# to one box.
- qbx_center_to_target_box[itarget_user] = \
- box_to_target_box[ibox] {id=tgt_write,if=in_bounds}
+ if in_bounds
+ qbx_center_to_target_box[itarget_user] = \
+ box_to_target_box[ibox] {id=tgt_write}
+ end
end
end
""",
@@ -263,104 +199,6 @@ class QBXFMMGeometryCodeGetter(object):
from boxtree.area_query import LeavesToBallsLookupBuilder
return LeavesToBallsLookupBuilder(self.cl_context)
- # {{{ check if a center may be used with global QBX
-
- @memoize_method
- def qbx_center_for_global_tester(self,
- coord_dtype, box_id_dtype, particle_id_dtype):
- from pyopencl.elementwise import ElementwiseTemplate
- return ElementwiseTemplate(
- arguments="""//CL:mako//
- /* input */
- %for iaxis in range(ambient_dim):
- coord_t *center_${iaxis}, /* [ncenters] */
- %endfor
- coord_t *radii, /* [ncenters] */
- box_id_t *qbx_center_to_target_box, /* [ncenters] */
-
- box_id_t *neighbor_source_boxes_starts,
- box_id_t *neighbor_source_boxes_lists,
- particle_id_t *box_point_source_starts,
- particle_id_t *box_point_source_counts_cumul,
- %for iaxis in range(ambient_dim):
- coord_t *point_sources_${iaxis},
- %endfor
-
- /* output */
- char *center_may_use_global_qbx,
- """,
- operation=r"""//CL:mako//
- particle_id_t icenter = i;
-
- %for iaxis in range(ambient_dim):
- coord_t my_center_${iaxis} = center_${iaxis}[icenter];
- %endfor
- coord_t radius = radii[icenter];
-
- // {{{ see if there are sources close enough to require local QBX
-
- bool found_too_close = false;
-
- coord_t radius_squared = radius * radius;
-
- box_id_t itgt_box = qbx_center_to_target_box[icenter];
-
- box_id_t nb_src_start = neighbor_source_boxes_starts[itgt_box];
- box_id_t nb_src_stop = neighbor_source_boxes_starts[itgt_box+1];
-
- for (
- box_id_t ibox_list = nb_src_start;
- ibox_list < nb_src_stop && !found_too_close;
- ++ibox_list)
- {
- box_id_t neighbor_box_id =
- neighbor_source_boxes_lists[ibox_list];
-
- box_id_t bps_start = box_point_source_starts[neighbor_box_id];
- box_id_t bps_stop =
- bps_start + box_point_source_counts_cumul[neighbor_box_id];
- for (
- box_id_t ipsrc = bps_start;
- ipsrc < bps_stop;
- ++ipsrc)
- {
- %for iaxis in range(ambient_dim):
- coord_t psource_${iaxis} = point_sources_${iaxis}[ipsrc];
- %endfor
-
- coord_t dist_squared = 0
- %for iaxis in range(ambient_dim):
- + (my_center_${iaxis} - psource_${iaxis})
- * (my_center_${iaxis} - psource_${iaxis})
- %endfor
- ;
-
- const coord_t too_close_slack = (1+1e-2)*(1+1e-2);
- found_too_close = found_too_close
- || (dist_squared*too_close_slack < radius_squared);
-
- if (found_too_close)
- break;
- }
- }
-
- // }}}
-
- center_may_use_global_qbx[icenter] = !found_too_close;
- """,
- name="qbx_test_center_for_global").build(
- self.cl_context,
- type_aliases=(
- ("box_id_t", box_id_dtype),
- ("particle_id_t", particle_id_dtype),
- ("coord_t", coord_dtype),
- ),
- var_values=(
- ("ambient_dim", self.ambient_dim),
- ))
-
- # }}}
-
@property
@memoize_method
def key_value_sort(self):
@@ -379,7 +217,10 @@ class QBXFMMGeometryCodeGetter(object):
VectorArg(particle_id_dtype, "filtered_target_id"),
VectorArg(particle_id_dtype, "count"),
],
+
+ # "Does this target have a QBX center?"
input_expr="(target_to_center[i] >= 0) ? 1 : 0",
+
scan_expr="a+b", neutral="0",
output_statement="""
if (prev_item != item)
@@ -419,32 +260,6 @@ class TargetInfo(DeviceDataRecord):
"""
-class CenterInfo(DeviceDataRecord):
- """Information on location of QBX centers.
- Returned from :meth:`QBXFMMGeometryData.center_info`.
-
- .. attribute:: centers
-
- Shape: ``[dim][ncenters]``
-
- .. attribute:: sides
-
- Shape: ``[ncenters]``
-
- -1 for inside, +1 for outside (relative to normal)
-
- .. attribute:: radii
-
- Shape: ``[ncenters]``
-
- .. attribute:: ncenters
- """
-
- @property
- def ncenters(self):
- return len(self.radii)
-
-
class CenterToTargetList(DeviceDataRecord):
"""A lookup table of targets covered by each QBX disk. Indexed by global
number of QBX center, ``lists[start[i]:start[i+1]]`` indicates numbers
@@ -516,19 +331,22 @@ class QBXFMMGeometryData(object):
.. attribute:: coord_dtype
+ .. rubric :: Expansion centers
+
+ .. attribute:: ncenters
+ .. automethod:: centers()
+ .. automethod:: radii()
+
.. rubric :: Methods
- .. automethod:: center_info()
.. automethod:: target_info()
.. automethod:: tree()
.. automethod:: traversal()
- .. automethod:: leaf_to_center_lookup
.. automethod:: qbx_center_to_target_box()
.. automethod:: global_qbx_flags()
.. automethod:: global_qbx_centers()
.. automethod:: user_target_to_center()
.. automethod:: center_to_tree_targets()
- .. automethod:: global_qbx_centers_box_target_lists()
.. automethod:: non_qbx_box_target_lists()
.. automethod:: plot()
"""
@@ -561,35 +379,36 @@ class QBXFMMGeometryData(object):
def coord_dtype(self):
return self.lpot_source.fine_density_discr.nodes().dtype
- # {{{ center info
+ # {{{ centers/radii
- @memoize_method
- def kept_center_indices(self, el_group):
- """Return indices of unit nodes (of the target discretization)
- that will give rise to QBX centers.
- """
-
- # FIXME: Be more careful about which nodes to keep
- return np.arange(0, el_group.nunit_nodes)
+ @property
+ def ncenters(self):
+ return len(self.centers()[0])
@memoize_method
- def center_info(self):
- """ Return a :class:`CenterInfo`. |cached|
- """
+ def centers(self):
+ """ Return an object array of (interleaved) center coordinates.
- lpot_source = self.lpot_source
+ ``coord_t [ambient_dim][ncenters]``
+ """
with cl.CommandQueue(self.cl_context) as queue:
from pytential.qbx.utils import get_interleaved_centers
- centers = get_interleaved_centers(queue, lpot_source)
- sides = cl.array.arange(queue, len(centers[0]), dtype=np.int32)
- sides = 2 * (sides & 1) - 1
- radii = lpot_source.panel_sizes("ncenters").with_queue(queue) / 2
+ from pytools.obj_array import make_obj_array
+ return make_obj_array([
+ ccomp.with_queue(None)
+ for ccomp in get_interleaved_centers(queue, self.lpot_source)])
- return CenterInfo(
- sides=sides,
- radii=radii,
- centers=centers).with_queue(None)
+ @memoize_method
+ def expansion_radii(self):
+ """Return an array of radii associated with the (interleaved)
+ expansion centers.
+
+ ``coord_t [ncenters]``
+ """
+ with cl.CommandQueue(self.cl_context) as queue:
+ from pytential.qbx.utils import get_interleaved_radii
+ return get_interleaved_radii(queue, self.lpot_source)
# }}}
@@ -603,9 +422,7 @@ class QBXFMMGeometryData(object):
lpot_src = self.lpot_source
with cl.CommandQueue(self.cl_context) as queue:
- center_info = self.center_info()
-
- ntargets = center_info.ncenters
+ ntargets = self.ncenters
target_discr_starts = []
for target_discr, qbx_side in self.target_discrs_and_qbx_sides:
@@ -619,8 +436,8 @@ class QBXFMMGeometryData(object):
self.coord_dtype)
code_getter.copy_targets_kernel()(
queue,
- targets=targets[:, :center_info.ncenters],
- points=center_info.centers)
+ targets=targets[:, :self.ncenters],
+ points=self.centers())
for start, (target_discr, _) in zip(
target_discr_starts, self.target_discrs_and_qbx_sides):
@@ -634,23 +451,6 @@ class QBXFMMGeometryData(object):
target_discr_starts=target_discr_starts,
ntargets=ntargets).with_queue(None)
- def target_radii(self):
- """Shape: ``[ntargets]``
-
- A list of target radii for FMM tree construction. In this list, the QBX
- centers have the radii determined by :meth:`center_info`, and all other
- targets have radius zero.
- """
-
- tgt_info = self.target_info()
- center_info = self.center_info()
-
- with cl.CommandQueue(self.cl_context) as queue:
- target_radii = cl.array.zeros(queue, tgt_info.ntargets, self.coord_dtype)
- target_radii[:center_info.ncenters] = center_info.radii
-
- return target_radii.with_queue(None)
-
def target_side_preferences(self):
"""Return one big array combining all the data from
the *side* part of :attr:`TargetInfo.target_discrs_and_qbx_sides`.
@@ -658,12 +458,11 @@ class QBXFMMGeometryData(object):
Shape: ``[ntargets]``, dtype: int8"""
tgt_info = self.target_info()
- center_info = self.center_info()
with cl.CommandQueue(self.cl_context) as queue:
target_side_preferences = cl.array.empty(
queue, tgt_info.ntargets, np.int8)
- target_side_preferences[:center_info.ncenters] = 0
+ target_side_preferences[:self.ncenters] = 0
for tdstart, (target_discr, qbx_side) in \
zip(tgt_info.target_discr_starts,
@@ -679,7 +478,7 @@ class QBXFMMGeometryData(object):
@memoize_method
def tree(self):
- """Build and return a :class:`boxtree.tree.TreeWithLinkedPointSources`
+ """Build and return a :class:`boxtree.tree.Tree`
for this source with these targets.
|cached|
@@ -693,6 +492,12 @@ class QBXFMMGeometryData(object):
nsources = lpot_src.fine_density_discr.nnodes
nparticles = nsources + target_info.ntargets
+ target_radii = None
+ if self.lpot_source.expansion_disks_in_tree_have_extent:
+ target_radii = cl.array.zeros(queue, target_info.ntargets,
+ self.coord_dtype)
+ target_radii[:self.ncenters] = self.expansion_radii()
+
refine_weights = cl.array.zeros(queue, nparticles, dtype=np.int32)
refine_weights[:nsources] = 1
refine_weights.finish()
@@ -707,13 +512,22 @@ class QBXFMMGeometryData(object):
# to avoid having too many disks overlap more than one box.
#
# FIXME: Should investigate this further.
+
tree, _ = code_getter.build_tree(queue,
particles=lpot_src.fine_density_discr.nodes(),
targets=target_info.targets,
+ target_radii=target_radii,
max_leaf_refine_weight=256,
refine_weights=refine_weights,
debug=self.debug,
- kind="adaptive-level-restricted")
+ stick_out_factor=0,
+ kind="adaptive")
+
+ if self.debug:
+ tgt_count_2 = cl.array.sum(
+ tree.box_target_counts_nonchild, queue=queue).get()
+
+ assert (tree.ntargets == tgt_count_2), (tree.ntargets, tgt_count_2)
return tree
@@ -731,31 +545,21 @@ class QBXFMMGeometryData(object):
trav, _ = self.code_getter.build_traversal(queue, self.tree(),
debug=self.debug)
- return trav
-
- def leaf_to_center_lookup(self):
- """Return a :class:`boxtree.area_query.LeavesToBallsLookup` to look up
- which which QBX disks overlap each leaf box.
- """
-
- center_info = self.center_info()
+ if self.lpot_source.expansion_disks_in_tree_have_extent:
+ trav = trav.merge_close_lists(queue)
- with cl.CommandQueue(self.cl_context) as queue:
- lbl, _ = self.code_getter.build_leaf_to_ball_lookup(queue,
- self.tree(), center_info.centers, center_info.radii)
- return lbl
+ return trav
@memoize_method
def qbx_center_to_target_box(self):
- """Return a lookup table of length :attr:`CenterInfo.ncenters`
- (see :meth:`center_info`) indicating the target box in which each
+ """Return a lookup table of length :attr:`ncenters`
+ indicating the target box in which each
QBX disk is located.
|cached|
"""
tree = self.tree()
trav = self.traversal()
- center_info = self.center_info()
qbx_center_to_target_box_lookup = \
self.code_getter.qbx_center_to_target_box_lookup(
@@ -782,24 +586,36 @@ class QBXFMMGeometryData(object):
queue, len(sorted_target_ids),
user_target_from_tree_target.dtype)
- evt, (qbx_center_to_target_box,) = qbx_center_to_target_box_lookup(
+ qbx_center_to_target_box = cl.array.empty(
+ queue, self.ncenters, tree.box_id_dtype)
+
+ if self.lpot_source.debug:
+ qbx_center_to_target_box.fill(-1)
+
+ evt, _ = qbx_center_to_target_box_lookup(
queue,
+ qbx_center_to_target_box=qbx_center_to_target_box,
box_to_target_box=box_to_target_box,
box_target_starts=tree.box_target_starts,
box_target_counts_nonchild=tree.box_target_counts_nonchild,
user_target_from_tree_target=user_target_from_tree_target,
- ncenters=center_info.ncenters)
+ ncenters=self.ncenters)
+
+ if self.lpot_source.debug:
+ assert 0 <= cl.array.min(qbx_center_to_target_box).get()
+ assert (
+ cl.array.max(qbx_center_to_target_box).get()
+ < len(trav.target_boxes))
return qbx_center_to_target_box.with_queue(None)
@memoize_method
def global_qbx_flags(self):
"""Return an array of :class:`numpy.int8` of length
- :attr:`CenterInfo.ncenters` (see :meth:`center_info`) indicating
- whether each center can use gloal QBX, i.e. whether a single expansion
- can mediate interactions from *all* sources to all targets for which it
- is valid. If global QBX can be used, the center's entry will be 1,
- otherwise it will be 0.
+ :attr:`ncenters` indicating whether each center can use gloal QBX, i.e.
+ whether a single expansion can mediate interactions from *all* sources
+ to all targets for which it is valid. If global QBX can be used, the
+ center's entry will be 1, otherwise it will be 0.
(If not, local QBX is needed, and the center may only be
able to mediate some of the interactions to a given target.)
@@ -807,10 +623,8 @@ class QBXFMMGeometryData(object):
|cached|
"""
- center_info = self.center_info()
-
with cl.CommandQueue(self.cl_context) as queue:
- result = cl.array.empty(queue, center_info.ncenters, np.int8)
+ result = cl.array.empty(queue, self.ncenters, np.int8)
result.fill(1)
return result.with_queue(None)
@@ -829,7 +643,7 @@ class QBXFMMGeometryData(object):
logger.info("find global qbx centers: start")
result, count, _ = copy_if(
- cl.array.arange(queue, self.center_info().ncenters,
+ cl.array.arange(queue, self.ncenters,
tree.particle_id_dtype),
"global_qbx_flags[i] != 0",
extra_args=[
@@ -857,16 +671,15 @@ class QBXFMMGeometryData(object):
tgt_assoc = QBXTargetAssociator(self.cl_context)
tgt_info = self.target_info()
- center_info = self.center_info()
from pytential.target import PointsTarget
with cl.CommandQueue(self.cl_context) as queue:
target_side_prefs = (self
- .target_side_preferences()[center_info.ncenters:].get(queue=queue))
+ .target_side_preferences()[self.ncenters:].get(queue=queue))
target_discrs_and_qbx_sides = [(
- PointsTarget(tgt_info.targets[:, center_info.ncenters:]),
+ PointsTarget(tgt_info.targets[:, self.ncenters:]),
target_side_prefs.astype(np.int32))]
# FIXME: try block...
@@ -878,8 +691,8 @@ class QBXFMMGeometryData(object):
with cl.CommandQueue(self.cl_context) as queue:
result = cl.array.empty(queue, tgt_info.ntargets, tree.particle_id_dtype)
- result[:center_info.ncenters].fill(target_state.NO_QBX_NEEDED)
- result[center_info.ncenters:] = tgt_assoc_result.target_to_center
+ result[:self.ncenters].fill(target_state.NO_QBX_NEEDED)
+ result[self.ncenters:] = tgt_assoc_result.target_to_center
return result
@@ -891,7 +704,6 @@ class QBXFMMGeometryData(object):
|cached|
"""
- center_info = self.center_info()
user_ttc = self.user_target_to_center()
with cl.CommandQueue(self.cl_context) as queue:
@@ -917,7 +729,7 @@ class QBXFMMGeometryData(object):
center_target_starts, targets_sorted_by_center, _ = \
self.code_getter.key_value_sort(queue,
filtered_tree_ttc, filtered_target_ids,
- center_info.ncenters, tree_ttc.dtype)
+ self.ncenters, tree_ttc.dtype)
logger.info("build center -> targets lookup table: done")
@@ -927,31 +739,6 @@ class QBXFMMGeometryData(object):
return result
- @memoize_method
- def global_qbx_centers_box_target_lists(self):
- """Build a list of targets per box consisting only of global QBX centers.
- Returns a :class:`boxtree.tree.FilteredTargetListsInUserOrder`.
- (I.e. no new target order is created for these targets, as we expect
- there to be (relatively) not many of them.)
-
- |cached|
- """
-
- center_info = self.center_info()
- with cl.CommandQueue(self.cl_context) as queue:
- logger.info("find global qbx centers box target list: start")
-
- flags = cl.array.zeros(queue, self.tree().ntargets, np.int8)
-
- flags[:center_info.ncenters] = self.global_qbx_flags()
-
- from boxtree.tree import filter_target_lists_in_user_order
- result = filter_target_lists_in_user_order(queue, self.tree(), flags)
-
- logger.info("find global qbx centers box target list: done")
-
- return result.with_queue(None)
-
@memoize_method
def non_qbx_box_target_lists(self):
"""Build a list of targets per box that don't need to bother with QBX.
@@ -973,7 +760,7 @@ class QBXFMMGeometryData(object):
# 'flags' is in user order, and should be.
- nqbx_centers = self.center_info().ncenters
+ nqbx_centers = self.ncenters
flags[:nqbx_centers] = 0
from boxtree.tree import filter_target_lists_in_tree_order
@@ -985,10 +772,14 @@ class QBXFMMGeometryData(object):
# {{{ plotting (for debugging)
- def plot(self):
+ def plot(self, draw_circles=False, draw_center_numbers=False,
+ highlight_centers=None):
"""Plot most of the information contained in a :class:`QBXFMMGeometryData`
object, for debugging.
+ :arg highlight_centers: If not *None*, an object with which the array of
+ centers can be indexed to find the highlighted centers.
+
.. note::
This only works for two-dimensional geometries.
@@ -1000,8 +791,8 @@ class QBXFMMGeometryData(object):
with cl.CommandQueue(self.cl_context) as queue:
import matplotlib.pyplot as pt
- nodes_host = self.lpot_source.fine_density_discr.nodes().get(queue)
- pt.plot(nodes_host[0], nodes_host[1], "x-")
+ from meshmode.discretization.visualization import draw_curve
+ draw_curve(self.lpot_source.fine_density_discr)
global_flags = self.global_qbx_flags().get(queue=queue)
@@ -1012,22 +803,35 @@ class QBXFMMGeometryData(object):
# {{{ draw centers and circles
- center_info = self.center_info()
+ centers = self.centers()
centers = [
- center_info.centers[0].get(queue),
- center_info.centers[1].get(queue)]
+ centers[0].get(queue),
+ centers[1].get(queue)]
pt.plot(centers[0][global_flags == 0],
centers[1][global_flags == 0], "oc",
label="centers needing local qbx")
+
+ if highlight_centers is not None:
+ pt.plot(centers[0][highlight_centers],
+ centers[1][highlight_centers], "oc",
+ label="highlighted centers",
+ markersize=15)
+
ax = pt.gca()
- for icenter, (cx, cy, r) in enumerate(zip(
- centers[0], centers[1], center_info.radii.get(queue))):
- ax.add_artist(
- pt.Circle((cx, cy), r, fill=False, ls="dotted", lw=1))
- pt.text(cx, cy,
- str(icenter), fontsize=8,
- ha="left", va="center",
- bbox=dict(facecolor='white', alpha=0.5, lw=0))
+
+ if draw_circles:
+ for icenter, (cx, cy, r) in enumerate(zip(
+ centers[0], centers[1],
+ self.expansion_radii().get(queue))):
+ ax.add_artist(
+ pt.Circle((cx, cy), r, fill=False, ls="dotted", lw=1))
+
+ if draw_center_numbers:
+ for icenter, (cx, cy, r) in enumerate(zip(centers[0], centers[1])):
+ pt.text(cx, cy,
+ str(icenter), fontsize=8,
+ ha="left", va="center",
+ bbox=dict(facecolor='white', alpha=0.5, lw=0))
# }}}
@@ -1054,9 +858,9 @@ class QBXFMMGeometryData(object):
tccount = 0
checked = 0
for tx, ty, tcenter in zip(
- targets[0][center_info.ncenters:],
- targets[1][center_info.ncenters:],
- ttc[center_info.ncenters:]):
+ targets[0][self.ncenters:],
+ targets[1][self.ncenters:],
+ ttc[self.ncenters:]):
checked += 1
if tcenter >= 0:
tccount += 1
diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py
index c942b45b64b9a7eb275eea8f965de189b7c3e491..8efd9a8800f614f9d54a5d30285f108c42072ba1 100644
--- a/pytential/qbx/refinement.py
+++ b/pytential/qbx/refinement.py
@@ -3,7 +3,7 @@ from __future__ import division, absolute_import, print_function
__copyright__ = """
Copyright (C) 2013 Andreas Kloeckner
-Copyright (C) 2016 Matt Wala
+Copyright (C) 2016, 2017 Matt Wala
"""
__license__ = """
@@ -32,98 +32,56 @@ import numpy as np
import pyopencl as cl
from pytools import memoize_method
-from pytential.qbx.utils import DiscrPlotterMixin
from boxtree.area_query import AreaQueryElementwiseTemplate
-from pyopencl.elementwise import ElementwiseTemplate
from boxtree.tools import InlineBinarySearch
from pytential.qbx.utils import QBX_TREE_C_PREAMBLE, QBX_TREE_MAKO_DEFS
-unwrap_args = AreaQueryElementwiseTemplate.unwrap_args
-
import logging
logger = logging.getLogger(__name__)
-__doc__ = """
-Refinement
-^^^^^^^^^^
-
-.. autoclass:: QBXLayerPotentialSourceRefiner
-"""
-
-# {{{ kernels
-
-TUNNEL_QUERY_DISTANCE_FINDER_TEMPLATE = ElementwiseTemplate(
- arguments=r"""//CL:mako//
- /* input */
- particle_id_t source_offset,
- particle_id_t panel_offset,
- int npanels,
- particle_id_t *panel_to_source_starts,
- particle_id_t *sorted_target_ids,
- coord_t *panel_sizes,
-
- /* output */
- float *tunnel_query_dists,
- /* input, dim-dependent size */
- %for ax in AXIS_NAMES[:dimensions]:
- coord_t *particles_${ax},
- %endfor
- """,
- operation=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r"""//CL:mako//
- /* Find my panel. */
- particle_id_t panel = bsearch(panel_to_source_starts, npanels + 1, i);
+# max_levels granularity for the stack used by the tree descent code in the
+# area query kernel.
+MAX_LEVELS_INCREMENT = 10
- /* Compute dist(tunnel region, panel center) */
- coord_vec_t center_of_mass;
- ${load_particle("INDEX_FOR_PANEL_PARTICLE(panel)", "center_of_mass")}
+__doc__ = """
+The refiner takes a layer potential source and refines it until it satisfies
+three global QBX refinement criteria:
- coord_vec_t center;
- ${load_particle("INDEX_FOR_SOURCE_PARTICLE(i)", "center")}
+ * *Condition 1* (Expansion disk undisturbed by sources)
+ A center must be closest to its own source.
- coord_t panel_size = panel_sizes[panel];
+ * *Condition 2* (Sufficient quadrature sampling from all source panels)
+ The quadrature contribution from each panel is as accurate
+ as from the center's own source panel.
- coord_t max_dist = 0;
+ * *Condition 3* (Panel size bounded based on kernel length scale)
+ The panel size is bounded by a kernel length scale. This
+ applies only to Helmholtz kernels.
- %for ax in AXIS_NAMES[:dimensions]:
- {
- max_dist = fmax(max_dist,
- distance(center_of_mass.${ax}, center.${ax} + panel_size / 2));
- max_dist = fmax(max_dist,
- distance(center_of_mass.${ax}, center.${ax} - panel_size / 2));
- }
- %endfor
+.. autoclass:: RefinerCodeContainer
+.. autoclass:: RefinerWrangler
+.. autoclass:: RefinerNotConvergedWarning
- // The atomic max operation supports only integer types.
- // However, max_dist is of a floating point type.
- // For comparison purposes we reinterpret the bits of max_dist
- // as an integer. The comparison result is the same as for positive
- // IEEE floating point numbers, so long as the float/int endianness
- // matches (fingers crossed).
- atomic_max(
- (volatile __global int *)
- &tunnel_query_dists[panel],
- as_int((float) max_dist));
- """,
- name="find_tunnel_query_distance",
- preamble=str(InlineBinarySearch("particle_id_t")))
+.. autofunction:: make_empty_refine_flags
+.. autofunction:: refine_for_global_qbx
+"""
+# {{{ kernels
-# Implements "Algorithm for triggering refinement based on Condition 1"
-#
-# Does not use r_max as we do not use Newton for checking center-panel closeness.
-CENTER_IS_CLOSEST_TO_ORIG_PANEL_REFINER = AreaQueryElementwiseTemplate(
+# Refinement checker for Condition 1.
+EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER = AreaQueryElementwiseTemplate(
extra_args=r"""
/* input */
- particle_id_t *box_to_panel_starts,
- particle_id_t *box_to_panel_lists,
+ particle_id_t *box_to_source_starts,
+ particle_id_t *box_to_source_lists,
particle_id_t *panel_to_source_starts,
particle_id_t *panel_to_center_starts,
particle_id_t source_offset,
particle_id_t center_offset,
particle_id_t *sorted_target_ids,
- coord_t *panel_sizes,
+ coord_t *center_danger_zone_radii_by_panel,
int npanels,
/* output */
@@ -136,39 +94,35 @@ CENTER_IS_CLOSEST_TO_ORIG_PANEL_REFINER = AreaQueryElementwiseTemplate(
%endfor
""",
ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r"""
+ /* Find the panel associated with this center. */
particle_id_t my_panel = bsearch(panel_to_center_starts, npanels + 1, i);
${load_particle("INDEX_FOR_CENTER_PARTICLE(i)", ball_center)}
- ${ball_radius} = panel_sizes[my_panel] / 2;
+ ${ball_radius} = center_danger_zone_radii_by_panel[my_panel];
""",
leaf_found_op=QBX_TREE_MAKO_DEFS + r"""
- for (particle_id_t panel_idx = box_to_panel_starts[${leaf_box_id}];
- panel_idx < box_to_panel_starts[${leaf_box_id} + 1];
- ++panel_idx)
+ /* Check that each source in the leaf box is sufficiently far from the
+ center; if not, mark the panel for refinement. */
+
+ for (particle_id_t source_idx = box_to_source_starts[${leaf_box_id}];
+ source_idx < box_to_source_starts[${leaf_box_id} + 1];
+ ++source_idx)
{
- particle_id_t panel = box_to_panel_lists[panel_idx];
+ particle_id_t source = box_to_source_lists[source_idx];
+ particle_id_t panel = bsearch(
+ panel_to_source_starts, npanels + 1, source);
- // Skip self.
if (my_panel == panel)
{
continue;
}
- bool is_close = false;
-
- for (particle_id_t source = panel_to_source_starts[panel];
- source < panel_to_source_starts[panel + 1];
- ++source)
- {
- coord_vec_t source_coords;
+ coord_vec_t source_coords;
+ ${load_particle("INDEX_FOR_SOURCE_PARTICLE(source)", "source_coords")}
- ${load_particle(
- "INDEX_FOR_SOURCE_PARTICLE(source)", "source_coords")}
-
- is_close |= (
- distance(${ball_center}, source_coords)
- <= panel_sizes[my_panel] / 2);
- }
+ bool is_close = (
+ distance(${ball_center}, source_coords)
+ <= center_danger_zone_radii_by_panel[my_panel]);
if (is_close)
{
@@ -178,27 +132,22 @@ CENTER_IS_CLOSEST_TO_ORIG_PANEL_REFINER = AreaQueryElementwiseTemplate(
}
}
""",
- name="refine_center_closest_to_orig_panel",
+ name="check_center_closest_to_orig_panel",
preamble=str(InlineBinarySearch("particle_id_t")))
-# Implements "Algorithm for triggering refinement based on Condition 2"
-CENTER_IS_FAR_FROM_NONNEIGHBOR_PANEL_REFINER = AreaQueryElementwiseTemplate(
+# Refinement checker for Condition 2.
+SUFFICIENT_SOURCE_QUADRATURE_RESOLUTION_CHECKER = AreaQueryElementwiseTemplate(
extra_args=r"""
/* input */
- particle_id_t *box_to_panel_starts,
- particle_id_t *box_to_panel_lists,
+ particle_id_t *box_to_center_starts,
+ particle_id_t *box_to_center_lists,
particle_id_t *panel_to_source_starts,
- particle_id_t *panel_to_center_starts,
particle_id_t source_offset,
particle_id_t center_offset,
- particle_id_t panel_offset,
particle_id_t *sorted_target_ids,
- coord_t *panel_sizes,
- particle_id_t *panel_adjacency_starts,
- particle_id_t *panel_adjacency_lists,
+ coord_t *source_danger_zone_radii_by_panel,
int npanels,
- coord_t *tunnel_query_dists,
/* output */
int *panel_refine_flags,
@@ -210,51 +159,29 @@ CENTER_IS_FAR_FROM_NONNEIGHBOR_PANEL_REFINER = AreaQueryElementwiseTemplate(
%endfor
""",
ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r"""
- particle_id_t my_panel = bsearch(panel_to_center_starts, npanels + 1, i);
- coord_vec_t my_center_coords;
+ /* Find the panel associated with this source. */
+ particle_id_t my_panel = bsearch(panel_to_source_starts, npanels + 1, i);
- ${load_particle("INDEX_FOR_CENTER_PARTICLE(i)", "my_center_coords")}
- ${load_particle("INDEX_FOR_PANEL_PARTICLE(my_panel)", ball_center)}
- ${ball_radius} = tunnel_query_dists[my_panel];
+ ${load_particle("INDEX_FOR_SOURCE_PARTICLE(i)", ball_center)}
+ ${ball_radius} = source_danger_zone_radii_by_panel[my_panel];
""",
leaf_found_op=QBX_TREE_MAKO_DEFS + r"""
- for (particle_id_t panel_idx = box_to_panel_starts[${leaf_box_id}];
- panel_idx < box_to_panel_starts[${leaf_box_id} + 1];
- ++panel_idx)
- {
- particle_id_t panel = box_to_panel_lists[panel_idx];
-
- bool is_self_or_adjacent = (my_panel == panel);
+ /* Check that each center in the leaf box is sufficiently far from the
+ panel; if not, mark the panel for refinement. */
- for (particle_id_t adj_panel_idx = panel_adjacency_starts[my_panel];
- adj_panel_idx < panel_adjacency_starts[my_panel + 1];
- ++adj_panel_idx)
- {
- is_self_or_adjacent |= (
- panel_adjacency_lists[adj_panel_idx] == panel);
- }
-
- // Skip self and adjacent panels.
- if (is_self_or_adjacent)
- {
- continue;
- }
-
- bool is_close = false;
-
- for (particle_id_t source = panel_to_source_starts[panel];
- source < panel_to_source_starts[panel + 1];
- ++source)
- {
- coord_vec_t source_coords;
+ for (particle_id_t center_idx = box_to_center_starts[${leaf_box_id}];
+ center_idx < box_to_center_starts[${leaf_box_id} + 1];
+ ++center_idx)
+ {
+ particle_id_t center = box_to_center_lists[center_idx];
- ${load_particle(
- "INDEX_FOR_SOURCE_PARTICLE(source)", "source_coords")}
+ coord_vec_t center_coords;
+ ${load_particle(
+ "INDEX_FOR_CENTER_PARTICLE(center)", "center_coords")}
- is_close |= (
- distance(my_center_coords, source_coords)
- <= panel_sizes[panel] / 2);
- }
+ bool is_close = (
+ distance(${ball_center}, center_coords)
+ <= source_danger_zone_radii_by_panel[my_panel]);
if (is_close)
{
@@ -264,123 +191,46 @@ CENTER_IS_FAR_FROM_NONNEIGHBOR_PANEL_REFINER = AreaQueryElementwiseTemplate(
}
}
""",
- name="refine_center_far_from_nonneighbor_panels",
+ name="check_source_quadrature_resolution",
preamble=str(InlineBinarySearch("particle_id_t")))
# }}}
-# {{{ lpot source refiner
+# {{{ code container
-class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin):
- """
- Driver for refining the QBX source grid. Follows [1]_.
+class RefinerCodeContainer(object):
- .. [1] Rachh, Manas, Andreas Klöckner, and Michael O'Neil. "Fast
- algorithms for Quadrature by Expansion I: Globally valid expansions."
-
- .. automethod:: get_refine_flags
- .. automethod:: __call__
- """
-
- def __init__(self, context):
- self.cl_context = context
- from pytential.qbx.utils import TreeWithQBXMetadataBuilder
- self.tree_builder = TreeWithQBXMetadataBuilder(self.cl_context)
- from boxtree.area_query import PeerListFinder
- self.peer_list_finder = PeerListFinder(self.cl_context)
-
- # {{{ kernels
+ def __init__(self, cl_context):
+ self.cl_context = cl_context
@memoize_method
- def get_tunnel_query_distance_finder(self, dimensions, coord_dtype,
- particle_id_dtype):
- from pyopencl.tools import dtype_to_ctype
- from boxtree.tools import AXIS_NAMES
- logger.info("refiner: building tunnel query distance finder kernel")
-
- knl = TUNNEL_QUERY_DISTANCE_FINDER_TEMPLATE.build(
+ def expansion_disk_undisturbed_by_sources_checker(
+ self, dimensions, coord_dtype, box_id_dtype, peer_list_idx_dtype,
+ particle_id_dtype, max_levels):
+ return EXPANSION_DISK_UNDISTURBED_BY_SOURCES_CHECKER.generate(
self.cl_context,
- type_aliases=(
- ("particle_id_t", particle_id_dtype),
- ("coord_t", coord_dtype),
- ),
- var_values=(
- ("dimensions", dimensions),
- ("AXIS_NAMES", AXIS_NAMES),
- ("coord_dtype", coord_dtype),
- ("dtype_to_ctype", dtype_to_ctype),
- ("vec_types", tuple(cl.array.vec.types.items())),
- ))
-
- logger.info("refiner: done building tunnel query distance finder kernel")
- return knl
-
- @memoize_method
- def get_center_is_closest_to_orig_panel_refiner(self, dimensions,
- coord_dtype, box_id_dtype,
- peer_list_idx_dtype,
- particle_id_dtype,
- max_levels):
- return CENTER_IS_CLOSEST_TO_ORIG_PANEL_REFINER.generate(self.cl_context,
dimensions, coord_dtype, box_id_dtype, peer_list_idx_dtype,
max_levels,
extra_type_aliases=(("particle_id_t", particle_id_dtype),))
@memoize_method
- def get_center_is_far_from_nonneighbor_panel_refiner(self, dimensions,
- coord_dtype,
- box_id_dtype,
- peer_list_idx_dtype,
- particle_id_dtype,
- max_levels):
- return CENTER_IS_FAR_FROM_NONNEIGHBOR_PANEL_REFINER.generate(self.cl_context,
+ def sufficient_source_quadrature_resolution_checker(
+ self, dimensions, coord_dtype, box_id_dtype, peer_list_idx_dtype,
+ particle_id_dtype, max_levels):
+ return SUFFICIENT_SOURCE_QUADRATURE_RESOLUTION_CHECKER.generate(
+ self.cl_context,
dimensions, coord_dtype, box_id_dtype, peer_list_idx_dtype,
max_levels,
extra_type_aliases=(("particle_id_t", particle_id_dtype),))
@memoize_method
- def get_2_to_1_panel_ratio_refiner(self):
- knl = lp.make_kernel([
- "{[panel]: 0<=panel<npanels}",
- "{[ineighbor]: neighbor_start<=ineighbor<neighbor_stop}"
- ],
- """
- for panel
- <> neighbor_start = panel_adjacency_starts[panel]
- <> neighbor_stop = panel_adjacency_starts[panel + 1]
- for ineighbor
- <> neighbor = panel_adjacency_lists[ineighbor]
-
- <> oversize = (refine_flags_prev[panel] == 0
- and (
- (panel_sizes[panel] > 2 * panel_sizes[neighbor]) or
- (panel_sizes[panel] > panel_sizes[neighbor] and
- refine_flags_prev[neighbor] == 1)))
-
- if oversize
- refine_flags[panel] = 1
- refine_flags_updated = 1 {id=write_refine_flags_updated}
- end
- end
- end
- """, [
- lp.GlobalArg("panel_adjacency_lists", shape=None),
- "..."
- ],
- options="return_dict",
- silenced_warnings="write_race(write_refine_flags_updated)",
- name="refine_2_to_1_adj_panel_size_ratio")
- knl = lp.split_iname(knl, "panel", 128, inner_tag="l.0", outer_tag="g.0")
- return knl
-
- @memoize_method
- def get_helmholtz_k_to_panel_ratio_refiner(self):
+ def kernel_length_scale_to_panel_size_ratio_checker(self):
knl = lp.make_kernel(
"{[panel]: 0<=panel<npanels}",
"""
for panel
- <> oversize = panel_sizes[panel] * helmholtz_k > 5
+ <> oversize = panel_sizes[panel] > kernel_length_scale
if oversize
refine_flags[panel] = 1
refine_flags_updated = 1 {id=write_refine_flags_updated}
@@ -389,21 +239,49 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin):
""",
options="return_dict",
silenced_warnings="write_race(write_refine_flags_updated)",
- name="refine_helmholtz_k_to_panel_size_ratio")
+ name="refine_kernel_length_scale_to_panel_size_ratio")
knl = lp.split_iname(knl, "panel", 128, inner_tag="l.0", outer_tag="g.0")
return knl
- # }}}
+ @memoize_method
+ def tree_builder(self):
+ from boxtree.tree_build import TreeBuilder
+ return TreeBuilder(self.cl_context)
+
+ @memoize_method
+ def peer_list_finder(self):
+ from boxtree.area_query import PeerListFinder
+ return PeerListFinder(self.cl_context)
+
+ def get_wrangler(self, queue):
+ """
+ :arg queue:
+ """
+ return RefinerWrangler(self, queue)
+
+# }}}
+
+
+# {{{ wrangler
- # {{{ refinement triggering
+class RefinerWrangler(object):
+
+ def __init__(self, code_container, queue):
+ self.code_container = code_container
+ self.queue = queue
+
+ # {{{ check subroutines for conditions 1-3
+
+ def check_expansion_disks_undisturbed_by_sources(self,
+ lpot_source, tree, peer_lists, refine_flags,
+ debug, wait_for=None):
- def refinement_check_center_is_closest_to_orig_panel(self, queue, tree,
- lpot_source, peer_lists, refine_flags, debug, wait_for=None):
# Avoid generating too many kernels.
from pytools import div_ceil
- max_levels = 10 * div_ceil(tree.nlevels, 10)
+ max_levels = MAX_LEVELS_INCREMENT * div_ceil(
+ tree.nlevels, MAX_LEVELS_INCREMENT)
- knl = self.get_center_is_closest_to_orig_panel_refiner(
+ knl = self.code_container.expansion_disk_undisturbed_by_sources_checker(
tree.dimensions,
tree.coord_dtype, tree.box_id_dtype,
peer_lists.peer_list_starts.dtype,
@@ -415,26 +293,32 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin):
if debug:
npanels_to_refine_prev = cl.array.sum(refine_flags).get()
- found_panel_to_refine = cl.array.zeros(queue, 1, np.int32)
+ found_panel_to_refine = cl.array.zeros(self.queue, 1, np.int32)
found_panel_to_refine.finish()
+ unwrap_args = AreaQueryElementwiseTemplate.unwrap_args
+
+ # FIXME: This shouldn't be by panel
+ center_danger_zone_radii_by_panel = \
+ lpot_source._expansion_radii("npanels")
evt = knl(
*unwrap_args(
tree, peer_lists,
- tree.box_to_qbx_panel_starts,
- tree.box_to_qbx_panel_lists,
+ tree.box_to_qbx_source_starts,
+ tree.box_to_qbx_source_lists,
tree.qbx_panel_to_source_starts,
tree.qbx_panel_to_center_starts,
tree.qbx_user_source_slice.start,
tree.qbx_user_center_slice.start,
tree.sorted_target_ids,
- lpot_source.panel_sizes("npanels"),
+ center_danger_zone_radii_by_panel,
tree.nqbxpanels,
refine_flags,
found_panel_to_refine,
*tree.sources),
range=slice(tree.nqbxcenters),
- queue=queue)
+ queue=self.queue,
+ wait_for=wait_for)
cl.wait_for_events([evt])
@@ -448,51 +332,51 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin):
return found_panel_to_refine.get()[0] == 1
- def refinement_check_center_is_far_from_nonneighbor_panels(self, queue,
- tree, lpot_source, peer_lists, tq_dists, refine_flags, debug,
- wait_for=None):
+ def check_sufficient_source_quadrature_resolution(
+ self, lpot_source, tree, peer_lists, refine_flags, debug,
+ wait_for=None):
+
# Avoid generating too many kernels.
from pytools import div_ceil
- max_levels = 10 * div_ceil(tree.nlevels, 10)
+ max_levels = MAX_LEVELS_INCREMENT * div_ceil(
+ tree.nlevels, MAX_LEVELS_INCREMENT)
- knl = self.get_center_is_far_from_nonneighbor_panel_refiner(
+ knl = self.code_container.sufficient_source_quadrature_resolution_checker(
tree.dimensions,
tree.coord_dtype, tree.box_id_dtype,
peer_lists.peer_list_starts.dtype,
tree.particle_id_dtype,
max_levels)
-
- logger.info("refiner: checking center is far from nonneighbor panels")
-
if debug:
npanels_to_refine_prev = cl.array.sum(refine_flags).get()
- found_panel_to_refine = cl.array.zeros(queue, 1, np.int32)
+ logger.info("refiner: checking sufficient quadrature resolution")
+
+ found_panel_to_refine = cl.array.zeros(self.queue, 1, np.int32)
found_panel_to_refine.finish()
- adjacency = self.get_adjacency_on_device(queue, lpot_source)
+ source_danger_zone_radii_by_panel = (
+ lpot_source._fine_panel_sizes("npanels")
+ .with_queue(self.queue) / 4)
+ unwrap_args = AreaQueryElementwiseTemplate.unwrap_args
evt = knl(
*unwrap_args(
tree, peer_lists,
- tree.box_to_qbx_panel_starts,
- tree.box_to_qbx_panel_lists,
+ tree.box_to_qbx_center_starts,
+ tree.box_to_qbx_center_lists,
tree.qbx_panel_to_source_starts,
- tree.qbx_panel_to_center_starts,
tree.qbx_user_source_slice.start,
tree.qbx_user_center_slice.start,
- tree.qbx_user_panel_slice.start,
tree.sorted_target_ids,
- lpot_source.panel_sizes("npanels"),
- adjacency.adjacency_starts,
- adjacency.adjacency_lists,
+ source_danger_zone_radii_by_panel,
tree.nqbxpanels,
- tq_dists,
refine_flags,
found_panel_to_refine,
*tree.sources),
- range=slice(tree.nqbxcenters),
- queue=queue)
+ range=slice(tree.nqbxsources),
+ queue=self.queue,
+ wait_for=wait_for)
cl.wait_for_events([evt])
@@ -502,24 +386,24 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin):
logger.debug("refiner: found {} panel(s) to refine".format(
npanels_to_refine - npanels_to_refine_prev))
- logger.info("refiner: done checking center is far from nonneighbor panels")
+ logger.info("refiner: done checking sufficient quadrature resolution")
return found_panel_to_refine.get()[0] == 1
- def refinement_check_helmholtz_k_to_panel_size_ratio(self, queue, lpot_source,
- helmholtz_k, refine_flags, debug, wait_for=None):
- knl = self.get_helmholtz_k_to_panel_ratio_refiner()
+ def check_kernel_length_scale_to_panel_size_ratio(self, lpot_source,
+ kernel_length_scale, refine_flags, debug, wait_for=None):
+ knl = self.code_container.kernel_length_scale_to_panel_size_ratio_checker()
- logger.info("refiner: checking helmholtz k to panel size ratio")
+ logger.info("refiner: checking kernel length scale to panel size ratio")
if debug:
npanels_to_refine_prev = cl.array.sum(refine_flags).get()
- evt, out = knl(queue,
- panel_sizes=lpot_source.panel_sizes("npanels"),
+ evt, out = knl(self.queue,
+ panel_sizes=lpot_source._panel_sizes("npanels"),
refine_flags=refine_flags,
refine_flags_updated=np.array(0),
- helmholtz_k=np.array(helmholtz_k),
+ kernel_length_scale=np.array(kernel_length_scale),
wait_for=wait_for)
cl.wait_for_events([evt])
@@ -530,246 +414,229 @@ class QBXLayerPotentialSourceRefiner(DiscrPlotterMixin):
logger.debug("refiner: found {} panel(s) to refine".format(
npanels_to_refine - npanels_to_refine_prev))
- logger.info("refiner: done checking helmholtz k to panel size ratio")
+ logger.info("refiner: done checking kernel length scale to panel size ratio")
return (out["refine_flags_updated"].get() == 1).all()
- def refinement_check_2_to_1_panel_size_ratio(self, queue, lpot_source,
- refine_flags, debug, wait_for=None):
- knl = self.get_2_to_1_panel_ratio_refiner()
- adjacency = self.get_adjacency_on_device(queue, lpot_source)
+ # }}}
- refine_flags_updated = False
+ def build_tree(self, lpot_source, use_base_fine_discr=False):
+ tb = self.code_container.tree_builder()
+ from pytential.qbx.utils import build_tree_with_qbx_metadata
+ return build_tree_with_qbx_metadata(
+ self.queue, tb, lpot_source, use_base_fine_discr=use_base_fine_discr)
- logger.info("refiner: checking 2-to-1 panel size ratio")
+ def find_peer_lists(self, tree):
+ plf = self.code_container.peer_list_finder()
+ peer_lists, evt = plf(self.queue, tree)
+ cl.wait_for_events([evt])
+ return peer_lists
- if debug:
- npanels_to_refine_prev = cl.array.sum(refine_flags).get()
+ def refine(self, density_discr, refiner, refine_flags, factory, debug):
+ """
+ Refine the underlying mesh and discretization.
+ """
+ if isinstance(refine_flags, cl.array.Array):
+ refine_flags = refine_flags.get(self.queue)
+ refine_flags = refine_flags.astype(np.bool)
- # Iterative refinement until no more panels can be marked
- while True:
- evt, out = knl(queue,
- npanels=lpot_source.density_discr.mesh.nelements,
- panel_sizes=lpot_source.panel_sizes("npanels"),
- refine_flags=refine_flags,
- # It's safe to pass this here, as the resulting data
- # race won't affect the final result of the
- # computation.
- refine_flags_prev=refine_flags,
- refine_flags_updated=np.array(0),
- panel_adjacency_starts=adjacency.adjacency_starts,
- panel_adjacency_lists=adjacency.adjacency_lists,
- wait_for=wait_for)
-
- cl.wait_for_events([evt])
-
- if (out["refine_flags_updated"].get() == 1).all():
- refine_flags_updated = True
- else:
- break
+ logger.info("refiner: calling meshmode")
- if debug:
- npanels_to_refine = cl.array.sum(refine_flags).get()
- if npanels_to_refine > npanels_to_refine_prev:
- logger.debug("refiner: found {} panel(s) to refine".format(
- npanels_to_refine - npanels_to_refine_prev))
+ refiner.refine(refine_flags)
+ from meshmode.discretization.connection import make_refinement_connection
+ conn = make_refinement_connection(refiner, density_discr, factory)
- logger.info("refiner: done checking 2-to-1 panel size ratio")
+ logger.info("refiner: done calling meshmode")
- return refine_flags_updated
+ return conn
- # }}}
+# }}}
- # {{{ other utilities
- def get_tunnel_query_dists(self, queue, tree, lpot_source):
- """
- Compute radii for the tubular neighborhood around each panel center of mass.
- """
- nqbxpanels = lpot_source.density_discr.mesh.nelements
- # atomic_max only works on float32
- tq_dists = cl.array.zeros(queue, nqbxpanels, np.float32)
- tq_dists.finish()
-
- knl = self.get_tunnel_query_distance_finder(tree.dimensions,
- tree.coord_dtype, tree.particle_id_dtype)
-
- evt = knl(tree.qbx_user_source_slice.start,
- tree.qbx_user_panel_slice.start,
- nqbxpanels,
- tree.qbx_panel_to_source_starts,
- tree.sorted_target_ids,
- lpot_source.panel_sizes("npanels"),
- tq_dists,
- *tree.sources,
- queue=queue,
- range=slice(tree.nqbxsources))
+class RefinerNotConvergedWarning(UserWarning):
+ pass
- cl.wait_for_events([evt])
- if tree.coord_dtype != tq_dists.dtype:
- tq_dists = tq_dists.astype(tree.coord_dtype)
+def make_empty_refine_flags(queue, lpot_source, use_base_fine_discr=False):
+ """Return an array on the device suitable for use as element refine flags.
- return tq_dists, evt
+ :arg queue: An instance of :class:`pyopencl.CommandQueue`.
+ :arg lpot_source: An instance of :class:`QBXLayerPotentialSource`.
- def get_adjacency_on_device(self, queue, lpot_source):
- """
- Take adjacency information from the mesh and place it onto the device.
- """
- from boxtree.tools import DeviceDataRecord
- adjacency = lpot_source.density_discr.mesh.nodal_adjacency
- adjacency_starts = cl.array.to_device(queue, adjacency.neighbors_starts)
- adjacency_lists = cl.array.to_device(queue, adjacency.neighbors)
+ :returns: A :class:`pyopencl.array.Array` suitable for use as refine flags,
+ initialized to zero.
+ """
+ discr = (lpot_source.base_fine_density_discr
+ if use_base_fine_discr
+ else lpot_source.density_discr)
+ result = cl.array.zeros(queue, discr.mesh.nelements, np.int32)
+ result.finish()
+ return result
- return DeviceDataRecord(
- adjacency_starts=adjacency_starts,
- adjacency_lists=adjacency_lists)
- def get_refine_flags(self, queue, lpot_source):
- """
- Return an array on the device suitable for use as element refine flags.
+# {{{ main entry point
- :arg queue: An instance of :class:`pyopencl.CommandQueue`.
- :arg lpot_source: An instance of :class:`QBXLayerPotentialSource`.
+def refine_for_global_qbx(lpot_source, code_container,
+ group_factory, fine_group_factory, kernel_length_scale=None,
+ # FIXME: Set debug=False once everything works.
+ refine_flags=None, debug=True, maxiter=50):
+ """
+ Entry point for calling the refiner.
- :returns: An instance of :class:`pyopencl.array.Array` suitable for
- use as refine flags, initialized to zero.
- """
- result = cl.array.zeros(
- queue, lpot_source.density_discr.mesh.nelements, np.int32)
- return result, result.events[0]
+ :arg lpot_source: An instance of :class:`QBXLayerPotentialSource`.
- # }}}
+ :arg code_container: An instance of :class:`RefinerCodeContainer`.
- def refine(self, queue, lpot_source, refine_flags, refiner, factory, debug):
- """
- Refine the underlying mesh and discretization.
- """
- if isinstance(refine_flags, cl.array.Array):
- refine_flags = refine_flags.get(queue)
- refine_flags = refine_flags.astype(np.bool)
+ :arg group_factory: An instance of
+ :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for
+ discretizing the coarse refined mesh.
- logger.info("refiner: calling meshmode")
+ :arg fine_group_factory: An instance of
+ :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for
+ oversample the refined mesh at the finest level.
- refiner.refine(refine_flags)
- from meshmode.discretization.connection import make_refinement_connection
+ :arg kernel_length_scale: The kernel length scale, or *None* if not
+ applicable. All panels are refined to below this size.
- conn = make_refinement_connection(
- refiner, lpot_source.density_discr,
- factory)
+ :arg refine_flags: A :class:`pyopencl.array.Array` indicating which
+ panels should get refined initially, or `None` if no initial
+ refinement should be done. Should have size equal to the number of
+ panels. See also :func:`make_empty_refine_flags()`.
- logger.info("refiner: done calling meshmode")
+ :arg maxiter: The maximum number of refiner iterations.
- new_density_discr = conn.to_discr
+ :returns: A tuple ``(lpot_source, *conn*)`` where ``lpot_source`` is the
+ refined layer potential source, and ``conn`` is a
+ :class:`meshmode.discretization.connection.DiscretizationConnection`
+ going from the original mesh to the refined mesh.
+ """
- new_lpot_source = lpot_source.copy(
- density_discr=new_density_discr,
- refined_for_global_qbx=True)
+ # TODO: Stop doing redundant checks by avoiding panels which no longer need
+ # refinement.
- return new_lpot_source, conn
+ from meshmode.mesh.refinement import Refiner
+ from meshmode.discretization.connection import (
+ ChainedDiscretizationConnection, make_same_mesh_connection)
- def __call__(self, lpot_source, discr_factory, helmholtz_k=None,
- # FIXME: Set debug=False once everything works.
- refine_flags=None, debug=True, maxiter=50):
- """
- Entry point for calling the refiner.
+ with cl.CommandQueue(lpot_source.cl_context) as queue:
+ wrangler = code_container.get_wrangler(queue)
- :arg lpot_source: An instance of :class:`QBXLayerPotentialSource`.
+ refiner = Refiner(lpot_source.density_discr.mesh)
+ connections = []
- :arg group_factory: An instance of
- :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for
- discretizing the refined mesh.
+ # Do initial refinement.
+ if refine_flags is not None:
+ conn = wrangler.refine(
+ lpot_source.density_discr, refiner, refine_flags, group_factory,
+ debug)
+ connections.append(conn)
+ lpot_source = lpot_source.copy(density_discr=conn.to_discr)
- :arg helmholtz_k: The Helmholtz parameter, or `None` if not applicable.
+ # {{{ first stage refinement
- :arg refine_flags: A :class:`pyopencl.array.Array` indicating which
- panels should get refined initially, or `None` if no initial
- refinement should be done. Should have size equal to the number of
- panels. See also :meth:`get_refine_flags()`.
+ must_refine = True
+ niter = 0
- :returns: A tuple ``(lpot_source, conns)`` where ``lpot_source`` is the
- refined layer potential source, and ``conns`` is a list of
- :class:`meshmode.discretization.connection.DiscretizationConnection`
- objects going from the original mesh to the refined mesh.
- """
- from meshmode.mesh.refinement import Refiner
- refiner = Refiner(lpot_source.density_discr.mesh)
- connections = []
+ while must_refine:
+ must_refine = False
+ niter += 1
- lpot_source = lpot_source.copy(
- debug=debug,
- refined_for_global_qbx=True)
+ if niter > maxiter:
+ from warnings import warn
+ warn(
+ "Max iteration count reached in QBX layer potential source"
+ " refiner.",
+ RefinerNotConvergedWarning)
+ break
- with cl.CommandQueue(self.cl_context) as queue:
- if refine_flags:
- lpot_source, conn = self.refine(
- queue, lpot_source, refine_flags, refiner, discr_factory,
- debug)
+ # Build tree and auxiliary data.
+ # FIXME: The tree should not have to be rebuilt at each iteration.
+ tree = wrangler.build_tree(lpot_source)
+ peer_lists = wrangler.find_peer_lists(tree)
+ refine_flags = make_empty_refine_flags(queue, lpot_source)
+
+ # Check condition 1.
+ must_refine |= wrangler.check_expansion_disks_undisturbed_by_sources(
+ lpot_source, tree, peer_lists, refine_flags, debug)
+
+ # Check condition 3.
+ if kernel_length_scale is not None:
+ must_refine |= (
+ wrangler.check_kernel_length_scale_to_panel_size_ratio(
+ lpot_source, kernel_length_scale, refine_flags, debug))
+
+ if must_refine:
+ conn = wrangler.refine(
+ lpot_source.density_discr, refiner, refine_flags,
+ group_factory, debug)
connections.append(conn)
+ lpot_source = lpot_source.copy(density_discr=conn.to_discr)
- done_refining = False
- niter = 0
+ del tree
+ del refine_flags
+ del peer_lists
- while not done_refining:
- niter += 1
+ # }}}
- if niter > maxiter:
- logger.warning(
- "Max iteration count reached in QBX layer potential source"
- " refiner.")
- break
+ # {{{ second stage refinement
- # Build tree and auxiliary data.
- # FIXME: The tree should not have to be rebuilt at each iteration.
- tree = self.tree_builder(queue, lpot_source)
- wait_for = []
+ must_refine = True
+ niter = 0
+ fine_connections = []
- peer_lists, evt = self.peer_list_finder(queue, tree, wait_for)
- wait_for = [evt]
+ base_fine_density_discr = lpot_source.density_discr
- refine_flags, evt = self.get_refine_flags(queue, lpot_source)
- wait_for.append(evt)
+ while must_refine:
+ must_refine = False
+ niter += 1
- tq_dists, evt = self.get_tunnel_query_dists(queue, tree, lpot_source)
- wait_for.append(evt)
+ if niter > maxiter:
+ from warnings import warn
+ warn(
+ "Max iteration count reached in QBX layer potential source"
+ " refiner.",
+ RefinerNotConvergedWarning)
+ break
- # Run refinement checkers.
- must_refine = False
+ # Build tree and auxiliary data.
+ # FIXME: The tree should not have to be rebuilt at each iteration.
+ tree = wrangler.build_tree(lpot_source, use_base_fine_discr=True)
+ peer_lists = wrangler.find_peer_lists(tree)
+ refine_flags = make_empty_refine_flags(
+ queue, lpot_source, use_base_fine_discr=True)
- must_refine |= \
- self.refinement_check_center_is_closest_to_orig_panel(
- queue, tree, lpot_source, peer_lists, refine_flags,
- debug, wait_for)
+ must_refine |= wrangler.check_sufficient_source_quadrature_resolution(
+ lpot_source, tree, peer_lists, refine_flags, debug)
- must_refine |= \
- self.refinement_check_center_is_far_from_nonneighbor_panels(
- queue, tree, lpot_source, peer_lists, tq_dists,
- refine_flags, debug, wait_for)
+ if must_refine:
+ conn = wrangler.refine(
+ base_fine_density_discr,
+ refiner, refine_flags, group_factory, debug)
+ base_fine_density_discr = conn.to_discr
+ fine_connections.append(conn)
+ lpot_source = lpot_source.copy(
+ base_resampler=ChainedDiscretizationConnection(
+ fine_connections))
- must_refine |= \
- self.refinement_check_2_to_1_panel_size_ratio(
- queue, lpot_source, refine_flags, debug, wait_for)
+ del tree
+ del refine_flags
+ del peer_lists
- if helmholtz_k:
- must_refine |= \
- self.refinement_check_helmholtz_k_to_panel_size_ratio(
- queue, lpot_source, helmholtz_k, refine_flags, debug,
- wait_for)
+ # }}}
- if must_refine:
- lpot_source, conn = self.refine(
- queue, lpot_source, refine_flags, refiner, discr_factory,
- debug)
- connections.append(conn)
+ lpot_source = lpot_source.copy(debug=debug, refined_for_global_qbx=True)
- del tree
- del peer_lists
- del tq_dists
- del refine_flags
- done_refining = not must_refine
+ if len(connections) == 0:
+ # FIXME: This is inefficient
+ connection = make_same_mesh_connection(
+ lpot_source.density_discr,
+ lpot_source.density_discr)
+ else:
+ connection = ChainedDiscretizationConnection(connections)
- return lpot_source, connections
+ return lpot_source, connection
# }}}
+
# vim: foldmethod=marker:filetype=pyopencl
diff --git a/pytential/qbx/target_assoc.py b/pytential/qbx/target_assoc.py
index 8890d992106e09509ceee424d73f8fa13ec9fa4f..4255d865f753da087e1cca1e7942248d14b684b0 100644
--- a/pytential/qbx/target_assoc.py
+++ b/pytential/qbx/target_assoc.py
@@ -36,7 +36,7 @@ from boxtree.tools import DeviceDataRecord
from boxtree.area_query import AreaQueryElementwiseTemplate
from boxtree.tools import InlineBinarySearch
from pytential.qbx.utils import (
- QBX_TREE_C_PREAMBLE, QBX_TREE_MAKO_DEFS, DiscrPlotterMixin)
+ QBX_TREE_C_PREAMBLE, QBX_TREE_MAKO_DEFS)
unwrap_args = AreaQueryElementwiseTemplate.unwrap_args
@@ -128,7 +128,7 @@ QBX_TARGET_MARKER = AreaQueryElementwiseTemplate(
particle_id_t source_offset,
particle_id_t target_offset,
particle_id_t *sorted_target_ids,
- coord_t *panel_sizes,
+ coord_t *tunnel_radius_by_source,
coord_t *box_to_search_dist,
/* output */
@@ -140,7 +140,7 @@ QBX_TARGET_MARKER = AreaQueryElementwiseTemplate(
coord_t *particles_${ax},
%endfor
""",
- ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r"""
+ ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r"""//CL//
coord_vec_t tgt_coords;
${load_particle("INDEX_FOR_TARGET_PARTICLE(i)", "tgt_coords")}
{
@@ -149,7 +149,7 @@ QBX_TARGET_MARKER = AreaQueryElementwiseTemplate(
${ball_radius} = box_to_search_dist[my_box];
}
""",
- leaf_found_op=QBX_TREE_MAKO_DEFS + r"""
+ leaf_found_op=QBX_TREE_MAKO_DEFS + r"""//CL//
for (particle_id_t source_idx = box_to_source_starts[${leaf_box_id}];
source_idx < box_to_source_starts[${leaf_box_id} + 1];
++source_idx)
@@ -158,7 +158,8 @@ QBX_TARGET_MARKER = AreaQueryElementwiseTemplate(
coord_vec_t source_coords;
${load_particle("INDEX_FOR_SOURCE_PARTICLE(source)", "source_coords")}
- if (distance(source_coords, tgt_coords) <= panel_sizes[source] / 2)
+ if (distance(source_coords, tgt_coords)
+ <= tunnel_radius_by_source[source])
{
target_status[i] = MARKED_QBX_CENTER_PENDING;
*found_target_close_to_panel = 1;
@@ -177,9 +178,8 @@ QBX_CENTER_FINDER = AreaQueryElementwiseTemplate(
particle_id_t center_offset,
particle_id_t target_offset,
particle_id_t *sorted_target_ids,
- coord_t *panel_sizes,
+ coord_t *expansion_radii_by_center_with_stick_out,
coord_t *box_to_search_dist,
- coord_t stick_out_factor,
int *target_flags,
/* input/output */
@@ -222,7 +222,7 @@ QBX_CENTER_FINDER = AreaQueryElementwiseTemplate(
coord_t my_dist_to_center = distance(tgt_coords, center_coords);
if (my_dist_to_center
- <= (panel_sizes[center] / 2) * (1 + stick_out_factor)
+ <= expansion_radii_by_center_with_stick_out[center]
&& my_dist_to_center < min_dist_to_center[i])
{
target_status[i] = MARKED_QBX_CENTER_FOUND;
@@ -245,7 +245,7 @@ QBX_FAILED_TARGET_ASSOCIATION_REFINER = AreaQueryElementwiseTemplate(
particle_id_t target_offset,
int npanels,
particle_id_t *sorted_target_ids,
- coord_t *panel_sizes,
+ coord_t *tunnel_radius_by_source,
int *target_status,
coord_t *box_to_search_dist,
@@ -258,7 +258,7 @@ QBX_FAILED_TARGET_ASSOCIATION_REFINER = AreaQueryElementwiseTemplate(
coord_t *particles_${ax},
%endfor
""",
- ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r"""
+ ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r"""//CL//
coord_vec_t target_coords;
${load_particle("INDEX_FOR_TARGET_PARTICLE(i)", "target_coords")}
{
@@ -267,7 +267,7 @@ QBX_FAILED_TARGET_ASSOCIATION_REFINER = AreaQueryElementwiseTemplate(
${ball_radius} = box_to_search_dist[my_box];
}
""",
- leaf_found_op=QBX_TREE_MAKO_DEFS + r"""
+ leaf_found_op=QBX_TREE_MAKO_DEFS + r"""//CL//
for (particle_id_t source_idx = box_to_source_starts[${leaf_box_id}];
source_idx < box_to_source_starts[${leaf_box_id} + 1];
++source_idx)
@@ -279,7 +279,7 @@ QBX_FAILED_TARGET_ASSOCIATION_REFINER = AreaQueryElementwiseTemplate(
bool is_close =
distance(target_coords, source_coords)
- <= panel_sizes[source] / 2;
+ <= tunnel_radius_by_source[source];
if (is_close && target_status[i] == MARKED_QBX_CENTER_PENDING)
{
@@ -318,11 +318,11 @@ class QBXTargetAssociation(DeviceDataRecord):
pass
-class QBXTargetAssociator(DiscrPlotterMixin):
+class QBXTargetAssociator(object):
def __init__(self, cl_context):
- from pytential.qbx.utils import TreeWithQBXMetadataBuilder
- self.tree_builder = TreeWithQBXMetadataBuilder(cl_context)
+ from boxtree.tree_build import TreeBuilder
+ self.tree_builder = TreeBuilder(cl_context)
self.cl_context = cl_context
from boxtree.area_query import PeerListFinder, SpaceInvaderQueryBuilder
self.peer_list_finder = PeerListFinder(cl_context)
@@ -398,19 +398,37 @@ class QBXTargetAssociator(DiscrPlotterMixin):
# Perform a space invader query over the sources.
source_slice = tree.user_source_ids[tree.qbx_user_source_slice]
sources = [axis.with_queue(queue)[source_slice] for axis in tree.sources]
- panel_sizes = lpot_source.panel_sizes("nsources").with_queue(queue)
+ tunnel_radius_by_source = \
+ lpot_source._close_target_tunnel_radius("nsources").with_queue(queue)
+
+ # Target-marking algorithm (TGTMARK):
+ #
+ # (1) Use a space invader query to tag each leaf box that intersects with the
+ # "near-source-detection tunnel" with the distance to the closest source.
+ #
+ # (2) Do an area query around all targets with the radius resulting
+ # from the space invader query, enumerate sources in that vicinity.
+ # If a source is found whose distance to the target is less than the
+ # source's tunnel radius, mark that target as pending.
+ # (or below: mark the source for refinement)
+
+ # Note that this comment is referred to below by "TGTMARK". If you
+ # remove this comment or change the algorithm here, make sure that
+ # the reference below is still accurate.
box_to_search_dist, evt = self.space_invader_query(
queue,
tree,
sources,
- panel_sizes / 2,
+ tunnel_radius_by_source,
peer_lists,
wait_for=wait_for)
wait_for = [evt]
logger.info("target association: marking targets close to panels")
+ tunnel_radius_by_source = lpot_source._close_target_tunnel_radius("nsources")
+
evt = knl(
*unwrap_args(
tree, peer_lists,
@@ -419,7 +437,7 @@ class QBXTargetAssociator(DiscrPlotterMixin):
tree.qbx_user_source_slice.start,
tree.qbx_user_target_slice.start,
tree.sorted_target_ids,
- panel_sizes,
+ tunnel_radius_by_source,
box_to_search_dist,
target_status,
found_target_close_to_panel,
@@ -463,13 +481,22 @@ class QBXTargetAssociator(DiscrPlotterMixin):
center_slice = \
tree.sorted_target_ids[tree.qbx_user_center_slice].with_queue(queue)
centers = [axis.with_queue(queue)[center_slice] for axis in tree.sources]
- panel_sizes = lpot_source.panel_sizes("ncenters").with_queue(queue)
+ expansion_radii_by_center = \
+ lpot_source._expansion_radii("ncenters").with_queue(queue)
+ expansion_radii_by_center_with_stick_out = \
+ expansion_radii_by_center * (1 + stick_out_factor)
+
+ # Idea:
+ #
+ # (1) Tag leaf boxes around centers with max distance to usable center.
+ # (2) Area query from targets with those radii to find closest eligible
+ # center.
box_to_search_dist, evt = self.space_invader_query(
queue,
tree,
centers,
- panel_sizes * ((1 + stick_out_factor) / 2),
+ expansion_radii_by_center_with_stick_out,
peer_lists,
wait_for=wait_for)
wait_for = [evt]
@@ -490,9 +517,8 @@ class QBXTargetAssociator(DiscrPlotterMixin):
tree.qbx_user_center_slice.start,
tree.qbx_user_target_slice.start,
tree.sorted_target_ids,
- panel_sizes,
+ expansion_radii_by_center_with_stick_out,
box_to_search_dist,
- stick_out_factor,
target_flags,
target_status,
target_assoc.target_to_center,
@@ -535,13 +561,16 @@ class QBXTargetAssociator(DiscrPlotterMixin):
# Perform a space invader query over the sources.
source_slice = tree.user_source_ids[tree.qbx_user_source_slice]
sources = [axis.with_queue(queue)[source_slice] for axis in tree.sources]
- panel_sizes = lpot_source.panel_sizes("nsources").with_queue(queue)
+ tunnel_radius_by_source = \
+ lpot_source._close_target_tunnel_radius("nsources").with_queue(queue)
+
+ # See (TGTMARK) above for algorithm.
box_to_search_dist, evt = self.space_invader_query(
queue,
tree,
sources,
- panel_sizes / 2,
+ tunnel_radius_by_source,
peer_lists,
wait_for=wait_for)
wait_for = [evt]
@@ -558,7 +587,7 @@ class QBXTargetAssociator(DiscrPlotterMixin):
tree.qbx_user_target_slice.start,
tree.nqbxpanels,
tree.sorted_target_ids,
- lpot_source.panel_sizes("nsources"),
+ lpot_source._close_target_tunnel_radius("nsources"),
target_status,
box_to_search_dist,
refine_flags,
@@ -654,8 +683,14 @@ class QBXTargetAssociator(DiscrPlotterMixin):
"""
with cl.CommandQueue(self.cl_context) as queue:
- tree = self.tree_builder(queue, lpot_source,
+ from pytential.qbx.utils import build_tree_with_qbx_metadata
+
+ tree = build_tree_with_qbx_metadata(
+ queue,
+ self.tree_builder,
+ lpot_source,
[discr for discr, _ in target_discrs_and_qbx_sides])
+
peer_lists, evt = self.peer_list_finder(queue, tree, wait_for)
wait_for = [evt]
diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py
index 60da4b5b2a4e914696ba0ea1083c0cc1a6f825d0..02a19e007eae86e90066224f97b0b261a02a7e6a 100644
--- a/pytential/qbx/utils.py
+++ b/pytential/qbx/utils.py
@@ -99,8 +99,8 @@ def get_interleaved_centers(queue, lpot_source):
next to corresponding exterior centers.
"""
knl = get_interleaver_kernel(lpot_source.density_discr.real_dtype)
- int_centers = lpot_source.centers(-1)
- ext_centers = lpot_source.centers(+1)
+ int_centers = get_centers_on_side(lpot_source, -1)
+ ext_centers = get_centers_on_side(lpot_source, +1)
result = []
wait_for = []
@@ -118,203 +118,433 @@ def get_interleaved_centers(queue, lpot_source):
# }}}
-# {{{ discr plotter mixin
-
-class DiscrPlotterMixin(object):
-
- def plot_discr(self, lpot_source, outfilename="discr.pdf"):
- with cl.CommandQueue(self.cl_context) as queue:
- tree = self.tree_builder(queue, lpot_source).get(queue=queue)
- from boxtree.visualization import TreePlotter
-
- import matplotlib
- matplotlib.use('Agg')
- import matplotlib.pyplot as plt
- tp = TreePlotter(tree)
- tp.draw_tree()
- sources = (tree.sources[0], tree.sources[1])
- sti = tree.sorted_target_ids
- plt.plot(sources[0][sti[tree.qbx_user_source_slice]],
- sources[1][sti[tree.qbx_user_source_slice]],
- lw=0, marker=".", markersize=1, label="sources")
- plt.plot(sources[0][sti[tree.qbx_user_center_slice]],
- sources[1][sti[tree.qbx_user_center_slice]],
- lw=0, marker=".", markersize=1, label="centers")
- plt.plot(sources[0][sti[tree.qbx_user_target_slice]],
- sources[1][sti[tree.qbx_user_target_slice]],
- lw=0, marker=".", markersize=1, label="targets")
- plt.axis("equal")
- plt.legend()
- plt.savefig(outfilename)
+# {{{ make interleaved radii
+
+def get_interleaved_radii(queue, lpot_source):
+ """
+ Return an array of shape (dim, ncenters) in which interior centers are placed
+ next to corresponding exterior centers.
+ """
+ knl = get_interleaver_kernel(lpot_source.density_discr.real_dtype)
+ radii = lpot_source._expansion_radii("nsources")
+
+ result = cl.array.empty(queue, len(radii) * 2, radii.dtype)
+ evt, _ = knl(queue, src1=radii, src2=radii, dst=result)
+ evt.wait()
+
+ return result
# }}}
-# {{{ tree creation
+# {{{ panel sizes
+
+def panel_sizes(discr, last_dim_length):
+ if last_dim_length not in ("nsources", "ncenters", "npanels"):
+ raise ValueError(
+ "invalid value of last_dim_length: %s" % last_dim_length)
+
+ # To get the panel size this does the equivalent of (∫ 1 ds)**(1/dim).
+ # FIXME: Kernel optimizations
+
+ if last_dim_length == "nsources" or last_dim_length == "ncenters":
+ knl = lp.make_kernel(
+ "{[i,j,k]: 0<=i<nelements and 0<=j,k<nunit_nodes}",
+ "panel_sizes[i,j] = sum(k, ds[i,k])**(1/dim)",
+ name="compute_size")
+
+ def panel_size_view(discr, group_nr):
+ return discr.groups[group_nr].view
+
+ elif last_dim_length == "npanels":
+ knl = lp.make_kernel(
+ "{[i,j]: 0<=i<nelements and 0<=j<nunit_nodes}",
+ "panel_sizes[i] = sum(j, ds[i,j])**(1/dim)",
+ name="compute_size")
+ from functools import partial
+
+ def panel_size_view(discr, group_nr):
+ return partial(el_view, discr, group_nr)
+
+ else:
+ raise ValueError("unknown dim length specified")
+
+ knl = lp.fix_parameters(knl, dim=discr.dim)
+
+ with cl.CommandQueue(discr.cl_context) as queue:
+ from pytential import bind, sym
+ ds = bind(
+ discr,
+ sym.area_element(ambient_dim=discr.ambient_dim, dim=discr.dim)
+ * sym.QWeight()
+ )(queue)
+ panel_sizes = cl.array.empty(
+ queue, discr.nnodes
+ if last_dim_length in ("nsources", "ncenters")
+ else discr.mesh.nelements, discr.real_dtype)
+ for group_nr, group in enumerate(discr.groups):
+ _, (result,) = knl(queue,
+ nelements=group.nelements,
+ nunit_nodes=group.nunit_nodes,
+ ds=group.view(ds),
+ panel_sizes=panel_size_view(
+ discr, group_nr)(panel_sizes))
+ panel_sizes.finish()
+ if last_dim_length == "ncenters":
+ from pytential.qbx.utils import get_interleaver_kernel
+ knl = get_interleaver_kernel(discr.real_dtype)
+ _, (panel_sizes,) = knl(queue, dstlen=2*discr.nnodes,
+ src1=panel_sizes, src2=panel_sizes)
+ return panel_sizes.with_queue(None)
+
+# }}}
-class TreeWithQBXMetadataBuilder(object):
- MAX_REFINE_WEIGHT = 64
+# {{{ element centers of mass
- class TreeWithQBXMetadata(Tree):
+def element_centers_of_mass(discr):
+ knl = lp.make_kernel(
+ """{[dim,k,i]:
+ 0<=dim<ndims and
+ 0<=k<nelements and
+ 0<=i<nunit_nodes}""",
"""
- .. attribute:: nqbxpanels
- .. attribuet:: nqbxsources
- .. attribute:: nqbxcenters
- .. attribute:: nqbxtargets
+ panels[dim, k] = sum(i, nodes[dim, k, i])/nunit_nodes
+ """,
+ default_offset=lp.auto, name="find_panel_centers_of_mass")
+
+ knl = lp.fix_parameters(knl, ndims=discr.ambient_dim)
+
+ knl = lp.split_iname(knl, "k", 128, inner_tag="l.0", outer_tag="g.0")
+ knl = lp.tag_inames(knl, dict(dim="ilp"))
+
+ with cl.CommandQueue(discr.cl_context) as queue:
+ mesh = discr.mesh
+ panels = cl.array.empty(queue, (mesh.ambient_dim, mesh.nelements),
+ dtype=discr.real_dtype)
+ for group_nr, group in enumerate(discr.groups):
+ _, (result,) = knl(queue,
+ nelements=group.nelements,
+ nunit_nodes=group.nunit_nodes,
+ nodes=group.view(discr.nodes()),
+ panels=el_view(discr, group_nr, panels))
+ panels.finish()
+ panels = panels.with_queue(None)
+ return tuple(panels[d, :] for d in range(mesh.ambient_dim))
- .. attribute:: box_to_qbx_panel_starts
- .. attribute:: box_to_qbx_panel_lists
+# }}}
- .. attribute:: box_to_qbx_source_starts
- .. attribute:: box_to_qbx_source_lists
- .. attribute:: box_to_qbx_center_starts
- .. attribute:: box_to_qbx_center_lists
+# {{{ compute center array
- .. attribute:: box_to_qbx_target_starts
- .. attribute:: box_to_qbx_target_lists
+def get_centers_on_side(lpot_src, sign):
+ adim = lpot_src.density_discr.ambient_dim
+ dim = lpot_src.density_discr.dim
- .. attribute:: qbx_panel_to_source_starts
- .. attribute:: qbx_panel_to_center_starts
+ from pytential import sym, bind
+ with cl.CommandQueue(lpot_src.cl_context) as queue:
+ nodes = bind(lpot_src.density_discr, sym.nodes(adim))(queue)
+ normals = bind(lpot_src.density_discr, sym.normal(adim, dim=dim))(queue)
+ expansion_radii = lpot_src._expansion_radii("nsources").with_queue(queue)
+ return (nodes + normals * sign * expansion_radii).as_vector(np.object)
- .. attribute:: qbx_user_source_slice
- .. attribute:: qbx_user_center_slice
- .. attribute:: qbx_user_panel_slice
- .. attribute:: qbx_user_target_slice
- """
- pass
+# }}}
- def __init__(self, context):
- self.context = context
- from boxtree.tree_build import TreeBuilder
- self.tree_builder = TreeBuilder(self.context)
- def __call__(self, queue, lpot_source, targets_list=()):
- """
- Return a :class:`TreeWithQBXMetadata` built from the given layer
- potential source.
+# {{{ el_view
- :arg queue: An instance of :class:`pyopencl.CommandQueue`
+def el_view(discr, group_nr, global_array):
+ """Return a view of *global_array* of shape
+ ``(..., discr.groups[group_nr].nelements)``
+ where *global_array* is of shape ``(..., nelements)``,
+ where *nelements* is the global (per-discretization) node count.
+ """
- :arg lpot_source: An instance of
- :class:`pytential.qbx.NewQBXLayerPotentialSource`.
+ group = discr.groups[group_nr]
+ el_nr_base = sum(group.nelements for group in discr.groups[:group_nr])
- :arg targets_list: A list of :class:`pytential.target.TargetBase`
- """
- # The ordering of particles is as follows:
- # - sources go first
- # - then centers
- # - then panels (=centers of mass)
- # - then targets
-
- logger.info("start building tree with qbx metadata")
-
- sources = lpot_source.density_discr.nodes()
- centers = get_interleaved_centers(queue, lpot_source)
- centers_of_mass = lpot_source.panel_centers_of_mass()
- targets = (tgt.nodes() for tgt in targets_list)
-
- particles = tuple(
- cl.array.concatenate(dim_coords, queue=queue)
- for dim_coords in zip(sources, centers, centers_of_mass, *targets))
-
- # Counts
- nparticles = len(particles[0])
- npanels = len(centers_of_mass[0])
- nsources = len(sources[0])
- ncenters = len(centers[0])
- # Each source gets an interior / exterior center.
- assert 2 * nsources == ncenters
- ntargets = sum(tgt.nnodes for tgt in targets_list)
-
- # Slices
- qbx_user_source_slice = slice(0, nsources)
- panel_slice_start = 3 * nsources
- qbx_user_center_slice = slice(nsources, panel_slice_start)
- target_slice_start = panel_slice_start + npanels
- qbx_user_panel_slice = slice(panel_slice_start, panel_slice_start + npanels)
- qbx_user_target_slice = slice(target_slice_start,
- target_slice_start + ntargets)
-
- # Build tree with sources, centers, and centers of mass. Split boxes
- # only because of sources.
- refine_weights = cl.array.zeros(queue, nparticles, np.int32)
- refine_weights[:nsources].fill(1)
-
- refine_weights.finish()
-
- tree, evt = self.tree_builder(queue, particles,
- max_leaf_refine_weight=self.MAX_REFINE_WEIGHT,
- refine_weights=refine_weights)
-
- # Compute box => particle class relations
- flags = refine_weights
- del refine_weights
- particle_classes = {}
-
- for class_name, particle_slice, fixup in (
- ("box_to_qbx_source", qbx_user_source_slice, 0),
- ("box_to_qbx_target", qbx_user_target_slice, -target_slice_start),
- ("box_to_qbx_center", qbx_user_center_slice, -nsources),
- ("box_to_qbx_panel", qbx_user_panel_slice, -panel_slice_start)):
- flags.fill(0)
- flags[particle_slice].fill(1)
- flags.finish()
-
- from boxtree.tree import filter_target_lists_in_user_order
- box_to_class = (
- filter_target_lists_in_user_order(queue, tree, flags)
- .with_queue(queue))
-
- if fixup:
- box_to_class.target_lists += fixup
- particle_classes[class_name + "_starts"] = box_to_class.target_starts
- particle_classes[class_name + "_lists"] = box_to_class.target_lists
-
- del flags
- del box_to_class
-
- # Compute panel => source relation
- qbx_panel_to_source_starts = cl.array.empty(
+ return global_array[
+ ..., el_nr_base:el_nr_base + group.nelements] \
+ .reshape(
+ global_array.shape[:-1]
+ + (group.nelements,))
+
+# }}}
+
+
+# {{{ discr plotter
+
+def plot_discr(lpot_source, outfilename="discr.pdf"):
+ with cl.CommandQueue(lpot_source.cl_context) as queue:
+ from boxtree.tree_builder import TreeBuilder
+ tree_builder = TreeBuilder(lpot_source.cl_context)
+ tree = tree_builder(queue, lpot_source).get(queue=queue)
+ from boxtree.visualization import TreePlotter
+
+ import matplotlib
+ matplotlib.use('Agg')
+ import matplotlib.pyplot as plt
+ tp = TreePlotter(tree)
+ tp.draw_tree()
+ sources = (tree.sources[0], tree.sources[1])
+ sti = tree.sorted_target_ids
+ plt.plot(sources[0][sti[tree.qbx_user_source_slice]],
+ sources[1][sti[tree.qbx_user_source_slice]],
+ lw=0, marker=".", markersize=1, label="sources")
+ plt.plot(sources[0][sti[tree.qbx_user_center_slice]],
+ sources[1][sti[tree.qbx_user_center_slice]],
+ lw=0, marker=".", markersize=1, label="centers")
+ plt.plot(sources[0][sti[tree.qbx_user_target_slice]],
+ sources[1][sti[tree.qbx_user_target_slice]],
+ lw=0, marker=".", markersize=1, label="targets")
+ plt.axis("equal")
+ plt.legend()
+ plt.savefig(outfilename)
+
+# }}}
+
+
+# {{{ tree-with-metadata: data structure
+
+class TreeWithQBXMetadata(Tree):
+ """A subclass of :class:`boxtree.tree.Tree`. Has all of that class's
+ attributes, along with the following:
+
+ .. attribute:: nqbxpanels
+ .. attribuet:: nqbxsources
+ .. attribute:: nqbxcenters
+ .. attribute:: nqbxtargets
+
+ .. ------------------------------------------------------------------------
+ .. rubric:: Box properties
+ .. ------------------------------------------------------------------------
+
+ .. rubric:: Box to QBX panels
+
+ .. attribute:: box_to_qbx_panel_starts
+
+ ``box_id_t [nboxes + 1]``
+
+ .. attribute:: box_to_qbx_panel_lists
+
+ ``particle_id_t [*]``
+
+ .. rubric:: Box to QBX sources
+
+ .. attribute:: box_to_qbx_source_starts
+
+ ``box_id_t [nboxes + 1]``
+
+ .. attribute:: box_to_qbx_source_lists
+
+ ``particle_id_t [*]``
+
+ .. rubric:: Box to QBX centers
+
+ .. attribute:: box_to_qbx_center_starts
+
+ ``box_id_t [nboxes + 1]``
+
+ .. attribute:: box_to_qbx_center_lists
+
+ ``particle_id_t [*]``
+
+ .. rubric:: Box to QBX targets
+
+ .. attribute:: box_to_qbx_target_starts
+
+ ``box_id_t [nboxes + 1]``
+
+ .. attribute:: box_to_qbx_target_lists
+
+ ``particle_id_t [*]``
+
+ .. ------------------------------------------------------------------------
+ .. rubric:: Panel properties
+ .. ------------------------------------------------------------------------
+
+ .. attribute:: qbx_panel_to_source_starts
+
+ ``particle_id_t [nqbxpanels + 1]``
+
+ .. attribute:: qbx_panel_to_center_starts
+
+ ``particle_id_t [nqbxpanels + 1]``
+
+ .. ------------------------------------------------------------------------
+ .. rubric:: Particle order indices
+ .. ------------------------------------------------------------------------
+
+ .. attribute:: qbx_user_source_slice
+ .. attribute:: qbx_user_center_slice
+ .. attribute:: qbx_user_panel_slice
+ .. attribute:: qbx_user_target_slice
+ """
+ pass
+
+# }}}
+
+
+# {{{ tree-with-metadata: creation
+
+MAX_REFINE_WEIGHT = 64
+
+
+def build_tree_with_qbx_metadata(
+ queue, tree_builder, lpot_source, targets_list=(),
+ use_base_fine_discr=False):
+ """Return a :class:`TreeWithQBXMetadata` built from the given layer
+ potential source. This contains particles of four different types:
+
+ * source particles and panel centers of mass either from
+ ``lpot_source.density_discr`` or ``lpot_source.base_fine_density_discr``
+ * centers from ``lpot_source.centers()``
+ * targets from ``targets_list``.
+
+ :arg queue: An instance of :class:`pyopencl.CommandQueue`
+
+ :arg lpot_source: An instance of
+ :class:`pytential.qbx.NewQBXLayerPotentialSource`.
+
+ :arg targets_list: A list of :class:`pytential.target.TargetBase`
+
+ :arg use_base_fine_discr: If *True*, builds a tree with sources/centers of
+ mass from ``lpot_source.base_fine_density_discr``. If *False* (default),
+ they are from ``lpot_source.density_discr``.
+ """
+ # The ordering of particles is as follows:
+ # - sources go first
+ # - then centers
+ # - then panels (=centers of mass)
+ # - then targets
+
+ logger.info("start building tree with qbx metadata")
+
+ sources = (
+ lpot_source.density_discr.nodes()
+ if not use_base_fine_discr
+ else lpot_source.fine_density_discr.nodes())
+
+ centers = get_interleaved_centers(queue, lpot_source)
+
+ centers_of_mass = (
+ lpot_source._panel_centers_of_mass()
+ if not use_base_fine_discr
+ else lpot_source._fine_panel_centers_of_mass())
+
+ targets = (tgt.nodes() for tgt in targets_list)
+
+ particles = tuple(
+ cl.array.concatenate(dim_coords, queue=queue)
+ for dim_coords in zip(sources, centers, centers_of_mass, *targets))
+
+ # Counts
+ nparticles = len(particles[0])
+ npanels = len(centers_of_mass[0])
+ nsources = len(sources[0])
+ ncenters = len(centers[0])
+ # Each source gets an interior / exterior center.
+ assert 2 * nsources == ncenters or use_base_fine_discr
+ ntargets = sum(tgt.nnodes for tgt in targets_list)
+
+ # Slices
+ qbx_user_source_slice = slice(0, nsources)
+
+ center_slice_start = nsources
+ qbx_user_center_slice = slice(center_slice_start, center_slice_start + ncenters)
+
+ panel_slice_start = center_slice_start + ncenters
+ qbx_user_panel_slice = slice(panel_slice_start, panel_slice_start + npanels)
+
+ target_slice_start = panel_slice_start + npanels
+ qbx_user_target_slice = slice(target_slice_start, target_slice_start + ntargets)
+
+ # Build tree with sources, centers, and centers of mass. Split boxes
+ # only because of sources.
+ refine_weights = cl.array.zeros(queue, nparticles, np.int32)
+ refine_weights[:nsources].fill(1)
+
+ refine_weights.finish()
+
+ tree, evt = tree_builder(queue, particles,
+ max_leaf_refine_weight=MAX_REFINE_WEIGHT,
+ refine_weights=refine_weights)
+
+ # Compute box => particle class relations
+ flags = refine_weights
+ del refine_weights
+ particle_classes = {}
+
+ for class_name, particle_slice, fixup in (
+ ("box_to_qbx_source", qbx_user_source_slice, 0),
+ ("box_to_qbx_target", qbx_user_target_slice, -target_slice_start),
+ ("box_to_qbx_center", qbx_user_center_slice, -center_slice_start),
+ ("box_to_qbx_panel", qbx_user_panel_slice, -panel_slice_start)):
+ flags.fill(0)
+ flags[particle_slice].fill(1)
+ flags.finish()
+
+ from boxtree.tree import filter_target_lists_in_user_order
+ box_to_class = (
+ filter_target_lists_in_user_order(queue, tree, flags)
+ .with_queue(queue))
+
+ if fixup:
+ box_to_class.target_lists += fixup
+ particle_classes[class_name + "_starts"] = box_to_class.target_starts
+ particle_classes[class_name + "_lists"] = box_to_class.target_lists
+
+ del flags
+ del box_to_class
+
+ # Compute panel => source relation
+ if use_base_fine_discr:
+ density_discr = lpot_source.fine_density_discr
+ else:
+ density_discr = lpot_source.density_discr
+
+ qbx_panel_to_source_starts = cl.array.empty(
queue, npanels + 1, dtype=tree.particle_id_dtype)
- el_offset = 0
- for group in lpot_source.density_discr.groups:
- qbx_panel_to_source_starts[el_offset:el_offset + group.nelements] = \
- cl.array.arange(queue, group.node_nr_base,
- group.node_nr_base + group.nnodes,
- group.nunit_nodes,
- dtype=tree.particle_id_dtype)
- el_offset += group.nelements
- qbx_panel_to_source_starts[-1] = nsources
-
- # Compute panel => center relation
- qbx_panel_to_center_starts = 2 * qbx_panel_to_source_starts
-
- # Transfer all tree attributes.
- tree_attrs = {}
- for attr_name in tree.__class__.fields:
- try:
- tree_attrs[attr_name] = getattr(tree, attr_name)
- except AttributeError:
- pass
-
- tree_attrs.update(particle_classes)
-
- logger.info("done building tree with qbx metadata")
-
- return self.TreeWithQBXMetadata(
- qbx_panel_to_source_starts=qbx_panel_to_source_starts,
- qbx_panel_to_center_starts=qbx_panel_to_center_starts,
- qbx_user_source_slice=qbx_user_source_slice,
- qbx_user_panel_slice=qbx_user_panel_slice,
- qbx_user_center_slice=qbx_user_center_slice,
- qbx_user_target_slice=qbx_user_target_slice,
- nqbxpanels=npanels,
- nqbxsources=nsources,
- nqbxcenters=ncenters,
- nqbxtargets=ntargets,
- **tree_attrs).with_queue(None)
+ el_offset = 0
+ for group in density_discr.groups:
+ qbx_panel_to_source_starts[el_offset:el_offset + group.nelements] = \
+ cl.array.arange(queue, group.node_nr_base,
+ group.node_nr_base + group.nnodes,
+ group.nunit_nodes,
+ dtype=tree.particle_id_dtype)
+ el_offset += group.nelements
+ qbx_panel_to_source_starts[-1] = nsources
+
+ # Compute panel => center relation
+ qbx_panel_to_center_starts = (
+ 2 * qbx_panel_to_source_starts
+ if not use_base_fine_discr
+ else None)
+
+ # Transfer all tree attributes.
+ tree_attrs = {}
+ for attr_name in tree.__class__.fields:
+ try:
+ tree_attrs[attr_name] = getattr(tree, attr_name)
+ except AttributeError:
+ pass
+
+ tree_attrs.update(particle_classes)
+
+ logger.info("done building tree with qbx metadata")
+
+ return TreeWithQBXMetadata(
+ qbx_panel_to_source_starts=qbx_panel_to_source_starts,
+ qbx_panel_to_center_starts=qbx_panel_to_center_starts,
+ qbx_user_source_slice=qbx_user_source_slice,
+ qbx_user_panel_slice=qbx_user_panel_slice,
+ qbx_user_center_slice=qbx_user_center_slice,
+ qbx_user_target_slice=qbx_user_target_slice,
+ nqbxpanels=npanels,
+ nqbxsources=nsources,
+ nqbxcenters=ncenters,
+ nqbxtargets=ntargets,
+ **tree_attrs).with_queue(None)
# }}}
diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py
index 0c30cb4bea611766c2ac50ec74ab180c080c3ce4..16679ac20f6cbfa4ca52469b3441a99e2ad49353 100644
--- a/pytential/symbolic/mappers.py
+++ b/pytential/symbolic/mappers.py
@@ -282,6 +282,30 @@ class DerivativeTaker(Mapper):
def __init__(self, ambient_axis):
self.ambient_axis = ambient_axis
+ def map_sum(self, expr):
+ from pymbolic.primitives import flattened_sum
+ return flattened_sum(tuple(self.rec(child) for child in expr.children))
+
+ def map_product(self, expr):
+ from pymbolic.primitives import is_constant
+ const = []
+ nonconst = []
+ for subexpr in expr.children:
+ if is_constant(subexpr):
+ const.append(subexpr)
+ else:
+ nonconst.append(subexpr)
+
+ if len(nonconst) > 1:
+ raise RuntimeError("DerivativeTaker doesn't support products with "
+ "more than one non-constant")
+
+ if not nonconst:
+ nonconst = [1]
+
+ from pytools import product
+ return product(const) * self.rec(nonconst[0])
+
def map_int_g(self, expr):
from sumpy.kernel import AxisTargetDerivative
return expr.copy(kernel=AxisTargetDerivative(self.ambient_axis, expr.kernel))
diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py
index 9b4622fde5f4bf8cc0527b05b09c38433fdf4267..2012b9928c0dad7e10947ba110c5ca65e9041faa 100644
--- a/pytential/symbolic/matrix.py
+++ b/pytential/symbolic/matrix.py
@@ -183,11 +183,13 @@ class MatrixBuilder(EvaluationMapperBase):
assert target_discr is source.density_discr
+ from pytential.qbx.utils import get_centers_on_side
+
assert abs(expr.qbx_forced_limit) > 0
_, (mat,) = mat_gen(self.queue,
target_discr.nodes(),
source.fine_density_discr.nodes(),
- source.centers(expr.qbx_forced_limit),
+ get_centers_on_side(source, expr.qbx_forced_limit),
**kernel_args)
mat = mat.get()
diff --git a/pytential/symbolic/pde/cahn_hilliard.py b/pytential/symbolic/pde/cahn_hilliard.py
new file mode 100644
index 0000000000000000000000000000000000000000..fd6ec4824cc107fb5d4ea4de8520997f59557ca5
--- /dev/null
+++ b/pytential/symbolic/pde/cahn_hilliard.py
@@ -0,0 +1,106 @@
+from __future__ import division, absolute_import, print_function
+
+__copyright__ = "Copyright (C) 2017 Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+__doc__ = """
+.. autoclass:: CahnHilliardOperator
+"""
+
+import numpy as np
+from pytential.symbolic.pde.scalar import L2WeightedPDEOperator
+from pytential import sym
+from functools import partial
+
+
+class CahnHilliardOperator(L2WeightedPDEOperator):
+ def __init__(self, lambda1, lambda2, c):
+ self.lambdas = (lambda1, lambda2)
+ self.c = c
+
+ def make_unknown(self, name):
+ return sym.make_sym_vector(name, 2)
+
+ def S_G(self, i, density, qbx_forced_limit, op_map=None): # noqa: N802
+ if op_map is None:
+ op_map = lambda x: x # noqa: E731
+
+ from sumpy.kernel import HelmholtzKernel
+ hhk = HelmholtzKernel(2, allow_evanescent=True)
+ hhk_scaling = 1j/4
+
+ if i == 0:
+ lam1, lam2 = self.lambdas
+ return (
+ # FIXME: Verify scaling
+ -1/(2*np.pi*(lam1**2-lam2**2)) / hhk_scaling
+ *
+ (
+ op_map(sym.S(hhk, density, k=1j*lam1,
+ qbx_forced_limit=qbx_forced_limit))
+ -
+ op_map(sym.S(hhk, density, k=1j*lam2,
+ qbx_forced_limit=qbx_forced_limit))))
+ else:
+ return (
+ # FIXME: Verify scaling
+
+ -1/(2*np.pi) / hhk_scaling
+ * op_map(sym.S(hhk, density, k=1j*self.lambdas[i-1],
+ qbx_forced_limit=qbx_forced_limit)))
+
+ def representation(self, unknown):
+ sig1, sig2 = unknown
+ S_G = partial(self.S_G, qbx_forced_limit=None) # noqa: N806
+
+ return S_G(1, sig1) + S_G(0, sig2)
+
+ def operator(self, unknown):
+ sig1, sig2 = unknown
+ lam1, lam2 = self.lambdas
+ S_G = partial(self.S_G, qbx_forced_limit=1) # noqa: N806
+
+ c = self.c
+
+ def Sn_G(i, density): # noqa
+ return self.S_G(i, density,
+ qbx_forced_limit="avg",
+ op_map=partial(sym.normal_derivative, 2))
+
+ d = sym.make_obj_array([
+ 0.5*sig1,
+ 0.5*lam2**2*sig1 - 0.5*sig2
+ ])
+ a = sym.make_obj_array([
+ # A11
+ Sn_G(1, sig1) + c*S_G(1, sig1)
+ # A12
+ + Sn_G(0, sig2) + c*S_G(0, sig2),
+
+ # A21
+ lam2**2*Sn_G(1, sig1)
+ # A22
+ - Sn_G(1, sig2) + lam1**2*Sn_G(0, sig2)
+ ])
+
+ return d+a
diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py
index 57b3ff97cc0dcdb9469deb46b59a5b8ae154f14a..bcf9400885c3d45d4c6bc02b538c27902d4b5685 100644
--- a/pytential/symbolic/primitives.py
+++ b/pytential/symbolic/primitives.py
@@ -84,6 +84,12 @@ Elementary numerics
Calculus
^^^^^^^^
.. autoclass:: Derivative
+.. autofunction:: dd_axis
+.. autofunction:: d_dx
+.. autofunction:: d_dy
+.. autofunction:: d_dz
+.. autofunction:: grad_mv
+.. autofunction:: grad
Layer potentials
^^^^^^^^^^^^^^^^
@@ -452,6 +458,9 @@ class Derivative(DerivativeBase):
def dd_axis(axis, ambient_dim, operand):
+ """Return the derivative along (XYZ) axis *axis*
+ (in *ambient_dim*-dimensional space) of *operand*.
+ """
d = Derivative()
unit_vector = np.zeros(ambient_dim)
@@ -469,6 +478,23 @@ d_dy = partial(dd_axis, 1)
d_dz = partial(dd_axis, 2)
+def grad_mv(ambient_dim, operand):
+ """Return the *ambient_dim*-dimensional gradient of
+ *operand* as a :class:`pymbolic.geometric_algebra.MultiVector`.
+ """
+
+ d = Derivative()
+ return d.resolve(
+ d.dnabla(ambient_dim) * d(operand))
+
+
+def grad(ambient_dim, operand):
+ """Return the *ambient_dim*-dimensional gradient of
+ *operand* as a :class:`numpy.ndarray`.
+ """
+ return grad_mv(ambient_dim, operand).as_vector()
+
+
# {{{ potentials
def hashable_kernel_args(kernel_arguments):
diff --git a/requirements.txt b/requirements.txt
index 0c3cc71e0819605fe98f5d4bc63c44ad3fcfaff5..b8bf37d218305231a886bb54047db0d56d9d4bf3 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,10 +1,10 @@
numpy
git+git://github.com/inducer/pymbolic
sympy==1.0
-git+git://github.com/inducer/modepy
-git+git://github.com/pyopencl/pyopencl
-git+git://github.com/inducer/islpy
-git+git://github.com/inducer/loopy
-git+git://github.com/inducer/boxtree
-git+git://github.com/inducer/meshmode
-git+git://github.com/inducer/sumpy
+git+https://github.com/inducer/modepy
+git+https://github.com/pyopencl/pyopencl
+git+https://github.com/inducer/islpy
+git+https://github.com/inducer/loopy
+git+https://github.com/inducer/boxtree
+git+https://github.com/inducer/meshmode
+git+https://github.com/inducer/sumpy
diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py
index d7e514510d320602496a8b13a562a9da99733eba..7326efc13d212be30c83adecedc7a3e579722605 100644
--- a/test/test_global_qbx.py
+++ b/test/test_global_qbx.py
@@ -39,7 +39,7 @@ from pytential.qbx import QBXLayerPotentialSource
from functools import partial
from meshmode.mesh.generation import ( # noqa
ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut,
- make_curve_mesh)
+ make_curve_mesh, generate_icosphere, generate_torus)
from extra_curve_data import horseshoe
@@ -53,134 +53,161 @@ RNG_SEED = 10
FAR_TARGET_DIST_FROM_SOURCE = 10
-# {{{ utilities for iterating over panels
+# {{{ source refinement checker
class ElementInfo(RecordWithoutPickling):
"""
.. attribute:: element_nr
- .. attribute:: neighbors
.. attribute:: discr_slice
- .. attribute:: mesh_slice
- .. attribute:: element_group
- .. attribute:: mesh_element_group
"""
__slots__ = ["element_nr",
- "neighbors",
"discr_slice"]
def iter_elements(discr):
discr_nodes_idx = 0
element_nr = 0
- adjacency = discr.mesh.nodal_adjacency
for discr_group in discr.groups:
start = element_nr
for element_nr in range(start, start + discr_group.nelements):
yield ElementInfo(
element_nr=element_nr,
- neighbors=list(adjacency.neighbors[
- slice(*adjacency.neighbors_starts[
- element_nr:element_nr+2])]),
discr_slice=slice(discr_nodes_idx,
discr_nodes_idx + discr_group.nunit_nodes))
discr_nodes_idx += discr_group.nunit_nodes
-# }}}
-
-@pytest.mark.parametrize(("curve_name", "curve_f", "nelements"), [
- ("20-to-1 ellipse", partial(ellipse, 20), 100),
- ("horseshoe", horseshoe, 64),
- ])
-def test_source_refinement(ctx_getter, curve_name, curve_f, nelements):
+def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None):
cl_ctx = ctx_getter()
queue = cl.CommandQueue(cl_ctx)
- # {{{ generate lpot source, run refiner
-
- order = 16
- helmholtz_k = 10
-
- mesh = make_curve_mesh(curve_f, np.linspace(0, 1, nelements+1), order)
-
from meshmode.discretization import Discretization
- from meshmode.discretization.poly_element import \
- InterpolatoryQuadratureSimplexGroupFactory
+ from meshmode.discretization.poly_element import (
+ InterpolatoryQuadratureSimplexGroupFactory,
+ QuadratureSimplexGroupFactory)
+
factory = InterpolatoryQuadratureSimplexGroupFactory(order)
+ fine_factory = QuadratureSimplexGroupFactory(4 * order)
discr = Discretization(cl_ctx, mesh, factory)
- from pytential.qbx.refinement import QBXLayerPotentialSourceRefiner
+ from pytential.qbx.refinement import (
+ RefinerCodeContainer, refine_for_global_qbx)
lpot_source = QBXLayerPotentialSource(discr, order)
del discr
- refiner = QBXLayerPotentialSourceRefiner(cl_ctx)
- lpot_source, conn = refiner(lpot_source, factory, helmholtz_k)
+ refiner_extra_kwargs = {}
+ if helmholtz_k is not None:
+ refiner_extra_kwargs["kernel_length_scale"] = 5/helmholtz_k
+
+ lpot_source, conn = refine_for_global_qbx(
+ lpot_source, RefinerCodeContainer(cl_ctx),
+ factory, fine_factory, **refiner_extra_kwargs)
+
+ from pytential.qbx.utils import get_centers_on_side
discr_nodes = lpot_source.density_discr.nodes().get(queue)
- int_centers = lpot_source.centers(-1)
+ fine_discr_nodes = lpot_source.fine_density_discr.nodes().get(queue)
+ int_centers = get_centers_on_side(lpot_source, -1)
int_centers = np.array([axis.get(queue) for axis in int_centers])
- ext_centers = lpot_source.centers(+1)
+ ext_centers = get_centers_on_side(lpot_source, +1)
ext_centers = np.array([axis.get(queue) for axis in ext_centers])
- panel_sizes = lpot_source.panel_sizes("npanels").get(queue)
-
- # }}}
+ expansion_radii = lpot_source._expansion_radii("npanels").get(queue)
+ panel_sizes = lpot_source._panel_sizes("npanels").get(queue)
+ fine_panel_sizes = lpot_source._fine_panel_sizes("npanels").get(queue)
# {{{ check if satisfying criteria
- def check_panel(panel):
- # Check 2-to-1 panel to neighbor size ratio.
- for neighbor in panel.neighbors:
- assert panel_sizes[panel.element_nr] / panel_sizes[neighbor] <= 2, \
- (panel_sizes[panel.element_nr], panel_sizes[neighbor])
-
- # Check wavenumber to panel size ratio.
- assert panel_sizes[panel.element_nr] * helmholtz_k <= 5
-
- def check_panel_pair(panel_1, panel_2):
- h_1 = panel_sizes[panel_1.element_nr]
- h_2 = panel_sizes[panel_2.element_nr]
-
- if panel_1.element_nr == panel_2.element_nr:
+ def check_disk_undisturbed_by_sources(centers_panel, sources_panel):
+ if centers_panel.element_nr == sources_panel.element_nr:
# Same panel
return
- panel_1_centers = int_centers[:, panel_1.discr_slice]
- panel_2_nodes = discr_nodes[:, panel_2.discr_slice]
+ my_int_centers = int_centers[:, centers_panel.discr_slice]
+ my_ext_centers = ext_centers[:, centers_panel.discr_slice]
+ all_centers = np.append(my_int_centers, my_ext_centers, axis=-1)
+
+ nodes = discr_nodes[:, sources_panel.discr_slice]
# =distance(centers of panel 1, panel 2)
dist = (
la.norm((
- panel_1_centers[..., np.newaxis] -
- panel_2_nodes[:, np.newaxis, ...]).T,
+ all_centers[..., np.newaxis] -
+ nodes[:, np.newaxis, ...]).T,
axis=-1)
.min())
- # Criterion 1:
+ # Criterion:
# A center cannot be closer to another panel than to its originating
# panel.
- assert dist >= h_1 / 2, (dist, h_1, panel_1.element_nr, panel_2.element_nr)
+ rad = expansion_radii[centers_panel.element_nr]
+ assert dist >= rad, \
+ (dist, rad, centers_panel.element_nr, sources_panel.element_nr)
- # Criterion 2:
- # A center cannot be closer to another panel than that panel's
- # centers - unless the panels are adjacent, to allow for refinement.
+ def check_sufficient_quadrature_resolution(centers_panel, sources_panel):
+ h = fine_panel_sizes[sources_panel.element_nr]
- if panel_2.element_nr in panel_1.neighbors:
- return
+ my_int_centers = int_centers[:, centers_panel.discr_slice]
+ my_ext_centers = ext_centers[:, centers_panel.discr_slice]
+ all_centers = np.append(my_int_centers, my_ext_centers, axis=-1)
- assert dist >= h_2 / 2, (dist, h_2, panel_1.element_nr, panel_2.element_nr)
+ nodes = fine_discr_nodes[:, sources_panel.discr_slice]
- for panel_1 in iter_elements(lpot_source.density_discr):
- check_panel(panel_1)
+ # =distance(centers of panel 1, panel 2)
+ dist = (
+ la.norm((
+ all_centers[..., np.newaxis] -
+ nodes[:, np.newaxis, ...]).T,
+ axis=-1)
+ .min())
+
+ # Criterion:
+ # The quadrature contribution from each panel is as accurate
+ # as from the center's own source panel.
+ assert dist >= h / 4, \
+ (dist, h, centers_panel.element_nr, sources_panel.element_nr)
+
+ def check_panel_size_to_helmholtz_k_ratio(panel):
+ # Check wavenumber to panel size ratio.
+ assert panel_sizes[panel.element_nr] * helmholtz_k <= 5
+
+ for i, panel_1 in enumerate(iter_elements(lpot_source.density_discr)):
for panel_2 in iter_elements(lpot_source.density_discr):
- check_panel_pair(panel_1, panel_2)
+ check_disk_undisturbed_by_sources(panel_1, panel_2)
+ for panel_2 in iter_elements(lpot_source.fine_density_discr):
+ check_sufficient_quadrature_resolution(panel_1, panel_2)
+ if helmholtz_k is not None:
+ check_panel_size_to_helmholtz_k_ratio(panel_1)
# }}}
+# }}}
+
+
+@pytest.mark.parametrize(("curve_name", "curve_f", "nelements"), [
+ ("20-to-1 ellipse", partial(ellipse, 20), 100),
+ ("horseshoe", horseshoe, 64),
+ ])
+def test_source_refinement_2d(ctx_getter, curve_name, curve_f, nelements):
+ helmholtz_k = 10
+ order = 8
+
+ mesh = make_curve_mesh(curve_f, np.linspace(0, 1, nelements+1), order)
+ run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k)
+
+
+@pytest.mark.parametrize(("surface_name", "surface_f", "order"), [
+ ("sphere", partial(generate_icosphere, 1), 4),
+ ("torus", partial(generate_torus, 3, 1, n_inner=10, n_outer=7), 6),
+ ])
+def test_source_refinement_3d(ctx_getter, surface_name, surface_f, order):
+ mesh = surface_f(order=order)
+ run_source_refinement_test(ctx_getter, mesh, order)
+
@pytest.mark.parametrize(("curve_name", "curve_f", "nelements"), [
("20-to-1 ellipse", partial(ellipse, 20), 100),
@@ -204,16 +231,13 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements):
discr = Discretization(cl_ctx, mesh, factory)
- from pytential.qbx.refinement import QBXLayerPotentialSourceRefiner
-
- lpot_source = QBXLayerPotentialSource(discr, order)
+ lpot_source, conn = QBXLayerPotentialSource(discr, order).with_refinement()
del discr
- refiner = QBXLayerPotentialSourceRefiner(cl_ctx)
- lpot_source, conn = refiner(lpot_source, factory)
+ from pytential.qbx.utils import get_centers_on_side
- int_centers = lpot_source.centers(-1)
- ext_centers = lpot_source.centers(+1)
+ int_centers = get_centers_on_side(lpot_source, -1)
+ ext_centers = get_centers_on_side(lpot_source, +1)
# }}}
@@ -223,7 +247,8 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements):
rng = PhiloxGenerator(cl_ctx, seed=RNG_SEED)
nsources = lpot_source.density_discr.nnodes
noise = rng.uniform(queue, nsources, dtype=np.float, a=0.01, b=1.0)
- panel_sizes = lpot_source.panel_sizes("nsources").with_queue(queue)
+ tunnel_radius = \
+ lpot_source._close_target_tunnel_radius("nsources").with_queue(queue)
def targets_from_sources(sign, dist):
from pytential import sym, bind
@@ -234,8 +259,8 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements):
from pytential.target import PointsTarget
- int_targets = PointsTarget(targets_from_sources(-1, noise * panel_sizes / 2))
- ext_targets = PointsTarget(targets_from_sources(+1, noise * panel_sizes / 2))
+ int_targets = PointsTarget(targets_from_sources(-1, noise * tunnel_radius))
+ ext_targets = PointsTarget(targets_from_sources(+1, noise * tunnel_radius))
far_targets = PointsTarget(targets_from_sources(+1, FAR_TARGET_DIST_FROM_SOURCE))
# Create target discretizations.
@@ -270,7 +295,7 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements):
QBXTargetAssociator(cl_ctx)(lpot_source, target_discrs)
.get(queue=queue))
- panel_sizes = lpot_source.panel_sizes("nsources").get(queue)
+ expansion_radii = lpot_source._expansion_radii("nsources").get(queue)
int_centers = np.array([axis.get(queue) for axis in int_centers])
ext_centers = np.array([axis.get(queue) for axis in ext_centers])
@@ -290,7 +315,7 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements):
target_to_source_result, target_to_side_result):
assert (target_to_side_result == true_side).all()
dists = la.norm((targets.T - centers.T[target_to_source_result]), axis=1)
- assert (dists <= panel_sizes[target_to_source_result] / 2).all()
+ assert (dists <= expansion_radii[target_to_source_result]).all()
# Checks that far targets are not assigned a center.
def check_far_targets(target_to_source_result):
diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py
index c67b2e5c5012d53242fa245e1128c1798cebe8c3..6e8a81000b31958f228615c17c9d0919221a7d41 100644
--- a/test/test_layer_pot.py
+++ b/test/test_layer_pot.py
@@ -34,7 +34,7 @@ from pyopencl.tools import ( # noqa
from functools import partial
from meshmode.mesh.generation import ( # noqa
- ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut,
+ ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, WobblyCircle,
make_curve_mesh)
from sumpy.visualization import FieldPlotter
from pytential import bind, sym, norm
@@ -160,7 +160,10 @@ def test_ellipse_eigenvalues(ctx_getter, ellipse_aspect, mode_nr, qbx_order):
if 0:
# plot geometry, centers, normals
- centers = qbx.centers(density_discr, 1)
+
+ from pytential.qbx.utils import get_centers_on_side
+ centers = get_centers_on_side(qbx, 1)
+
nodes_h = nodes.get()
centers_h = [centers[0].get(), centers[1].get()]
pt.plot(nodes_h[0], nodes_h[1], "x-")
@@ -309,10 +312,15 @@ def run_int_eq_test(
if source_order is None:
source_order = 4*target_order
+ refiner_extra_kwargs = {}
+
+ if k != 0:
+ refiner_extra_kwargs["kernel_length_scale"] = 5/k
+
qbx, _ = QBXLayerPotentialSource(
pre_density_discr, fine_order=source_order, qbx_order=qbx_order,
# Don't use FMM for now
- fmm_order=False).with_refinement()
+ fmm_order=False).with_refinement(**refiner_extra_kwargs)
density_discr = qbx.density_discr
@@ -694,6 +702,13 @@ def get_starfish_mesh(refinement_increment, target_order):
target_order)
+def get_wobbly_circle_mesh(refinement_increment, target_order):
+ nelements = [3000, 5000, 7000][refinement_increment]
+ return make_curve_mesh(WobblyCircle.random(30, seed=30),
+ np.linspace(0, 1, nelements+1),
+ target_order)
+
+
def get_sphere_mesh(refinement_increment, target_order):
from meshmode.mesh.generation import generate_icosphere
mesh = generate_icosphere(1, target_order)
@@ -801,11 +816,17 @@ def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order,
# FIXME: FMM kernel generation slow
direct_eval = (k != 0)
+ refiner_extra_kwargs = {}
+
+ if k != 0:
+ refiner_extra_kwargs["kernel_length_scale"] = 5/k
+
qbx, _ = QBXLayerPotentialSource(
pre_density_discr, 4*target_order,
qbx_order, fmm_order=(
False if direct_eval else qbx_order + order_bump)
- ).with_refinement()
+ ).with_refinement(**refiner_extra_kwargs)
+
density_discr = qbx.density_discr
# {{{ compute values of a solution to the PDE
@@ -903,8 +924,12 @@ def test_off_surface_eval(ctx_getter, use_fmm, do_plot=False):
pre_density_discr = Discretization(
cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order))
- qbx, _ = QBXLayerPotentialSource(pre_density_discr, 4*target_order, qbx_order,
- fmm_order=fmm_order).with_refinement()
+ qbx, _ = QBXLayerPotentialSource(
+ pre_density_discr,
+ 4*target_order,
+ qbx_order,
+ fmm_order=fmm_order,
+ ).with_refinement()
density_discr = qbx.density_discr
@@ -919,10 +944,11 @@ def test_off_surface_eval(ctx_getter, use_fmm, do_plot=False):
(qbx, PointsTarget(fplot.points)),
op)(queue, sigma=sigma)
- print(fld_in_vol)
-
err = cl.clmath.fabs(fld_in_vol - (-1))
+ linf_err = cl.array.max(err).get()
+ print("l_inf error:", linf_err)
+
if do_plot:
fplot.show_scalar_in_matplotlib(fld_in_vol.get())
import matplotlib.pyplot as pt
@@ -930,7 +956,86 @@ def test_off_surface_eval(ctx_getter, use_fmm, do_plot=False):
pt.show()
# FIXME: Why does the FMM only meet this sloppy tolerance?
- assert (err < 1e-2).get().all()
+ assert linf_err < 1e-2
+
+# }}}
+
+
+# {{{ test off-surface eval vs direct
+
+def test_off_surface_eval_vs_direct(ctx_getter, do_plot=False):
+ logging.basicConfig(level=logging.INFO)
+
+ cl_ctx = ctx_getter()
+ queue = cl.CommandQueue(cl_ctx)
+
+ # prevent cache 'splosion
+ from sympy.core.cache import clear_cache
+ clear_cache()
+
+ nelements = 300
+ target_order = 8
+ qbx_order = 3
+
+ mesh = make_curve_mesh(WobblyCircle.random(8, seed=30),
+ 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))
+ direct_qbx, _ = QBXLayerPotentialSource(
+ pre_density_discr, 4*target_order, qbx_order,
+ fmm_order=False,
+ target_stick_out_factor=0.05,
+ ).with_refinement()
+ fmm_qbx, _ = QBXLayerPotentialSource(
+ pre_density_discr, 4*target_order, qbx_order,
+ fmm_order=qbx_order + 3,
+ expansion_disks_in_tree_have_extent=True,
+ target_stick_out_factor=0.05,
+ ).with_refinement()
+
+ fplot = FieldPlotter(np.zeros(2), extent=5, npoints=1000)
+ from pytential.target import PointsTarget
+ ptarget = PointsTarget(fplot.points)
+ from sumpy.kernel import LaplaceKernel
+
+ op = sym.D(LaplaceKernel(2), sym.var("sigma"), qbx_forced_limit=None)
+
+ from pytential.qbx import QBXTargetAssociationFailedException
+ try:
+ direct_density_discr = direct_qbx.density_discr
+ direct_sigma = direct_density_discr.zeros(queue) + 1
+ direct_fld_in_vol = bind((direct_qbx, ptarget), op)(
+ queue, sigma=direct_sigma)
+
+ except QBXTargetAssociationFailedException as e:
+ fplot.show_scalar_in_matplotlib(e.failed_target_flags.get(queue))
+ import matplotlib.pyplot as pt
+ pt.show()
+ raise
+
+ fmm_density_discr = fmm_qbx.density_discr
+ fmm_sigma = fmm_density_discr.zeros(queue) + 1
+ fmm_fld_in_vol = bind((fmm_qbx, ptarget), op)(queue, sigma=fmm_sigma)
+
+ err = cl.clmath.fabs(fmm_fld_in_vol - direct_fld_in_vol)
+
+ linf_err = cl.array.max(err).get()
+ print("l_inf error:", linf_err)
+
+ if do_plot:
+ #fplot.show_scalar_in_mayavi(0.1*cl.clmath.log10(1e-15 + err).get(queue))
+ fplot.show_scalar_in_mayavi(fmm_fld_in_vol.get(queue))
+ import mayavi.mlab as mlab
+ mlab.show()
+
+ assert linf_err < 1e-3
# }}}
diff --git a/test/test_muller.py b/test/test_muller.py
index 30d1a47d8564f607a125bdfd0df5f6e67f811606..fb23e911d32c1c720a3284004014cb83088d75d7 100644
--- a/test/test_muller.py
+++ b/test/test_muller.py
@@ -37,10 +37,14 @@ def test_muller(true_roots):
:arg n: number of zeros sought
:return: (roots, niter)
"""
+ np.random.seed(15)
+ z_start = np.random.rand(3) + 1j*np.random.rand(3)
+
eps = 1e-12
from pytential.muller import muller_deflate
roots, niter = muller_deflate(
- lambda z: poly_with_roots(z, true_roots), len(true_roots))
+ lambda z: poly_with_roots(z, true_roots), len(true_roots),
+ z_start=z_start)
for r_i in roots:
min_dist, true_root = min(