From 232104954cb1d36b34db1b15ca43310ca46d9ee9 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 24 Jun 2017 09:26:42 -0500 Subject: [PATCH 01/24] Move fmmlib P2QBXL to parallel implementation --- pytential/qbx/fmmlib.py | 98 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 87 insertions(+), 11 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index fd1dc4db..e474897f 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -258,33 +258,109 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): qbx_center_to_target_box = geo_data.qbx_center_to_target_box() qbx_centers = geo_data.centers() - formta = self.get_routine("%ddformta") + # {{{ parallel + print("par data prep") + + center_source_counts = [0] 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] + source_count = sum( + self.tree.box_source_counts_nonchild[lists[isrc_box]] + for isrc_box in range(isrc_box_start, isrc_box_stop)) + + center_source_counts.append(source_count) + + center_source_counts = np.array(center_source_counts) + center_source_starts = np.cumsum(center_source_counts) + nsources_total = center_source_starts[-1] + + sources = np.empty((3, nsources_total), dtype=np.float64) + charges = np.empty(nsources_total, dtype=np.complex128) - ctr_coeffs = 0 + isource = 0 + for itgt_center, tgt_icenter in enumerate(geo_data.global_qbx_centers()): + assert isource == center_source_starts[itgt_center] + itgt_box = qbx_center_to_target_box[tgt_icenter] + + isrc_box_start = starts[itgt_box] + isrc_box_stop = starts[itgt_box+1] for isrc_box in range(isrc_box_start, isrc_box_stop): src_ibox = lists[isrc_box] src_pslice = self._get_source_slice(src_ibox) + ns = self.tree.box_source_counts_nonchild[src_ibox] + sources[:, isource:isource+ns] = self._get_sources(src_pslice) + charges[isource:isource+ns] = src_weights[src_pslice] + isource += ns + + centers = qbx_centers[:, geo_data.global_qbx_centers()] + + rscale_vec = np.empty(len(center_source_counts) - 1, dtype=np.float64) + rscale_vec.fill(rscale) # FIXME + + formta_vec = self.get_vec_routine("%ddformta") + print("par data prep done") + + ier, loc_exp_pre = formta_vec( + self.helmholtz_k, rscale_vec, + sources, center_source_starts[:-1], + charges, center_source_starts[:-1], + center_source_counts[1:], + centers, self.nterms) + + if np.any(ier != 0): + raise RuntimeError("formta returned an error") + + loc_exp_pre = np.moveaxis(loc_exp_pre, -1, 0) + local_exps[geo_data.global_qbx_centers()] = loc_exp_pre + + if 0: + # {{{ sequential - 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") + formta = self.get_routine("%ddformta") - ctr_coeffs += coeffs + local_exps_1 = np.zeros_like(local_exps) - local_exps[tgt_icenter] += ctr_coeffs + 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_1[tgt_icenter] += ctr_coeffs + + diff = local_exps - local_exps_1 + + # At least in 3D, this will have two "triangles" of non-zeros. + print(diff) + + # }}} + + # }}} return local_exps -- GitLab From 8c98d43536266841d03a0b394c76b72530f336fe Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 25 Jun 2017 17:50:51 -0500 Subject: [PATCH 02/24] Track boxtree switch to 'transposed' expansion shape, start on M2QBXL parallelization --- pytential/qbx/fmmlib.py | 95 +++++++++++++++++++++++++++++------------ 1 file changed, 67 insertions(+), 28 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index e474897f..14130b45 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -317,8 +317,7 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): if np.any(ier != 0): raise RuntimeError("formta returned an error") - loc_exp_pre = np.moveaxis(loc_exp_pre, -1, 0) - local_exps[geo_data.global_qbx_centers()] = loc_exp_pre + local_exps[geo_data.global_qbx_centers()] = loc_exp_pre.T if 0: # {{{ sequential @@ -372,40 +371,80 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): qbx_centers = geo_data.centers() centers = self.tree.box_centers - mploc = self.get_translation_routine("%ddmploc") + if 0: + mploc = self.get_translation_routine("%ddmploc", vec_suffix="_csr") - 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) + for isrc_level, ssn in enumerate(geo_data.traversal().sep_smaller_by_level): - # FIXME - rscale = 1 + # FIXME + rscale = 1 - kwargs = {} - if self.dim == 3: - # FIXME Is this right? - kwargs["radius"] = self.tree.root_extent * 2**(-isrc_level) + 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 + print("par data prep lev %d" % isrc_level) + itgt_center = np.arange(len(geo_data.global_qbx_centers())) + tgt_icenter = geo_data.global_qbx_centers() icontaining_tgt_box = qbx_center_to_target_box[tgt_icenter] - tgt_center = qbx_centers[:, tgt_icenter] + tgt_centers = qbx_centers[:, tgt_icenter] + + print("end par data prep") + + ctr_loc = ctr_loc + mploc( + self.helmholtz_k, + rscale, src_center, multipole_exps[src_ibox].T, + rscale, tgt_center, self.nterms, **kwargs)[..., 0].T + + + + if 1: + # {{{ sequential + + mploc = self.get_translation_routine("%ddmploc") + + local_exps_1 = self.qbx_local_expansion_zeros() + + 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 isrc_box in range( - ssn.starts[icontaining_tgt_box], - ssn.starts[icontaining_tgt_box+1]): + for itgt_center, tgt_icenter in enumerate( + geo_data.global_qbx_centers()): + ctr_loc = 0 - src_ibox = ssn.lists[isrc_box] - src_center = centers[:, src_ibox] + icontaining_tgt_box = qbx_center_to_target_box[tgt_icenter] - ctr_loc = ctr_loc + mploc( - self.helmholtz_k, - rscale, src_center, multipole_exps[src_ibox], - rscale, tgt_center, self.nterms, **kwargs)[..., 0] + 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].T, + rscale, tgt_center, self.nterms, **kwargs)[..., 0].T + + local_exps_1[tgt_icenter] += ctr_loc + + # }}} - local_exps[tgt_icenter] += ctr_loc + return local_exps_1 return local_exps @@ -455,8 +494,8 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): 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] + rscale, src_center, local_exps[src_ibox].T, + rscale, tgt_center, local_order, **kwargs)[..., 0].T qbx_expansions[tgt_icenter] += tmp_loc_exp @@ -486,7 +525,7 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): center = qbx_centers[:, src_icenter] pot, grad = taeval(self.helmholtz_k, rscale, - center, qbx_expansions[src_icenter], + center, qbx_expansions[src_icenter].T, all_targets[:, center_itgt]) self.add_potgrad_onto_output(output, center_itgt, pot, grad) -- GitLab From 85e68cb71ba2fdbeb6e0163efc2f812d52cd6256 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 27 Jun 2017 15:19:27 -0500 Subject: [PATCH 03/24] Parallelize fmmlib m2qbxl --- pytential/qbx/fmmlib.py | 120 +++++++++++++++++++++++++++++++--------- 1 file changed, 94 insertions(+), 26 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index 14130b45..cdde0c58 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -308,11 +308,17 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): print("par data prep done") ier, loc_exp_pre = formta_vec( - self.helmholtz_k, rscale_vec, - sources, center_source_starts[:-1], - charges, center_source_starts[:-1], - center_source_counts[1:], - centers, self.nterms) + zk=self.helmholtz_k, + rscale=rscale_vec, + + sources=sources, + sources_offsets=center_source_starts[:-1], + charges=charges, + charges_offsets=center_source_starts[:-1], + + nsources=center_source_counts[1:], + center=centers, + nterms=self.nterms) if np.any(ier != 0): raise RuntimeError("formta returned an error") @@ -371,46 +377,107 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): qbx_centers = geo_data.centers() centers = self.tree.box_centers - if 0: - mploc = self.get_translation_routine("%ddmploc", vec_suffix="_csr") + if 1: + # {{{ parallel - for isrc_level, ssn in enumerate(geo_data.traversal().sep_smaller_by_level): + mploc = self.get_translation_routine("%ddmploc", vec_suffix="_imany") + + 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 + + print("par data prep lev %d" % isrc_level) + + ngqbx_centers = len(geo_data.global_qbx_centers()) + tgt_icenter_vec = geo_data.global_qbx_centers() + icontaining_tgt_box_vec = qbx_center_to_target_box[tgt_icenter_vec] + + tgt_centers = qbx_centers[:, tgt_icenter_vec] + + # FIXME + rscale2 = np.ones(ngqbx_centers, np.float64) kwargs = {} if self.dim == 3: # FIXME Is this right? - kwargs["radius"] = self.tree.root_extent * 2**(-isrc_level) + kwargs["radius"] = ( + np.ones(ngqbx_centers) + * self.tree.root_extent * 2**(-isrc_level)) - print("par data prep lev %d" % isrc_level) + nsrc_boxes_per_gqbx_center = ( + ssn.starts[icontaining_tgt_box_vec+1] + - ssn.starts[icontaining_tgt_box_vec]) + nsrc_boxes = np.sum(nsrc_boxes_per_gqbx_center) + + src_boxes_starts = np.empty(ngqbx_centers+1, dtype=np.int32) + src_boxes_starts[0] = 0 + src_boxes_starts[1:] = np.cumsum(nsrc_boxes_per_gqbx_center) - itgt_center = np.arange(len(geo_data.global_qbx_centers())) - tgt_icenter = geo_data.global_qbx_centers() - icontaining_tgt_box = qbx_center_to_target_box[tgt_icenter] + # FIXME + rscale1 = np.ones(nsrc_boxes) + rscale1_offsets = np.arange(nsrc_boxes) + + src_ibox = np.empty(nsrc_boxes, dtype=np.int32) + for itgt_center, tgt_icenter in enumerate( + geo_data.global_qbx_centers()): + icontaining_tgt_box = qbx_center_to_target_box[tgt_icenter] + src_ibox[ + src_boxes_starts[itgt_center]: + src_boxes_starts[itgt_center+1]] = ( + ssn.lists[ + ssn.starts[icontaining_tgt_box]: + ssn.starts[icontaining_tgt_box+1]]) - tgt_centers = qbx_centers[:, tgt_icenter] + del itgt_center + del tgt_icenter + del icontaining_tgt_box print("end par data prep") - ctr_loc = ctr_loc + mploc( - self.helmholtz_k, - rscale, src_center, multipole_exps[src_ibox].T, - rscale, tgt_center, self.nterms, **kwargs)[..., 0].T + # These get max'd/added onto: pass initialized versions. + ier = np.zeros(ngqbx_centers, dtype=np.int32) + expn2 = np.zeros( + (ngqbx_centers,) + self.expansion_shape(self.qbx_order), + dtype=self.dtype) + expn2 = mploc( + zk=self.helmholtz_k, + rscale1=rscale1, + rscale1_offsets=rscale1_offsets, + rscale1_starts=src_boxes_starts, - if 1: + center1=centers, + center1_offsets=src_ibox, + center1_starts=src_boxes_starts, + + expn1=source_mpoles_view.T, + expn1_offsets=src_ibox - source_level_start_ibox, + expn1_starts=src_boxes_starts, + + rscale2=rscale2, + center2=tgt_centers, # FIXME: wrong layout, will copy + expn2=expn2.T, + ier=ier, **kwargs).T + + local_exps[geo_data.global_qbx_centers()] += expn2 + + # }}} + + if 0: # {{{ sequential mploc = self.get_translation_routine("%ddmploc") local_exps_1 = self.qbx_local_expansion_zeros() - 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) + 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 @@ -437,14 +504,15 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): ctr_loc = ctr_loc + mploc( self.helmholtz_k, - rscale, src_center, multipole_exps[src_ibox].T, + rscale, src_center, + source_mpoles_view[src_ibox - source_level_start_ibox].T, rscale, tgt_center, self.nterms, **kwargs)[..., 0].T local_exps_1[tgt_icenter] += ctr_loc - # }}} + # diff = local_exps - local_exps_1 - return local_exps_1 + # }}} return local_exps -- GitLab From e965e6df55013338c27fba5b4bed28d9f96a948a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 27 Jun 2017 16:30:21 -0500 Subject: [PATCH 04/24] fmmlib: Minor p2qbxl simplification [ci skip] --- pytential/qbx/fmmlib.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index cdde0c58..5723a886 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -395,8 +395,6 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): tgt_icenter_vec = geo_data.global_qbx_centers() icontaining_tgt_box_vec = qbx_center_to_target_box[tgt_icenter_vec] - tgt_centers = qbx_centers[:, tgt_icenter_vec] - # FIXME rscale2 = np.ones(ngqbx_centers, np.float64) @@ -459,7 +457,8 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): expn1_starts=src_boxes_starts, rscale2=rscale2, - center2=tgt_centers, # FIXME: wrong layout, will copy + # FIXME: center2 has wrong layout, will copy + center2=qbx_centers[:, tgt_icenter_vec], expn2=expn2.T, ier=ier, **kwargs).T -- GitLab From 28f4d9e368aa349735d856bb9ba671cb79070f5f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 06:15:57 -0500 Subject: [PATCH 05/24] Refactor p2qbxl driver towards making interaction list cacheable --- pytential/qbx/fmm.py | 12 +++-- pytential/qbx/fmmlib.py | 110 +++++++++++++++++----------------------- 2 files changed, 54 insertions(+), 68 deletions(-) diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index c7d6c029..68060da2 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -191,13 +191,18 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), # {{{ qbx-related - def form_global_qbx_locals(self, starts, lists, src_weights): + def form_global_qbx_locals(self, src_weights): local_exps = self.qbx_local_expansion_zeros() geo_data = self.geo_data if len(geo_data.global_qbx_centers()) == 0: return local_exps + traversal = geo_data.traversal() + + starts = traversal.neighbor_source_boxes_starts + lists = traversal.neighbor_source_boxes_lists + kwargs = self.extra_kwargs.copy() kwargs.update(self.box_source_list_kwargs()) @@ -472,10 +477,7 @@ def drive_fmm(expansion_wrangler, src_weights): # {{{ wrangle qbx expansions 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) + qbx_expansions = wrangler.form_global_qbx_locals(src_weights) logger.info("translate from list 3 multipoles to qbx local expansions") qbx_expansions = qbx_expansions + \ diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index 5723a886..c97f88f0 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -23,13 +23,21 @@ THE SOFTWARE. """ import numpy as np -from pytools import memoize_method +from pytools import memoize_method, Record import pyopencl as cl # noqa import pyopencl.array # noqa: F401 from boxtree.pyfmmlib_integration import HelmholtzExpansionWrangler from sumpy.kernel import HelmholtzKernel +import logging +logger = logging.getLogger(__name__) + + +class P2QBXLInfo(Record): + pass + + class QBXFMMLibExpansionWranglerCodeContainer(object): def __init__(self, cl_context, multipole_expansion_factory, local_expansion_factory, @@ -239,29 +247,28 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): # }}} - # {{{ 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() + # {{{ p2qbxl + + #@memoize_method + def _info_for_form_global_qbx_locals(self, src_weights): + logger.info("preparing interaction list for p2qbxl: start") rscale = 1 # FIXME + geo_data = self.geo_data + traversal = geo_data.traversal() - if len(geo_data.global_qbx_centers()) == 0: - return local_exps + starts = traversal.neighbor_source_boxes_starts + lists = traversal.neighbor_source_boxes_lists qbx_center_to_target_box = geo_data.qbx_center_to_target_box() qbx_centers = geo_data.centers() - # {{{ parallel - - print("par data prep") - center_source_counts = [0] for itgt_center, tgt_icenter in enumerate(geo_data.global_qbx_centers()): itgt_box = qbx_center_to_target_box[tgt_icenter] @@ -304,71 +311,50 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): rscale_vec = np.empty(len(center_source_counts) - 1, dtype=np.float64) rscale_vec.fill(rscale) # FIXME - formta_vec = self.get_vec_routine("%ddformta") - print("par data prep done") - - ier, loc_exp_pre = formta_vec( - zk=self.helmholtz_k, - rscale=rscale_vec, + logger.info("preparing interaction list for p2qbxl: done") + return P2QBXLInfo( sources=sources, - sources_offsets=center_source_starts[:-1], + centers=centers, charges=charges, - charges_offsets=center_source_starts[:-1], - - nsources=center_source_counts[1:], - center=centers, - nterms=self.nterms) - - if np.any(ier != 0): - raise RuntimeError("formta returned an error") - - local_exps[geo_data.global_qbx_centers()] = loc_exp_pre.T + center_source_starts=center_source_starts, + center_source_counts=center_source_counts, + rscale_vec=rscale_vec, + ) - if 0: - # {{{ sequential - - formta = self.get_routine("%ddformta") - - local_exps_1 = np.zeros_like(local_exps) - - 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] + def form_global_qbx_locals(self, src_weights): + geo_data = self.geo_data - src_pslice = self._get_source_slice(src_ibox) + local_exps = self.qbx_local_expansion_zeros(src_weights) - 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") + if len(geo_data.global_qbx_centers()) == 0: + return local_exps - ctr_coeffs += coeffs + formta_vec = self.get_vec_routine("%ddformta") + info = self._info_for_form_global_qbx_locals() - local_exps_1[tgt_icenter] += ctr_coeffs + ier, loc_exp_pre = formta_vec( + zk=self.helmholtz_k, + rscale=info.rscale_vec, - diff = local_exps - local_exps_1 + sources=info.sources, + sources_offsets=info.center_source_starts[:-1], + charges=info.charges, + charges_offsets=info.center_source_starts[:-1], - # At least in 3D, this will have two "triangles" of non-zeros. - print(diff) + nsources=info.center_source_counts[1:], + center=info.centers, + nterms=self.nterms) - # }}} + if np.any(ier != 0): + raise RuntimeError("formta returned an error") - # }}} + local_exps[geo_data.global_qbx_centers()] = loc_exp_pre.T return local_exps + # }}} + def translate_box_multipoles_to_qbx_local(self, multipole_exps): local_exps = self.qbx_local_expansion_zeros() @@ -608,8 +594,6 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): return cl.array.to_device(self.queue, potential) * scale_factor - # }}} - # }}} # vim: foldmethod=marker -- GitLab From a79aaf70218378fd64aab322149e3b192641b708 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 06:19:31 -0500 Subject: [PATCH 06/24] Fix: Refactor p2qbxl driver towards making interaction list cacheable --- pytential/qbx/fmmlib.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index c97f88f0..6da05574 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -325,13 +325,13 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): def form_global_qbx_locals(self, src_weights): geo_data = self.geo_data - local_exps = self.qbx_local_expansion_zeros(src_weights) + local_exps = self.qbx_local_expansion_zeros() if len(geo_data.global_qbx_centers()) == 0: return local_exps formta_vec = self.get_vec_routine("%ddformta") - info = self._info_for_form_global_qbx_locals() + info = self._info_for_form_global_qbx_locals(src_weights) ier, loc_exp_pre = formta_vec( zk=self.helmholtz_k, -- GitLab From 3619ba482ab20d235b99a878d62f2d4fb6eb77ca Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 07:14:51 -0500 Subject: [PATCH 07/24] Pyfmmlib: Use kwargs throughout, move towards being PDE-generic --- pytential/qbx/fmmlib.py | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index 6da05574..c19dc3a4 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -334,17 +334,18 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): info = self._info_for_form_global_qbx_locals(src_weights) ier, loc_exp_pre = formta_vec( - zk=self.helmholtz_k, rscale=info.rscale_vec, sources=info.sources, sources_offsets=info.center_source_starts[:-1], - charges=info.charges, - charges_offsets=info.center_source_starts[:-1], + charge=info.charges, + charge_offsets=info.center_source_starts[:-1], nsources=info.center_source_counts[1:], center=info.centers, - nterms=self.nterms) + nterms=self.nterms, + + **self.kernel_kwargs) if np.any(ier != 0): raise RuntimeError("formta returned an error") @@ -428,8 +429,6 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): dtype=self.dtype) expn2 = mploc( - zk=self.helmholtz_k, - rscale1=rscale1, rscale1_offsets=rscale1_offsets, rscale1_starts=src_boxes_starts, @@ -446,7 +445,10 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): # FIXME: center2 has wrong layout, will copy center2=qbx_centers[:, tgt_icenter_vec], expn2=expn2.T, - ier=ier, **kwargs).T + ier=ier, + + **kwargs, + **self.kernel_kwargs).T local_exps[geo_data.global_qbx_centers()] += expn2 @@ -546,9 +548,16 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): 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].T, - rscale, tgt_center, local_order, **kwargs)[..., 0].T + rscale1=rscale, + center1=src_center, + expn1=local_exps[src_ibox].T, + + rscale2=rscale, + center2=tgt_center, + nterms2=local_order, + + **kwargs, + **self.kernel_kwargs)[..., 0].T qbx_expansions[tgt_icenter] += tmp_loc_exp @@ -577,9 +586,12 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): center = qbx_centers[:, src_icenter] - pot, grad = taeval(self.helmholtz_k, rscale, - center, qbx_expansions[src_icenter].T, - all_targets[:, center_itgt]) + pot, grad = taeval( + rscale=rscale, + center=center, + expn=qbx_expansions[src_icenter].T, + ztarg=all_targets[:, center_itgt], + **self.kernel_kwargs) self.add_potgrad_onto_output(output, center_itgt, pot, grad) -- GitLab From c2f7e2339dc517cfe29898a183ecd8a54dcb5e87 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 08:33:00 -0500 Subject: [PATCH 08/24] Fiddle with ellipsoid inteq test case --- test/test_layer_pot.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 7e86a96e..0bc22061 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -310,7 +310,10 @@ def run_int_eq_test(cl_ctx, queue, case, resolution): if case.fmm_backend is None: fmm_order = False else: - fmm_order = case.qbx_order + 5 + if hasattr(case, "fmm_order"): + fmm_order = case.fmm_order + else: + fmm_order = case.qbx_order + 5 qbx = QBXLayerPotentialSource( pre_density_discr, fine_order=source_order, qbx_order=case.qbx_order, @@ -757,7 +760,7 @@ class EllipseIntEqTestCase(CurveIntEqTestCase): class EllipsoidIntEqTestCase(IntEqTestCase): - resolutions = [2, 1] + resolutions = [2, 0.8] name = "ellipsoid" def get_mesh(self, resolution, target_order): @@ -773,13 +776,14 @@ class EllipsoidIntEqTestCase(IntEqTestCase): return perform_flips(mesh, np.ones(mesh.nelements)) fmm_backend = "fmmlib" + fmm_order = 25 use_refinement = False neumann_alpha = 0 # no double layers in FMMlib backend yet inner_radius = 0.4 outer_radius = 5 - qbx_order = 2 + qbx_order = 4 target_order = qbx_order check_tangential_deriv = False -- GitLab From 04b6a0a97848b6013ca5e066701b337cfa7c2247 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 08:38:20 -0500 Subject: [PATCH 09/24] Split layer pot test file into multiple files --- test/test_layer_pot.py | 780 +------------------------------- test/test_layer_pot_identity.py | 268 +++++++++++ test/test_scalar_int_eq.py | 640 ++++++++++++++++++++++++++ 3 files changed, 909 insertions(+), 779 deletions(-) create mode 100644 test/test_layer_pot_identity.py create mode 100644 test/test_scalar_int_eq.py diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index 0bc22061..01db67a7 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -24,11 +24,10 @@ THE SOFTWARE. import numpy as np -import numpy.linalg as la +import numpy.linalg as la # noqa import pyopencl as cl import pyopencl.clmath # noqa import pytest -from pytools import Record from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) @@ -50,15 +49,6 @@ except ImportError: pass -def make_circular_point_group(ambient_dim, npoints, radius, - center=np.array([0., 0.]), func=lambda x: x): - t = func(np.linspace(0, 1, npoints, endpoint=False)) * (2 * np.pi) - center = np.asarray(center) - result = np.zeros((ambient_dim, npoints)) - result[:2, :] = center[:, np.newaxis] + radius*np.vstack((np.cos(t), np.sin(t))) - return result - - # {{{ geometry test def test_geometry(ctx_getter): @@ -286,774 +276,6 @@ def test_ellipse_eigenvalues(ctx_getter, ellipse_aspect, mode_nr, qbx_order): # }}} -# {{{ integral equation test backend - -def run_int_eq_test(cl_ctx, queue, case, resolution): - mesh = case.get_mesh(resolution, case.target_order) - print("%d elements" % mesh.nelements) - - 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(case.target_order)) - - source_order = 4*case.target_order - - refiner_extra_kwargs = {} - - if case.k != 0: - refiner_extra_kwargs["kernel_length_scale"] = 5/case.k - - if case.fmm_backend is None: - fmm_order = False - else: - if hasattr(case, "fmm_order"): - fmm_order = case.fmm_order - else: - fmm_order = case.qbx_order + 5 - - qbx = QBXLayerPotentialSource( - pre_density_discr, fine_order=source_order, qbx_order=case.qbx_order, - fmm_order=fmm_order, fmm_backend=case.fmm_backend) - - if case.use_refinement: - qbx, _ = qbx.with_refinement(**refiner_extra_kwargs) - - density_discr = qbx.density_discr - - # {{{ plot geometry - - if 0: - if mesh.ambient_dim == 2: - # show geometry, centers, normals - nodes_h = density_discr.nodes().get(queue=queue) - pt.plot(nodes_h[0], nodes_h[1], "x-") - normal = bind(density_discr, sym.normal(2))(queue).as_vector(np.object) - pt.quiver(nodes_h[0], nodes_h[1], - normal[0].get(queue), normal[1].get(queue)) - pt.gca().set_aspect("equal") - pt.show() - - elif mesh.ambient_dim == 3: - from meshmode.discretization.visualization import make_visualizer - bdry_vis = make_visualizer(queue, density_discr, case.target_order) - - bdry_normals = bind(density_discr, sym.normal(3))(queue)\ - .as_vector(dtype=object) - - bdry_vis.write_vtk_file("source-%s.vtu" % resolution, [ - ("bdry_normals", bdry_normals), - ]) - - else: - raise ValueError("invalid mesh dim") - - # }}} - - # {{{ set up operator - - from pytential.symbolic.pde.scalar import ( - DirichletOperator, - NeumannOperator) - - from sumpy.kernel import LaplaceKernel, HelmholtzKernel - if case.k: - knl = HelmholtzKernel(mesh.ambient_dim) - knl_kwargs = {"k": sym.var("k")} - concrete_knl_kwargs = {"k": case.k} - else: - knl = LaplaceKernel(mesh.ambient_dim) - knl_kwargs = {} - concrete_knl_kwargs = {} - - if knl.is_complex_valued: - dtype = np.complex128 - else: - dtype = np.float64 - - if case.bc_type == "dirichlet": - op = DirichletOperator(knl, case.loc_sign, use_l2_weighting=True, - kernel_arguments=knl_kwargs) - elif case.bc_type == "neumann": - op = NeumannOperator(knl, case.loc_sign, use_l2_weighting=True, - use_improved_operator=False, kernel_arguments=knl_kwargs, - alpha=case.neumann_alpha) - else: - assert False - - op_u = op.operator(sym.var("u")) - - # }}} - - # {{{ set up test data - - if case.loc_sign < 0: - test_src_geo_radius = case.outer_radius - test_tgt_geo_radius = case.inner_radius - else: - test_src_geo_radius = case.inner_radius - test_tgt_geo_radius = case.outer_radius - - point_sources = make_circular_point_group( - mesh.ambient_dim, 10, test_src_geo_radius, - func=lambda x: x**1.5) - test_targets = make_circular_point_group( - mesh.ambient_dim, 20, test_tgt_geo_radius) - - np.random.seed(22) - source_charges = np.random.randn(point_sources.shape[1]) - source_charges[-1] = -np.sum(source_charges[:-1]) - source_charges = source_charges.astype(dtype) - assert np.sum(source_charges) < 1e-15 - - source_charges_dev = cl.array.to_device(queue, source_charges) - - # }}} - - # {{{ establish BCs - - from pytential.source import PointPotentialSource - from pytential.target import PointsTarget - - point_source = PointPotentialSource(cl_ctx, point_sources) - - pot_src = sym.IntG( - # FIXME: qbx_forced_limit--really? - knl, sym.var("charges"), qbx_forced_limit=None, **knl_kwargs) - - test_direct = bind((point_source, PointsTarget(test_targets)), pot_src)( - queue, charges=source_charges_dev, **concrete_knl_kwargs) - - if case.bc_type == "dirichlet": - bc = bind((point_source, density_discr), pot_src)( - queue, charges=source_charges_dev, **concrete_knl_kwargs) - - elif case.bc_type == "neumann": - bc = bind( - (point_source, density_discr), - sym.normal_derivative( - qbx.ambient_dim, pot_src, where=sym.DEFAULT_TARGET) - )(queue, charges=source_charges_dev, **concrete_knl_kwargs) - - # }}} - - # {{{ solve - - bound_op = bind(qbx, op_u) - - rhs = bind(density_discr, op.prepare_rhs(sym.var("bc")))(queue, bc=bc) - - from pytential.solve import gmres - gmres_result = gmres( - bound_op.scipy_op(queue, "u", dtype, **concrete_knl_kwargs), - rhs, - tol=case.gmres_tol, - progress=True, - hard_failure=True) - - print("gmres state:", gmres_result.state) - u = gmres_result.solution - - # }}} - - # {{{ build matrix for spectrum check - - if 0: - from sumpy.tools import build_matrix - mat = build_matrix(bound_op.scipy_op("u", dtype=dtype, k=case.k)) - w, v = la.eig(mat) - if 0: - pt.imshow(np.log10(1e-20+np.abs(mat))) - pt.colorbar() - pt.show() - - #assert abs(s[-1]) < 1e-13, "h - #assert abs(s[-2]) > 1e-7 - #from pudb import set_trace; set_trace() - - # }}} - - # {{{ error check - - bound_tgt_op = bind((qbx, PointsTarget(test_targets)), - op.representation(sym.var("u"))) - - test_via_bdry = bound_tgt_op(queue, u=u, k=case.k) - - err = test_direct-test_via_bdry - - err = err.get() - test_direct = test_direct.get() - test_via_bdry = test_via_bdry.get() - - # {{{ remove effect of net source charge - - if case.k == 0 and case.bc_type == "neumann" and case.loc_sign == -1: - - # remove constant offset in interior Laplace Neumann error - tgt_ones = np.ones_like(test_direct) - tgt_ones = tgt_ones/la.norm(tgt_ones) - err = err - np.vdot(tgt_ones, err)*tgt_ones - - # }}} - - rel_err_2 = la.norm(err)/la.norm(test_direct) - rel_err_inf = la.norm(err, np.inf)/la.norm(test_direct, np.inf) - - # }}} - - print("rel_err_2: %g rel_err_inf: %g" % (rel_err_2, rel_err_inf)) - - # {{{ test tangential derivative - - if case.check_tangential_deriv: - bound_t_deriv_op = bind(qbx, - op.representation( - sym.var("u"), - map_potentials=lambda pot: sym.tangential_derivative(2, pot), - qbx_forced_limit=case.loc_sign)) - - #print(bound_t_deriv_op.code) - - tang_deriv_from_src = bound_t_deriv_op( - queue, u=u, **concrete_knl_kwargs).as_scalar().get() - - tang_deriv_ref = (bind( - (point_source, density_discr), - sym.tangential_derivative(2, pot_src) - )(queue, charges=source_charges_dev, **concrete_knl_kwargs) - .as_scalar().get()) - - if 0: - pt.plot(tang_deriv_ref.real) - pt.plot(tang_deriv_from_src.real) - pt.show() - - td_err = (tang_deriv_from_src - tang_deriv_ref) - - rel_td_err_inf = la.norm(td_err, np.inf)/la.norm(tang_deriv_ref, np.inf) - - print("rel_td_err_inf: %g" % rel_td_err_inf) - - else: - rel_td_err_inf = None - - # }}} - - # {{{ 3D plotting - - if 0: - from meshmode.discretization.visualization import make_visualizer - bdry_vis = make_visualizer(queue, density_discr, case.target_order) - - bdry_normals = bind(density_discr, sym.normal(3))(queue)\ - .as_vector(dtype=object) - - bdry_vis.write_vtk_file("source-%s.vtu" % resolution, [ - ("u", u), - ("bc", bc), - ("bdry_normals", bdry_normals), - ]) - - 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 - fplot = FieldPlotter( - bbox_center, extent=2*2*bbox_size, npoints=(150, 150, 1)) - - qbx_stick_out = qbx.copy(target_stick_out_factor=0.15) - from pytential.target import PointsTarget - from pytential.qbx import QBXTargetAssociationFailedException - - try: - solved_pot = bind( - (qbx_stick_out, PointsTarget(fplot.points)), - op.representation(sym.var("u")) - )(queue, u=u, k=case.k) - except QBXTargetAssociationFailedException as e: - fplot.write_vtk_file( - "failed-targets.vts", - [ - ("failed_targets", e.failed_target_flags.get(queue)) - ]) - raise - - solved_pot = solved_pot.get() - - true_pot = bind((point_source, PointsTarget(fplot.points)), pot_src)( - queue, charges=source_charges_dev, **concrete_knl_kwargs).get() - - #fplot.show_scalar_in_mayavi(solved_pot.real, max_val=5) - fplot.write_vtk_file( - "potential.vts", - [ - ("solved_pot", solved_pot), - ("true_pot", true_pot), - ("pot_diff", solved_pot-true_pot), - ] - ) - - # }}} - - # {{{ 2D plotting - - if 0: - fplot = FieldPlotter(np.zeros(2), - extent=1.25*2*max(test_src_geo_radius, test_tgt_geo_radius), - npoints=200) - - #pt.plot(u) - #pt.show() - - fld_from_src = bind((point_source, PointsTarget(fplot.points)), - pot_src)(queue, charges=source_charges_dev, **concrete_knl_kwargs) - fld_from_bdry = bind( - (qbx, PointsTarget(fplot.points)), - op.representation(sym.var("u")) - )(queue, u=u, k=case.k) - fld_from_src = fld_from_src.get() - fld_from_bdry = fld_from_bdry.get() - - nodes = density_discr.nodes().get(queue=queue) - - def prep(): - pt.plot(point_sources[0], point_sources[1], "o", - label="Monopole 'Point Charges'") - pt.plot(test_targets[0], test_targets[1], "v", - label="Observation Points") - pt.plot(nodes[0], nodes[1], "k-", - label=r"$\Gamma$") - - from matplotlib.cm import get_cmap - cmap = get_cmap() - cmap._init() - if 0: - cmap._lut[(cmap.N*99)//100:, -1] = 0 # make last percent transparent? - - prep() - if 1: - pt.subplot(131) - pt.title("Field error (loc_sign=%s)" % case.loc_sign) - log_err = np.log10(1e-20+np.abs(fld_from_src-fld_from_bdry)) - log_err = np.minimum(-3, log_err) - fplot.show_scalar_in_matplotlib(log_err, cmap=cmap) - - #from matplotlib.colors import Normalize - #im.set_norm(Normalize(vmin=-6, vmax=1)) - - cb = pt.colorbar(shrink=0.9) - cb.set_label(r"$\log_{10}(\mathdefault{Error})$") - - if 1: - pt.subplot(132) - prep() - pt.title("Source Field") - fplot.show_scalar_in_matplotlib( - fld_from_src.real, max_val=3) - - pt.colorbar(shrink=0.9) - if 1: - pt.subplot(133) - prep() - pt.title("Solved Field") - fplot.show_scalar_in_matplotlib( - fld_from_bdry.real, max_val=3) - - pt.colorbar(shrink=0.9) - - # total field - #fplot.show_scalar_in_matplotlib( - #fld_from_src.real+fld_from_bdry.real, max_val=0.1) - - #pt.colorbar() - - pt.legend(loc="best", prop=dict(size=15)) - from matplotlib.ticker import NullFormatter - pt.gca().xaxis.set_major_formatter(NullFormatter()) - pt.gca().yaxis.set_major_formatter(NullFormatter()) - - pt.gca().set_aspect("equal") - - if 0: - border_factor_top = 0.9 - border_factor = 0.3 - - xl, xh = pt.xlim() - xhsize = 0.5*(xh-xl) - pt.xlim(xl-border_factor*xhsize, xh+border_factor*xhsize) - - yl, yh = pt.ylim() - yhsize = 0.5*(yh-yl) - pt.ylim(yl-border_factor_top*yhsize, yh+border_factor*yhsize) - - #pt.savefig("helmholtz.pdf", dpi=600) - pt.show() - - # }}} - - class Result(Record): - pass - - return Result( - h_max=qbx.h_max, - rel_err_2=rel_err_2, - rel_err_inf=rel_err_inf, - rel_td_err_inf=rel_td_err_inf, - gmres_result=gmres_result) - -# }}} - - -# {{{ integral equation test frontend - -class IntEqTestCase: - def __init__(self, helmholtz_k, bc_type, loc_sign): - self.helmholtz_k = helmholtz_k - self.bc_type = bc_type - self.loc_sign = loc_sign - - @property - def k(self): - return self.helmholtz_k - - def __str__(self): - return ("name: %s, bc_type: %s, loc_sign: %s, " - "helmholtz_k: %s, qbx_order: %d, target_order: %d" - % (self.name, self.bc_type, self.loc_sign, self.helmholtz_k, - self.qbx_order, self.target_order)) - - fmm_backend = "sumpy" - gmres_tol = 1e-14 - - -class CurveIntEqTestCase(IntEqTestCase): - resolutions = [30, 40, 50] - - def get_mesh(self, resolution, target_order): - return make_curve_mesh( - self.curve_func, - np.linspace(0, 1, resolution+1), - target_order) - - fmm_backend = None - use_refinement = True - neumann_alpha = None # default - - inner_radius = 0.1 - outer_radius = 2 - - qbx_order = 5 - target_order = qbx_order - - check_tangential_deriv = True - - -class EllipseIntEqTestCase(CurveIntEqTestCase): - name = "3-to-1 ellipse" - - def curve_func(self, x): - return ellipse(3, x) - - -class EllipsoidIntEqTestCase(IntEqTestCase): - resolutions = [2, 0.8] - name = "ellipsoid" - - def get_mesh(self, resolution, target_order): - from meshmode.mesh.io import generate_gmsh, FileSource - mesh = generate_gmsh( - FileSource("ellipsoid.step"), 2, order=2, - other_options=[ - "-string", - "Mesh.CharacteristicLengthMax = %g;" % resolution]) - - from meshmode.mesh.processing import perform_flips - # Flip elements--gmsh generates inside-out geometry. - return perform_flips(mesh, np.ones(mesh.nelements)) - - fmm_backend = "fmmlib" - fmm_order = 25 - use_refinement = False - neumann_alpha = 0 # no double layers in FMMlib backend yet - - inner_radius = 0.4 - outer_radius = 5 - - qbx_order = 4 - target_order = qbx_order - check_tangential_deriv = False - - # We're only expecting three digits based on FMM settings. Who are we - # kidding? - gmres_tol = 1e-5 - - -@pytest.mark.parametrize("case", [ - EllipseIntEqTestCase(helmholtz_k=helmholtz_k, bc_type=bc_type, - loc_sign=loc_sign) - 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 -def test_integral_equation(ctx_getter, case): - logging.basicConfig(level=logging.INFO) - - cl_ctx = ctx_getter() - queue = cl.CommandQueue(cl_ctx) - - if case.fmm_backend == "fmmlib": - pytest.importorskip("pyfmmlib") - - # prevent cache 'splosion - from sympy.core.cache import clear_cache - clear_cache() - - from pytools.convergence import EOCRecorder - print("qbx_order: %d, %s" % (case.qbx_order, case)) - - eoc_rec_target = EOCRecorder() - eoc_rec_td = EOCRecorder() - - for resolution in case.resolutions: - result = run_int_eq_test(cl_ctx, queue, case, resolution) - - eoc_rec_target.add_data_point(result.h_max, result.rel_err_2) - - if result.rel_td_err_inf is not None: - eoc_rec_td.add_data_point(result.h_max, result.rel_td_err_inf) - - if case.bc_type == "dirichlet": - tgt_order = case.qbx_order - elif case.bc_type == "neumann": - tgt_order = case.qbx_order-1 - else: - assert False - - print("TARGET ERROR:") - print(eoc_rec_target) - assert eoc_rec_target.order_estimate() > tgt_order - 1.3 - - if case.check_tangential_deriv: - print("TANGENTIAL DERIVATIVE ERROR:") - print(eoc_rec_td) - assert eoc_rec_td.order_estimate() > tgt_order - 2.3 - -# }}} - - -# {{{ integral identity tester - -d1 = sym.Derivative() -d2 = sym.Derivative() - - -def get_starfish_mesh(refinement_increment, target_order): - nelements = [30, 50, 70][refinement_increment] - return make_curve_mesh(starfish, - np.linspace(0, 1, nelements+1), - 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) - from meshmode.mesh.refinement import Refiner - - refiner = Refiner(mesh) - for i in range(refinement_increment): - flags = np.ones(mesh.nelements, dtype=bool) - refiner.refine(flags) - mesh = refiner.get_current_mesh() - - return mesh - - -@pytest.mark.parametrize(("mesh_name", "mesh_getter", "qbx_order"), [ - #("circle", partial(ellipse, 1)), - #("3-to-1 ellipse", partial(ellipse, 3)), - ("starfish", get_starfish_mesh, 5), - ("sphere", get_sphere_mesh, 3), - ]) -@pytest.mark.parametrize(("zero_op_name", "k"), [ - ("green", 0), - ("green", 1.2), - ("green_grad", 0), - ("green_grad", 1.2), - ("zero_calderon", 0), - ]) -# sample invocation to copy and paste: -# 'test_identities(cl._csc, "green", "starfish", get_starfish_mesh, 4, 0)' -def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, k): - 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() - - if mesh_name == "sphere" and k != 0: - pytest.skip("both direct eval and generating the FMM kernels are too slow") - - if mesh_name == "sphere" and zero_op_name == "green_grad": - pytest.skip("does not achieve sufficient precision") - - target_order = 8 - - order_table = { - "green": qbx_order, - "green_grad": qbx_order-1, - "zero_calderon": qbx_order-1, - } - - from pytools.convergence import EOCRecorder - eoc_rec = EOCRecorder() - - for refinement_increment in [0, 1, 2]: - mesh = mesh_getter(refinement_increment, target_order) - if mesh is None: - break - - d = mesh.ambient_dim - - u_sym = sym.var("u") - grad_u_sym = sym.make_sym_mv("grad_u", d) - dn_u_sym = sym.var("dn_u") - - from sumpy.kernel import LaplaceKernel, HelmholtzKernel - lap_k_sym = LaplaceKernel(d) - if k == 0: - k_sym = lap_k_sym - knl_kwargs = {} - else: - k_sym = HelmholtzKernel(d) - knl_kwargs = {"k": sym.var("k")} - - zero_op_table = { - "green": - sym.S(k_sym, dn_u_sym, qbx_forced_limit=-1, **knl_kwargs) - - sym.D(k_sym, u_sym, qbx_forced_limit="avg", **knl_kwargs) - - 0.5*u_sym, - - "green_grad": - d1.resolve(d1.dnabla(d) * d1(sym.S(k_sym, dn_u_sym, - qbx_forced_limit="avg", **knl_kwargs))) - - d2.resolve(d2.dnabla(d) * d2(sym.D(k_sym, u_sym, - qbx_forced_limit="avg", **knl_kwargs))) - - 0.5*grad_u_sym, - - # only for k==0: - "zero_calderon": - -sym.Dp(lap_k_sym, sym.S(lap_k_sym, u_sym)) - - 0.25*u_sym + sym.Sp(lap_k_sym, sym.Sp(lap_k_sym, u_sym)) - } - - zero_op = zero_op_table[zero_op_name] - - from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import \ - InterpolatoryQuadratureSimplexGroupFactory - from pytential.qbx import QBXLayerPotentialSource - pre_density_discr = Discretization( - cl_ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(target_order)) - - if d == 2: - order_bump = 15 - elif d == 3: - order_bump = 8 - - 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=qbx_order + order_bump - ).with_refinement(**refiner_extra_kwargs) - - density_discr = qbx.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)] - - if k != 0: - if d == 2: - angle = 0.3 - wave_vec = np.array([np.cos(angle), np.sin(angle)]) - u = np.exp(1j*k*np.tensordot(wave_vec, nodes_host, axes=1)) - grad_u = 1j*k*wave_vec[:, np.newaxis]*u - else: - center = np.array([3, 1, 2]) - diff = nodes_host - center[:, np.newaxis] - r = la.norm(diff, axis=0) - u = np.exp(1j*k*r) / r - grad_u = diff * (1j*k*u/r - u/r**2) - else: - 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) - grad_u_dev = cl.array.to_device(queue, grad_u) - - key = (qbx_order, mesh_name, refinement_increment, zero_op_name) - - bound_op = bind(qbx, zero_op) - error = bound_op( - queue, u=u_dev, dn_u=dn_u_dev, grad_u=grad_u_dev, k=k) - if 0: - pt.plot(error) - pt.show() - - l2_error_norm = norm(density_discr, queue, error) - print(key, l2_error_norm) - - eoc_rec.add_data_point(qbx.h_max, l2_error_norm) - - print(eoc_rec) - tgt_order = order_table[zero_op_name] - assert eoc_rec.order_estimate() > tgt_order - 1.3 - -# }}} - - # {{{ test off-surface eval @pytest.mark.parametrize("use_fmm", [True, False]) diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py new file mode 100644 index 00000000..96f46ec3 --- /dev/null +++ b/test/test_layer_pot_identity.py @@ -0,0 +1,268 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2013-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 +import numpy.linalg as la +import pyopencl as cl +import pyopencl.clmath # noqa +import pytest +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +from functools import partial +from meshmode.mesh.generation import ( # noqa + ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, WobblyCircle, + make_curve_mesh) +# from sumpy.visualization import FieldPlotter +from pytential import bind, sym, norm + +import logging +logger = logging.getLogger(__name__) + +circle = partial(ellipse, 1) + +try: + import matplotlib.pyplot as pt +except ImportError: + pass + + +# {{{ integral identity tester + +d1 = sym.Derivative() +d2 = sym.Derivative() + + +def get_starfish_mesh(refinement_increment, target_order): + nelements = [30, 50, 70][refinement_increment] + return make_curve_mesh(starfish, + np.linspace(0, 1, nelements+1), + 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) + from meshmode.mesh.refinement import Refiner + + refiner = Refiner(mesh) + for i in range(refinement_increment): + flags = np.ones(mesh.nelements, dtype=bool) + refiner.refine(flags) + mesh = refiner.get_current_mesh() + + return mesh + + +@pytest.mark.parametrize(("mesh_name", "mesh_getter", "qbx_order"), [ + #("circle", partial(ellipse, 1)), + #("3-to-1 ellipse", partial(ellipse, 3)), + ("starfish", get_starfish_mesh, 5), + ("sphere", get_sphere_mesh, 3), + ]) +@pytest.mark.parametrize(("zero_op_name", "k"), [ + ("green", 0), + ("green", 1.2), + ("green_grad", 0), + ("green_grad", 1.2), + ("zero_calderon", 0), + ]) +# sample invocation to copy and paste: +# 'test_identities(cl._csc, "green", "starfish", get_starfish_mesh, 4, 0)' +def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, k): + 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() + + if mesh_name == "sphere" and k != 0: + pytest.skip("both direct eval and generating the FMM kernels are too slow") + + if mesh_name == "sphere" and zero_op_name == "green_grad": + pytest.skip("does not achieve sufficient precision") + + target_order = 8 + + order_table = { + "green": qbx_order, + "green_grad": qbx_order-1, + "zero_calderon": qbx_order-1, + } + + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + for refinement_increment in [0, 1, 2]: + mesh = mesh_getter(refinement_increment, target_order) + if mesh is None: + break + + d = mesh.ambient_dim + + u_sym = sym.var("u") + grad_u_sym = sym.make_sym_mv("grad_u", d) + dn_u_sym = sym.var("dn_u") + + from sumpy.kernel import LaplaceKernel, HelmholtzKernel + lap_k_sym = LaplaceKernel(d) + if k == 0: + k_sym = lap_k_sym + knl_kwargs = {} + else: + k_sym = HelmholtzKernel(d) + knl_kwargs = {"k": sym.var("k")} + + zero_op_table = { + "green": + sym.S(k_sym, dn_u_sym, qbx_forced_limit=-1, **knl_kwargs) + - sym.D(k_sym, u_sym, qbx_forced_limit="avg", **knl_kwargs) + - 0.5*u_sym, + + "green_grad": + d1.resolve(d1.dnabla(d) * d1(sym.S(k_sym, dn_u_sym, + qbx_forced_limit="avg", **knl_kwargs))) + - d2.resolve(d2.dnabla(d) * d2(sym.D(k_sym, u_sym, + qbx_forced_limit="avg", **knl_kwargs))) + - 0.5*grad_u_sym, + + # only for k==0: + "zero_calderon": + -sym.Dp(lap_k_sym, sym.S(lap_k_sym, u_sym)) + - 0.25*u_sym + sym.Sp(lap_k_sym, sym.Sp(lap_k_sym, u_sym)) + } + + zero_op = zero_op_table[zero_op_name] + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + from pytential.qbx import QBXLayerPotentialSource + pre_density_discr = Discretization( + cl_ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + if d == 2: + order_bump = 15 + elif d == 3: + order_bump = 8 + + 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=qbx_order + order_bump + ).with_refinement(**refiner_extra_kwargs) + + density_discr = qbx.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)] + + if k != 0: + if d == 2: + angle = 0.3 + wave_vec = np.array([np.cos(angle), np.sin(angle)]) + u = np.exp(1j*k*np.tensordot(wave_vec, nodes_host, axes=1)) + grad_u = 1j*k*wave_vec[:, np.newaxis]*u + else: + center = np.array([3, 1, 2]) + diff = nodes_host - center[:, np.newaxis] + r = la.norm(diff, axis=0) + u = np.exp(1j*k*r) / r + grad_u = diff * (1j*k*u/r - u/r**2) + else: + 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) + grad_u_dev = cl.array.to_device(queue, grad_u) + + key = (qbx_order, mesh_name, refinement_increment, zero_op_name) + + bound_op = bind(qbx, zero_op) + error = bound_op( + queue, u=u_dev, dn_u=dn_u_dev, grad_u=grad_u_dev, k=k) + if 0: + pt.plot(error) + pt.show() + + l2_error_norm = norm(density_discr, queue, error) + print(key, l2_error_norm) + + eoc_rec.add_data_point(qbx.h_max, l2_error_norm) + + print(eoc_rec) + tgt_order = order_table[zero_op_name] + assert eoc_rec.order_estimate() > tgt_order - 1.3 + +# }}} + + +# You can test individual routines by typing +# $ python test_layer_pot_identity.py 'test_routine()' + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from py.test.cmdline import main + main([__file__]) + +# vim: fdm=marker diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py new file mode 100644 index 00000000..10e20519 --- /dev/null +++ b/test/test_scalar_int_eq.py @@ -0,0 +1,640 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2014 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 +import numpy.linalg as la +import pyopencl as cl +import pyopencl.clmath # noqa +import pytest +from pytools import Record +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) + +from functools import partial +from meshmode.mesh.generation import ( # noqa + ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, WobblyCircle, + make_curve_mesh) +from sumpy.visualization import FieldPlotter +from pytential import bind, sym + +import logging +logger = logging.getLogger(__name__) + +circle = partial(ellipse, 1) + +try: + import matplotlib.pyplot as pt +except ImportError: + pass + + +def make_circular_point_group(ambient_dim, npoints, radius, + center=np.array([0., 0.]), func=lambda x: x): + t = func(np.linspace(0, 1, npoints, endpoint=False)) * (2 * np.pi) + center = np.asarray(center) + result = np.zeros((ambient_dim, npoints)) + result[:2, :] = center[:, np.newaxis] + radius*np.vstack((np.cos(t), np.sin(t))) + return result + + +# {{{ test backend + +def run_int_eq_test(cl_ctx, queue, case, resolution): + mesh = case.get_mesh(resolution, case.target_order) + print("%d elements" % mesh.nelements) + + 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(case.target_order)) + + source_order = 4*case.target_order + + refiner_extra_kwargs = {} + + if case.k != 0: + refiner_extra_kwargs["kernel_length_scale"] = 5/case.k + + if case.fmm_backend is None: + fmm_order = False + else: + if hasattr(case, "fmm_order"): + fmm_order = case.fmm_order + else: + fmm_order = case.qbx_order + 5 + + qbx = QBXLayerPotentialSource( + pre_density_discr, fine_order=source_order, qbx_order=case.qbx_order, + fmm_order=fmm_order, fmm_backend=case.fmm_backend) + + if case.use_refinement: + qbx, _ = qbx.with_refinement(**refiner_extra_kwargs) + + density_discr = qbx.density_discr + + # {{{ plot geometry + + if 0: + if mesh.ambient_dim == 2: + # show geometry, centers, normals + nodes_h = density_discr.nodes().get(queue=queue) + pt.plot(nodes_h[0], nodes_h[1], "x-") + normal = bind(density_discr, sym.normal(2))(queue).as_vector(np.object) + pt.quiver(nodes_h[0], nodes_h[1], + normal[0].get(queue), normal[1].get(queue)) + pt.gca().set_aspect("equal") + pt.show() + + elif mesh.ambient_dim == 3: + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, density_discr, case.target_order) + + bdry_normals = bind(density_discr, sym.normal(3))(queue)\ + .as_vector(dtype=object) + + bdry_vis.write_vtk_file("source-%s.vtu" % resolution, [ + ("bdry_normals", bdry_normals), + ]) + + else: + raise ValueError("invalid mesh dim") + + # }}} + + # {{{ set up operator + + from pytential.symbolic.pde.scalar import ( + DirichletOperator, + NeumannOperator) + + from sumpy.kernel import LaplaceKernel, HelmholtzKernel + if case.k: + knl = HelmholtzKernel(mesh.ambient_dim) + knl_kwargs = {"k": sym.var("k")} + concrete_knl_kwargs = {"k": case.k} + else: + knl = LaplaceKernel(mesh.ambient_dim) + knl_kwargs = {} + concrete_knl_kwargs = {} + + if knl.is_complex_valued: + dtype = np.complex128 + else: + dtype = np.float64 + + if case.bc_type == "dirichlet": + op = DirichletOperator(knl, case.loc_sign, use_l2_weighting=True, + kernel_arguments=knl_kwargs) + elif case.bc_type == "neumann": + op = NeumannOperator(knl, case.loc_sign, use_l2_weighting=True, + use_improved_operator=False, kernel_arguments=knl_kwargs, + alpha=case.neumann_alpha) + else: + assert False + + op_u = op.operator(sym.var("u")) + + # }}} + + # {{{ set up test data + + if case.loc_sign < 0: + test_src_geo_radius = case.outer_radius + test_tgt_geo_radius = case.inner_radius + else: + test_src_geo_radius = case.inner_radius + test_tgt_geo_radius = case.outer_radius + + point_sources = make_circular_point_group( + mesh.ambient_dim, 10, test_src_geo_radius, + func=lambda x: x**1.5) + test_targets = make_circular_point_group( + mesh.ambient_dim, 20, test_tgt_geo_radius) + + np.random.seed(22) + source_charges = np.random.randn(point_sources.shape[1]) + source_charges[-1] = -np.sum(source_charges[:-1]) + source_charges = source_charges.astype(dtype) + assert np.sum(source_charges) < 1e-15 + + source_charges_dev = cl.array.to_device(queue, source_charges) + + # }}} + + # {{{ establish BCs + + from pytential.source import PointPotentialSource + from pytential.target import PointsTarget + + point_source = PointPotentialSource(cl_ctx, point_sources) + + pot_src = sym.IntG( + # FIXME: qbx_forced_limit--really? + knl, sym.var("charges"), qbx_forced_limit=None, **knl_kwargs) + + test_direct = bind((point_source, PointsTarget(test_targets)), pot_src)( + queue, charges=source_charges_dev, **concrete_knl_kwargs) + + if case.bc_type == "dirichlet": + bc = bind((point_source, density_discr), pot_src)( + queue, charges=source_charges_dev, **concrete_knl_kwargs) + + elif case.bc_type == "neumann": + bc = bind( + (point_source, density_discr), + sym.normal_derivative( + qbx.ambient_dim, pot_src, where=sym.DEFAULT_TARGET) + )(queue, charges=source_charges_dev, **concrete_knl_kwargs) + + # }}} + + # {{{ solve + + bound_op = bind(qbx, op_u) + + rhs = bind(density_discr, op.prepare_rhs(sym.var("bc")))(queue, bc=bc) + + from pytential.solve import gmres + gmres_result = gmres( + bound_op.scipy_op(queue, "u", dtype, **concrete_knl_kwargs), + rhs, + tol=case.gmres_tol, + progress=True, + hard_failure=True) + + print("gmres state:", gmres_result.state) + u = gmres_result.solution + + # }}} + + # {{{ build matrix for spectrum check + + if 0: + from sumpy.tools import build_matrix + mat = build_matrix(bound_op.scipy_op("u", dtype=dtype, k=case.k)) + w, v = la.eig(mat) + if 0: + pt.imshow(np.log10(1e-20+np.abs(mat))) + pt.colorbar() + pt.show() + + #assert abs(s[-1]) < 1e-13, "h + #assert abs(s[-2]) > 1e-7 + #from pudb import set_trace; set_trace() + + # }}} + + # {{{ error check + + bound_tgt_op = bind((qbx, PointsTarget(test_targets)), + op.representation(sym.var("u"))) + + test_via_bdry = bound_tgt_op(queue, u=u, k=case.k) + + err = test_direct-test_via_bdry + + err = err.get() + test_direct = test_direct.get() + test_via_bdry = test_via_bdry.get() + + # {{{ remove effect of net source charge + + if case.k == 0 and case.bc_type == "neumann" and case.loc_sign == -1: + + # remove constant offset in interior Laplace Neumann error + tgt_ones = np.ones_like(test_direct) + tgt_ones = tgt_ones/la.norm(tgt_ones) + err = err - np.vdot(tgt_ones, err)*tgt_ones + + # }}} + + rel_err_2 = la.norm(err)/la.norm(test_direct) + rel_err_inf = la.norm(err, np.inf)/la.norm(test_direct, np.inf) + + # }}} + + print("rel_err_2: %g rel_err_inf: %g" % (rel_err_2, rel_err_inf)) + + # {{{ test tangential derivative + + if case.check_tangential_deriv: + bound_t_deriv_op = bind(qbx, + op.representation( + sym.var("u"), + map_potentials=lambda pot: sym.tangential_derivative(2, pot), + qbx_forced_limit=case.loc_sign)) + + #print(bound_t_deriv_op.code) + + tang_deriv_from_src = bound_t_deriv_op( + queue, u=u, **concrete_knl_kwargs).as_scalar().get() + + tang_deriv_ref = (bind( + (point_source, density_discr), + sym.tangential_derivative(2, pot_src) + )(queue, charges=source_charges_dev, **concrete_knl_kwargs) + .as_scalar().get()) + + if 0: + pt.plot(tang_deriv_ref.real) + pt.plot(tang_deriv_from_src.real) + pt.show() + + td_err = (tang_deriv_from_src - tang_deriv_ref) + + rel_td_err_inf = la.norm(td_err, np.inf)/la.norm(tang_deriv_ref, np.inf) + + print("rel_td_err_inf: %g" % rel_td_err_inf) + + else: + rel_td_err_inf = None + + # }}} + + # {{{ 3D plotting + + if 0: + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, density_discr, case.target_order) + + bdry_normals = bind(density_discr, sym.normal(3))(queue)\ + .as_vector(dtype=object) + + bdry_vis.write_vtk_file("source-%s.vtu" % resolution, [ + ("u", u), + ("bc", bc), + ("bdry_normals", bdry_normals), + ]) + + 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 + fplot = FieldPlotter( + bbox_center, extent=2*2*bbox_size, npoints=(150, 150, 1)) + + qbx_stick_out = qbx.copy(target_stick_out_factor=0.15) + from pytential.target import PointsTarget + from pytential.qbx import QBXTargetAssociationFailedException + + try: + solved_pot = bind( + (qbx_stick_out, PointsTarget(fplot.points)), + op.representation(sym.var("u")) + )(queue, u=u, k=case.k) + except QBXTargetAssociationFailedException as e: + fplot.write_vtk_file( + "failed-targets.vts", + [ + ("failed_targets", e.failed_target_flags.get(queue)) + ]) + raise + + solved_pot = solved_pot.get() + + true_pot = bind((point_source, PointsTarget(fplot.points)), pot_src)( + queue, charges=source_charges_dev, **concrete_knl_kwargs).get() + + #fplot.show_scalar_in_mayavi(solved_pot.real, max_val=5) + fplot.write_vtk_file( + "potential.vts", + [ + ("solved_pot", solved_pot), + ("true_pot", true_pot), + ("pot_diff", solved_pot-true_pot), + ] + ) + + # }}} + + # {{{ 2D plotting + + if 0: + fplot = FieldPlotter(np.zeros(2), + extent=1.25*2*max(test_src_geo_radius, test_tgt_geo_radius), + npoints=200) + + #pt.plot(u) + #pt.show() + + fld_from_src = bind((point_source, PointsTarget(fplot.points)), + pot_src)(queue, charges=source_charges_dev, **concrete_knl_kwargs) + fld_from_bdry = bind( + (qbx, PointsTarget(fplot.points)), + op.representation(sym.var("u")) + )(queue, u=u, k=case.k) + fld_from_src = fld_from_src.get() + fld_from_bdry = fld_from_bdry.get() + + nodes = density_discr.nodes().get(queue=queue) + + def prep(): + pt.plot(point_sources[0], point_sources[1], "o", + label="Monopole 'Point Charges'") + pt.plot(test_targets[0], test_targets[1], "v", + label="Observation Points") + pt.plot(nodes[0], nodes[1], "k-", + label=r"$\Gamma$") + + from matplotlib.cm import get_cmap + cmap = get_cmap() + cmap._init() + if 0: + cmap._lut[(cmap.N*99)//100:, -1] = 0 # make last percent transparent? + + prep() + if 1: + pt.subplot(131) + pt.title("Field error (loc_sign=%s)" % case.loc_sign) + log_err = np.log10(1e-20+np.abs(fld_from_src-fld_from_bdry)) + log_err = np.minimum(-3, log_err) + fplot.show_scalar_in_matplotlib(log_err, cmap=cmap) + + #from matplotlib.colors import Normalize + #im.set_norm(Normalize(vmin=-6, vmax=1)) + + cb = pt.colorbar(shrink=0.9) + cb.set_label(r"$\log_{10}(\mathdefault{Error})$") + + if 1: + pt.subplot(132) + prep() + pt.title("Source Field") + fplot.show_scalar_in_matplotlib( + fld_from_src.real, max_val=3) + + pt.colorbar(shrink=0.9) + if 1: + pt.subplot(133) + prep() + pt.title("Solved Field") + fplot.show_scalar_in_matplotlib( + fld_from_bdry.real, max_val=3) + + pt.colorbar(shrink=0.9) + + # total field + #fplot.show_scalar_in_matplotlib( + #fld_from_src.real+fld_from_bdry.real, max_val=0.1) + + #pt.colorbar() + + pt.legend(loc="best", prop=dict(size=15)) + from matplotlib.ticker import NullFormatter + pt.gca().xaxis.set_major_formatter(NullFormatter()) + pt.gca().yaxis.set_major_formatter(NullFormatter()) + + pt.gca().set_aspect("equal") + + if 0: + border_factor_top = 0.9 + border_factor = 0.3 + + xl, xh = pt.xlim() + xhsize = 0.5*(xh-xl) + pt.xlim(xl-border_factor*xhsize, xh+border_factor*xhsize) + + yl, yh = pt.ylim() + yhsize = 0.5*(yh-yl) + pt.ylim(yl-border_factor_top*yhsize, yh+border_factor*yhsize) + + #pt.savefig("helmholtz.pdf", dpi=600) + pt.show() + + # }}} + + class Result(Record): + pass + + return Result( + h_max=qbx.h_max, + rel_err_2=rel_err_2, + rel_err_inf=rel_err_inf, + rel_td_err_inf=rel_td_err_inf, + gmres_result=gmres_result) + +# }}} + + +# {{{ test cases + +class IntEqTestCase: + def __init__(self, helmholtz_k, bc_type, loc_sign): + self.helmholtz_k = helmholtz_k + self.bc_type = bc_type + self.loc_sign = loc_sign + + @property + def k(self): + return self.helmholtz_k + + def __str__(self): + return ("name: %s, bc_type: %s, loc_sign: %s, " + "helmholtz_k: %s, qbx_order: %d, target_order: %d" + % (self.name, self.bc_type, self.loc_sign, self.helmholtz_k, + self.qbx_order, self.target_order)) + + fmm_backend = "sumpy" + gmres_tol = 1e-14 + + +class CurveIntEqTestCase(IntEqTestCase): + resolutions = [30, 40, 50] + + def get_mesh(self, resolution, target_order): + return make_curve_mesh( + self.curve_func, + np.linspace(0, 1, resolution+1), + target_order) + + fmm_backend = None + use_refinement = True + neumann_alpha = None # default + + inner_radius = 0.1 + outer_radius = 2 + + qbx_order = 5 + target_order = qbx_order + + check_tangential_deriv = True + + +class EllipseIntEqTestCase(CurveIntEqTestCase): + name = "3-to-1 ellipse" + + def curve_func(self, x): + return ellipse(3, x) + + +class EllipsoidIntEqTestCase(IntEqTestCase): + resolutions = [2, 0.8] + name = "ellipsoid" + + def get_mesh(self, resolution, target_order): + from meshmode.mesh.io import generate_gmsh, FileSource + mesh = generate_gmsh( + FileSource("ellipsoid.step"), 2, order=2, + other_options=[ + "-string", + "Mesh.CharacteristicLengthMax = %g;" % resolution]) + + from meshmode.mesh.processing import perform_flips + # Flip elements--gmsh generates inside-out geometry. + return perform_flips(mesh, np.ones(mesh.nelements)) + + fmm_backend = "fmmlib" + fmm_order = 25 + use_refinement = False + neumann_alpha = 0 # no double layers in FMMlib backend yet + + inner_radius = 0.4 + outer_radius = 5 + + qbx_order = 4 + target_order = qbx_order + check_tangential_deriv = False + + # We're only expecting three digits based on FMM settings. Who are we + # kidding? + gmres_tol = 1e-5 + +# }}} + + +# {{{ test frontend + +@pytest.mark.parametrize("case", [ + EllipseIntEqTestCase(helmholtz_k=helmholtz_k, bc_type=bc_type, + loc_sign=loc_sign) + 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 +def test_integral_equation(ctx_getter, case): + logging.basicConfig(level=logging.INFO) + + cl_ctx = ctx_getter() + queue = cl.CommandQueue(cl_ctx) + + if case.fmm_backend == "fmmlib": + pytest.importorskip("pyfmmlib") + + # prevent cache 'splosion + from sympy.core.cache import clear_cache + clear_cache() + + from pytools.convergence import EOCRecorder + print("qbx_order: %d, %s" % (case.qbx_order, case)) + + eoc_rec_target = EOCRecorder() + eoc_rec_td = EOCRecorder() + + for resolution in case.resolutions: + result = run_int_eq_test(cl_ctx, queue, case, resolution) + + eoc_rec_target.add_data_point(result.h_max, result.rel_err_2) + + if result.rel_td_err_inf is not None: + eoc_rec_td.add_data_point(result.h_max, result.rel_td_err_inf) + + if case.bc_type == "dirichlet": + tgt_order = case.qbx_order + elif case.bc_type == "neumann": + tgt_order = case.qbx_order-1 + else: + assert False + + print("TARGET ERROR:") + print(eoc_rec_target) + assert eoc_rec_target.order_estimate() > tgt_order - 1.3 + + if case.check_tangential_deriv: + print("TANGENTIAL DERIVATIVE ERROR:") + print(eoc_rec_td) + assert eoc_rec_td.order_estimate() > tgt_order - 2.3 + +# }}} + + +# You can test individual routines by typing +# $ python test_scalar_int_eq.py 'test_routine()' + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from py.test.cmdline import main + main([__file__]) + +# vim: fdm=marker -- GitLab From 7683a23569bac110af0ee2ec59887239540d5199 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 12:10:25 -0500 Subject: [PATCH 10/24] Refactor integral identity test case for object-based test cases --- test/test_layer_pot_identity.py | 228 +++++++++++++++++++++----------- 1 file changed, 150 insertions(+), 78 deletions(-) diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index 96f46ec3..8486cd70 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -37,6 +37,7 @@ from meshmode.mesh.generation import ( # noqa make_curve_mesh) # from sumpy.visualization import FieldPlotter from pytential import bind, sym, norm +from sumpy.kernel import LaplaceKernel, HelmholtzKernel import logging logger = logging.getLogger(__name__) @@ -49,19 +50,10 @@ except ImportError: pass -# {{{ integral identity tester - d1 = sym.Derivative() d2 = sym.Derivative() -def get_starfish_mesh(refinement_increment, target_order): - nelements = [30, 50, 70][refinement_increment] - return make_curve_mesh(starfish, - np.linspace(0, 1, nelements+1), - 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), @@ -83,24 +75,142 @@ def get_sphere_mesh(refinement_increment, target_order): return mesh -@pytest.mark.parametrize(("mesh_name", "mesh_getter", "qbx_order"), [ - #("circle", partial(ellipse, 1)), - #("3-to-1 ellipse", partial(ellipse, 3)), - ("starfish", get_starfish_mesh, 5), - ("sphere", get_sphere_mesh, 3), - ]) -@pytest.mark.parametrize(("zero_op_name", "k"), [ - ("green", 0), - ("green", 1.2), - ("green_grad", 0), - ("green_grad", 1.2), - ("zero_calderon", 0), - ]) -# sample invocation to copy and paste: -# 'test_identities(cl._csc, "green", "starfish", get_starfish_mesh, 4, 0)' -def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, k): +class StarfishGeometry(object): + mesh_name = "starfish" + dim = 2 + + resolutions = [30, 50, 70] + + def get_mesh(self, nelements, target_order): + return make_curve_mesh(starfish, + np.linspace(0, 1, nelements+1), + target_order) + + +class SphereGeometry(object): + mesh_name = "sphere" + dim = 3 + + resolutions = [0, 1, 2] + + def get_mesh(self, resolution, tgt_order): + return get_sphere_mesh(resolution, tgt_order) + + +class GreenExpr(object): + zero_op_name = "green" + + def get_zero_op(self, kernel, **knl_kwargs): + + u_sym = sym.var("u") + dn_u_sym = sym.var("dn_u") + + return ( + sym.S(kernel, dn_u_sym, qbx_forced_limit=-1, **knl_kwargs) + - sym.D(kernel, u_sym, qbx_forced_limit="avg", **knl_kwargs) + - 0.5*u_sym) + + order_drop = 0 + + +class GradGreenExpr(object): + zero_op_name = "grad_green" + + def get_zero_op(self, kernel, **knl_kwargs): + d = kernel.dim + u_sym = sym.var("u") + grad_u_sym = sym.make_sym_mv("grad_u", d) + dn_u_sym = sym.var("dn_u") + + return ( + d1.resolve(d1.dnabla(d) * d1(sym.S(kernel, dn_u_sym, + qbx_forced_limit="avg", **knl_kwargs))) + - d2.resolve(d2.dnabla(d) * d2(sym.D(kernel, u_sym, + qbx_forced_limit="avg", **knl_kwargs))) + - 0.5*grad_u_sym + ) + + order_drop = 1 + + +class ZeroCalderonExpr(object): + zero_op_name = "calderon" + + def get_zero_op(self, kernel, **knl_kwargs): + assert isinstance(kernel, LaplaceKernel) + assert not knl_kwargs + + u_sym = sym.var("u") + + return ( + -sym.Dp(kernel, sym.S(kernel, u_sym)) + - 0.25*u_sym + sym.Sp(kernel, sym.Sp(kernel, u_sym)) + ) + + order_drop = 1 + + +class StaticTestCase(object): + def check(self): + pass + + +class SphereGreenTest(StaticTestCase): + expr = GreenExpr() + geometry = SphereGeometry() + k = 1.2 + qbx_order = 3 + fmm_order = 15 + + fmm_backend = "fmmlib" + + +class DynamicTestCase(object): + fmm_backend = "sumpy" + + def __init__(self, geometry, expr, k): + self.geometry = geometry + self.expr = expr + self.k = k + self.qbx_order = 5 if geometry.dim == 2 else 3 + + if geometry.dim == 2: + order_bump = 15 + elif geometry.dim == 3: + order_bump = 8 + + self.fmm_order = self.qbx_order + order_bump + + def check(self): + if self.geometry.mesh_name == "sphere" and self.k != 0: + pytest.skip("both direct eval and generating the FMM kernels " + "are too slow") + + if (self.geometry.mesh_name == "sphere" + and self.expr.zero_op_name == "green_grad"): + pytest.skip("does not achieve sufficient precision") + + +# {{{ integral identity tester + +@pytest.mark.parametrize("case", [ + tc + for geom in [ + StarfishGeometry(), + SphereGeometry(), + ] + for tc in [ + DynamicTestCase(geom, GreenExpr(), 0), + DynamicTestCase(geom, GreenExpr(), 1.2), + DynamicTestCase(geom, GradGreenExpr(), 0), + DynamicTestCase(geom, GradGreenExpr(), 1.2), + DynamicTestCase(geom, ZeroCalderonExpr(), 0), + ]]) +def test_identity_convergence(ctx_getter, case): logging.basicConfig(level=logging.INFO) + case.check() + cl_ctx = ctx_getter() queue = cl.CommandQueue(cl_ctx) @@ -108,35 +218,19 @@ def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, from sympy.core.cache import clear_cache clear_cache() - if mesh_name == "sphere" and k != 0: - pytest.skip("both direct eval and generating the FMM kernels are too slow") - - if mesh_name == "sphere" and zero_op_name == "green_grad": - pytest.skip("does not achieve sufficient precision") - target_order = 8 - order_table = { - "green": qbx_order, - "green_grad": qbx_order-1, - "zero_calderon": qbx_order-1, - } - from pytools.convergence import EOCRecorder eoc_rec = EOCRecorder() - for refinement_increment in [0, 1, 2]: - mesh = mesh_getter(refinement_increment, target_order) + for resolution in case.geometry.resolutions: + mesh = case.geometry.get_mesh(resolution, target_order) if mesh is None: break d = mesh.ambient_dim + k = case.k - u_sym = sym.var("u") - grad_u_sym = sym.make_sym_mv("grad_u", d) - dn_u_sym = sym.var("dn_u") - - from sumpy.kernel import LaplaceKernel, HelmholtzKernel lap_k_sym = LaplaceKernel(d) if k == 0: k_sym = lap_k_sym @@ -145,27 +239,6 @@ def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, k_sym = HelmholtzKernel(d) knl_kwargs = {"k": sym.var("k")} - zero_op_table = { - "green": - sym.S(k_sym, dn_u_sym, qbx_forced_limit=-1, **knl_kwargs) - - sym.D(k_sym, u_sym, qbx_forced_limit="avg", **knl_kwargs) - - 0.5*u_sym, - - "green_grad": - d1.resolve(d1.dnabla(d) * d1(sym.S(k_sym, dn_u_sym, - qbx_forced_limit="avg", **knl_kwargs))) - - d2.resolve(d2.dnabla(d) * d2(sym.D(k_sym, u_sym, - qbx_forced_limit="avg", **knl_kwargs))) - - 0.5*grad_u_sym, - - # only for k==0: - "zero_calderon": - -sym.Dp(lap_k_sym, sym.S(lap_k_sym, u_sym)) - - 0.25*u_sym + sym.Sp(lap_k_sym, sym.Sp(lap_k_sym, u_sym)) - } - - zero_op = zero_op_table[zero_op_name] - from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ InterpolatoryQuadratureSimplexGroupFactory @@ -174,19 +247,15 @@ def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, cl_ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(target_order)) - if d == 2: - order_bump = 15 - elif d == 3: - order_bump = 8 - refiner_extra_kwargs = {} - if k != 0: - refiner_extra_kwargs["kernel_length_scale"] = 5/k + if case.k != 0: + refiner_extra_kwargs["kernel_length_scale"] = 5/case.k qbx, _ = QBXLayerPotentialSource( pre_density_discr, 4*target_order, - qbx_order, fmm_order=qbx_order + order_bump + case.qbx_order, fmm_order=case.fmm_order, + fmm_backend=case.fmm_backend, ).with_refinement(**refiner_extra_kwargs) density_discr = qbx.density_discr @@ -203,12 +272,14 @@ def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, wave_vec = np.array([np.cos(angle), np.sin(angle)]) u = np.exp(1j*k*np.tensordot(wave_vec, nodes_host, axes=1)) grad_u = 1j*k*wave_vec[:, np.newaxis]*u - else: + elif d == 3: center = np.array([3, 1, 2]) diff = nodes_host - center[:, np.newaxis] r = la.norm(diff, axis=0) u = np.exp(1j*k*r) / r grad_u = diff * (1j*k*u/r - u/r**2) + else: + raise ValueError("invalid dim") else: center = np.array([3, 1, 2])[:d] diff = nodes_host - center[:, np.newaxis] @@ -233,11 +304,12 @@ def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, dn_u_dev = cl.array.to_device(queue, dn_u) grad_u_dev = cl.array.to_device(queue, grad_u) - key = (qbx_order, mesh_name, refinement_increment, zero_op_name) + key = (case.qbx_order, case.geometry.mesh_name, resolution, + case.expr.zero_op_name) - bound_op = bind(qbx, zero_op) + bound_op = bind(qbx, case.expr.get_zero_op(k_sym, **knl_kwargs)) error = bound_op( - queue, u=u_dev, dn_u=dn_u_dev, grad_u=grad_u_dev, k=k) + queue, u=u_dev, dn_u=dn_u_dev, grad_u=grad_u_dev, k=case.k) if 0: pt.plot(error) pt.show() @@ -248,7 +320,7 @@ def test_identities(ctx_getter, zero_op_name, mesh_name, mesh_getter, qbx_order, eoc_rec.add_data_point(qbx.h_max, l2_error_norm) print(eoc_rec) - tgt_order = order_table[zero_op_name] + tgt_order = case.qbx_order - case.expr.order_drop assert eoc_rec.order_estimate() > tgt_order - 1.3 # }}} -- GitLab From 4157448c6789ce5e18ae64d90e20e574a7502ac6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 12:14:25 -0500 Subject: [PATCH 11/24] Track boxtree's FMMLibExpansionWrangler rename --- pytential/qbx/fmmlib.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index c19dc3a4..0542cb5b 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -26,7 +26,7 @@ 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 HelmholtzExpansionWrangler +from boxtree.pyfmmlib_integration import FMMLibExpansionWrangler from sumpy.kernel import HelmholtzKernel @@ -54,7 +54,7 @@ class QBXFMMLibExpansionWranglerCodeContainer(object): source_extra_kwargs={}, kernel_extra_kwargs=None): - return QBXFMMLibHelmholtzExpansionWrangler(self, queue, geo_data, dtype, + return QBXFMMLibExpansionWrangler(self, queue, geo_data, dtype, qbx_order, fmm_level_to_order, source_extra_kwargs, kernel_extra_kwargs) @@ -113,7 +113,7 @@ class ToHostTransferredGeoDataWrapper(object): # {{{ fmmlib expansion wrangler -class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): +class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): def __init__(self, code, queue, geo_data, dtype, qbx_order, fmm_level_to_order, source_extra_kwargs, @@ -170,7 +170,7 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): from pytools import single_valued assert single_valued(self.level_orders) - super(QBXFMMLibHelmholtzExpansionWrangler, self).__init__( + super(QBXFMMLibExpansionWrangler, self).__init__( self.geo_data.tree(), helmholtz_k=helmholtz_k, @@ -204,8 +204,10 @@ class QBXFMMLibHelmholtzExpansionWrangler(HelmholtzExpansionWrangler): for k in self.outputs]) def reorder_sources(self, source_array): - source_array = source_array.get(queue=self.queue) - return (super(QBXFMMLibHelmholtzExpansionWrangler, self) + if isinstance(source_array, cl.array.Array): + source_array = source_array.get(queue=self.queue) + + return (super(QBXFMMLibExpansionWrangler, self) .reorder_sources(source_array)) def reorder_potentials(self, potentials): -- GitLab From e5df281d3875e8281391214053dee953babec7e0 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 12:19:54 -0500 Subject: [PATCH 12/24] Ditch non-parallel m2qbxl --- pytential/qbx/fmmlib.py | 196 +++++++++++++++------------------------- 1 file changed, 71 insertions(+), 125 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index 0542cb5b..a14d8921 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -366,142 +366,88 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): qbx_centers = geo_data.centers() centers = self.tree.box_centers - if 1: - # {{{ parallel + mploc = self.get_translation_routine("%ddmploc", vec_suffix="_imany") - mploc = self.get_translation_routine("%ddmploc", vec_suffix="_imany") + 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) - 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) + print("par data prep lev %d" % isrc_level) - # FIXME - - print("par data prep lev %d" % isrc_level) - - ngqbx_centers = len(geo_data.global_qbx_centers()) - 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) - - kwargs = {} - if self.dim == 3: - # FIXME Is this right? - kwargs["radius"] = ( - np.ones(ngqbx_centers) - * self.tree.root_extent * 2**(-isrc_level)) - - nsrc_boxes_per_gqbx_center = ( - ssn.starts[icontaining_tgt_box_vec+1] - - ssn.starts[icontaining_tgt_box_vec]) - nsrc_boxes = np.sum(nsrc_boxes_per_gqbx_center) - - src_boxes_starts = np.empty(ngqbx_centers+1, dtype=np.int32) - src_boxes_starts[0] = 0 - src_boxes_starts[1:] = np.cumsum(nsrc_boxes_per_gqbx_center) - - # FIXME - rscale1 = np.ones(nsrc_boxes) - rscale1_offsets = np.arange(nsrc_boxes) - - src_ibox = np.empty(nsrc_boxes, dtype=np.int32) - for itgt_center, tgt_icenter in enumerate( - geo_data.global_qbx_centers()): - icontaining_tgt_box = qbx_center_to_target_box[tgt_icenter] - src_ibox[ - src_boxes_starts[itgt_center]: - src_boxes_starts[itgt_center+1]] = ( - ssn.lists[ - ssn.starts[icontaining_tgt_box]: - ssn.starts[icontaining_tgt_box+1]]) - - del itgt_center - del tgt_icenter - del icontaining_tgt_box - - print("end par data prep") - - # These get max'd/added onto: pass initialized versions. - ier = np.zeros(ngqbx_centers, dtype=np.int32) - expn2 = np.zeros( - (ngqbx_centers,) + self.expansion_shape(self.qbx_order), - dtype=self.dtype) - - expn2 = mploc( - rscale1=rscale1, - rscale1_offsets=rscale1_offsets, - rscale1_starts=src_boxes_starts, - - center1=centers, - center1_offsets=src_ibox, - center1_starts=src_boxes_starts, - - expn1=source_mpoles_view.T, - expn1_offsets=src_ibox - source_level_start_ibox, - expn1_starts=src_boxes_starts, + ngqbx_centers = len(geo_data.global_qbx_centers()) + tgt_icenter_vec = geo_data.global_qbx_centers() + icontaining_tgt_box_vec = qbx_center_to_target_box[tgt_icenter_vec] - rscale2=rscale2, - # FIXME: center2 has wrong layout, will copy - center2=qbx_centers[:, tgt_icenter_vec], - expn2=expn2.T, - ier=ier, + # FIXME + rscale2 = np.ones(ngqbx_centers, np.float64) - **kwargs, - **self.kernel_kwargs).T - - local_exps[geo_data.global_qbx_centers()] += expn2 - - # }}} - - if 0: - # {{{ sequential - - mploc = self.get_translation_routine("%ddmploc") - - local_exps_1 = self.qbx_local_expansion_zeros() - - 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] + kwargs = {} + if self.dim == 3: + # FIXME Is this right? + kwargs["radius"] = ( + np.ones(ngqbx_centers) + * self.tree.root_extent * 2**(-isrc_level)) + + nsrc_boxes_per_gqbx_center = ( + ssn.starts[icontaining_tgt_box_vec+1] + - ssn.starts[icontaining_tgt_box_vec]) + nsrc_boxes = np.sum(nsrc_boxes_per_gqbx_center) + + src_boxes_starts = np.empty(ngqbx_centers+1, dtype=np.int32) + src_boxes_starts[0] = 0 + src_boxes_starts[1:] = np.cumsum(nsrc_boxes_per_gqbx_center) + + # FIXME + rscale1 = np.ones(nsrc_boxes) + rscale1_offsets = np.arange(nsrc_boxes) + + src_ibox = np.empty(nsrc_boxes, dtype=np.int32) + for itgt_center, tgt_icenter in enumerate( + geo_data.global_qbx_centers()): + icontaining_tgt_box = qbx_center_to_target_box[tgt_icenter] + src_ibox[ + src_boxes_starts[itgt_center]: + src_boxes_starts[itgt_center+1]] = ( + ssn.lists[ + ssn.starts[icontaining_tgt_box]: + ssn.starts[icontaining_tgt_box+1]]) + + del itgt_center + del tgt_icenter + del icontaining_tgt_box + + print("end par data prep") + + # These get max'd/added onto: pass initialized versions. + ier = np.zeros(ngqbx_centers, dtype=np.int32) + expn2 = np.zeros( + (ngqbx_centers,) + self.expansion_shape(self.qbx_order), + dtype=self.dtype) - for isrc_box in range( - ssn.starts[icontaining_tgt_box], - ssn.starts[icontaining_tgt_box+1]): + expn2 = mploc( + rscale1=rscale1, + rscale1_offsets=rscale1_offsets, + rscale1_starts=src_boxes_starts, - src_ibox = ssn.lists[isrc_box] - src_center = centers[:, src_ibox] + center1=centers, + center1_offsets=src_ibox, + center1_starts=src_boxes_starts, - ctr_loc = ctr_loc + mploc( - self.helmholtz_k, - rscale, src_center, - source_mpoles_view[src_ibox - source_level_start_ibox].T, - rscale, tgt_center, self.nterms, **kwargs)[..., 0].T + expn1=source_mpoles_view.T, + expn1_offsets=src_ibox - source_level_start_ibox, + expn1_starts=src_boxes_starts, - local_exps_1[tgt_icenter] += ctr_loc + rscale2=rscale2, + # FIXME: center2 has wrong layout, will copy + center2=qbx_centers[:, tgt_icenter_vec], + expn2=expn2.T, + ier=ier, - # diff = local_exps - local_exps_1 + **kwargs, + **self.kernel_kwargs).T - # }}} + local_exps[geo_data.global_qbx_centers()] += expn2 return local_exps -- GitLab From 5e40f891c24352c5779ad3e4b285f19f6dac2412 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 12:32:46 -0500 Subject: [PATCH 13/24] Initial (non-working) support for Fmmlib double layers --- pytential/qbx/fmmlib.py | 61 +++++++++++++++++++++++++++++++++++------ 1 file changed, 52 insertions(+), 9 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index a14d8921..83823b82 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -130,11 +130,18 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): # {{{ digest out_kernels - from sumpy.kernel import AxisTargetDerivative + from sumpy.kernel import AxisTargetDerivative, DirectionalSourceDerivative k_names = [] + source_deriv_names = [] def is_supported_helmknl(knl): + if isinstance(knl, DirectionalSourceDerivative): + source_deriv_names.append(knl.dir_vec_name) + knl = knl.inner_kernel + else: + source_deriv_names.append(None) + result = isinstance(knl, HelmholtzKernel) and knl.dim == 3 if result: k_names.append(knl.helmholtz_k_name) @@ -154,6 +161,12 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): "only the 3D Helmholtz kernel and its target derivatives " "are supported for now") + from pytools import is_single_valued + if not is_single_valued(source_deriv_names): + raise ValueError("not all kernels passed are the same in" + "whether they represent a source derivative") + + source_deriv_name = source_deriv_names[0] self.outputs = outputs # }}} @@ -170,10 +183,18 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): from pytools import single_valued assert single_valued(self.level_orders) + dipole_vec = None + if source_deriv_name is not None: + dipole_vec = np.array([ + d_i.get(queue=queue) + for d_i in source_extra_kwargs[source_deriv_name]], + order="F") + super(QBXFMMLibExpansionWrangler, self).__init__( self.geo_data.tree(), helmholtz_k=helmholtz_k, + dipole_vec=dipole_vec, # FIXME nterms=fmm_level_to_order(0), @@ -288,8 +309,13 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): center_source_starts = np.cumsum(center_source_counts) nsources_total = center_source_starts[-1] - sources = np.empty((3, nsources_total), dtype=np.float64) - charges = np.empty(nsources_total, dtype=np.complex128) + sources = np.empty((self.dim, nsources_total), dtype=np.float64) + charge = np.empty(nsources_total, dtype=np.complex128) + + full_dipvec = self.dipole_vec + if full_dipvec is not None: + dipvec = np.empty((self.dim, nsources_total), dtype=np.float64, + order="F") isource = 0 for itgt_center, tgt_icenter in enumerate(geo_data.global_qbx_centers()): @@ -305,7 +331,11 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): src_pslice = self._get_source_slice(src_ibox) ns = self.tree.box_source_counts_nonchild[src_ibox] sources[:, isource:isource+ns] = self._get_sources(src_pslice) - charges[isource:isource+ns] = src_weights[src_pslice] + charge[isource:isource+ns] = src_weights[src_pslice] + + if full_dipvec is not None: + dipvec[:, isource:isource+ns] = full_dipvec[:, src_pslice] + isource += ns centers = qbx_centers[:, geo_data.global_qbx_centers()] @@ -313,15 +343,26 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): rscale_vec = np.empty(len(center_source_counts) - 1, dtype=np.float64) rscale_vec.fill(rscale) # FIXME + source_kwargs = {} + if self.dipole_vec is None: + source_kwargs["charge"] = charge + source_kwargs["charge_offsets"] = center_source_starts[:-1] + + else: + source_kwargs["dipstr"] = charge + source_kwargs["dipstr_offsets"] = center_source_starts[:-1] + source_kwargs["dipvec"] = dipvec + source_kwargs["dipvec_offsets"] = center_source_starts[:-1] + logger.info("preparing interaction list for p2qbxl: done") return P2QBXLInfo( sources=sources, centers=centers, - charges=charges, center_source_starts=center_source_starts, center_source_counts=center_source_counts, rscale_vec=rscale_vec, + source_kwargs=source_kwargs, ) def form_global_qbx_locals(self, src_weights): @@ -332,22 +373,24 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): if len(geo_data.global_qbx_centers()) == 0: return local_exps - formta_vec = self.get_vec_routine("%ddformta") + formta_vec = self.get_vec_routine("%ddformta" + self.dp_suffix) info = self._info_for_form_global_qbx_locals(src_weights) + kwargs = {} + kwargs.update(info.source_kwargs) + kwargs.update(self.kernel_kwargs) + ier, loc_exp_pre = formta_vec( rscale=info.rscale_vec, sources=info.sources, sources_offsets=info.center_source_starts[:-1], - charge=info.charges, - charge_offsets=info.center_source_starts[:-1], nsources=info.center_source_counts[1:], center=info.centers, nterms=self.nterms, - **self.kernel_kwargs) + **kwargs) if np.any(ier != 0): raise RuntimeError("formta returned an error") -- GitLab From 165440b40e3314ea9f796532f8d372293ad24a0f Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 13:12:11 -0500 Subject: [PATCH 14/24] Add manually-moved lp test changes from merge --- test/test_layer_pot_identity.py | 4 +- test/test_scalar_int_eq.py | 114 ++++++++++++++++++++++++++------ 2 files changed, 97 insertions(+), 21 deletions(-) diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index 8486cd70..9b2fdb1b 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -254,8 +254,10 @@ def test_identity_convergence(ctx_getter, case): qbx, _ = QBXLayerPotentialSource( pre_density_discr, 4*target_order, - case.qbx_order, fmm_order=case.fmm_order, + case.qbx_order, + fmm_order=case.fmm_order, fmm_backend=case.fmm_backend, + _expansion_disks_in_tree_have_extent=True, ).with_refinement(**refiner_extra_kwargs) density_discr = qbx.density_discr diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py index 10e20519..187cf5f9 100644 --- a/test/test_scalar_int_eq.py +++ b/test/test_scalar_int_eq.py @@ -77,9 +77,6 @@ def run_int_eq_test(cl_ctx, queue, case, resolution): refiner_extra_kwargs = {} - if case.k != 0: - refiner_extra_kwargs["kernel_length_scale"] = 5/case.k - if case.fmm_backend is None: fmm_order = False else: @@ -93,7 +90,12 @@ def run_int_eq_test(cl_ctx, queue, case, resolution): fmm_order=fmm_order, fmm_backend=case.fmm_backend) if case.use_refinement: + if case.k != 0: + refiner_extra_kwargs["kernel_length_scale"] = 5/case.k + + print("%d elements before refinement" % pre_density_discr.mesh.nelements) qbx, _ = qbx.with_refinement(**refiner_extra_kwargs) + print("%d elements after refinement" % qbx.density_discr.mesh.nelements) density_discr = qbx.density_discr @@ -112,12 +114,12 @@ def run_int_eq_test(cl_ctx, queue, case, resolution): elif mesh.ambient_dim == 3: from meshmode.discretization.visualization import make_visualizer - bdry_vis = make_visualizer(queue, density_discr, case.target_order) + bdry_vis = make_visualizer(queue, density_discr, case.target_order+3) bdry_normals = bind(density_discr, sym.normal(3))(queue)\ .as_vector(dtype=object) - bdry_vis.write_vtk_file("source-%s.vtu" % resolution, [ + bdry_vis.write_vtk_file("pre-solve-source-%s.vtu" % resolution, [ ("bdry_normals", bdry_normals), ]) @@ -320,7 +322,7 @@ def run_int_eq_test(cl_ctx, queue, case, resolution): if 0: from meshmode.discretization.visualization import make_visualizer - bdry_vis = make_visualizer(queue, density_discr, case.target_order) + bdry_vis = make_visualizer(queue, density_discr, case.target_order+3) bdry_normals = bind(density_discr, sym.normal(3))(queue)\ .as_vector(dtype=object) @@ -338,13 +340,13 @@ def run_int_eq_test(cl_ctx, queue, case, resolution): fplot = FieldPlotter( bbox_center, extent=2*2*bbox_size, npoints=(150, 150, 1)) - qbx_stick_out = qbx.copy(target_stick_out_factor=0.15) + qbx_tgt_tol = qbx.copy(target_association_tolerance=0.15) from pytential.target import PointsTarget from pytential.qbx import QBXTargetAssociationFailedException try: solved_pot = bind( - (qbx_stick_out, PointsTarget(fplot.points)), + (qbx_tgt_tol, PointsTarget(fplot.points)), op.representation(sym.var("u")) )(queue, u=u, k=case.k) except QBXTargetAssociationFailedException as e: @@ -512,7 +514,6 @@ class CurveIntEqTestCase(IntEqTestCase): np.linspace(0, 1, resolution+1), target_order) - fmm_backend = None use_refinement = True neumann_alpha = None # default @@ -521,6 +522,7 @@ class CurveIntEqTestCase(IntEqTestCase): qbx_order = 5 target_order = qbx_order + fmm_backend = None check_tangential_deriv = True @@ -532,7 +534,22 @@ class EllipseIntEqTestCase(CurveIntEqTestCase): return ellipse(3, x) -class EllipsoidIntEqTestCase(IntEqTestCase): +class Helmholtz3DIntEqTestCase(IntEqTestCase): + fmm_backend = "fmmlib" + use_refinement = False + neumann_alpha = 0 # no double layers in FMMlib backend yet + + qbx_order = 2 + fmm_order = 7 + target_order = qbx_order + check_tangential_deriv = False + + # We're only expecting three digits based on FMM settings. Who are we + # kidding? + gmres_tol = 1e-5 + + +class EllipsoidIntEqTestCase(Helmholtz3DIntEqTestCase): resolutions = [2, 0.8] name = "ellipsoid" @@ -548,21 +565,78 @@ class EllipsoidIntEqTestCase(IntEqTestCase): # Flip elements--gmsh generates inside-out geometry. return perform_flips(mesh, np.ones(mesh.nelements)) - fmm_backend = "fmmlib" - fmm_order = 25 - use_refinement = False - neumann_alpha = 0 # no double layers in FMMlib backend yet + fmm_order = 13 inner_radius = 0.4 outer_radius = 5 - qbx_order = 4 - target_order = qbx_order - check_tangential_deriv = False - # We're only expecting three digits based on FMM settings. Who are we - # kidding? - gmres_tol = 1e-5 +class MergedCubesIntEqTestCase(Helmholtz3DIntEqTestCase): + resolutions = [1.4] + name = "merged-cubes" + + def get_mesh(self, resolution, target_order): + from meshmode.mesh.io import generate_gmsh, FileSource + mesh = generate_gmsh( + FileSource("merged-cubes.step"), 2, order=2, + other_options=[ + "-string", + "Mesh.CharacteristicLengthMax = %g;" % resolution]) + + from meshmode.mesh.processing import perform_flips + # Flip elements--gmsh generates inside-out geometry. + mesh = perform_flips(mesh, np.ones(mesh.nelements)) + + return mesh + + use_refinement = True + + inner_radius = 0.4 + outer_radius = 12 + + +class ManyEllipsoidIntEqTestCase(Helmholtz3DIntEqTestCase): + resolutions = [2, 1] + name = "ellipsoid" + + nx = 2 + ny = 2 + nz = 2 + + def get_mesh(self, resolution, target_order): + from meshmode.mesh.io import generate_gmsh, FileSource + base_mesh = generate_gmsh( + FileSource("ellipsoid.step"), 2, order=2, + other_options=[ + "-string", + "Mesh.CharacteristicLengthMax = %g;" % resolution]) + + from meshmode.mesh.processing import perform_flips + # Flip elements--gmsh generates inside-out geometry. + base_mesh = perform_flips(base_mesh, np.ones(base_mesh.nelements)) + + from meshmode.mesh.processing import affine_map, merge_disjoint_meshes + from meshmode.mesh.tools import rand_rotation_matrix + pitch = 10 + meshes = [ + affine_map( + base_mesh, + A=rand_rotation_matrix(3), + b=pitch*np.array([ + (ix-self.nx//2), + (iy-self.ny//2), + (iz-self.ny//2)])) + for ix in range(self.nx) + for iy in range(self.ny) + for iz in range(self.nz) + ] + + mesh = merge_disjoint_meshes(meshes, single_group=True) + return mesh + + inner_radius = 0.4 + # This should sit in the area just outside the middle ellipsoid + outer_radius = 5 # }}} -- GitLab From a18bededa5588c0ef01b3efce9f46226ed67a292 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 14:19:18 -0500 Subject: [PATCH 15/24] Add optional visualization to identities test --- test/test_layer_pot_identity.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index 9b2fdb1b..98a5a085 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -206,7 +206,7 @@ class DynamicTestCase(object): DynamicTestCase(geom, GradGreenExpr(), 1.2), DynamicTestCase(geom, ZeroCalderonExpr(), 0), ]]) -def test_identity_convergence(ctx_getter, case): +def test_identity_convergence(ctx_getter, case, visualize=False): logging.basicConfig(level=logging.INFO) case.check() @@ -321,6 +321,19 @@ def test_identity_convergence(ctx_getter, case): eoc_rec.add_data_point(qbx.h_max, l2_error_norm) + if visualize: + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, density_discr, target_order+3) + + bdry_normals = bind(density_discr, sym.normal(3))(queue)\ + .as_vector(dtype=object) + + bdry_vis.write_vtk_file("source-%s.vtu" % resolution, [ + ("u", u_dev), + ("bdry_normals", bdry_normals), + ("error", error), + ]) + print(eoc_rec) tgt_order = case.qbx_order - case.expr.order_drop assert eoc_rec.order_estimate() > tgt_order - 1.3 -- GitLab From 77d6c6c81cf5950ad74e9d4f418512730c32f292 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 17:52:04 -0500 Subject: [PATCH 16/24] Fmmlib: Py2.7-compatible kwarg passing --- pytential/qbx/fmmlib.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index 83823b82..5e32d492 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -468,6 +468,8 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): (ngqbx_centers,) + self.expansion_shape(self.qbx_order), dtype=self.dtype) + kwargs.update(self.kernel_kwargs) + expn2 = mploc( rscale1=rscale1, rscale1_offsets=rscale1_offsets, @@ -487,8 +489,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): expn2=expn2.T, ier=ier, - **kwargs, - **self.kernel_kwargs).T + **kwargs).T local_exps[geo_data.global_qbx_centers()] += expn2 -- GitLab From ba99200be6eee46ebbc4c957eaed6868b15e1d29 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 17:52:59 -0500 Subject: [PATCH 17/24] Fmmlib: Dipoles coming in through the QBX FMM are already reordered --- pytential/qbx/fmmlib.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index 5e32d492..3d040817 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -195,6 +195,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): helmholtz_k=helmholtz_k, dipole_vec=dipole_vec, + dipoles_already_reordered=True, # FIXME nterms=fmm_level_to_order(0), -- GitLab From 29fc4178f4d0ecf28a22debd0dc078a2f206b53b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 17:54:32 -0500 Subject: [PATCH 18/24] FMMLib: Create formta sources temp array in F order for a suspected minor perf gain --- pytential/qbx/fmmlib.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index 3d040817..bfa201e3 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -310,7 +310,8 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): center_source_starts = np.cumsum(center_source_counts) nsources_total = center_source_starts[-1] - sources = np.empty((self.dim, nsources_total), dtype=np.float64) + sources = np.empty((self.dim, nsources_total), dtype=np.float64, + order="F") charge = np.empty(nsources_total, dtype=np.complex128) full_dipvec = self.dipole_vec -- GitLab From 20e99d437182a1c1d5adfd6ce32bfe8c3ae5f1b6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 18:10:13 -0500 Subject: [PATCH 19/24] API terminology: don't talk about expansion *disks* (no such thing in 3D) --- pytential/qbx/__init__.py | 30 +++++++++++++++--------------- pytential/qbx/geometry.py | 6 +++--- test/test_layer_pot.py | 4 ++-- test/test_layer_pot_identity.py | 2 +- test/test_stokes.py | 2 +- 5 files changed, 22 insertions(+), 22 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 16814fec..3568601f 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -75,8 +75,8 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): # FIXME default debug=False once everything works debug=True, _refined_for_global_qbx=False, - _expansion_disks_in_tree_have_extent=False, - _expansion_disk_stick_out_factor=0, + _expansions_in_tree_have_extent=False, + _expansion_stick_out_factor=0, performance_data_file=None, fmm_backend="sumpy", target_stick_out_factor=_not_provided): @@ -155,9 +155,9 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): 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._expansion_disk_stick_out_factor = _expansion_disk_stick_out_factor + self._expansions_in_tree_have_extent = \ + _expansions_in_tree_have_extent + self._expansion_stick_out_factor = _expansion_stick_out_factor self.performance_data_file = performance_data_file def copy( @@ -168,8 +168,8 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): fmm_level_to_order=None, base_resampler=None, target_association_tolerance=_not_provided, - _expansion_disks_in_tree_have_extent=_not_provided, - _expansion_disk_stick_out_factor=_not_provided, + _expansions_in_tree_have_extent=_not_provided, + _expansion_stick_out_factor=_not_provided, performance_data_file=None, debug=_not_provided, @@ -220,16 +220,16 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): _refined_for_global_qbx if _refined_for_global_qbx is not _not_provided else self._refined_for_global_qbx), - _expansion_disks_in_tree_have_extent=( + _expansions_in_tree_have_extent=( # False is a valid value here - _expansion_disks_in_tree_have_extent - if _expansion_disks_in_tree_have_extent is not _not_provided - else self._expansion_disks_in_tree_have_extent), - _expansion_disk_stick_out_factor=( + _expansions_in_tree_have_extent + if _expansions_in_tree_have_extent is not _not_provided + else self._expansions_in_tree_have_extent), + _expansion_stick_out_factor=( # 0 is a valid value here - _expansion_disk_stick_out_factor - if _expansion_disk_stick_out_factor is not _not_provided - else self._expansion_disk_stick_out_factor), + _expansion_stick_out_factor + if _expansion_stick_out_factor is not _not_provided + else self._expansion_stick_out_factor), performance_data_file=( performance_data_file or self.performance_data_file), fmm_backend=self.fmm_backend) diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index b08370df..b43dbd4a 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -493,7 +493,7 @@ class QBXFMMGeometryData(object): nparticles = nsources + target_info.ntargets target_radii = None - if self.lpot_source._expansion_disks_in_tree_have_extent: + if self.lpot_source._expansions_in_tree_have_extent: target_radii = cl.array.zeros(queue, target_info.ntargets, self.coord_dtype) target_radii[:self.ncenters] = self.expansion_radii() @@ -520,7 +520,7 @@ class QBXFMMGeometryData(object): max_leaf_refine_weight=32, refine_weights=refine_weights, debug=self.debug, - stick_out_factor=lpot_src._expansion_disk_stick_out_factor, + stick_out_factor=lpot_src._expansion_stick_out_factor, kind="adaptive") if self.debug: @@ -545,7 +545,7 @@ class QBXFMMGeometryData(object): trav, _ = self.code_getter.build_traversal(queue, self.tree(), debug=self.debug) - if self.lpot_source._expansion_disks_in_tree_have_extent: + if self.lpot_source._expansions_in_tree_have_extent: trav = trav.merge_close_lists(queue) return trav diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py index bb104993..280a918c 100644 --- a/test/test_layer_pot.py +++ b/test/test_layer_pot.py @@ -146,7 +146,7 @@ def test_ellipse_eigenvalues(ctx_getter, ellipse_aspect, mode_nr, qbx_order): qbx, _ = QBXLayerPotentialSource( pre_density_discr, 4*target_order, qbx_order, fmm_order=fmm_order, - _expansion_disks_in_tree_have_extent=True, + _expansions_in_tree_have_extent=True, ).with_refinement() density_discr = qbx.density_discr @@ -383,7 +383,7 @@ def test_off_surface_eval_vs_direct(ctx_getter, do_plot=False): fmm_qbx, _ = QBXLayerPotentialSource( pre_density_discr, 4*target_order, qbx_order, fmm_order=qbx_order + 3, - _expansion_disks_in_tree_have_extent=True, + _expansions_in_tree_have_extent=True, target_association_tolerance=0.05, ).with_refinement() diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index 98a5a085..285fe576 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -257,7 +257,7 @@ def test_identity_convergence(ctx_getter, case, visualize=False): case.qbx_order, fmm_order=case.fmm_order, fmm_backend=case.fmm_backend, - _expansion_disks_in_tree_have_extent=True, + _expansions_in_tree_have_extent=True, ).with_refinement(**refiner_extra_kwargs) density_discr = qbx.density_discr diff --git a/test/test_stokes.py b/test/test_stokes.py index a966e5f6..e6fb8537 100644 --- a/test/test_stokes.py +++ b/test/test_stokes.py @@ -74,7 +74,7 @@ def run_exterior_stokes_2d(ctx_factory, nelements, coarse_density_discr, fine_order=ovsmp_target_order, qbx_order=qbx_order, fmm_order=fmm_order, target_association_tolerance=target_association_tolerance, - _expansion_disks_in_tree_have_extent=True, + _expansions_in_tree_have_extent=True, ).with_refinement() density_discr = qbx.density_discr -- GitLab From 9002c558f08e7bc6cf15fdeee64389c7d19c4825 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 28 Jun 2017 18:28:56 -0500 Subject: [PATCH 20/24] Tweak Helmholtz sphere Green test --- test/test_layer_pot_identity.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index 285fe576..af073cc3 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -160,7 +160,11 @@ class SphereGreenTest(StaticTestCase): geometry = SphereGeometry() k = 1.2 qbx_order = 3 - fmm_order = 15 + fmm_order = 10 + + resolutions = [0, 1] + + _expansion_stick_out_factor = 0.75 fmm_backend = "fmmlib" @@ -223,7 +227,10 @@ def test_identity_convergence(ctx_getter, case, visualize=False): from pytools.convergence import EOCRecorder eoc_rec = EOCRecorder() - for resolution in case.geometry.resolutions: + for resolution in ( + getattr(case, "resolutions", None) + or case.geometry.resolutions + ): mesh = case.geometry.get_mesh(resolution, target_order) if mesh is None: break @@ -258,6 +265,8 @@ def test_identity_convergence(ctx_getter, case, visualize=False): fmm_order=case.fmm_order, fmm_backend=case.fmm_backend, _expansions_in_tree_have_extent=True, + _expansion_stick_out_factor=getattr( + case, "_expansion_stick_out_factor", 0), ).with_refinement(**refiner_extra_kwargs) density_discr = qbx.density_discr @@ -323,7 +332,7 @@ def test_identity_convergence(ctx_getter, case, visualize=False): if visualize: from meshmode.discretization.visualization import make_visualizer - bdry_vis = make_visualizer(queue, density_discr, target_order+3) + bdry_vis = make_visualizer(queue, density_discr, target_order) bdry_normals = bind(density_discr, sym.normal(3))(queue)\ .as_vector(dtype=object) -- GitLab From 8f208305e31b0b6567d96b3f7a81805d3ae58b3c Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 29 Jun 2017 05:05:46 -0500 Subject: [PATCH 21/24] Fmmlib: Py2.7-compatible kwarg passing --- pytential/qbx/fmmlib.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index bfa201e3..dc3d5924 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -525,6 +525,8 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): # FIXME Is this right? kwargs["radius"] = self.tree.root_extent * 2**(-isrc_level) + kwargs.update(self.kernel_kwargs) + for tgt_icenter in range(geo_data.ncenters): isrc_box = qbx_center_to_target_box[tgt_icenter] @@ -550,8 +552,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): center2=tgt_center, nterms2=local_order, - **kwargs, - **self.kernel_kwargs)[..., 0].T + **kwargs)[..., 0].T qbx_expansions[tgt_icenter] += tmp_loc_exp -- GitLab From a8c896560a00bee6d7220fd48a1630957b05ef16 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 29 Jun 2017 06:14:36 -0500 Subject: [PATCH 22/24] Fmmlib parallelization: Switch to indirect-addressed p2qbxl --- pytential/qbx/fmmlib.py | 88 +++++++++++++++++++++++------------------ 1 file changed, 49 insertions(+), 39 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index dc3d5924..c948a1be 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -278,12 +278,10 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): # {{{ p2qbxl - #@memoize_method - def _info_for_form_global_qbx_locals(self, src_weights): + @memoize_method + def _info_for_form_global_qbx_locals(self): logger.info("preparing interaction list for p2qbxl: start") - rscale = 1 # FIXME - geo_data = self.geo_data traversal = geo_data.traversal() @@ -309,15 +307,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): center_source_counts = np.array(center_source_counts) center_source_starts = np.cumsum(center_source_counts) nsources_total = center_source_starts[-1] - - sources = np.empty((self.dim, nsources_total), dtype=np.float64, - order="F") - charge = np.empty(nsources_total, dtype=np.complex128) - - full_dipvec = self.dipole_vec - if full_dipvec is not None: - dipvec = np.empty((self.dim, nsources_total), dtype=np.float64, - order="F") + center_source_offsets = np.empty(nsources_total, np.int32) isource = 0 for itgt_center, tgt_icenter in enumerate(geo_data.global_qbx_centers()): @@ -332,39 +322,28 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): src_pslice = self._get_source_slice(src_ibox) ns = self.tree.box_source_counts_nonchild[src_ibox] - sources[:, isource:isource+ns] = self._get_sources(src_pslice) - charge[isource:isource+ns] = src_weights[src_pslice] - - if full_dipvec is not None: - dipvec[:, isource:isource+ns] = full_dipvec[:, src_pslice] + center_source_offsets[isource:isource+ns] = np.arange( + src_pslice.start, src_pslice.stop) 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 - source_kwargs = {} - if self.dipole_vec is None: - source_kwargs["charge"] = charge - source_kwargs["charge_offsets"] = center_source_starts[:-1] - - else: - source_kwargs["dipstr"] = charge - source_kwargs["dipstr_offsets"] = center_source_starts[:-1] - source_kwargs["dipvec"] = dipvec - source_kwargs["dipvec_offsets"] = center_source_starts[:-1] + nsources_vec = np.ones(self.tree.nsources, np.int32) logger.info("preparing interaction list for p2qbxl: done") return P2QBXLInfo( - sources=sources, centers=centers, center_source_starts=center_source_starts, - center_source_counts=center_source_counts, + center_source_offsets=center_source_offsets, + nsources_vec=nsources_vec, rscale_vec=rscale_vec, - source_kwargs=source_kwargs, + ngqbx_centers=centers.shape[1], ) def form_global_qbx_locals(self, src_weights): @@ -375,34 +354,63 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): if len(geo_data.global_qbx_centers()) == 0: return local_exps - formta_vec = self.get_vec_routine("%ddformta" + self.dp_suffix) - info = self._info_for_form_global_qbx_locals(src_weights) + formta_imany = self.get_routine("%ddformta" + self.dp_suffix, + suffix="_imany") + info = self._info_for_form_global_qbx_locals() kwargs = {} - kwargs.update(info.source_kwargs) kwargs.update(self.kernel_kwargs) - ier, loc_exp_pre = formta_vec( + if self.dipole_vec is None: + kwargs["charge"] = src_weights + kwargs["charge_offsets"] = info.center_source_offsets + 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 + + # These get max'd/added onto: pass initialized versions. + ier = np.zeros(info.ngqbx_centers, dtype=np.int32) + expn = np.zeros( + (info.ngqbx_centers,) + self.expansion_shape(self.qbx_order), + dtype=self.dtype) + + ier, expn = formta_imany( rscale=info.rscale_vec, - sources=info.sources, - sources_offsets=info.center_source_starts[:-1], + sources=self._get_single_sources_array(), + sources_offsets=info.center_source_offsets, + sources_starts=info.center_source_starts, + + nsources=info.nsources_vec, + nsources_offsets=info.center_source_offsets, + nsources_starts=info.center_source_starts, - nsources=info.center_source_counts[1:], center=info.centers, nterms=self.nterms, + ier=ier, + expn=expn.T, + **kwargs) if np.any(ier != 0): raise RuntimeError("formta returned an error") - local_exps[geo_data.global_qbx_centers()] = loc_exp_pre.T + local_exps[geo_data.global_qbx_centers()] = expn.T return local_exps # }}} + # {{{ m2qbxl + def translate_box_multipoles_to_qbx_local(self, multipole_exps): local_exps = self.qbx_local_expansion_zeros() @@ -497,6 +505,8 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): return local_exps + # }}} + def translate_box_local_to_qbx_local(self, local_exps): qbx_expansions = self.qbx_local_expansion_zeros() -- GitLab From 8ca203b67bf57276304ee251bb84bb1a3f698cf8 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 29 Jun 2017 07:43:47 -0500 Subject: [PATCH 23/24] Hack gradient checking into the inteq test script --- test/test_scalar_int_eq.py | 36 +++++++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py index 187cf5f9..9443261a 100644 --- a/test/test_scalar_int_eq.py +++ b/test/test_scalar_int_eq.py @@ -253,7 +253,8 @@ def run_int_eq_test(cl_ctx, queue, case, resolution): # {{{ error check - bound_tgt_op = bind((qbx, PointsTarget(test_targets)), + points_target = PointsTarget(test_targets) + bound_tgt_op = bind((qbx, points_target), op.representation(sym.var("u"))) test_via_bdry = bound_tgt_op(queue, u=u, k=case.k) @@ -282,6 +283,36 @@ def run_int_eq_test(cl_ctx, queue, case, resolution): print("rel_err_2: %g rel_err_inf: %g" % (rel_err_2, rel_err_inf)) + # {{{ test gradient + + if case.check_gradient: + bound_grad_op = bind((qbx, points_target), + op.representation( + sym.var("u"), + map_potentials=lambda pot: sym.grad(mesh.ambient_dim, pot), + qbx_forced_limit=None)) + + #print(bound_t_deriv_op.code) + + grad_from_src = bound_grad_op( + queue, u=u, **concrete_knl_kwargs) + + grad_ref = (bind( + (point_source, points_target), + sym.grad(mesh.ambient_dim, pot_src) + )(queue, charges=source_charges_dev, **concrete_knl_kwargs) + ) + + grad_err = (grad_from_src - grad_ref) + + rel_grad_err_inf = ( + la.norm(grad_err[0].get(), np.inf) + / + la.norm(grad_ref[0].get(), np.inf)) + + print("rel_grad_err_inf: %g" % rel_grad_err_inf) + + # }}} # {{{ test tangential derivative if case.check_tangential_deriv: @@ -525,6 +556,7 @@ class CurveIntEqTestCase(IntEqTestCase): fmm_backend = None check_tangential_deriv = True + check_gradient = False class EllipseIntEqTestCase(CurveIntEqTestCase): @@ -570,6 +602,8 @@ class EllipsoidIntEqTestCase(Helmholtz3DIntEqTestCase): inner_radius = 0.4 outer_radius = 5 + check_gradient = True + class MergedCubesIntEqTestCase(Helmholtz3DIntEqTestCase): resolutions = [1.4] -- GitLab From 16e052c8e1da78daea1f313781e834a65b562b2e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 29 Jun 2017 07:57:46 -0500 Subject: [PATCH 24/24] Fix fmmlib kernel processing [ci skip] --- pytential/qbx/fmmlib.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index c948a1be..9ec0194c 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -137,14 +137,15 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): def is_supported_helmknl(knl): if isinstance(knl, DirectionalSourceDerivative): - source_deriv_names.append(knl.dir_vec_name) + source_deriv_name = knl.dir_vec_name knl = knl.inner_kernel else: - source_deriv_names.append(None) + source_deriv_name = None result = isinstance(knl, HelmholtzKernel) and knl.dim == 3 if result: k_names.append(knl.helmholtz_k_name) + source_deriv_names.append(source_deriv_name) return result ifgrad = False @@ -163,7 +164,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): from pytools import is_single_valued if not is_single_valued(source_deriv_names): - raise ValueError("not all kernels passed are the same in" + raise ValueError("not all kernels passed are the same in " "whether they represent a source derivative") source_deriv_name = source_deriv_names[0] -- GitLab