From dbba7115ccf9cc0a291bc04e78723fd79af9f4d9 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 5 Jun 2017 19:29:40 -0400 Subject: [PATCH 01/13] Add an explanatory comment to L2QBXL --- pytential/qbx/interactions.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pytential/qbx/interactions.py b/pytential/qbx/interactions.py index 38271cdf..172953f3 100644 --- a/pytential/qbx/interactions.py +++ b/pytential/qbx/interactions.py @@ -241,6 +241,8 @@ class L2QBXL(E2EBase): <> src_ibox = target_boxes[isrc_box] \ {id=read_src_ibox} + # Is the box number on the level currently under + # consideration? <> in_range = (target_base_ibox <= src_ibox and src_ibox < target_base_ibox + nboxes) -- GitLab From 6764e63456b905faa8cb3fcd89e5008ee096e33b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 5 Jun 2017 19:30:21 -0400 Subject: [PATCH 02/13] Work towards multiple FMM backends, FMMLib backend --- pytential/qbx/__init__.py | 31 +++- pytential/qbx/fmm.py | 16 +- pytential/qbx/fmmlib.py | 382 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 412 insertions(+), 17 deletions(-) create mode 100644 pytential/qbx/fmmlib.py diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index ad657877..5292f7bb 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -78,7 +78,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): debug=True, refined_for_global_qbx=False, expansion_disks_in_tree_have_extent=False, - performance_data_file=None): + performance_data_file=None, + fmm_backend="sumpy"): """ :arg fine_order: The total degree to which the (upsampled) underlying quadrature is exact. @@ -115,6 +116,7 @@ class QBXLayerPotentialSource(LayerPotentialSource): self.density_discr = density_discr self.fmm_level_to_order = fmm_level_to_order self.target_stick_out_factor = target_stick_out_factor + self.fmm_backend = fmm_backend # Default values are lazily provided if these are None self._base_resampler = base_resampler @@ -164,7 +166,8 @@ class QBXLayerPotentialSource(LayerPotentialSource): 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) + performance_data_file=self.performance_data_file, + fmm_backend=self.fmm_backend) # }}} @@ -418,12 +421,24 @@ class QBXLayerPotentialSource(LayerPotentialSource): fmm_local_factory = partial(local_expn_class, base_kernel) qbx_local_factory = partial(local_expn_class, base_kernel) - from pytential.qbx.fmm import \ - QBXExpansionWranglerCodeContainer - return QBXExpansionWranglerCodeContainer( - self.cl_context, - fmm_mpole_factory, fmm_local_factory, qbx_local_factory, - out_kernels) + if self.fmm_backend == "sumpy": + from pytential.qbx.fmm import \ + QBXSumpyExpansionWranglerCodeContainer + return QBXSumpyExpansionWranglerCodeContainer( + self.cl_context, + fmm_mpole_factory, fmm_local_factory, qbx_local_factory, + out_kernels) + + elif self.fmm_backend == "fmmlib": + from pytential.qbx.fmmlib import \ + QBXFMMLibExpansionWranglerCodeContainer + return QBXFMMLibExpansionWranglerCodeContainer( + self.cl_context, + fmm_mpole_factory, fmm_local_factory, qbx_local_factory, + out_kernels) + + else: + raise ValueError("invalid FMM backend: %s" % self.fmm_backend) def exec_compute_potential_insn_fmm(self, queue, insn, bound_expr, evaluate): # {{{ build list of unique target discretizations used diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index 426de45d..ad07bb5c 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -37,7 +37,7 @@ logger = logging.getLogger(__name__) __doc__ = """ -.. autoclass:: QBXExpansionWranglerCodeContainer +.. autoclass:: QBXSumpyExpansionWranglerCodeContainer .. autoclass:: QBXExpansionWrangler @@ -45,9 +45,9 @@ __doc__ = """ """ -# {{{ expansion wrangler +# {{{ sumpy expansion wrangler -class QBXExpansionWranglerCodeContainer(SumpyExpansionWranglerCodeContainer): +class QBXSumpyExpansionWranglerCodeContainer(SumpyExpansionWranglerCodeContainer): def __init__(self, cl_context, multipole_expansion_factory, local_expansion_factory, qbx_local_expansion_factory, out_kernels): @@ -175,12 +175,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), # {{{ source/target dispatch - def box_source_list_kwargs(self): - return dict( - box_source_starts=self.tree.box_source_starts, - box_source_counts_nonchild=( - self.tree.box_source_counts_nonchild), - sources=self.tree.sources) + # box_source_list_kwargs inherited from superclass def box_target_list_kwargs(self): # This only covers the non-QBX targets. @@ -529,6 +524,8 @@ def drive_fmm(expansion_wrangler, src_weights): # }}} +# {{{ performance model + def write_performance_model(outf, geo_data): from pymbolic import var p_fmm = var("p_fmm") @@ -714,5 +711,6 @@ def write_performance_model(outf, geo_data): # }}} +# }}} # vim: foldmethod=marker diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py new file mode 100644 index 00000000..f7497fd9 --- /dev/null +++ b/pytential/qbx/fmmlib.py @@ -0,0 +1,382 @@ +from __future__ import division, absolute_import + +__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. +""" + +import numpy as np +from pytools import memoize_method +import pyopencl as cl # noqa +import pyopencl.array # noqa: F401 +from boxtree.pyfmmlib_integration import HelmholtzExpansionWrangler + + +class QBXFMMLibExpansionWranglerCodeContainer(object): + def __init__(self, cl_context, + multipole_expansion_factory, local_expansion_factory, + qbx_local_expansion_factory, out_kernels): + self.cl_context = cl_context + self.multipole_expansion_factory = multipole_expansion_factory + self.local_expansion_factory = local_expansion_factory + self.qbx_local_expansion_factory = qbx_local_expansion_factory + + self.out_kernels = out_kernels + + def get_wrangler(self, queue, geo_data, dtype, + qbx_order, fmm_level_to_order, + source_extra_kwargs={}, + kernel_extra_kwargs=None): + + from sumpy.kernel import HelmholtzKernel + for out_knl in self.out_kernels: + if not isinstance(out_knl, HelmholtzKernel): + raise NotImplementedError( + "only the Helmholtz kernel is supported for now") + + return QBXFMMLibHelmholtzExpansionWrangler(self, queue, geo_data, dtype, + qbx_order, fmm_level_to_order, + source_extra_kwargs, + kernel_extra_kwargs, + self.out_kernels) + +# }}} + + +# {{{ host geo data wrapper + +class ToHostTransferredGeoDataWrapper(object): + def __init__(self, queue, geo_data): + self.queue = queue + self.geo_data = geo_data + + @memoize_method + def tree(self): + return self.traversal().tree + + @memoize_method + def traversal(self): + return self.geo_data.traversal().get(queue=self.queue) + + @property + def ncenters(self): + return self.geo_data.ncenters + + @memoize_method + def centers(self): + return np.array([ + ci.get(queue=self.queue) + for ci in self.geo_data.centers()]) + + @memoize_method + def global_qbx_centers(self): + return self.geo_data.global_qbx_centers().get(queue=self.queue) + + @memoize_method + def qbx_center_to_target_box(self): + return self.geo_data.qbx_center_to_target_box().get(queue=self.queue) + + @memoize_method + def non_qbx_box_target_lists(self): + return self.geo_data.non_qbx_box_target_lists().get(queue=self.queue) + + @memoize_method + def center_to_tree_targets(self): + return self.geo_data.center_to_tree_targets().get(queue=self.queue) + +# }}} + + +# {{{ fmmlib expansion wrangler + +class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): + def __init__(self, code_container, queue, geo_data, dtype, + qbx_order, fmm_level_to_order, + source_extra_kwargs, + kernel_extra_kwargs, + out_kernels): + + from pytools import single_valued + k_name = single_valued(out_knl.helmholtz_k_name for out_knl in out_kernels) + helmholtz_k = kernel_extra_kwargs[k_name] + + self.code_container = code_container + self.queue = queue + self.geo_data = ToHostTransferredGeoDataWrapper(queue, geo_data) + self.qbx_order = qbx_order + + self.level_orders = [ + fmm_level_to_order(level) + for level in range(self.geo_data.tree().nlevels)] + + # FIXME: For now + from pytools import single_valued + assert single_valued(self.level_orders) + + super(QBXFMMLibHelmholtzExpansionWrangler, self).__init__( + # FMMLib is CPU-only--get the tree out of OpenCL-land + self.geo_data.tree(), + + helmholtz_k=helmholtz_k, + + # FIXME + nterms=fmm_level_to_order(0)) + + def potential_zeros(self): + """This ought to be called ``non_qbx_potential_zeros``, but since + it has to override the superclass's behavior to integrate seamlessly, + it needs to be called just :meth:`potential_zeros`. + """ + + nqbtl = self.geo_data.non_qbx_box_target_lists() + + # from pytools.obj_array import make_obj_array + # return make_obj_array([ + # cl.array.zeros( + # self.queue, + # nqbtl.nfiltered_targets, + # dtype=self.dtype) + # for k in self.code.out_kernels]) + + return np.zeros(nqbtl.nfiltered_targets, self.dtype) + + def full_potential_zeros(self): + # The superclass generates a full field of zeros, for all + # (not just non-QBX) targets. + return super(QBXFMMLibHelmholtzExpansionWrangler, self).potential_zeros() + + def reorder_sources(self, source_array): + source_array = source_array.get(queue=self.queue) + return (super(QBXFMMLibHelmholtzExpansionWrangler, self) + .reorder_sources(source_array)) + + def reorder_potentials(self, potentials): + raise NotImplementedError("reorder_potentials should not " + "be called on a QBXFMMLibHelmholtzExpansionWrangler") + + # Because this is a multi-stage, more complicated process that combines + # potentials from non-QBX targets and QBX targets. + + # {{{ override target lists to only hit non-QBX targets + + def box_target_starts(self): + nqbtl = self.geo_data.non_qbx_box_target_lists() + return nqbtl.box_target_starts + + def box_target_counts_nonchild(self): + nqbtl = self.geo_data.non_qbx_box_target_lists() + return nqbtl.box_target_counts_nonchild + + def targets(self): + nqbtl = self.geo_data.non_qbx_box_target_lists() + return nqbtl.targets + + # }}} + + # {{{ qbx-related + + def qbx_local_expansion_zeros(self): + return np.zeros( + (self.geo_data.ncenters,) + self.expansion_shape(self.qbx_order), + dtype=self.dtype) + + def form_global_qbx_locals(self, starts, lists, src_weights): + local_exps = self.qbx_local_expansion_zeros() + + rscale = 1 # FIXME + geo_data = self.geo_data + + if len(geo_data.global_qbx_centers()) == 0: + return local_exps + + qbx_center_to_target_box = geo_data.qbx_center_to_target_box() + qbx_centers = geo_data.centers() + + formta = self.get_routine("%ddformta") + + for itgt_center, tgt_icenter in enumerate(geo_data.global_qbx_centers()): + itgt_box = qbx_center_to_target_box[tgt_icenter] + + isrc_box_start = starts[itgt_box] + isrc_box_stop = starts[itgt_box+1] + + tgt_center = qbx_centers[:, tgt_icenter] + + ctr_coeffs = 0 + + for isrc_box in range(isrc_box_start, isrc_box_stop): + src_ibox = lists[isrc_box] + + src_pslice = self._get_source_slice(src_ibox) + + ier, coeffs = formta( + self.helmholtz_k, rscale, + self._get_sources(src_pslice), src_weights[src_pslice], + tgt_center, self.nterms) + if ier: + raise RuntimeError("formta failed") + + ctr_coeffs += coeffs + + local_exps[tgt_icenter] += ctr_coeffs + + return local_exps + + def translate_box_multipoles_to_qbx_local(self, multipole_exps): + local_exps = self.qbx_local_expansion_zeros() + + geo_data = self.geo_data + qbx_center_to_target_box = geo_data.qbx_center_to_target_box() + qbx_centers = geo_data.centers() + centers = self.tree.box_centers + + mploc = self.get_translation_routine("%ddmploc") + + for isrc_level, ssn in enumerate(geo_data.traversal().sep_smaller_by_level): + source_level_start_ibox, source_mpoles_view = \ + self.multipole_expansions_view(multipole_exps, isrc_level) + + # FIXME + rscale = 1 + + kwargs = {} + if self.dim == 3: + # FIXME Is this right? + kwargs["radius"] = self.tree.root_extent * 2**(-isrc_level) + + for itgt_center, tgt_icenter in enumerate(geo_data.global_qbx_centers()): + ctr_loc = 0 + + icontaining_tgt_box = qbx_center_to_target_box[tgt_icenter] + + tgt_center = qbx_centers[:, tgt_icenter] + + for isrc_box in range( + ssn.starts[icontaining_tgt_box], + ssn.starts[icontaining_tgt_box+1]): + + src_ibox = ssn.lists[isrc_box] + src_center = centers[:, src_ibox] + + ctr_loc = ctr_loc + mploc( + self.helmholtz_k, + rscale, src_center, multipole_exps[src_ibox], + rscale, tgt_center, self.nterms, **kwargs)[..., 0] + + local_exps[tgt_icenter] += ctr_loc + + return local_exps + + def translate_box_local_to_qbx_local(self, local_exps): + qbx_expansions = self.qbx_local_expansion_zeros() + + geo_data = self.geo_data + if geo_data.ncenters == 0: + return qbx_expansions + trav = geo_data.traversal() + qbx_center_to_target_box = geo_data.qbx_center_to_target_box() + qbx_centers = geo_data.centers() + + rscale = 1 # FIXME + + locloc = self.get_translation_routine("%ddlocloc") + + for isrc_level in range(geo_data.tree().nlevels): + local_order = self.level_orders[isrc_level] + + lev_box_start, lev_box_stop = self.tree.level_start_box_nrs[ + isrc_level:isrc_level+2] + target_level_start_ibox, target_locals_view = \ + self.local_expansions_view(local_exps, isrc_level) + assert target_level_start_ibox == lev_box_start + + kwargs = {} + if self.dim == 3: + # FIXME Is this right? + kwargs["radius"] = self.tree.root_extent * 2**(-isrc_level) + + for tgt_icenter in range(geo_data.ncenters): + isrc_box = qbx_center_to_target_box[tgt_icenter] + + tgt_center = qbx_centers[:, tgt_icenter] + + # The box's expansions which we're translating here + # (our source) is, globally speaking, a target box. + + src_ibox = trav.target_boxes[isrc_box] + + # Is the box number on the level currently under + # consideration? + in_range = (lev_box_start <= src_ibox and src_ibox < lev_box_stop) + + if in_range: + src_center = self.tree.box_centers[:, src_ibox] + tmp_loc_exp = locloc( + self.helmholtz_k, + rscale, src_center, local_exps[src_ibox], + rscale, tgt_center, local_order, **kwargs)[..., 0] + + qbx_expansions[tgt_icenter] += tmp_loc_exp + + return qbx_expansions + + def eval_qbx_expansions(self, qbx_expansions): + pot = self.full_potential_zeros() + + geo_data = self.geo_data + ctt = geo_data.center_to_tree_targets() + + rscale = 1 # FIXME + + taeval = self.get_expn_eval_routine("ta") + + for iglobal_center + src_icenter = global_qbx_centers[iglobal_center] + + icenter_tgt_start = center_to_targets_starts[src_icenter] + icenter_tgt_end = center_to_targets_starts[src_icenter+1] + + for icenter_tgt + + center_itgt = center_to_targets_lists[icenter_tgt] + + center = qbx_centers[:, src_icenter] + b[idim] = targets[idim, center_itgt] - center[idim] + + """] + [""" + <> coeff{i} = qbx_expansions[src_icenter, {i}] + """.format(i=i) for i in range(ncoeffs)] + [ + + ] + loopy_insns + [""" + + result[{i},center_itgt] = kernel_scaling * result_{i}_p \ + {{id_prefix=write_result}} + """.format(i=i) for i in range(len(result_names))] + [""" + tmp_pot = taeval(self.helmholtz_k, rscale, + center, qbx_expansions[src_icenter], + self._get_targets(tgt_pslice)) + + end + end + # }}} + +# }}} + +# vim: foldmethod=marker -- GitLab From a3767abc84037d4bce6c6024083b6cae77a83ad9 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 5 Jun 2017 19:32:52 -0400 Subject: [PATCH 03/13] Work towards reviving the C-C Maxwell solver --- examples/maxwell.py | 112 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 112 insertions(+) create mode 100644 examples/maxwell.py diff --git a/examples/maxwell.py b/examples/maxwell.py new file mode 100644 index 00000000..7bde17a6 --- /dev/null +++ b/examples/maxwell.py @@ -0,0 +1,112 @@ +from __future__ import division +import numpy as np +import pyopencl as cl +from sumpy.visualization import FieldPlotter +#from mayavi import mlab +from sumpy.kernel import one_kernel_2d, LaplaceKernel, HelmholtzKernel # noqa + +import faulthandler +from six.moves import range +faulthandler.enable() + +import logging +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger(__name__) + +cl_ctx = cl.create_some_context() +queue = cl.CommandQueue(cl_ctx) + +target_order = 5 +qbx_order = 2 +nelements = 60 +mode_nr = 0 + +h = .5 + +k = 2 + + +def main(): + from meshmode.mesh.io import generate_gmsh, FileSource + mesh = generate_gmsh( + FileSource("ellipsoid.step"), 2, order=2, + other_options=["-string", "Mesh.CharacteristicLengthMax = %g;" % h]) + + from meshmode.mesh.processing import find_bounding_box + bbox_min, bbox_max = find_bounding_box(mesh) + bbox_center = 0.5*(bbox_min+bbox_max) + bbox_size = max(bbox_max-bbox_min) / 2 + + logger.info("%d elements" % mesh.nelements) + + from pytential.qbx import QBXLayerPotentialSource + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + density_discr = Discretization( + cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + qbx = QBXLayerPotentialSource(density_discr, 4*target_order, qbx_order, + fmm_order=qbx_order, fmm_backend="fmmlib") + + from pytential.symbolic.pde.maxwell import PECAugmentedMFIEOperator + pde_op = PECAugmentedMFIEOperator() + from pytential import bind, sym + + jt_sym = sym.make_sym_vector("jt", 2) + rho_sym = sym.var("rho") + repr_op = pde_op.scattered_volume_field(jt_sym, rho_sym) + + # {{{ make a density + + nodes = density_discr.nodes().with_queue(queue) + + angle = cl.clmath.atan2(nodes[1], nodes[0]) + + sigma = cl.clmath.cos(mode_nr*angle) + if 0: + sigma = 0*angle + from random import randrange + for i in range(5): + sigma[randrange(len(sigma))] = 1 + + sigma = sigma.astype(np.complex128) + + jt = sym.make_obj_array([sigma, sigma]) + rho = sigma + + # }}} + + #mlab.figure(bgcolor=(1, 1, 1)) + if 1: + fplot = FieldPlotter(bbox_center, extent=1.5*bbox_size, npoints=30) + + qbx_stick_out = qbx.copy(target_stick_out_factor=0.1) + from pytential.target import PointsTarget + fld_in_vol = bind( + (qbx_stick_out, PointsTarget(fplot.points)), + repr_op)(queue, jt=jt, rho=rho, k=k).get() + + #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) + fplot.write_vtk_file( + "potential.vts", + [ + ("potential", fld_in_vol) + ] + ) + + bdry_normals = bind(density_discr, sym.normal())(queue)\ + .as_vector(dtype=object) + + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, density_discr, target_order) + + bdry_vis.write_vtk_file("source.vtu", [ + ("sigma", sigma), + ("bdry_normals", bdry_normals), + ]) + + +if __name__ == "__main__": + main() -- GitLab From dd808c2c9846fd38b1569270b34518417818fc37 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 5 Jun 2017 19:38:08 -0400 Subject: [PATCH 04/13] Work towards reviving the C-C Maxwell solver (non-driver bits) --- pytential/symbolic/pde/maxwell/__init__.py | 21 +++++----- pytential/symbolic/primitives.py | 49 ++++++++++++++++++++++ 2 files changed, 59 insertions(+), 11 deletions(-) diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py index 68cceb5e..9096998a 100644 --- a/pytential/symbolic/pde/maxwell/__init__.py +++ b/pytential/symbolic/pde/maxwell/__init__.py @@ -37,7 +37,9 @@ class PECAugmentedMFIEOperator: see notes/mfie.tm """ - def __init__(self, k): + def __init__(self, k=sym.var("k")): + from sumpy.kernel import HelmholtzKernel + self.kernel = HelmholtzKernel(3) self.k = k def j_operator(self, loc, Jt): @@ -58,9 +60,8 @@ class PECAugmentedMFIEOperator: def scattered_boundary_field(self, Jt, rho, loc): Jxyz = cse(tangential_to_xyz(Jt), "Jxyz") - k = self.k - A = S(k, Jxyz) + A = S(self.kernel, Jxyz, k=self.k) grad_phi = grad_S(k, rho, 3) # use - n x n x v = v_tangential @@ -72,18 +73,16 @@ class PECAugmentedMFIEOperator: return join_fields(E_scat, H_scat) def scattered_volume_field(self, Jt, rho): - Jxyz = cse(tangential_to_xyz(Jt), "Jxyz") - - k = self.k + Jxyz = sym.cse(sym.tangential_to_xyz(Jt), "Jxyz") - A = S(k, Jxyz) - grad_phi = grad_S(k, rho, 3) + A = sym.S(self.kernel, Jxyz, k=self.k, qbx_forced_limit=None) + #grad_phi = grad_S(self.kernel, rho, 3, k=self.k) - E_scat = 1j*k*A - grad_phi - H_scat = curl_S_volume(k, Jxyz) + E_scat = 1j*self.k*A# - grad_phi + #H_scat = curl_S_volume(k, Jxyz) from pytools.obj_array import join_fields - return join_fields(E_scat, H_scat) + return E_scat #join_fields(E_scat, H_scat) # }}} diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py index b4d1c615..c4c63b99 100644 --- a/pytential/symbolic/primitives.py +++ b/pytential/symbolic/primitives.py @@ -894,6 +894,55 @@ def Dp(kernel, *args, **kwargs): # noqa # }}} +# {{{ geometric operations + +def tangential_onb(ambient_dim, dim=None, where=None): + if dim is None: + dim = ambient_dim - 1 + + pd_mat = cse(parametrization_derivative_matrix(ambient_dim, dim, where), + "pd_matrix", cse_scope.DISCRETIZATION) + + # {{{ Gram-Schmidt + + orth_pd_mat = np.zeros_like(pd_mat) + for k in range(pd_mat.shape[1]): + avec = pd_mat[:, k] + q = avec + for j in range(k): + q = q - np.dot(avec, orth_pd_mat[:, j])*orth_pd_mat[:, j] + + orth_pd_mat[:, k] = cse(q/sqrt(np.sum(q**2)), "orth_pd_vec") + + # }}} + + return orth_pd_mat + + +def xyz_to_tangential(xyz_vec, where=None): + ambient_dim = len(xyz_vec) + tonb = tangential_onb(ambient_dim, where=where) + return make_obj_array([ + np.dot(tonb[:, i], xyz_vec) + for i in range(ambient_dim - 1) + ]) + + +def tangential_to_xyz(tangential_vec, where=None): + ambient_dim = len(tangential_vec) + 1 + tonb = tangential_onb(ambient_dim, where=where) + return sum( + tonb[:, i] * tangential_vec[i] + for i in range(ambient_dim - 1)) + + +def project_to_tangential(xyz_vec, which=None): + return tangential_to_xyz( + cse(xyz_to_tangential(xyz_vec, which), which)) + +# }}} + + def pretty(expr): # Doesn't quite belong here, but this is exposed to the user as # "pytential.sym", so in here it goes. -- GitLab From f0f1e4eb61c27c8d1c52499f5984c772409cb575 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 5 Jun 2017 21:52:41 -0400 Subject: [PATCH 05/13] Fix 'sym.normal()' invocation in layerpot-3d example --- examples/layerpot-3d.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/layerpot-3d.py b/examples/layerpot-3d.py index da1f50e7..68d9e33e 100644 --- a/examples/layerpot-3d.py +++ b/examples/layerpot-3d.py @@ -90,7 +90,9 @@ if 1: ] ) - bdry_normals = bind(density_discr, sym.normal())(queue).as_vector(dtype=object) + bdry_normals = bind( + density_discr, + sym.normal(density_discr.ambient_dim))(queue).as_vector(dtype=object) from meshmode.discretization.visualization import make_visualizer bdry_vis = make_visualizer(queue, density_discr, target_order) -- GitLab From 65af801053d1dfe5bcdc487f96c3ff20407030b6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 5 Jun 2017 21:57:08 -0400 Subject: [PATCH 06/13] Upgrade FMM progress messages to 'info' level --- pytential/qbx/fmm.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index ad07bb5c..4ad50fc3 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -366,12 +366,12 @@ def drive_fmm(expansion_wrangler, src_weights): logger.info("start qbx fmm") - logger.debug("reorder source weights") + logger.info("reorder source weights") src_weights = wrangler.reorder_sources(src_weights) # {{{ construct local multipoles - logger.debug("construct local multipoles") + logger.info("construct local multipoles") mpole_exps = wrangler.form_multipoles( traversal.level_start_source_box_nrs, traversal.source_boxes, @@ -381,7 +381,7 @@ def drive_fmm(expansion_wrangler, src_weights): # {{{ propagate multipoles upward - logger.debug("propagate multipoles upward") + logger.info("propagate multipoles upward") wrangler.coarsen_multipoles( traversal.level_start_source_parent_box_nrs, traversal.source_parent_boxes, @@ -391,7 +391,7 @@ def drive_fmm(expansion_wrangler, src_weights): # {{{ direct evaluation from neighbor source boxes ("list 1") - logger.debug("direct evaluation from neighbor source boxes ('list 1')") + logger.info("direct evaluation from neighbor source boxes ('list 1')") non_qbx_potentials = wrangler.eval_direct( traversal.target_boxes, traversal.neighbor_source_boxes_starts, @@ -402,7 +402,7 @@ def drive_fmm(expansion_wrangler, src_weights): # {{{ translate separated siblings' ("list 2") mpoles to local - logger.debug("translate separated siblings' ('list 2') mpoles to local") + logger.info("translate separated siblings' ('list 2') mpoles to local") local_exps = wrangler.multipole_to_local( traversal.level_start_target_or_target_parent_box_nrs, traversal.target_or_target_parent_boxes, @@ -414,7 +414,7 @@ def drive_fmm(expansion_wrangler, src_weights): # {{{ evaluate sep. smaller mpoles ("list 3") at particles - logger.debug("evaluate sep. smaller mpoles at particles ('list 3 far')") + logger.info("evaluate sep. smaller mpoles at particles ('list 3 far')") # (the point of aiming this stage at particles is specifically to keep its # contribution *out* of the downward-propagating local expansions) @@ -432,7 +432,7 @@ def drive_fmm(expansion_wrangler, src_weights): # {{{ form locals for separated bigger mpoles ("list 4") - logger.debug("form locals for separated bigger mpoles ('list 4 far')") + logger.info("form locals for separated bigger mpoles ('list 4 far')") local_exps = local_exps + wrangler.form_locals( traversal.level_start_target_or_target_parent_box_nrs, @@ -448,7 +448,7 @@ def drive_fmm(expansion_wrangler, src_weights): # {{{ propagate local_exps downward - logger.debug("propagate local_exps downward") + logger.info("propagate local_exps downward") wrangler.refine_locals( traversal.level_start_target_or_target_parent_box_nrs, traversal.target_or_target_parent_boxes, @@ -458,7 +458,7 @@ def drive_fmm(expansion_wrangler, src_weights): # {{{ evaluate locals - logger.debug("evaluate locals") + logger.info("evaluate locals") non_qbx_potentials = non_qbx_potentials + wrangler.eval_locals( traversal.level_start_target_box_nrs, traversal.target_boxes, @@ -468,22 +468,22 @@ def drive_fmm(expansion_wrangler, src_weights): # {{{ wrangle qbx expansions - logger.debug("form global qbx expansions from list 1") + logger.info("form global qbx expansions from list 1") qbx_expansions = wrangler.form_global_qbx_locals( traversal.neighbor_source_boxes_starts, traversal.neighbor_source_boxes_lists, src_weights) - logger.debug("translate from list 3 multipoles to qbx local expansions") + logger.info("translate from list 3 multipoles to qbx local expansions") qbx_expansions = qbx_expansions + \ wrangler.translate_box_multipoles_to_qbx_local(mpole_exps) - logger.debug("translate from box local expansions to contained " + logger.info("translate from box local expansions to contained " "qbx local expansions") qbx_expansions = qbx_expansions + \ wrangler.translate_box_local_to_qbx_local(local_exps) - logger.debug("evaluate qbx local expansions") + logger.info("evaluate qbx local expansions") qbx_potentials = wrangler.eval_qbx_expansions( qbx_expansions) @@ -491,7 +491,7 @@ def drive_fmm(expansion_wrangler, src_weights): # {{{ reorder potentials - logger.debug("reorder potentials") + logger.info("reorder potentials") nqbtl = geo_data.non_qbx_box_target_lists() -- GitLab From 8ced5cba0ca0c6396fe098a1fa0fc40f4743d9c2 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 5 Jun 2017 21:58:34 -0400 Subject: [PATCH 07/13] Add finalize_potentials to QBX wrangler interface --- pytential/qbx/fmm.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index 4ad50fc3..858c3e2d 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -337,6 +337,9 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), # }}} + def finalize_potential(self, potential): + return potential + # }}} @@ -509,7 +512,9 @@ def drive_fmm(expansion_wrangler, src_weights): all_potentials_in_tree_order += qbx_potentials def reorder_potentials(x): - return x[tree.sorted_target_ids] + # "finalize" gives host FMMs (like FMMlib) a chance to turn the + # potential back into a CL array. + return wrangler.finalize_potential(x[tree.sorted_target_ids]) from pytools.obj_array import with_object_array_or_scalar result = with_object_array_or_scalar( -- GitLab From 87615ef78e1750a1a7273214336a2fb2d38257d1 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 5 Jun 2017 21:59:12 -0400 Subject: [PATCH 08/13] Tweak Maxwell example to evaluate only a single SLP for testing --- examples/maxwell.py | 54 +++++++++++++++++++++++++++++---------------- 1 file changed, 35 insertions(+), 19 deletions(-) diff --git a/examples/maxwell.py b/examples/maxwell.py index 7bde17a6..033ff902 100644 --- a/examples/maxwell.py +++ b/examples/maxwell.py @@ -17,13 +17,13 @@ cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) target_order = 5 -qbx_order = 2 +qbx_order = 3 nelements = 60 mode_nr = 0 -h = .5 +h = 1 -k = 2 +k = 4 def main(): @@ -31,6 +31,7 @@ def main(): mesh = generate_gmsh( FileSource("ellipsoid.step"), 2, order=2, other_options=["-string", "Mesh.CharacteristicLengthMax = %g;" % h]) + print("%d elements" % mesh.nelements) from meshmode.mesh.processing import find_bounding_box bbox_min, bbox_max = find_bounding_box(mesh) @@ -48,7 +49,7 @@ def main(): cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) qbx = QBXLayerPotentialSource(density_discr, 4*target_order, qbx_order, - fmm_order=qbx_order, fmm_backend="fmmlib") + fmm_order=qbx_order + 3, fmm_backend="fmmlib") from pytential.symbolic.pde.maxwell import PECAugmentedMFIEOperator pde_op = PECAugmentedMFIEOperator() @@ -80,13 +81,39 @@ def main(): #mlab.figure(bgcolor=(1, 1, 1)) if 1: - fplot = FieldPlotter(bbox_center, extent=1.5*bbox_size, npoints=30) + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, density_discr, target_order) + + bdry_normals = bind(density_discr, sym.normal(3))(queue)\ + .as_vector(dtype=object) + + bdry_vis.write_vtk_file("source.vtu", [ + ("sigma", sigma), + ("bdry_normals", bdry_normals), + ]) + + fplot = FieldPlotter(bbox_center, extent=2*bbox_size, npoints=(150, 150, 1)) qbx_stick_out = qbx.copy(target_stick_out_factor=0.1) from pytential.target import PointsTarget - fld_in_vol = bind( - (qbx_stick_out, PointsTarget(fplot.points)), - repr_op)(queue, jt=jt, rho=rho, k=k).get() + from pytential.qbx import QBXTargetAssociationFailedException + try: + fld_in_vol = bind( + (qbx_stick_out, PointsTarget(fplot.points)), + sym.S(pde_op.kernel, rho_sym, k=sym.var("k"), + qbx_forced_limit=None) + )(queue, jt=jt, rho=rho, k=k) + except QBXTargetAssociationFailedException as e: + fplot.write_vtk_file( + "failed-targets.vts", + [ + ("failed_targets", e.failed_target_flags.get(queue)) + ]) + raise + + fld_in_vol = fld_in_vol.get() + # fld_in_vol = sym.make_obj_array( + # [fiv.get() for fiv in fld_in_vol]) #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) fplot.write_vtk_file( @@ -96,17 +123,6 @@ def main(): ] ) - bdry_normals = bind(density_discr, sym.normal())(queue)\ - .as_vector(dtype=object) - - from meshmode.discretization.visualization import make_visualizer - bdry_vis = make_visualizer(queue, density_discr, target_order) - - bdry_vis.write_vtk_file("source.vtu", [ - ("sigma", sigma), - ("bdry_normals", bdry_normals), - ]) - if __name__ == "__main__": main() -- GitLab From d3d89778ad995ea81349dffd456aca075c64d236 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 5 Jun 2017 22:00:12 -0400 Subject: [PATCH 09/13] Get FMMlib QBX FMM to the point of evaluating a credible-looking 3D Helmholtz potential --- pytential/qbx/fmm.py | 14 +++++++++----- pytential/qbx/fmmlib.py | 41 ++++++++++++++++++++--------------------- 2 files changed, 29 insertions(+), 26 deletions(-) diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index 858c3e2d..503aeeb3 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -498,13 +498,17 @@ def drive_fmm(expansion_wrangler, src_weights): nqbtl = geo_data.non_qbx_box_target_lists() + # FIXME from pytools.obj_array import make_obj_array all_potentials_in_tree_order = make_obj_array([ - cl.array.zeros( - wrangler.queue, - tree.ntargets, - dtype=wrangler.dtype) - for k in wrangler.code.out_kernels]) + wrangler.full_potential_zeros() + ]) + qbx_potentials = make_obj_array([ + qbx_potentials + ]) + non_qbx_potentials = make_obj_array([ + non_qbx_potentials + ]) for ap_i, nqp_i in zip(all_potentials_in_tree_order, non_qbx_potentials): ap_i[nqbtl.unfiltered_from_filtered_target_indices] = nqp_i diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index f7497fd9..cdbcf48a 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -101,6 +101,11 @@ class ToHostTransferredGeoDataWrapper(object): def center_to_tree_targets(self): return self.geo_data.center_to_tree_targets().get(queue=self.queue) + @memoize_method + def all_targets(self): + """All (not just non-QBX) targets packaged into a single array.""" + return np.array(list(self.tree().targets)) + # }}} @@ -342,39 +347,33 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): geo_data = self.geo_data ctt = geo_data.center_to_tree_targets() + global_qbx_centers = geo_data.global_qbx_centers() + qbx_centers = geo_data.centers() + + all_targets = geo_data.all_targets() rscale = 1 # FIXME taeval = self.get_expn_eval_routine("ta") - for iglobal_center - src_icenter = global_qbx_centers[iglobal_center] + for isrc_center, src_icenter in enumerate(global_qbx_centers): + for icenter_tgt in range( + ctt.starts[src_icenter], + ctt.starts[src_icenter+1]): - icenter_tgt_start = center_to_targets_starts[src_icenter] - icenter_tgt_end = center_to_targets_starts[src_icenter+1] - - for icenter_tgt - - center_itgt = center_to_targets_lists[icenter_tgt] + center_itgt = ctt.lists[icenter_tgt] center = qbx_centers[:, src_icenter] - b[idim] = targets[idim, center_itgt] - center[idim] - """] + [""" - <> coeff{i} = qbx_expansions[src_icenter, {i}] - """.format(i=i) for i in range(ncoeffs)] + [ + pot[center_itgt] += taeval(self.helmholtz_k, rscale, + center, qbx_expansions[src_icenter], + all_targets[:, center_itgt]) - ] + loopy_insns + [""" + return pot - result[{i},center_itgt] = kernel_scaling * result_{i}_p \ - {{id_prefix=write_result}} - """.format(i=i) for i in range(len(result_names))] + [""" - tmp_pot = taeval(self.helmholtz_k, rscale, - center, qbx_expansions[src_icenter], - self._get_targets(tgt_pslice)) + def finalize_potential(self, potential): + return cl.array.to_device(self.queue, potential) - end - end # }}} # }}} -- GitLab From 83d6b866835762c0717ad95c19bfa7d60a47ec61 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 5 Jun 2017 22:03:57 -0400 Subject: [PATCH 10/13] Add ellipsoid test geometry --- examples/ellipsoid.step | 138 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100644 examples/ellipsoid.step diff --git a/examples/ellipsoid.step b/examples/ellipsoid.step new file mode 100644 index 00000000..883cd89f --- /dev/null +++ b/examples/ellipsoid.step @@ -0,0 +1,138 @@ +ISO-10303-21; +HEADER; +FILE_DESCRIPTION(('FreeCAD Model'),'2;1'); +FILE_NAME('/home/andreas/ellipsoid.step','2017-06-05T10:39:15',('Author' + ),(''),'Open CASCADE STEP processor 6.8','FreeCAD','Unknown'); +FILE_SCHEMA(('AUTOMOTIVE_DESIGN_CC2 { 1 2 10303 214 -1 1 5 4 }')); +ENDSEC; +DATA; +#1 = APPLICATION_PROTOCOL_DEFINITION('committee draft', + 'automotive_design',1997,#2); +#2 = APPLICATION_CONTEXT( + 'core data for automotive mechanical design processes'); +#3 = SHAPE_DEFINITION_REPRESENTATION(#4,#10); +#4 = PRODUCT_DEFINITION_SHAPE('','',#5); +#5 = PRODUCT_DEFINITION('design','',#6,#9); +#6 = PRODUCT_DEFINITION_FORMATION('','',#7); +#7 = PRODUCT('ASSEMBLY','ASSEMBLY','',(#8)); +#8 = MECHANICAL_CONTEXT('',#2,'mechanical'); +#9 = PRODUCT_DEFINITION_CONTEXT('part definition',#2,'design'); +#10 = SHAPE_REPRESENTATION('',(#11,#15),#19); +#11 = AXIS2_PLACEMENT_3D('',#12,#13,#14); +#12 = CARTESIAN_POINT('',(0.,0.,0.)); +#13 = DIRECTION('',(0.,0.,1.)); +#14 = DIRECTION('',(1.,0.,-0.)); +#15 = AXIS2_PLACEMENT_3D('',#16,#17,#18); +#16 = CARTESIAN_POINT('',(0.,0.,0.)); +#17 = DIRECTION('',(0.,0.,1.)); +#18 = DIRECTION('',(1.,0.,0.)); +#19 = ( GEOMETRIC_REPRESENTATION_CONTEXT(3) +GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT((#23)) GLOBAL_UNIT_ASSIGNED_CONTEXT( +(#20,#21,#22)) REPRESENTATION_CONTEXT('Context #1', + '3D Context with UNIT and UNCERTAINTY') ); +#20 = ( LENGTH_UNIT() NAMED_UNIT(*) SI_UNIT(.MILLI.,.METRE.) ); +#21 = ( NAMED_UNIT(*) PLANE_ANGLE_UNIT() SI_UNIT($,.RADIAN.) ); +#22 = ( NAMED_UNIT(*) SI_UNIT($,.STERADIAN.) SOLID_ANGLE_UNIT() ); +#23 = UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-07),#20, + 'distance_accuracy_value','confusion accuracy'); +#24 = PRODUCT_TYPE('part',$,(#7)); +#25 = ADVANCED_BREP_SHAPE_REPRESENTATION('',(#11,#26),#69); +#26 = MANIFOLD_SOLID_BREP('',#27); +#27 = CLOSED_SHELL('',(#28)); +#28 = ADVANCED_FACE('',(#29),#33,.T.); +#29 = FACE_BOUND('',#30,.T.); +#30 = VERTEX_LOOP('',#31); +#31 = VERTEX_POINT('',#32); +#32 = CARTESIAN_POINT('',(2.449293598295E-16,-5.999039130647E-32,-2.)); +#33 = ( BOUNDED_SURFACE() B_SPLINE_SURFACE(2,2,( + (#34,#35,#36,#37,#38) + ,(#39,#40,#41,#42,#43) + ,(#44,#45,#46,#47,#48) + ,(#49,#50,#51,#52,#53) + ,(#54,#55,#56,#57,#58) + ,(#59,#60,#61,#62,#63) + ,(#64,#65,#66,#67,#68 +)),.UNSPECIFIED.,.T.,.F.,.F.) B_SPLINE_SURFACE_WITH_KNOTS((1,2,2,2,2,1), + (3,2,3),(-2.094395102393,0.,2.094395102393,4.188790204786, + 6.28318530718,8.377580409573),(-1.570796326795,0.,1.570796326795), +.UNSPECIFIED.) GEOMETRIC_REPRESENTATION_ITEM() RATIONAL_B_SPLINE_SURFACE +(( + (1.,0.707106781187,1.,0.707106781187,1.) + ,(0.5,0.353553390593,0.5,0.353553390593,0.5) + ,(1.,0.707106781187,1.,0.707106781187,1.) + ,(0.5,0.353553390593,0.5,0.353553390593,0.5) + ,(1.,0.707106781187,1.,0.707106781187,1.) + ,(0.5,0.353553390593,0.5,0.353553390593,0.5) +,(1.,0.707106781187,1.,0.707106781187,1. +))) REPRESENTATION_ITEM('') SURFACE() ); +#34 = CARTESIAN_POINT('',(2.449293598295E-16,0.,-2.)); +#35 = CARTESIAN_POINT('',(4.,0.,-2.)); +#36 = CARTESIAN_POINT('',(4.,0.,0.)); +#37 = CARTESIAN_POINT('',(4.,0.,2.)); +#38 = CARTESIAN_POINT('',(2.449293598295E-16,0.,2.)); +#39 = CARTESIAN_POINT('',(2.449293598295E-16,4.2423009549E-16,-2.)); +#40 = CARTESIAN_POINT('',(4.,6.928203230276,-2.)); +#41 = CARTESIAN_POINT('',(4.,6.928203230276,0.)); +#42 = CARTESIAN_POINT('',(4.,6.928203230276,2.)); +#43 = CARTESIAN_POINT('',(2.449293598295E-16,4.2423009549E-16,2.)); +#44 = CARTESIAN_POINT('',(-1.224646799147E-16,2.12115047745E-16,-2.)); +#45 = CARTESIAN_POINT('',(-2.,3.464101615138,-2.)); +#46 = CARTESIAN_POINT('',(-2.,3.464101615138,0.)); +#47 = CARTESIAN_POINT('',(-2.,3.464101615138,2.)); +#48 = CARTESIAN_POINT('',(-1.224646799147E-16,2.12115047745E-16,2.)); +#49 = CARTESIAN_POINT('',(-4.898587196589E-16,5.999039130647E-32,-2.)); +#50 = CARTESIAN_POINT('',(-8.,9.797174393179E-16,-2.)); +#51 = CARTESIAN_POINT('',(-8.,9.797174393179E-16,0.)); +#52 = CARTESIAN_POINT('',(-8.,9.797174393179E-16,2.)); +#53 = CARTESIAN_POINT('',(-4.898587196589E-16,5.999039130647E-32,2.)); +#54 = CARTESIAN_POINT('',(-1.224646799147E-16,-2.12115047745E-16,-2.)); +#55 = CARTESIAN_POINT('',(-2.,-3.464101615138,-2.)); +#56 = CARTESIAN_POINT('',(-2.,-3.464101615138,0.)); +#57 = CARTESIAN_POINT('',(-2.,-3.464101615138,2.)); +#58 = CARTESIAN_POINT('',(-1.224646799147E-16,-2.12115047745E-16,2.)); +#59 = CARTESIAN_POINT('',(2.449293598295E-16,-4.2423009549E-16,-2.)); +#60 = CARTESIAN_POINT('',(4.,-6.928203230276,-2.)); +#61 = CARTESIAN_POINT('',(4.,-6.928203230276,0.)); +#62 = CARTESIAN_POINT('',(4.,-6.928203230276,2.)); +#63 = CARTESIAN_POINT('',(2.449293598295E-16,-4.2423009549E-16,2.)); +#64 = CARTESIAN_POINT('',(2.449293598295E-16,0.,-2.)); +#65 = CARTESIAN_POINT('',(4.,0.,-2.)); +#66 = CARTESIAN_POINT('',(4.,0.,0.)); +#67 = CARTESIAN_POINT('',(4.,0.,2.)); +#68 = CARTESIAN_POINT('',(2.449293598295E-16,0.,2.)); +#69 = ( GEOMETRIC_REPRESENTATION_CONTEXT(3) +GLOBAL_UNCERTAINTY_ASSIGNED_CONTEXT((#73)) GLOBAL_UNIT_ASSIGNED_CONTEXT( +(#70,#71,#72)) REPRESENTATION_CONTEXT('Context #1', + '3D Context with UNIT and UNCERTAINTY') ); +#70 = ( LENGTH_UNIT() NAMED_UNIT(*) SI_UNIT(.MILLI.,.METRE.) ); +#71 = ( NAMED_UNIT(*) PLANE_ANGLE_UNIT() SI_UNIT($,.RADIAN.) ); +#72 = ( NAMED_UNIT(*) SI_UNIT($,.STERADIAN.) SOLID_ANGLE_UNIT() ); +#73 = UNCERTAINTY_MEASURE_WITH_UNIT(LENGTH_MEASURE(1.E-07),#70, + 'distance_accuracy_value','confusion accuracy'); +#74 = SHAPE_DEFINITION_REPRESENTATION(#75,#25); +#75 = PRODUCT_DEFINITION_SHAPE('','',#76); +#76 = PRODUCT_DEFINITION('design','',#77,#80); +#77 = PRODUCT_DEFINITION_FORMATION('','',#78); +#78 = PRODUCT('Ellipsoid','Ellipsoid','',(#79)); +#79 = MECHANICAL_CONTEXT('',#2,'mechanical'); +#80 = PRODUCT_DEFINITION_CONTEXT('part definition',#2,'design'); +#81 = CONTEXT_DEPENDENT_SHAPE_REPRESENTATION(#82,#84); +#82 = ( REPRESENTATION_RELATIONSHIP('','',#25,#10) +REPRESENTATION_RELATIONSHIP_WITH_TRANSFORMATION(#83) +SHAPE_REPRESENTATION_RELATIONSHIP() ); +#83 = ITEM_DEFINED_TRANSFORMATION('','',#11,#15); +#84 = PRODUCT_DEFINITION_SHAPE('Placement','Placement of an item',#85); +#85 = NEXT_ASSEMBLY_USAGE_OCCURRENCE('1','=>[0:1:1:2]','',#5,#76,$); +#86 = PRODUCT_TYPE('part',$,(#78)); +#87 = MECHANICAL_DESIGN_GEOMETRIC_PRESENTATION_REPRESENTATION('',(#88), + #69); +#88 = STYLED_ITEM('color',(#89),#28); +#89 = PRESENTATION_STYLE_ASSIGNMENT((#90)); +#90 = SURFACE_STYLE_USAGE(.BOTH.,#91); +#91 = SURFACE_SIDE_STYLE('',(#92)); +#92 = SURFACE_STYLE_FILL_AREA(#93); +#93 = FILL_AREA_STYLE('',(#94)); +#94 = FILL_AREA_STYLE_COLOUR('',#95); +#95 = COLOUR_RGB('',0.800000011921,0.800000011921,0.800000011921); +ENDSEC; +END-ISO-10303-21; -- GitLab From c32669fa3540a471493e9de8093e22a08b6a9d77 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 6 Jun 2017 12:01:38 -0400 Subject: [PATCH 11/13] Multi-output pyfmmlib QBX FMM working for a single output, repair top-level compatibility with sumpy FMM --- pytential/qbx/fmm.py | 16 ++------- pytential/qbx/fmmlib.py | 73 ++++++++++++++++++++++------------------- 2 files changed, 42 insertions(+), 47 deletions(-) diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index 503aeeb3..249cc9b2 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -126,7 +126,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), # {{{ data vector utilities - def potential_zeros(self): + def output_zeros(self): """This ought to be called ``non_qbx_potential_zeros``, but since it has to override the superclass's behavior to integrate seamlessly, it needs to be called just :meth:`potential_zeros`. @@ -142,7 +142,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), dtype=self.dtype) for k in self.code.out_kernels]) - def full_potential_zeros(self): + def full_output_zeros(self): # The superclass generates a full field of zeros, for all # (not just non-QBX) targets. return SumpyExpansionWrangler.potential_zeros(self) @@ -498,17 +498,7 @@ def drive_fmm(expansion_wrangler, src_weights): nqbtl = geo_data.non_qbx_box_target_lists() - # FIXME - from pytools.obj_array import make_obj_array - all_potentials_in_tree_order = make_obj_array([ - wrangler.full_potential_zeros() - ]) - qbx_potentials = make_obj_array([ - qbx_potentials - ]) - non_qbx_potentials = make_obj_array([ - non_qbx_potentials - ]) + all_potentials_in_tree_order = wrangler.full_output_zeros() for ap_i, nqp_i in zip(all_potentials_in_tree_order, non_qbx_potentials): ap_i[nqbtl.unfiltered_from_filtered_target_indices] = nqp_i diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index cdbcf48a..914e5389 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -27,6 +27,7 @@ from pytools import memoize_method import pyopencl as cl # noqa import pyopencl.array # noqa: F401 from boxtree.pyfmmlib_integration import HelmholtzExpansionWrangler +from sumpy.kernel import HelmholtzKernel class QBXFMMLibExpansionWranglerCodeContainer(object): @@ -45,17 +46,10 @@ class QBXFMMLibExpansionWranglerCodeContainer(object): source_extra_kwargs={}, kernel_extra_kwargs=None): - from sumpy.kernel import HelmholtzKernel - for out_knl in self.out_kernels: - if not isinstance(out_knl, HelmholtzKernel): - raise NotImplementedError( - "only the Helmholtz kernel is supported for now") - return QBXFMMLibHelmholtzExpansionWrangler(self, queue, geo_data, dtype, qbx_order, fmm_level_to_order, source_extra_kwargs, - kernel_extra_kwargs, - self.out_kernels) + kernel_extra_kwargs) # }}} @@ -112,21 +106,26 @@ class ToHostTransferredGeoDataWrapper(object): # {{{ fmmlib expansion wrangler class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): - def __init__(self, code_container, queue, geo_data, dtype, + def __init__(self, code, queue, geo_data, dtype, qbx_order, fmm_level_to_order, source_extra_kwargs, - kernel_extra_kwargs, - out_kernels): + kernel_extra_kwargs): - from pytools import single_valued - k_name = single_valued(out_knl.helmholtz_k_name for out_knl in out_kernels) - helmholtz_k = kernel_extra_kwargs[k_name] - - self.code_container = code_container + self.code = code self.queue = queue self.geo_data = ToHostTransferredGeoDataWrapper(queue, geo_data) self.qbx_order = qbx_order + for out_knl in self.code.out_kernels: + if not isinstance(out_knl, HelmholtzKernel) and out_knl.dim == 3: + raise NotImplementedError( + "only the 3D Helmholtz kernel is supported for now") + + from pytools import single_valued + k_name = single_valued(out_knl.helmholtz_k_name + for out_knl in code.out_kernels) + helmholtz_k = kernel_extra_kwargs[k_name] + self.level_orders = [ fmm_level_to_order(level) for level in range(self.geo_data.tree().nlevels)] @@ -144,28 +143,26 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): # FIXME nterms=fmm_level_to_order(0)) - def potential_zeros(self): - """This ought to be called ``non_qbx_potential_zeros``, but since + def output_zeros(self): + """This ought to be called ``non_qbx_output_zeros``, but since it has to override the superclass's behavior to integrate seamlessly, - it needs to be called just :meth:`potential_zeros`. + it needs to be called just :meth:`output_zeros`. """ nqbtl = self.geo_data.non_qbx_box_target_lists() - # from pytools.obj_array import make_obj_array - # return make_obj_array([ - # cl.array.zeros( - # self.queue, - # nqbtl.nfiltered_targets, - # dtype=self.dtype) - # for k in self.code.out_kernels]) + from pytools.obj_array import make_obj_array + return make_obj_array([ + np.zeros(nqbtl.nfiltered_targets, self.dtype) + for k in self.code.out_kernels]) - return np.zeros(nqbtl.nfiltered_targets, self.dtype) + def full_output_zeros(self): + """This includes QBX and non-QBX targets.""" - def full_potential_zeros(self): - # The superclass generates a full field of zeros, for all - # (not just non-QBX) targets. - return super(QBXFMMLibHelmholtzExpansionWrangler, self).potential_zeros() + from pytools.obj_array import make_obj_array + return make_obj_array([ + np.zeros(self.tree.ntargets, self.dtype) + for k in self.code.out_kernels]) def reorder_sources(self, source_array): source_array = source_array.get(queue=self.queue) @@ -179,6 +176,12 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): # Because this is a multi-stage, more complicated process that combines # potentials from non-QBX targets and QBX targets. + def add_potgrad_onto_output(self, output, output_slice, pot, grad): + assert (len(self.code.out_kernels) == 1 + and isinstance(self.code.out_kernels[0], HelmholtzKernel)) + + output[0][output_slice] += pot + # {{{ override target lists to only hit non-QBX targets def box_target_starts(self): @@ -343,7 +346,7 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): return qbx_expansions def eval_qbx_expansions(self, qbx_expansions): - pot = self.full_potential_zeros() + output = self.full_output_zeros() geo_data = self.geo_data ctt = geo_data.center_to_tree_targets() @@ -365,11 +368,13 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): center = qbx_centers[:, src_icenter] - pot[center_itgt] += taeval(self.helmholtz_k, rscale, + pot, grad = taeval(self.helmholtz_k, rscale, center, qbx_expansions[src_icenter], all_targets[:, center_itgt]) - return pot + self.add_potgrad_onto_output(output, center_itgt, pot, grad) + + return output def finalize_potential(self, potential): return cl.array.to_device(self.queue, potential) -- GitLab From 0b9253e2c4a80a883e3091d59f7d07219c3d2263 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 6 Jun 2017 14:39:42 -0400 Subject: [PATCH 12/13] Fix the sumpy-bassed QBX FMM after the output_zeros rename --- pytential/qbx/fmm.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index 249cc9b2..c7d6c029 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -127,9 +127,9 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), # {{{ data vector utilities def output_zeros(self): - """This ought to be called ``non_qbx_potential_zeros``, but since + """This ought to be called ``non_qbx_output_zeros``, but since it has to override the superclass's behavior to integrate seamlessly, - it needs to be called just :meth:`potential_zeros`. + it needs to be called just :meth:`output_zeros`. """ nqbtl = self.geo_data.non_qbx_box_target_lists() @@ -145,7 +145,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), def full_output_zeros(self): # The superclass generates a full field of zeros, for all # (not just non-QBX) targets. - return SumpyExpansionWrangler.potential_zeros(self) + return SumpyExpansionWrangler.output_zeros(self) def qbx_local_expansion_zeros(self): order = self.qbx_order @@ -306,7 +306,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), return qbx_expansions def eval_qbx_expansions(self, qbx_expansions): - pot = self.full_potential_zeros() + pot = self.full_output_zeros() geo_data = self.geo_data if len(geo_data.global_qbx_centers()) == 0: -- GitLab From a6fdfdd4532135780bf852d90cf707d525edb4c0 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 6 Jun 2017 14:40:30 -0400 Subject: [PATCH 13/13] out_vectors support for the fmmlib qbx FMM --- examples/maxwell.py | 31 ++++++++++++++------ pytential/qbx/fmmlib.py | 63 ++++++++++++++++++++++++++++++++++------- 2 files changed, 74 insertions(+), 20 deletions(-) diff --git a/examples/maxwell.py b/examples/maxwell.py index 033ff902..108a96c3 100644 --- a/examples/maxwell.py +++ b/examples/maxwell.py @@ -31,6 +31,11 @@ def main(): mesh = generate_gmsh( FileSource("ellipsoid.step"), 2, order=2, other_options=["-string", "Mesh.CharacteristicLengthMax = %g;" % h]) + + from meshmode.mesh.processing import perform_flips + # Flip elements--gmsh generates inside-out geometry. + mesh = perform_flips(mesh, np.ones(mesh.nelements)) + print("%d elements" % mesh.nelements) from meshmode.mesh.processing import find_bounding_box @@ -49,15 +54,15 @@ def main(): cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) qbx = QBXLayerPotentialSource(density_discr, 4*target_order, qbx_order, - fmm_order=qbx_order + 3, fmm_backend="fmmlib") + fmm_order=qbx_order + 10, fmm_backend="fmmlib") from pytential.symbolic.pde.maxwell import PECAugmentedMFIEOperator pde_op = PECAugmentedMFIEOperator() from pytential import bind, sym - jt_sym = sym.make_sym_vector("jt", 2) rho_sym = sym.var("rho") - repr_op = pde_op.scattered_volume_field(jt_sym, rho_sym) + # jt_sym = sym.make_sym_vector("jt", 2) + # repr_op = pde_op.scattered_volume_field(jt_sym, rho_sym) # {{{ make a density @@ -100,8 +105,16 @@ def main(): try: fld_in_vol = bind( (qbx_stick_out, PointsTarget(fplot.points)), - sym.S(pde_op.kernel, rho_sym, k=sym.var("k"), - qbx_forced_limit=None) + sym.make_obj_array([ + # sym.S(pde_op.kernel, rho_sym, k=sym.var("k"), + # qbx_forced_limit=None), + # sym.d_dx(3, sym.S(pde_op.kernel, rho_sym, k=sym.var("k"), + # qbx_forced_limit=None)), + sym.d_dy(3, sym.S(pde_op.kernel, rho_sym, k=sym.var("k"), + qbx_forced_limit=None)), + # sym.d_dz(3, sym.S(pde_op.kernel, rho_sym, k=sym.var("k"), + # qbx_forced_limit=None)), + ]) )(queue, jt=jt, rho=rho, k=k) except QBXTargetAssociationFailedException as e: fplot.write_vtk_file( @@ -111,15 +124,15 @@ def main(): ]) raise - fld_in_vol = fld_in_vol.get() - # fld_in_vol = sym.make_obj_array( - # [fiv.get() for fiv in fld_in_vol]) + fld_in_vol = sym.make_obj_array( + [fiv.get() for fiv in fld_in_vol]) #fplot.show_scalar_in_mayavi(fld_in_vol.real, max_val=5) fplot.write_vtk_file( "potential.vts", [ - ("potential", fld_in_vol) + ("potential", fld_in_vol[0]), + ("grad", fld_in_vol[1:]), ] ) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index 914e5389..14d523f1 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -113,17 +113,45 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): self.code = code self.queue = queue + + # FMMLib is CPU-only. This wrapper gets the geometry out of + # OpenCL-land. self.geo_data = ToHostTransferredGeoDataWrapper(queue, geo_data) + self.qbx_order = qbx_order + # {{{ digest out_kernels + + from sumpy.kernel import AxisTargetDerivative + + k_names = [] + + def is_supported_helmknl(knl): + result = isinstance(knl, HelmholtzKernel) and knl.dim == 3 + if result: + k_names.append(knl.helmholtz_k_name) + return result + + ifgrad = False + outputs = [] for out_knl in self.code.out_kernels: - if not isinstance(out_knl, HelmholtzKernel) and out_knl.dim == 3: + if is_supported_helmknl(out_knl): + outputs.append(()) + elif (isinstance(out_knl, AxisTargetDerivative) + and is_supported_helmknl(out_knl.inner_kernel)): + outputs.append((out_knl.axis,)) + ifgrad = True + else: raise NotImplementedError( - "only the 3D Helmholtz kernel is supported for now") + "only the 3D Helmholtz kernel and its target derivatives " + "are supported for now") + + self.outputs = outputs + + # }}} from pytools import single_valued - k_name = single_valued(out_knl.helmholtz_k_name - for out_knl in code.out_kernels) + k_name = single_valued(k_names) helmholtz_k = kernel_extra_kwargs[k_name] self.level_orders = [ @@ -135,13 +163,16 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): assert single_valued(self.level_orders) super(QBXFMMLibHelmholtzExpansionWrangler, self).__init__( - # FMMLib is CPU-only--get the tree out of OpenCL-land self.geo_data.tree(), helmholtz_k=helmholtz_k, # FIXME - nterms=fmm_level_to_order(0)) + nterms=fmm_level_to_order(0), + + ifgrad=ifgrad) + + # {{{ data vector helpers def output_zeros(self): """This ought to be called ``non_qbx_output_zeros``, but since @@ -154,7 +185,7 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): from pytools.obj_array import make_obj_array return make_obj_array([ np.zeros(nqbtl.nfiltered_targets, self.dtype) - for k in self.code.out_kernels]) + for k in self.outputs]) def full_output_zeros(self): """This includes QBX and non-QBX targets.""" @@ -162,7 +193,7 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): from pytools.obj_array import make_obj_array return make_obj_array([ np.zeros(self.tree.ntargets, self.dtype) - for k in self.code.out_kernels]) + for k in self.outputs]) def reorder_sources(self, source_array): source_array = source_array.get(queue=self.queue) @@ -177,10 +208,20 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): # potentials from non-QBX targets and QBX targets. def add_potgrad_onto_output(self, output, output_slice, pot, grad): - assert (len(self.code.out_kernels) == 1 - and isinstance(self.code.out_kernels[0], HelmholtzKernel)) + for i_out, out in enumerate(self.outputs): + if len(out) == 0: + output[i_out][output_slice] += pot + elif len(out) == 1: + axis, = out + if isinstance(grad, np.ndarray): + output[i_out][output_slice] += grad[axis] + else: + assert grad == 0 + else: + raise ValueError("element '%s' of outputs array not " + "understood" % out) - output[0][output_slice] += pot + # }}} # {{{ override target lists to only hit non-QBX targets -- GitLab