fmmlib unhappy about empty interaction list
When trying to evaluate a layer potential using fmm_backend='fmmlib', if the sibling interaction lists are empty I get an error.
For example, this error appears when evaluating a single layer potential over a unit sphere at target points on an enclosing cube, both of which have rather coarse meshes (see attached files enclosed_sphere.msh and enclosing_cube.msh). Here is the file and error produced (created using reprexpy):
import numpy as np
import pyopencl as cl
from meshmode.mesh.io import read_gmsh
from meshmode.discretization import Discretization
from meshmode.discretization.poly_element import \
InterpolatoryQuadratureSimplexGroupFactory
cl_ctx = cl.create_some_context()
queue = cl.CommandQueue(cl_ctx)
# Get target
from pytential.target import PointsTarget
target_mesh = read_gmsh('enclosing_cube.msh')
target = PointsTarget(target_mesh.vertices)
# Read mesh and make source discretization
source_mesh = read_gmsh('enclosed_sphere.msh')
ambient_dim = 3
bdry_quad_order = 1
fine_order = 1
qbx_order = 1
fmm_order = 3
pre_density_discr = Discretization(
cl_ctx,
source_mesh,
InterpolatoryQuadratureSimplexGroupFactory(bdry_quad_order))
from pytential.qbx import QBXLayerPotentialSource
# Make qbx -- nb that this bug only appears when fmm_backend='fmmlib' is set
qbx, _ = QBXLayerPotentialSource(
pre_density_discr, fine_order=fine_order, qbx_order=qbx_order,
fmm_order=fmm_order, fmm_backend='fmmlib',
).with_refinement()
#> /home/bensepan/py_proj/inducer/pytential/pytential/symbolic/compiler.py:540: UserWarning: mishandling CSE scope
#> warn("mishandling CSE scope")
density_discr = qbx.density_discr
from pytential import sym, bind
from sumpy.kernel import LaplaceKernel
# Bind layer potential operator
op = sym.S(LaplaceKernel(ambient_dim),
sym.var("sigma"),
qbx_forced_limit=None)
bound_op = bind((qbx, target), op)
sigma = cl.array.to_device(queue, np.zeros(density_discr.nnodes))
bound_op(queue, sigma=sigma)
#> Traceback (most recent call last):
#> ~/py_proj/firedrake-complex/lib/python3.6/site-packages/pytools/__init__.py in wrapper(obj, *args, **kwargs)
#> 581 try:
#> --> 582 return getattr(obj, cache_dict_name)[key]
#> 583 except AttributeError:
#> AttributeError: 'QBXFMMLibExpansionWrangler' object has no attribute '_memoize_dic_m2l_rotation_matrices'
#> During handling of the above exception, another exception occurred:
#> Traceback (most recent call last):
#> <ipython-input-28-b5ab809e6f4a> in <module>
#> ----> 1 bound_op(queue, sigma=sigma)
#> ~/py_proj/inducer/pytential/pytential/symbolic/execution.py in __call__(self, queue, **args)
#> 729 :class:`pyopencl.array.Array`, or an object array of these.
#> 730 """
#> --> 731 return self.eval(queue, args)
#> 732
#> 733
#> ~/py_proj/inducer/pytential/pytential/symbolic/execution.py in eval(self, queue, context, timing_data)
#> 719 exec_mapper = EvaluationMapper(
#> 720 self, queue, context, timing_data=timing_data)
#> --> 721 return self.code.execute(exec_mapper)
#> 722
#> 723 def __call__(self, queue, **args):
#> ~/py_proj/inducer/pytential/pytential/symbolic/compiler.py in execute(self, exec_mapper, pre_assign_check)
#> 390 self.get_exec_function(insn, exec_mapper)
#> 391 (exec_mapper.queue, insn, exec_mapper.bound_expr,
#> --> 392 exec_mapper))
#> 393
#> 394 assignees = insn.get_assignees()
#> ~/py_proj/inducer/pytential/pytential/symbolic/execution.py in exec_compute_potential_insn(self, queue, insn, bound_expr, evaluate)
#> 312 result, timing_data = (
#> 313 source.exec_compute_potential_insn(
#> --> 314 queue, insn, bound_expr, evaluate, return_timing_data))
#> 315
#> 316 if return_timing_data:
#> ~/py_proj/inducer/pytential/pytential/qbx/__init__.py in exec_compute_potential_insn(self, queue, insn, bound_expr, evaluate, return_timing_data)
#> 558
#> 559 return self._dispatch_compute_potential_insn(
#> --> 560 queue, insn, bound_expr, evaluate, func, extra_args)
#> 561
#> 562 def cost_model_compute_potential_insn(self, queue, insn, bound_expr, evaluate):
#> ~/py_proj/inducer/pytential/pytential/qbx/__init__.py in _dispatch_compute_potential_insn(self, queue, insn, bound_expr, evaluate, func, extra_args)
#> 602 extra_args = {}
#> 603
#> --> 604 return func(queue, insn, bound_expr, evaluate, **extra_args)
#> 605
#> 606 @property
#> ~/py_proj/inducer/pytential/pytential/qbx/__init__.py in exec_compute_potential_insn_fmm(self, queue, insn, bound_expr, evaluate, fmm_driver)
#> 736 all_potentials_on_every_target, extra_outputs = (
#> 737 fmm_driver(
#> --> 738 wrangler, strengths, geo_data, fmm_kernel, kernel_extra_kwargs))
#> 739
#> 740 result = []
#> ~/py_proj/inducer/pytential/pytential/qbx/__init__.py in drive_fmm(***failed resolving arguments***)
#> 553 else:
#> 554 timing_data = None
#> --> 555 return drive_fmm(wrangler, strengths, timing_data), timing_data
#> 556
#> 557 extra_args["fmm_driver"] = drive_fmm
#> ~/py_proj/inducer/pytential/pytential/qbx/fmm.py in drive_fmm(expansion_wrangler, src_weights, timing_data, traversal)
#> 463 traversal.from_sep_siblings_starts,
#> 464 traversal.from_sep_siblings_lists,
#> --> 465 mpole_exps)
#> 466
#> 467 recorder.add("multipole_to_local", timing_future)
#> ~/py_proj/firedrake-complex/lib/python3.6/site-packages/pytools/__init__.py in wrapper(*args, **kwargs)
#> 2266 self.logger,
#> 2267 self.description or wrapped.__name__):
#> -> 2268 return wrapped(*args, **kwargs)
#> 2269
#> 2270 from functools import update_wrapper
#> ~/py_proj/inducer/boxtree/boxtree/tools.py in wrapper(*args, **kwargs)
#> 561 def wrapper(*args, **kwargs):
#> 562 timer = ProcessTimer()
#> --> 563 retval = wrapped(*args, **kwargs)
#> 564 timer.done()
#> 565
#> ~/py_proj/inducer/boxtree/boxtree/pyfmmlib_integration.py in multipole_to_local(self, level_start_target_or_target_parent_box_nrs, target_or_target_parent_boxes, starts, lists, mpole_exps)
#> 727 # Precomputed rotation matrices (matrices of larger order can be used
#> 728 # for translations of smaller order)
#> --> 729 rotmatf, rotmatb, rotmat_order = self.m2l_rotation_matrices()
#> 730
#> 731 for lev in range(self.tree.nlevels):
#> ~/py_proj/firedrake-complex/lib/python3.6/site-packages/pytools/__init__.py in wrapper(obj, *args, **kwargs)
#> 582 return getattr(obj, cache_dict_name)[key]
#> 583 except AttributeError:
#> --> 584 result = function(obj, *args, **kwargs)
#> 585 setattr(obj, cache_dict_name, {key: result})
#> 586 return result
#> ~/py_proj/inducer/boxtree/boxtree/pyfmmlib_integration.py in m2l_rotation_matrices(self)
#> 705
#> 706 ier, rotmatf = (
#> --> 707 rotmat_builder(rotmat_order, m2l_rotation_angles))
#> 708 assert (0 == ier).all()
#> 709
#> ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (0,)Created on 2020-02-25 by the reprexpy package
I've attached the python script fmmlib_empty_interaction_list_reprex.py and traceback in case you find this format more helpful.