diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 7caff1a776467aba8abdca40f3d2ef9857ffbe74..3d98c18b9ba244574b14541f4c1fefbdc965f6c5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -31,6 +31,7 @@ Python 3.5 Conda: - 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" + allow_failure: true # takes very long tags: - linux except: diff --git a/.test-conda-env-py3-requirements.txt b/.test-conda-env-py3-requirements.txt index e3fac95dc0e8084cd5f7d50ed120d53a5655797a..9b143e4367f52301caa436f0d5a2304909a24c1a 100644 --- a/.test-conda-env-py3-requirements.txt +++ b/.test-conda-env-py3-requirements.txt @@ -1,5 +1,5 @@ git+https://gitlab.tiker.net/inducer/boxtree git+https://github.com/inducer/pymbolic git+https://github.com/inducer/loopy -git+https://github.com/inducer/sumpy +git+https://gitlab.tiker.net/inducer/sumpy@rscale git+https://github.com/inducer/meshmode diff --git a/examples/fmm-manyarm-starfish-nan.py b/examples/fmm-manyarm-starfish-nan.py deleted file mode 100644 index 7a8b1d3f58148c03f71782419bb792a22c4e9148..0000000000000000000000000000000000000000 --- a/examples/fmm-manyarm-starfish-nan.py +++ /dev/null @@ -1,109 +0,0 @@ -import numpy as np -import numpy.linalg as la # noqa -import pyopencl as cl -import pyopencl.clmath # noqa - -from meshmode.mesh.generation import ( # noqa - ellipse, cloverleaf, NArmedStarfish, drop, n_gon, qbx_peanut, - make_curve_mesh, starfish) -from pytential import bind, sym, norm - -import logging -logger = logging.getLogger(__name__) - - -def get_starfish_mesh(nelements, target_order): - return make_curve_mesh( - NArmedStarfish(20, 0.8), - np.linspace(0, 1, nelements+1), - target_order) - - -WITH_EXTENTS = True -EXPANSION_STICKOUT_FACTOR = 0.5 - - -def get_green_error( - queue, mesh_getter, nelements, fmm_order, qbx_order, k=0): - - target_order = 12 - - mesh = mesh_getter(nelements, target_order) - - d = mesh.ambient_dim - - # u_sym = sym.var("u") - dn_u_sym = sym.var("dn_u") - - from sumpy.kernel import LaplaceKernel - k_sym = LaplaceKernel(d) - zero_op = ( - sym.S(k_sym, dn_u_sym, qbx_forced_limit=-1) - # - sym.D(k_sym, u_sym, qbx_forced_limit="avg") - # - 0.5*u_sym - ) - - from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import ( - InterpolatoryQuadratureSimplexGroupFactory) - pre_density_discr = Discretization( - queue.context, mesh, - InterpolatoryQuadratureSimplexGroupFactory(target_order)) - - from pytential.qbx import QBXLayerPotentialSource - lpot_source = QBXLayerPotentialSource( - pre_density_discr, 4*target_order, - qbx_order, fmm_order=fmm_order, - _expansions_in_tree_have_extent=WITH_EXTENTS, - _expansion_stick_out_factor=EXPANSION_STICKOUT_FACTOR, - ) - - lpot_source, _ = lpot_source.with_refinement() - - density_discr = lpot_source.density_discr - - # {{{ compute values of a solution to the PDE - - nodes_host = density_discr.nodes().get(queue) - normal = bind(density_discr, sym.normal(d))(queue).as_vector(np.object) - normal_host = [normal[j].get() for j in range(d)] - - center = np.array([3, 1, 2])[:d] - diff = nodes_host - center[:, np.newaxis] - dist_squared = np.sum(diff**2, axis=0) - dist = np.sqrt(dist_squared) - if d == 2: - u = np.log(dist) - grad_u = diff/dist_squared - elif d == 3: - u = 1/dist - grad_u = -diff/dist**3 - else: - assert False - - dn_u = 0 - for i in range(d): - dn_u = dn_u + normal_host[i]*grad_u[i] - - # }}} - - u_dev = cl.array.to_device(queue, u) - dn_u_dev = cl.array.to_device(queue, dn_u) - - bound_op = bind(lpot_source, zero_op) - error = bound_op( - queue, u=u_dev, dn_u=dn_u_dev, k=k) - - print(len(error), np.where(np.isnan(error.get()))) - return norm(density_discr, queue, error) - - -def main(): - cl_ctx = cl.create_some_context() - queue = cl.CommandQueue(cl_ctx) - - get_green_error(queue, get_starfish_mesh, 500, 20, 2, k=0) - - -if __name__ == "__main__": - main() diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index ee716f53773e5372372931891268da839241fb20..a59d0fb6673989ff6737a9e184fc81d15e94b4ac 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -696,7 +696,9 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): queue, target_discr.nodes(), self.quad_stage2_density_discr.nodes(), utils.get_centers_on_side(self, o.qbx_forced_limit), - [strengths], **kernel_args) + [strengths], + expansion_radii=self._expansion_radii("nsources"), + **kernel_args) result.append((o.name, output_for_each_kernel[o.kernel_index])) else: # no on-disk kernel caching @@ -759,6 +761,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): targets=target_discr.nodes(), sources=self.quad_stage2_density_discr.nodes(), centers=geo_data.centers(), + expansion_radii=geo_data.expansion_radii(), strengths=[strengths], qbx_tgt_numbers=qbx_tgt_numbers, qbx_center_numbers=qbx_center_numbers, diff --git a/pytential/qbx/direct.py b/pytential/qbx/direct.py index b8ab7d249de9fd30cccd72aee0115318f052e426..6dc5cd9abbb7319d0cd7a4029a3a2b22b6a710e5 100644 --- a/pytential/qbx/direct.py +++ b/pytential/qbx/direct.py @@ -39,6 +39,7 @@ class LayerPotentialOnTargetAndCenterSubset(LayerPotentialBase): <> a[idim] = center[idim,icenter] - src[idim,isrc] {id=compute_a} <> b[idim] = tgt[idim,itgt_overall] - center[idim,icenter] \ {id=compute_b} + <> rscale = expansion_radii[icenter] end """ @@ -50,6 +51,7 @@ class LayerPotentialOnTargetAndCenterSubset(LayerPotentialBase): shape=(self.dim, "ntargets_total"), order="C"), lp.GlobalArg("center", None, shape=(self.dim, "ncenters_total"), order="C"), + lp.GlobalArg("expansion_radii", None, shape="ncenters_total"), lp.GlobalArg("qbx_tgt_numbers", None, shape="ntargets"), lp.GlobalArg("qbx_center_numbers", None, shape="ntargets"), lp.ValueArg("nsources", np.int32), diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index d03affe3504b0c7be81f9f739c0ce8c393c49f53..caa70c61a555fcfb44be55c2ed03990f12cf53a5 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -27,7 +27,8 @@ from six.moves import range, zip import numpy as np # noqa import pyopencl as cl # noqa import pyopencl.array # noqa -from sumpy.fmm import SumpyExpansionWranglerCodeContainer, SumpyExpansionWrangler +from sumpy.fmm import (SumpyExpansionWranglerCodeContainer, + SumpyExpansionWrangler, level_to_rscale) from pytools import memoize_method from pytential.qbx.interactions import P2QBXLFromCSR, M2QBXL, L2QBXL, QBXL2P @@ -59,7 +60,7 @@ class QBXSumpyExpansionWranglerCodeContainer(SumpyExpansionWranglerCodeContainer @memoize_method def qbx_local_expansion(self, order): - return self.qbx_local_expansion_factory(order) + return self.qbx_local_expansion_factory(order, self.use_rscale) @memoize_method def p2qbxl(self, order): @@ -213,6 +214,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), global_qbx_centers=geo_data.global_qbx_centers(), qbx_center_to_target_box=geo_data.qbx_center_to_target_box(), qbx_centers=geo_data.centers(), + qbx_expansion_radii=geo_data.expansion_radii(), source_box_starts=starts, source_box_lists=lists, @@ -250,6 +252,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), centers=self.tree.box_centers, qbx_centers=geo_data.centers(), + qbx_expansion_radii=geo_data.expansion_radii(), src_expansions=source_mpoles_view, src_base_ibox=source_level_start_ibox, @@ -258,6 +261,8 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), src_box_starts=ssn.starts, src_box_lists=ssn.lists, + src_rscale=level_to_rscale(self.tree, isrc_level), + wait_for=wait_for, **self.kernel_extra_kwargs) @@ -295,10 +300,13 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), centers=self.tree.box_centers, qbx_centers=geo_data.centers(), + qbx_expansion_radii=geo_data.expansion_radii(), expansions=target_locals_view, qbx_expansions=qbx_expansions, + src_rscale=level_to_rscale(self.tree, isrc_level), + wait_for=wait_for, **self.kernel_extra_kwargs) @@ -323,6 +331,8 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), evt, pot_res = qbxl2p(self.queue, qbx_centers=geo_data.centers(), + qbx_expansion_radii=geo_data.expansion_radii(), + global_qbx_centers=geo_data.global_qbx_centers(), center_to_targets_starts=ctt.starts, diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index 25ecbed94c54ea342a9bbf9e287a3ec4a82961a2..9f1d32407bb9309aa36b3eefc266f658818bee09 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -26,8 +26,9 @@ import numpy as np from pytools import memoize_method, Record import pyopencl as cl # noqa import pyopencl.array # noqa: F401 -from boxtree.pyfmmlib_integration import FMMLibExpansionWrangler -from sumpy.kernel import HelmholtzKernel +from boxtree.pyfmmlib_integration import ( + FMMLibExpansionWrangler, level_to_rscale) +from sumpy.kernel import LaplaceKernel, HelmholtzKernel import logging @@ -146,11 +147,16 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): else: source_deriv_name = None - result = isinstance(knl, HelmholtzKernel) and knl.dim == 3 - if result: + if isinstance(knl, HelmholtzKernel) and knl.dim in [2, 3]: k_names.append(knl.helmholtz_k_name) source_deriv_names.append(source_deriv_name) - return result + return True + elif isinstance(knl, LaplaceKernel) and knl.dim in [2, 3]: + k_names.append(None) + source_deriv_names.append(source_deriv_name) + return True + + return False ifgrad = False outputs = [] @@ -163,8 +169,8 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): ifgrad = True else: raise NotImplementedError( - "only the 3D Helmholtz kernel and its target derivatives " - "are supported for now") + "only the 2/3D Laplace and Helmholtz kernel " + "and their derivatives are supported") from pytools import is_single_valued if not is_single_valued(source_deriv_names): @@ -178,7 +184,10 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): from pytools import single_valued k_name = single_valued(k_names) - helmholtz_k = kernel_extra_kwargs[k_name] + if k_name is None: + helmholtz_k = 0 + else: + helmholtz_k = kernel_extra_kwargs[k_name] self.level_orders = [ fmm_level_to_order(level) @@ -333,10 +342,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): isource += ns centers = qbx_centers[:, geo_data.global_qbx_centers()] - - rscale = 1 # FIXME - rscale_vec = np.empty(len(center_source_counts) - 1, dtype=np.float64) - rscale_vec.fill(rscale) # FIXME + rscale_vec = geo_data.expansion_radii()[geo_data.global_qbx_centers()] nsources_vec = np.ones(self.tree.nsources, np.int32) @@ -372,13 +378,18 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): kwargs["charge_starts"] = info.center_source_starts else: - kwargs["dipstr"] = src_weights kwargs["dipstr_offsets"] = info.center_source_offsets kwargs["dipstr_starts"] = info.center_source_starts - kwargs["dipvec"] = self.dipole_vec - kwargs["dipvec_offsets"] = info.center_source_offsets - kwargs["dipvec_starts"] = info.center_source_starts + if self.dim == 2 and self.eqn_letter == "l": + kwargs["dipstr"] = -src_weights * ( + self.dipole_vec[0] + 1j*self.dipole_vec[1]) + else: + kwargs["dipstr"] = src_weights + + kwargs["dipvec"] = self.dipole_vec + kwargs["dipvec_offsets"] = info.center_source_offsets + kwargs["dipvec_starts"] = info.center_source_starts # These get max'd/added onto: pass initialized versions. ier = np.zeros(info.ngqbx_centers, dtype=np.int32) @@ -440,11 +451,10 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): tgt_icenter_vec = geo_data.global_qbx_centers() icontaining_tgt_box_vec = qbx_center_to_target_box[tgt_icenter_vec] - # FIXME - rscale2 = np.ones(ngqbx_centers, np.float64) + rscale2 = geo_data.expansion_radii()[geo_data.global_qbx_centers()] kwargs = {} - if self.dim == 3: + if self.dim == 3 and self.eqn_letter == "h": kwargs["radius"] = (0.5 * geo_data.expansion_radii()[geo_data.global_qbx_centers()]) @@ -457,8 +467,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): src_boxes_starts[0] = 0 src_boxes_starts[1:] = np.cumsum(nsrc_boxes_per_gqbx_center) - # FIXME - rscale1 = np.ones(nsrc_boxes) + rscale1 = np.ones(nsrc_boxes) * level_to_rscale(self.tree, isrc_level) rscale1_offsets = np.arange(nsrc_boxes) src_ibox = np.empty(nsrc_boxes, dtype=np.int32) @@ -478,8 +487,12 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): print("end par data prep") - # These get max'd/added onto: pass initialized versions. - ier = np.zeros(ngqbx_centers, dtype=np.int32) + if self.dim == 3: + # This gets max'd onto: pass initialized version. + ier = np.zeros(ngqbx_centers, dtype=np.int32) + kwargs["ier"] = ier + + # This gets added onto: pass initialized version. expn2 = np.zeros( (ngqbx_centers,) + self.expansion_shape(self.qbx_order), dtype=self.dtype) @@ -503,10 +516,13 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): # FIXME: center2 has wrong layout, will copy center2=qbx_centers[:, tgt_icenter_vec], expn2=expn2.T, - ier=ier, **kwargs).T + if self.dim == 3: + if ier.any(): + raise RuntimeError("m2qbxl failed") + local_exps[geo_data.global_qbx_centers()] += expn2 return local_exps @@ -522,8 +538,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): trav = geo_data.traversal() qbx_center_to_target_box = geo_data.qbx_center_to_target_box() qbx_centers = geo_data.centers() - - rscale = 1 # FIXME + qbx_radii = geo_data.expansion_radii() locloc = self.get_translation_routine("%ddlocloc") @@ -540,7 +555,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): kwargs.update(self.kernel_kwargs) for tgt_icenter in range(geo_data.ncenters): - if self.dim == 3: + if self.dim == 3 and self.eqn_letter == "h": # Yuck: This keeps overwriting 'radius' in the dict. kwargs["radius"] = 0.5 * ( geo_data.expansion_radii()[tgt_icenter]) @@ -561,11 +576,11 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): if in_range: src_center = self.tree.box_centers[:, src_ibox] tmp_loc_exp = locloc( - rscale1=rscale, + rscale1=level_to_rscale(self.tree, isrc_level), center1=src_center, expn1=local_exps[src_ibox].T, - rscale2=rscale, + rscale2=qbx_radii[tgt_icenter], center2=tgt_center, nterms2=local_order, @@ -582,11 +597,10 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): ctt = geo_data.center_to_tree_targets() global_qbx_centers = geo_data.global_qbx_centers() qbx_centers = geo_data.centers() + qbx_radii = geo_data.expansion_radii() all_targets = geo_data.all_targets() - rscale = 1 # FIXME - taeval = self.get_expn_eval_routine("ta") for isrc_center, src_icenter in enumerate(global_qbx_centers): @@ -599,7 +613,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): center = qbx_centers[:, src_icenter] pot, grad = taeval( - rscale=rscale, + rscale=qbx_radii[src_icenter], center=center, expn=qbx_expansions[src_icenter].T, ztarg=all_targets[:, center_itgt], diff --git a/pytential/qbx/interactions.py b/pytential/qbx/interactions.py index 172953f3e91297fc83a36117ccb51c938cb111f3..0d78dcdb2d4562fbc5f1fc515434bdb9a6eea4ca 100644 --- a/pytential/qbx/interactions.py +++ b/pytential/qbx/interactions.py @@ -32,11 +32,18 @@ from sumpy.e2e import E2EBase from sumpy.e2p import E2PBase +PYTENTIAL_KERNEL_VERSION = 5 + + # {{{ form qbx expansions from points class P2QBXLFromCSR(P2EBase): default_name = "p2qbxl_from_csr" + def get_cache_key(self): + return super(P2QBXLFromCSR, self).get_cache_key() + ( + PYTENTIAL_KERNEL_VERSION,) + def get_kernel(self): ncoeffs = len(self.expansion) @@ -54,6 +61,7 @@ class P2QBXLFromCSR(P2EBase): None, shape=None), lp.GlobalArg("qbx_centers", None, shape="dim, ncenters", dim_tags="sep,c"), + lp.GlobalArg("qbx_expansion_radii", None, shape="ncenters"), lp.GlobalArg("qbx_expansions", None, shape=("ncenters", ncoeffs)), lp.ValueArg("ncenters", np.int32), @@ -65,12 +73,16 @@ class P2QBXLFromCSR(P2EBase): [ "{[itgt_center]: 0<=itgt_center tgt_icenter = global_qbx_centers[itgt_center] + <> center[idim] = qbx_centers[idim, tgt_icenter] {dup=idim} + <> rscale = qbx_expansion_radii[tgt_icenter] + <> itgt_box = qbx_center_to_target_box[tgt_icenter] <> isrc_box_start = source_box_starts[itgt_box] @@ -81,8 +93,6 @@ class P2QBXLFromCSR(P2EBase): <> isrc_start = box_source_starts[src_ibox] <> isrc_end = isrc_start+box_source_counts_nonchild[src_ibox] - <> center[idim] = qbx_centers[idim, tgt_icenter] - for isrc <> a[idim] = center[idim] - sources[idim, isrc] \ {dup=idim} @@ -132,6 +142,9 @@ class M2QBXL(E2EBase): default_name = "m2qbxl_from_csr" + def get_cache_key(self): + return super(M2QBXL, self).get_cache_key() + (PYTENTIAL_KERNEL_VERSION,) + def get_kernel(self): ncoeff_src = len(self.src_expansion) ncoeff_tgt = len(self.tgt_expansion) @@ -149,6 +162,7 @@ class M2QBXL(E2EBase): <> tgt_center[idim] = qbx_centers[idim, icenter] \ {id=fetch_tgt_center} + <> tgt_rscale = qbx_expansion_radii[icenter] <> isrc_start = src_box_starts[icontaining_tgt_box] <> isrc_stop = src_box_starts[icontaining_tgt_box+1] @@ -179,10 +193,12 @@ class M2QBXL(E2EBase): """], [ lp.GlobalArg("centers", None, shape="dim, aligned_nboxes"), + lp.ValueArg("src_rscale", None), lp.GlobalArg("src_box_starts, src_box_lists", None, shape=None, strides=(1,)), lp.GlobalArg("qbx_centers", None, shape="dim, ncenters", dim_tags="sep,c"), + lp.GlobalArg("qbx_expansion_radii", None, shape="ncenters"), lp.ValueArg("aligned_nboxes,nsrc_level_boxes", np.int32), lp.ValueArg("src_base_ibox", np.int32), lp.GlobalArg("src_expansions", None, @@ -211,7 +227,15 @@ class M2QBXL(E2EBase): return knl def __call__(self, queue, **kwargs): - return self.get_cached_optimized_kernel()(queue, **kwargs) + centers = kwargs.pop("centers") + # "1" may be passed for rscale, which won't have its type + # meaningfully inferred. Make the type of rscale explicit. + src_rscale = centers.dtype.type(kwargs.pop("src_rscale")) + + return self.get_cached_optimized_kernel()(queue, + centers=centers, + src_rscale=src_rscale, + **kwargs) # }}} @@ -221,6 +245,9 @@ class M2QBXL(E2EBase): class L2QBXL(E2EBase): default_name = "l2qbxl" + def get_cache_key(self): + return super(L2QBXL, self).get_cache_key() + (PYTENTIAL_KERNEL_VERSION,) + def get_kernel(self): ncoeff_src = len(self.src_expansion) ncoeff_tgt = len(self.tgt_expansion) @@ -249,6 +276,9 @@ class L2QBXL(E2EBase): if in_range <> tgt_center[idim] = qbx_centers[idim, icenter] <> src_center[idim] = centers[idim, src_ibox] {dup=idim} + + <> tgt_rscale = qbx_expansion_radii[icenter] + <> d[idim] = tgt_center[idim] - src_center[idim] {dup=idim} """] + [""" @@ -268,8 +298,10 @@ class L2QBXL(E2EBase): lp.GlobalArg("target_boxes", None, shape=None, offset=lp.auto), lp.GlobalArg("centers", None, shape="dim, naligned_boxes"), + lp.ValueArg("src_rscale", None), lp.GlobalArg("qbx_centers", None, shape="dim, ncenters", dim_tags="sep,c"), + lp.GlobalArg("qbx_expansion_radii", None, shape="ncenters"), lp.ValueArg("naligned_boxes,target_base_ibox,nboxes", np.int32), lp.GlobalArg("expansions", None, shape=("nboxes", ncoeff_src), offset=lp.auto), @@ -298,7 +330,15 @@ class L2QBXL(E2EBase): return knl def __call__(self, queue, **kwargs): - return self.get_cached_optimized_kernel()(queue, **kwargs) + centers = kwargs.pop("centers") + # "1" may be passed for rscale, which won't have its type + # meaningfully inferred. Make the type of rscale explicit. + src_rscale = centers.dtype.type(kwargs.pop("src_rscale")) + + return self.get_cached_optimized_kernel()(queue, + centers=centers, + src_rscale=src_rscale, + **kwargs) # }}} @@ -308,6 +348,9 @@ class L2QBXL(E2EBase): class QBXL2P(E2PBase): default_name = "qbx_potential_from_local" + def get_cache_key(self): + return super(QBXL2P, self).get_cache_key() + (PYTENTIAL_KERNEL_VERSION,) + def get_kernel(self): ncoeffs = len(self.expansion) @@ -328,11 +371,13 @@ class QBXL2P(E2PBase): <> icenter_tgt_start = center_to_targets_starts[src_icenter] <> icenter_tgt_end = center_to_targets_starts[src_icenter+1] + <> center[idim] = qbx_centers[idim, src_icenter] {dup=idim} + <> rscale = qbx_expansion_radii[src_icenter] + for icenter_tgt <> center_itgt = center_to_targets_lists[icenter_tgt] - <> center[idim] = qbx_centers[idim, src_icenter] {dup=idim} <> b[idim] = targets[idim, center_itgt] - center[idim] """] + [""" diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py index 9f2f101ae64c8a654ec28b41b98fc2b92f395d3b..84cd46ce23990fb893f47f3e3c2d7672aa40129f 100644 --- a/pytential/symbolic/matrix.py +++ b/pytential/symbolic/matrix.py @@ -190,6 +190,7 @@ class MatrixBuilder(EvaluationMapperBase): target_discr.nodes(), source.quad_stage2_density_discr.nodes(), get_centers_on_side(source, expr.qbx_forced_limit), + expansion_radii=self.dep_source._expansion_radii("nsources"), **kernel_args) mat = mat.get() diff --git a/requirements.txt b/requirements.txt index 7360a7c78e2a50cc32a45024544cdfb468da0934..71e9d6a4a72be6036f6108d94cad50417ee58438 100644 --- a/requirements.txt +++ b/requirements.txt @@ -7,5 +7,5 @@ git+https://github.com/inducer/islpy git+https://github.com/inducer/loopy git+https://gitlab.tiker.net/inducer/boxtree git+https://github.com/inducer/meshmode -git+https://github.com/inducer/sumpy +git+https://gitlab.tiker.net/inducer/sumpy@rscale git+https://github.com/inducer/pyfmmlib diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index 1e1c92135e903ff36db851e2ace0a1862116be98..e8044d39b751ff1ad4a0705d86f910fcd4871950 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -108,7 +108,7 @@ class SphereGeometry(object): mesh_name = "sphere" dim = 3 - resolutions = [0, 1, 2] + resolutions = [0, 1] def get_mesh(self, resolution, tgt_order): return get_sphere_mesh(resolution, tgt_order) @@ -183,7 +183,7 @@ class StarfishGreenTest(StaticTestCase): _expansion_stick_out_factor = 0.5 - fmm_backend = "sumpy" + fmm_backend = "fmmlib" class WobblyCircleGreenTest(StaticTestCase): @@ -201,7 +201,7 @@ class WobblyCircleGreenTest(StaticTestCase): class SphereGreenTest(StaticTestCase): expr = GreenExpr() geometry = SphereGeometry() - k = 1.2 + k = 0 qbx_order = 3 fmm_order = 10 @@ -215,11 +215,12 @@ class SphereGreenTest(StaticTestCase): class DynamicTestCase(object): fmm_backend = "sumpy" - def __init__(self, geometry, expr, k): + def __init__(self, geometry, expr, k, fmm_backend="sumpy"): self.geometry = geometry self.expr = expr self.k = k self.qbx_order = 5 if geometry.dim == 2 else 3 + self.fmm_backend = fmm_backend if geometry.dim == 2: order_bump = 15 @@ -229,7 +230,9 @@ class DynamicTestCase(object): self.fmm_order = self.qbx_order + order_bump def check(self): - if self.geometry.mesh_name == "sphere" and self.k != 0: + if (self.geometry.mesh_name == "sphere" + and self.k != 0 + and self.fmm_backend == "sumpy"): pytest.skip("both direct eval and generating the FMM kernels " "are too slow") @@ -237,6 +240,9 @@ class DynamicTestCase(object): and self.expr.zero_op_name == "green_grad"): pytest.skip("does not achieve sufficient precision") + if self.fmm_backend == "fmmlib": + pytest.importorskip("pyfmmlib") + # {{{ integral identity tester @@ -252,6 +258,9 @@ class DynamicTestCase(object): DynamicTestCase(geom, GradGreenExpr(), 0), DynamicTestCase(geom, GradGreenExpr(), 1.2), DynamicTestCase(geom, ZeroCalderonExpr(), 0), + + DynamicTestCase(geom, GreenExpr(), 0, fmm_backend="fmmlib"), + DynamicTestCase(geom, GreenExpr(), 1.2, fmm_backend="fmmlib"), ]]) def test_identity_convergence(ctx_getter, case, visualize=False): logging.basicConfig(level=logging.INFO) diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py index 9443261aac95d329db2721d50b8fcb9f847b12d7..ec8a4c94c569843228e451948bf7b8231d599e25 100644 --- a/test/test_scalar_int_eq.py +++ b/test/test_scalar_int_eq.py @@ -683,8 +683,6 @@ class ManyEllipsoidIntEqTestCase(Helmholtz3DIntEqTestCase): for helmholtz_k in [0, 1.2] for bc_type in ["dirichlet", "neumann"] for loc_sign in [-1, +1] - ] + [ - EllipsoidIntEqTestCase(0.7, "neumann", +1) ]) # Sample test run: # 'test_integral_equation(cl._csc, EllipseIntEqTestCase(0, "dirichlet", +1), 5)' # noqa: E501