diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 16814fec0e722a2082ad82649c09e10136ee6114..3568601ff3e02c81cb78839663cd386657f54e81 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/fmm.py b/pytential/qbx/fmm.py index 6afa047b692e02a536d93def31f5c08097bcc587..e9eaa2b2d0251c89ea4e56adbeaf28e8dc9c4cd8 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 bb95cb30a904b8592ea7b25f6d6d3d6089e19df6..9ec0194cbd1c3d8494d484f5c14063ca46d32d56 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 FMMLibExpansionWrangler 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, @@ -46,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) @@ -105,7 +113,7 @@ class ToHostTransferredGeoDataWrapper(object): # {{{ fmmlib expansion wrangler -class QBXFMMLibHelmholtzExpansionWrangler(FMMLibExpansionWrangler): +class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): def __init__(self, code, queue, geo_data, dtype, qbx_order, fmm_level_to_order, source_extra_kwargs, @@ -122,14 +130,22 @@ class QBXFMMLibHelmholtzExpansionWrangler(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_name = knl.dir_vec_name + knl = knl.inner_kernel + else: + 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 @@ -146,6 +162,12 @@ class QBXFMMLibHelmholtzExpansionWrangler(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 # }}} @@ -162,10 +184,19 @@ class QBXFMMLibHelmholtzExpansionWrangler(FMMLibExpansionWrangler): from pytools import single_valued assert single_valued(self.level_orders) - super(QBXFMMLibHelmholtzExpansionWrangler, self).__init__( + 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, + dipoles_already_reordered=True, # FIXME nterms=fmm_level_to_order(0), @@ -196,8 +227,10 @@ class QBXFMMLibHelmholtzExpansionWrangler(FMMLibExpansionWrangler): 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): @@ -239,55 +272,146 @@ class QBXFMMLibHelmholtzExpansionWrangler(FMMLibExpansionWrangler): # }}} - # {{{ 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): + 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() - formta = self.get_routine("%ddformta") - + 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) - ctr_coeffs = 0 + center_source_counts = np.array(center_source_counts) + center_source_starts = np.cumsum(center_source_counts) + nsources_total = center_source_starts[-1] + center_source_offsets = np.empty(nsources_total, np.int32) + + 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] + 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 + + nsources_vec = np.ones(self.tree.nsources, np.int32) + + logger.info("preparing interaction list for p2qbxl: done") + + return P2QBXLInfo( + centers=centers, + center_source_starts=center_source_starts, + center_source_offsets=center_source_offsets, + nsources_vec=nsources_vec, + rscale_vec=rscale_vec, + ngqbx_centers=centers.shape[1], + ) + + def form_global_qbx_locals(self, src_weights): + geo_data = self.geo_data + + local_exps = self.qbx_local_expansion_zeros() - 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 + + formta_imany = self.get_routine("%ddformta" + self.dp_suffix, + suffix="_imany") + info = self._info_for_form_global_qbx_locals() + + kwargs = {} + kwargs.update(self.kernel_kwargs) + + 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) - ctr_coeffs += coeffs + ier, expn = formta_imany( + rscale=info.rscale_vec, - local_exps[tgt_icenter] += ctr_coeffs + 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, + + 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()] = expn.T return local_exps + # }}} + + # {{{ m2qbxl + def translate_box_multipoles_to_qbx_local(self, multipole_exps): local_exps = self.qbx_local_expansion_zeros() @@ -296,43 +420,94 @@ class QBXFMMLibHelmholtzExpansionWrangler(FMMLibExpansionWrangler): qbx_centers = geo_data.centers() centers = self.tree.box_centers - mploc = self.get_translation_routine("%ddmploc") + mploc = self.get_translation_routine("%ddmploc", vec_suffix="_imany") - for isrc_level, ssn in enumerate(geo_data.traversal().sep_smaller_by_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) + + 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 - rscale = 1 + 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)) - for itgt_center, tgt_icenter in enumerate(geo_data.global_qbx_centers()): - ctr_loc = 0 + 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) - tgt_center = qbx_centers[:, tgt_icenter] + kwargs.update(self.kernel_kwargs) + + expn2 = mploc( + rscale1=rscale1, + rscale1_offsets=rscale1_offsets, + rscale1_starts=src_boxes_starts, - for isrc_box in range( - ssn.starts[icontaining_tgt_box], - ssn.starts[icontaining_tgt_box+1]): + center1=centers, + center1_offsets=src_ibox, + center1_starts=src_boxes_starts, - src_ibox = ssn.lists[isrc_box] - src_center = centers[:, src_ibox] + expn1=source_mpoles_view.T, + expn1_offsets=src_ibox - source_level_start_ibox, + expn1_starts=src_boxes_starts, - ctr_loc = ctr_loc + mploc( - self.helmholtz_k, - rscale, src_center, multipole_exps[src_ibox], - rscale, tgt_center, self.nterms, **kwargs)[..., 0] + rscale2=rscale2, + # FIXME: center2 has wrong layout, will copy + center2=qbx_centers[:, tgt_icenter_vec], + expn2=expn2.T, + ier=ier, - local_exps[tgt_icenter] += ctr_loc + **kwargs).T + + local_exps[geo_data.global_qbx_centers()] += expn2 return local_exps + # }}} + def translate_box_local_to_qbx_local(self, local_exps): qbx_expansions = self.qbx_local_expansion_zeros() @@ -361,6 +536,8 @@ class QBXFMMLibHelmholtzExpansionWrangler(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] @@ -378,9 +555,15 @@ class QBXFMMLibHelmholtzExpansionWrangler(FMMLibExpansionWrangler): if in_range: src_center = self.tree.box_centers[:, src_ibox] tmp_loc_exp = locloc( - self.helmholtz_k, - rscale, src_center, local_exps[src_ibox], - rscale, tgt_center, local_order, **kwargs)[..., 0] + rscale1=rscale, + center1=src_center, + expn1=local_exps[src_ibox].T, + + rscale2=rscale, + center2=tgt_center, + nterms2=local_order, + + **kwargs)[..., 0].T qbx_expansions[tgt_icenter] += tmp_loc_exp @@ -409,9 +592,12 @@ class QBXFMMLibHelmholtzExpansionWrangler(FMMLibExpansionWrangler): center = qbx_centers[:, src_icenter] - pot, grad = taeval(self.helmholtz_k, rscale, - center, qbx_expansions[src_icenter], - 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) @@ -426,8 +612,6 @@ class QBXFMMLibHelmholtzExpansionWrangler(FMMLibExpansionWrangler): return cl.array.to_device(self.queue, potential) * scale_factor - # }}} - # }}} # vim: foldmethod=marker diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 8c4611c01b0e3821e954e7cd2949020259025197..af8da4ca928524fb0085c7744331981118806568 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -462,7 +462,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() @@ -489,7 +489,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: @@ -514,7 +514,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 34d11d3f2924be38c109467a68d87fbf4617dbe0..280a918ca3399a3caa8427e1334d6d69babd0699 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): @@ -156,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 @@ -289,847 +279,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.fmm_backend is None: - fmm_order = False - else: - fmm_order = case.fmm_order - - 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: - 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 - - # {{{ 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 + 5) - - bdry_normals = bind(density_discr, sym.normal(3))(queue)\ - .as_vector(dtype=object) - - bdry_vis.write_vtk_file("pre-solve-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+5) - - 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_tgt_tol = qbx.copy(target_association_tolerance=0.15) - from pytential.target import PointsTarget - from pytential.qbx import QBXTargetAssociationFailedException - - try: - solved_pot = bind( - (qbx_tgt_tol, 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) - - use_refinement = True - neumann_alpha = None # default - - inner_radius = 0.1 - outer_radius = 2 - - qbx_order = 5 - target_order = qbx_order - fmm_backend = None - - check_tangential_deriv = True - - -class EllipseIntEqTestCase(CurveIntEqTestCase): - name = "3-to-1 ellipse" - - def curve_func(self, x): - return ellipse(3, x) - - -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, 1] - 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_order = 13 - - inner_radius = 0.4 - outer_radius = 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 - - -@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, - _expansion_disks_in_tree_have_extent=True, - ).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]) @@ -1234,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 new file mode 100644 index 0000000000000000000000000000000000000000..af073cc3561c05bc8a1329d638f0877fc615b68e --- /dev/null +++ b/test/test_layer_pot_identity.py @@ -0,0 +1,364 @@ +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 +from sumpy.kernel import LaplaceKernel, HelmholtzKernel + +import logging +logger = logging.getLogger(__name__) + +circle = partial(ellipse, 1) + +try: + import matplotlib.pyplot as pt +except ImportError: + pass + + +d1 = sym.Derivative() +d2 = sym.Derivative() + + +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 + + +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 = 10 + + resolutions = [0, 1] + + _expansion_stick_out_factor = 0.75 + + 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, visualize=False): + logging.basicConfig(level=logging.INFO) + + case.check() + + cl_ctx = ctx_getter() + queue = cl.CommandQueue(cl_ctx) + + # prevent cache 'splosion + from sympy.core.cache import clear_cache + clear_cache() + + target_order = 8 + + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + for resolution in ( + getattr(case, "resolutions", None) + or case.geometry.resolutions + ): + mesh = case.geometry.get_mesh(resolution, target_order) + if mesh is None: + break + + d = mesh.ambient_dim + k = case.k + + 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")} + + 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)) + + refiner_extra_kwargs = {} + + if case.k != 0: + refiner_extra_kwargs["kernel_length_scale"] = 5/case.k + + qbx, _ = QBXLayerPotentialSource( + pre_density_discr, 4*target_order, + case.qbx_order, + 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 + + # {{{ 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 + 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] + 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 = (case.qbx_order, case.geometry.mesh_name, resolution, + case.expr.zero_op_name) + + 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=case.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) + + if visualize: + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, density_discr, target_order) + + bdry_normals = bind(density_discr, sym.normal(3))(queue)\ + .as_vector(dtype=object) + + bdry_vis.write_vtk_file("source-%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 + +# }}} + + +# 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 0000000000000000000000000000000000000000..9443261aac95d329db2721d50b8fcb9f847b12d7 --- /dev/null +++ b/test/test_scalar_int_eq.py @@ -0,0 +1,748 @@ +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.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: + 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 + + # {{{ 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+3) + + bdry_normals = bind(density_discr, sym.normal(3))(queue)\ + .as_vector(dtype=object) + + bdry_vis.write_vtk_file("pre-solve-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 + + 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) + + 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 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: + 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+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), + ("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_tgt_tol = qbx.copy(target_association_tolerance=0.15) + from pytential.target import PointsTarget + from pytential.qbx import QBXTargetAssociationFailedException + + try: + solved_pot = bind( + (qbx_tgt_tol, 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) + + use_refinement = True + neumann_alpha = None # default + + inner_radius = 0.1 + outer_radius = 2 + + qbx_order = 5 + target_order = qbx_order + fmm_backend = None + + check_tangential_deriv = True + check_gradient = False + + +class EllipseIntEqTestCase(CurveIntEqTestCase): + name = "3-to-1 ellipse" + + def curve_func(self, x): + return ellipse(3, x) + + +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" + + 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_order = 13 + + inner_radius = 0.4 + outer_radius = 5 + + check_gradient = True + + +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 + +# }}} + + +# {{{ 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 diff --git a/test/test_stokes.py b/test/test_stokes.py index a966e5f6a84e2b1572fec9c3175b576b49a352aa..e6fb853751268a921c32a7f51ede8c3a1ce7fb7f 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