From 77130f0a59d6606414d817b9be667c42203836d3 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 8 Aug 2017 01:16:39 -0500 Subject: [PATCH 01/15] Initial adaptation to the rigorous-bound 2-away FMM in boxtree --- pytential/qbx/fmm.py | 26 +++++++++++++------------- pytential/qbx/geometry.py | 2 +- requirements.txt | 2 +- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py index e9eaa2b2..3ad00338 100644 --- a/pytential/qbx/fmm.py +++ b/pytential/qbx/fmm.py @@ -237,7 +237,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`), wait_for = multipole_exps.events - for isrc_level, ssn in enumerate(traversal.sep_smaller_by_level): + for isrc_level, ssn in enumerate(traversal.from_sep_smaller_by_level): m2qbxl = self.code.m2qbxl( self.level_orders[isrc_level], self.qbx_order) @@ -414,8 +414,8 @@ def drive_fmm(expansion_wrangler, src_weights): local_exps = wrangler.multipole_to_local( traversal.level_start_target_or_target_parent_box_nrs, traversal.target_or_target_parent_boxes, - traversal.sep_siblings_starts, - traversal.sep_siblings_lists, + traversal.from_sep_siblings_starts, + traversal.from_sep_siblings_lists, mpole_exps) # }}} @@ -430,11 +430,11 @@ def drive_fmm(expansion_wrangler, src_weights): non_qbx_potentials = non_qbx_potentials + wrangler.eval_multipoles( traversal.level_start_target_box_nrs, traversal.target_boxes, - traversal.sep_smaller_by_level, + traversal.from_sep_smaller_by_level, mpole_exps) # assert that list 3 close has been merged into list 1 - assert traversal.sep_close_smaller_starts is None + assert traversal.from_sep_close_smaller_starts is None # }}} @@ -445,12 +445,12 @@ def drive_fmm(expansion_wrangler, src_weights): local_exps = local_exps + wrangler.form_locals( traversal.level_start_target_or_target_parent_box_nrs, traversal.target_or_target_parent_boxes, - traversal.sep_bigger_starts, - traversal.sep_bigger_lists, + traversal.from_sep_bigger_starts, + traversal.from_sep_bigger_lists, src_weights) # assert that list 4 close has been merged into list 1 - assert traversal.sep_close_bigger_starts is None + assert traversal.from_sep_close_bigger_starts is None # }}} @@ -590,7 +590,7 @@ def write_performance_model(outf, geo_data): def process_list2(): nm2l = 0 for itgt_box, tgt_ibox in enumerate(traversal.target_or_target_parent_boxes): - start, end = traversal.sep_siblings_starts[itgt_box:itgt_box+2] + start, end = traversal.from_sep_siblings_starts[itgt_box:itgt_box+2] nm2l += (end-start) @@ -608,7 +608,7 @@ def write_performance_model(outf, geo_data): for itgt_box, tgt_ibox in enumerate(traversal.target_boxes): ntargets = box_target_counts_nonchild[tgt_ibox] - for sep_smaller_list in traversal.sep_smaller_by_level: + for sep_smaller_list in traversal.from_sep_smaller_by_level: start, end = sep_smaller_list.starts[itgt_box:itgt_box+2] nmp_eval += ntargets * (end-start) @@ -625,9 +625,9 @@ def write_performance_model(outf, geo_data): nform_local = 0 for itgt_box, tgt_ibox in enumerate(traversal.target_or_target_parent_boxes): - start, end = traversal.sep_bigger_starts[itgt_box:itgt_box+2] + start, end = traversal.from_sep_bigger_starts[itgt_box:itgt_box+2] - for src_ibox in traversal.sep_bigger_lists[start:end]: + for src_ibox in traversal.from_sep_bigger_lists[start:end]: nsources = tree.box_source_counts_nonchild[src_ibox] nform_local += nsources @@ -696,7 +696,7 @@ def write_performance_model(outf, geo_data): def process_m2qbxl(): nqbx_m2l = 0 - for isrc_level, ssn in enumerate(traversal.sep_smaller_by_level): + for isrc_level, ssn in enumerate(traversal.from_sep_smaller_by_level): for itgt_center, tgt_icenter in enumerate(global_qbx_centers): icontaining_tgt_box = qbx_center_to_target_box[tgt_icenter] diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 96bacdf1..6a0f0b08 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -138,7 +138,7 @@ class QBXFMMGeometryCodeGetter(object): @memoize_method def build_traversal(self): from boxtree.traversal import FMMTraversalBuilder - return FMMTraversalBuilder(self.cl_context) + return FMMTraversalBuilder(self.cl_context, well_sep_is_n_away=2) @memoize_method def qbx_center_to_target_box_lookup(self, particle_id_dtype, box_id_dtype): diff --git a/requirements.txt b/requirements.txt index cdd2c65a..16a1273c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,7 +5,7 @@ git+https://github.com/inducer/modepy git+https://github.com/pyopencl/pyopencl git+https://github.com/inducer/islpy git+https://github.com/inducer/loopy -git+https://github.com/inducer/boxtree +git+https://github.com/inducer/boxtree@two-away git+https://github.com/inducer/meshmode git+https://github.com/inducer/sumpy git+https://github.com/inducer/pyfmmlib -- GitLab From 676dc67c5662add38408b5113a61d21525f63727 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 8 Aug 2017 10:12:09 -0500 Subject: [PATCH 02/15] Fix git address for boxtree two-away branch --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 16a1273c..39a134cc 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,7 +5,7 @@ git+https://github.com/inducer/modepy git+https://github.com/pyopencl/pyopencl git+https://github.com/inducer/islpy git+https://github.com/inducer/loopy -git+https://github.com/inducer/boxtree@two-away +git+https://gitlab.tiker.net/inducer/boxtree@two-away git+https://github.com/inducer/meshmode git+https://github.com/inducer/sumpy git+https://github.com/inducer/pyfmmlib -- GitLab From 5bd32c79c4f007385ce30c83bf577e7a887b74d6 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 8 Aug 2017 11:04:03 -0500 Subject: [PATCH 03/15] Fix up another from_ interaction list rename fallout bit --- pytential/qbx/fmmlib.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py index bbb37ab0..6965db33 100644 --- a/pytential/qbx/fmmlib.py +++ b/pytential/qbx/fmmlib.py @@ -427,7 +427,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler): mploc = self.get_translation_routine("%ddmploc", vec_suffix="_imany") for isrc_level, ssn in enumerate( - geo_data.traversal().sep_smaller_by_level): + geo_data.traversal().from_sep_smaller_by_level): source_level_start_ibox, source_mpoles_view = \ self.multipole_expansions_view(multipole_exps, isrc_level) -- GitLab From f6e98e99b14863ea15d33e2abd2d32166630a12e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 9 Aug 2017 14:33:54 -0500 Subject: [PATCH 04/15] Use correct boxtree branch from Conda, too --- .test-conda-env-py3-requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.test-conda-env-py3-requirements.txt b/.test-conda-env-py3-requirements.txt index 7cefac68..1acfd2ba 100644 --- a/.test-conda-env-py3-requirements.txt +++ b/.test-conda-env-py3-requirements.txt @@ -1,4 +1,4 @@ -git+https://github.com/inducer/boxtree +git+https://gitlab.tiker.net/inducer/boxtree@two-away git+https://github.com/inducer/pymbolic git+https://github.com/inducer/loopy git+https://github.com/inducer/sumpy -- GitLab From 228514965a713f554c180b7a3d269f41a13ef9ff Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 9 Aug 2017 16:33:18 -0500 Subject: [PATCH 05/15] Adjust stick-out factor for two-away --- test/test_layer_pot_identity.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index af073cc3..f0507aca 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -155,6 +155,20 @@ class StaticTestCase(object): pass +class StarfishGreenTest(StaticTestCase): + expr = GreenExpr() + geometry = StarfishGeometry() + k = 0 + qbx_order = 3 + fmm_order = 6 + + resolutions = [3000, 4000] + + _expansion_stick_out_factor = 0.5 + + fmm_backend = "sumpy" + + class SphereGreenTest(StaticTestCase): expr = GreenExpr() geometry = SphereGeometry() @@ -164,7 +178,7 @@ class SphereGreenTest(StaticTestCase): resolutions = [0, 1] - _expansion_stick_out_factor = 0.75 + _expansion_stick_out_factor = 0.5 fmm_backend = "fmmlib" -- GitLab From 3829c3cb2f7c787a94921ca1d08f5aa7d74e3623 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 9 Aug 2017 16:33:46 -0500 Subject: [PATCH 06/15] Make n-away externally controllabe from QBX source constructor --- pytential/qbx/__init__.py | 8 +++++--- pytential/qbx/geometry.py | 5 +++-- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 50fb8e58..0357a8e7 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -75,8 +75,9 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): # FIXME default debug=False once everything works debug=True, _refined_for_global_qbx=False, - _expansions_in_tree_have_extent=False, - _expansion_stick_out_factor=0, + _expansions_in_tree_have_extent=True, + _expansion_stick_out_factor=0.5, + _well_sep_is_n_away=2, performance_data_file=None, fmm_backend="sumpy", target_stick_out_factor=_not_provided): @@ -158,6 +159,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): self._expansions_in_tree_have_extent = \ _expansions_in_tree_have_extent self._expansion_stick_out_factor = _expansion_stick_out_factor + self._well_sep_is_n_away = _well_sep_is_n_away self.performance_data_file = performance_data_file def copy( @@ -436,7 +438,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): def qbx_fmm_code_getter(self): from pytential.qbx.geometry import QBXFMMGeometryCodeGetter return QBXFMMGeometryCodeGetter(self.cl_context, self.ambient_dim, - debug=self.debug) + debug=self.debug, _well_sep_is_n_away=self._well_sep_is_n_away) # {{{ fmm-based execution diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 6a0f0b08..3f560872 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -104,10 +104,11 @@ class target_state(Enum): # noqa class QBXFMMGeometryCodeGetter(object): - def __init__(self, cl_context, ambient_dim, debug): + def __init__(self, cl_context, ambient_dim, debug, _well_sep_is_n_away): self.cl_context = cl_context self.ambient_dim = ambient_dim self.debug = debug + self._well_sep_is_n_away = _well_sep_is_n_away @memoize_method def copy_targets_kernel(self): @@ -138,7 +139,7 @@ class QBXFMMGeometryCodeGetter(object): @memoize_method def build_traversal(self): from boxtree.traversal import FMMTraversalBuilder - return FMMTraversalBuilder(self.cl_context, well_sep_is_n_away=2) + return FMMTraversalBuilder(self.cl_context, well_sep_is_n_away=self._well_sep_is_n_away) @memoize_method def qbx_center_to_target_box_lookup(self, particle_id_dtype, box_id_dtype): -- GitLab From 74e77b2f0b69dc77a93d841cbbd0fe8c579d5436 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 9 Aug 2017 17:45:37 -0500 Subject: [PATCH 07/15] Fix overly long line to placate flake8 --- pytential/qbx/geometry.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index 3f560872..84d17c55 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -139,7 +139,9 @@ class QBXFMMGeometryCodeGetter(object): @memoize_method def build_traversal(self): from boxtree.traversal import FMMTraversalBuilder - return FMMTraversalBuilder(self.cl_context, well_sep_is_n_away=self._well_sep_is_n_away) + return FMMTraversalBuilder( + self.cl_context, + well_sep_is_n_away=self._well_sep_is_n_away) @memoize_method def qbx_center_to_target_box_lookup(self, particle_id_dtype, box_id_dtype): -- GitLab From a7097558a3744baeb54d4c1c7f5d430010ec8c7b Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 10 Aug 2017 15:51:01 -0500 Subject: [PATCH 08/15] Reset out-of-whack manual starfish green test parameters [ci skip] --- test/test_layer_pot_identity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index f0507aca..b58c4344 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -162,7 +162,7 @@ class StarfishGreenTest(StaticTestCase): qbx_order = 3 fmm_order = 6 - resolutions = [3000, 4000] + resolutions = [30, 50, 70] _expansion_stick_out_factor = 0.5 -- GitLab From ee49521aa43db2323cf3704ea9b04dfb5a7bf703 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 11 Aug 2017 11:59:34 -0500 Subject: [PATCH 09/15] Make starfish test in layer pot identity test use new n-armed starfish --- test/test_layer_pot_identity.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index b58c4344..67f8ffc3 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -34,6 +34,7 @@ from pyopencl.tools import ( # noqa from functools import partial from meshmode.mesh.generation import ( # noqa ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut, WobblyCircle, + NArmedStarfish, make_curve_mesh) # from sumpy.visualization import FieldPlotter from pytential import bind, sym, norm @@ -76,15 +77,25 @@ def get_sphere_mesh(refinement_increment, target_order): class StarfishGeometry(object): - mesh_name = "starfish" + def __init__(self, n_arms=5, amplitude=0.25): + self.n_arms = n_arms + self.amplitude = amplitude + + @property + def mesh_name(self): + return "%d-starfish-%s" % ( + self.n_arms, + self.amplitude) + 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) + return make_curve_mesh( + NArmedStarfish(self.n_arms, self.amplitude), + np.linspace(0, 1, nelements+1), + target_order) class SphereGeometry(object): -- GitLab From 04253e81f31033124dac803873ef8a99adf8c739 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 11 Aug 2017 12:00:16 -0500 Subject: [PATCH 10/15] Layer pot identity test: Add WobblyCircleGreenTest --- test/test_layer_pot_identity.py | 32 +++++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index 67f8ffc3..f9ea50c2 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -55,13 +55,6 @@ 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) @@ -98,6 +91,19 @@ class StarfishGeometry(object): target_order) +class WobblyCircleGeometry(object): + dim = 2 + mesh_name = "wobbly-circle" + + resolutions = [2000, 3000, 4000] + + def get_mesh(self, resolution, target_order): + return make_curve_mesh( + WobblyCircle.random(30, seed=30), + np.linspace(0, 1, resolution+1), + target_order) + + class SphereGeometry(object): mesh_name = "sphere" dim = 3 @@ -180,6 +186,18 @@ class StarfishGreenTest(StaticTestCase): fmm_backend = "sumpy" +class WobblyCircleGreenTest(StaticTestCase): + expr = GreenExpr() + geometry = WobblyCircleGeometry() + k = 0 + qbx_order = 3 + fmm_order = 10 + + _expansion_stick_out_factor = 0.5 + + fmm_backend = "sumpy" + + class SphereGreenTest(StaticTestCase): expr = GreenExpr() geometry = SphereGeometry() -- GitLab From 333437be27d1bb61153f1f769b2892f14f9aba9e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 11 Aug 2017 12:00:57 -0500 Subject: [PATCH 11/15] Layer pot identity test: fix plotting in <3D [ci skip] --- test/test_layer_pot_identity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py index f9ea50c2..e7b50f2e 100644 --- a/test/test_layer_pot_identity.py +++ b/test/test_layer_pot_identity.py @@ -377,7 +377,7 @@ def test_identity_convergence(ctx_getter, case, visualize=False): 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)\ + bdry_normals = bind(density_discr, sym.normal(mesh.ambient_dim))(queue)\ .as_vector(dtype=object) bdry_vis.write_vtk_file("source-%s.vtu" % resolution, [ -- GitLab From 4ebd54d9aec7842b323e6cc84f1c1ebc3fa14068 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 14 Aug 2017 18:24:34 -0500 Subject: [PATCH 12/15] Make max_leaf_refine_weight externally settable on the LP source, fix external setting of well_sep_is_n_away --- pytential/qbx/__init__.py | 12 +++++++++++- pytential/qbx/geometry.py | 13 +------------ 2 files changed, 12 insertions(+), 13 deletions(-) diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py index 8e6fd1ce..ee716f53 100644 --- a/pytential/qbx/__init__.py +++ b/pytential/qbx/__init__.py @@ -78,6 +78,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): _expansions_in_tree_have_extent=True, _expansion_stick_out_factor=0.5, _well_sep_is_n_away=2, + _max_leaf_refine_weight=32, geometry_data_inspector=None, fmm_backend="sumpy", target_stick_out_factor=_not_provided): @@ -160,8 +161,14 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): _expansions_in_tree_have_extent self._expansion_stick_out_factor = _expansion_stick_out_factor self._well_sep_is_n_away = _well_sep_is_n_away + self._max_leaf_refine_weight = _max_leaf_refine_weight self.geometry_data_inspector = geometry_data_inspector + # /!\ *All* parameters set here must also be set by copy() below, + # otherwise they will be reset to their default values behind your + # back if the layer potential source is ever copied. (such as + # during refinement) + def copy( self, density_discr=None, @@ -233,9 +240,12 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase): _expansion_stick_out_factor if _expansion_stick_out_factor is not _not_provided else self._expansion_stick_out_factor), + _well_sep_is_n_away=self._well_sep_is_n_away, + _max_leaf_refine_weight=self._max_leaf_refine_weight, geometry_data_inspector=( geometry_data_inspector or self.geometry_data_inspector), - fmm_backend=self.fmm_backend) + fmm_backend=self.fmm_backend, + ) # }}} diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py index b9ffb88f..4138fbb6 100644 --- a/pytential/qbx/geometry.py +++ b/pytential/qbx/geometry.py @@ -496,22 +496,11 @@ class QBXFMMGeometryData(object): refine_weights[:nsources] = 1 refine_weights.finish() - # NOTE: max_leaf_refine_weight has an impact on accuracy. - # For instance, if a panel contains 64*4 = 256 nodes, then - # a box with 128 sources will contain at most half a - # panel, meaning that its width will be on the order h/2, - # which means many QBX disks (diameter h) will be forced - # to cross boxes. - # So we set max_leaf_refine weight comfortably large - # to avoid having too many disks overlap more than one box. - # - # FIXME: Should investigate this further. - tree, _ = code_getter.build_tree(queue, particles=lpot_src.quad_stage2_density_discr.nodes(), targets=target_info.targets, target_radii=target_radii, - max_leaf_refine_weight=32, + max_leaf_refine_weight=lpot_src._max_leaf_refine_weight, refine_weights=refine_weights, debug=self.debug, stick_out_factor=lpot_src._expansion_stick_out_factor, -- GitLab From 52dece6869133380c85ab6b850809aea708cb0d2 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 15 Aug 2017 19:15:27 -0500 Subject: [PATCH 13/15] Add NaN reproducer [ci skip] --- examples/fmm-manyarm-starfish-nan.py | 109 +++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 examples/fmm-manyarm-starfish-nan.py diff --git a/examples/fmm-manyarm-starfish-nan.py b/examples/fmm-manyarm-starfish-nan.py new file mode 100644 index 00000000..7a8b1d3f --- /dev/null +++ b/examples/fmm-manyarm-starfish-nan.py @@ -0,0 +1,109 @@ +import numpy as np +import numpy.linalg as la # noqa +import pyopencl as cl +import pyopencl.clmath # noqa + +from meshmode.mesh.generation import ( # noqa + ellipse, cloverleaf, NArmedStarfish, drop, n_gon, qbx_peanut, + make_curve_mesh, starfish) +from pytential import bind, sym, norm + +import logging +logger = logging.getLogger(__name__) + + +def get_starfish_mesh(nelements, target_order): + return make_curve_mesh( + NArmedStarfish(20, 0.8), + np.linspace(0, 1, nelements+1), + target_order) + + +WITH_EXTENTS = True +EXPANSION_STICKOUT_FACTOR = 0.5 + + +def get_green_error( + queue, mesh_getter, nelements, fmm_order, qbx_order, k=0): + + target_order = 12 + + mesh = mesh_getter(nelements, target_order) + + d = mesh.ambient_dim + + # u_sym = sym.var("u") + dn_u_sym = sym.var("dn_u") + + from sumpy.kernel import LaplaceKernel + k_sym = LaplaceKernel(d) + zero_op = ( + sym.S(k_sym, dn_u_sym, qbx_forced_limit=-1) + # - sym.D(k_sym, u_sym, qbx_forced_limit="avg") + # - 0.5*u_sym + ) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory) + pre_density_discr = Discretization( + queue.context, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + from pytential.qbx import QBXLayerPotentialSource + lpot_source = QBXLayerPotentialSource( + pre_density_discr, 4*target_order, + qbx_order, fmm_order=fmm_order, + _expansions_in_tree_have_extent=WITH_EXTENTS, + _expansion_stick_out_factor=EXPANSION_STICKOUT_FACTOR, + ) + + lpot_source, _ = lpot_source.with_refinement() + + density_discr = lpot_source.density_discr + + # {{{ compute values of a solution to the PDE + + nodes_host = density_discr.nodes().get(queue) + normal = bind(density_discr, sym.normal(d))(queue).as_vector(np.object) + normal_host = [normal[j].get() for j in range(d)] + + center = np.array([3, 1, 2])[:d] + diff = nodes_host - center[:, np.newaxis] + dist_squared = np.sum(diff**2, axis=0) + dist = np.sqrt(dist_squared) + if d == 2: + u = np.log(dist) + grad_u = diff/dist_squared + elif d == 3: + u = 1/dist + grad_u = -diff/dist**3 + else: + assert False + + dn_u = 0 + for i in range(d): + dn_u = dn_u + normal_host[i]*grad_u[i] + + # }}} + + u_dev = cl.array.to_device(queue, u) + dn_u_dev = cl.array.to_device(queue, dn_u) + + bound_op = bind(lpot_source, zero_op) + error = bound_op( + queue, u=u_dev, dn_u=dn_u_dev, k=k) + + print(len(error), np.where(np.isnan(error.get()))) + return norm(density_discr, queue, error) + + +def main(): + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + + get_green_error(queue, get_starfish_mesh, 500, 20, 2, k=0) + + +if __name__ == "__main__": + main() -- GitLab From aeea38633849449b01db8d60c8e8d662abde7d65 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 19 Aug 2017 15:01:00 -0500 Subject: [PATCH 14/15] Track move of pyopencl back to inducer namespace --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 39a134cc..346d30f2 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,7 +2,7 @@ numpy git+git://github.com/inducer/pymbolic sympy==1.0 git+https://github.com/inducer/modepy -git+https://github.com/pyopencl/pyopencl +git+https://github.com/inducer/pyopencl git+https://github.com/inducer/islpy git+https://github.com/inducer/loopy git+https://gitlab.tiker.net/inducer/boxtree@two-away -- GitLab From 4079cd47431a97c9bf49a159a7c113a4c4eb0b7d Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 27 Aug 2017 22:45:48 -0500 Subject: [PATCH 15/15] Reset boxtree requirements to master branch --- .test-conda-env-py3-requirements.txt | 2 +- requirements.txt | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.test-conda-env-py3-requirements.txt b/.test-conda-env-py3-requirements.txt index 1acfd2ba..e3fac95d 100644 --- a/.test-conda-env-py3-requirements.txt +++ b/.test-conda-env-py3-requirements.txt @@ -1,4 +1,4 @@ -git+https://gitlab.tiker.net/inducer/boxtree@two-away +git+https://gitlab.tiker.net/inducer/boxtree git+https://github.com/inducer/pymbolic git+https://github.com/inducer/loopy git+https://github.com/inducer/sumpy diff --git a/requirements.txt b/requirements.txt index 346d30f2..7360a7c7 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,7 +5,7 @@ git+https://github.com/inducer/modepy git+https://github.com/inducer/pyopencl git+https://github.com/inducer/islpy git+https://github.com/inducer/loopy -git+https://gitlab.tiker.net/inducer/boxtree@two-away +git+https://gitlab.tiker.net/inducer/boxtree git+https://github.com/inducer/meshmode git+https://github.com/inducer/sumpy git+https://github.com/inducer/pyfmmlib -- GitLab