diff --git a/boxtree/fmm.py b/boxtree/fmm.py index 2d18a7baf4ebc982b9973200c5d4a356a48d97a5..63092716048b8b0caaa22f35788d1d6678b69790 100644 --- a/boxtree/fmm.py +++ b/boxtree/fmm.py @@ -185,6 +185,144 @@ def drive_fmm(traversal, expansion_wrangler, src_weights): return result +def drive_fmm_local_quad(intersection, intersection_starts, area_l1_diff, + area_l1_diff_starts, tree, traversal, expansion_wrangler, + src_weights): + """Top-level driver routine for a fast multipole calculation. + + In part, this is intended as a template for custom FMMs, in the sense that + you may copy and paste its + `source code `_ + as a starting point. + + Nonetheless, many common applications (such as point-to-point FMMs) can be + covered by supplying the right *expansion_wrangler* to this routine. + + :arg traversal: A :class:`boxtree.traversal.FMMTraversalInfo` instance. + :arg expansion_wrangler: An object exhibiting the + :class:`ExpansionWranglerInterface`. + :arg src_weights: Source 'density/weights/charges'. + Passed unmodified to *expansion_wrangler*. + + Returns the potentials computed by *expansion_wrangler*. + """ + wrangler = expansion_wrangler + + # Interface guidelines: Attributes of the tree are assumed to be known + # to the expansion wrangler and should not be passed. + + logger.info("start fmm") + + logger.debug("reorder source weights") + + src_weights_user = src_weights.copy() + src_weights = wrangler.reorder_sources(src_weights) + + # {{{ "Step 2.1:" Construct local multipoles + + logger.debug("construct local multipoles") + mpole_exps = wrangler.form_multipoles( + traversal.level_start_source_box_nrs, + traversal.source_boxes, + src_weights) + + # }}} + + # {{{ "Step 2.2:" Propagate multipoles upward + + logger.debug("propagate multipoles upward") + wrangler.coarsen_multipoles( + traversal.level_start_source_parent_box_nrs, + traversal.source_parent_boxes, + mpole_exps) + + # mpole_exps is called Phi in [1] + + # }}} + + # {{{ "Stage 3MOD:" Direct evaluation from neighbor source boxes ("list 1") + + logger.debug("Perform local quadrature") + + pot_local = wrangler.local_quadrature(intersection, intersection_starts) + pot_diff = wrangler.pot_diff_fmm_local(area_l1_diff, area_l1_diff_starts, + src_weights_user) + + # }}} + + # {{{ "Stage 4:" translate separated siblings' ("list 2") mpoles to local + + logger.debug("translate separated siblings' ('list 2') mpoles to local") + local_exps = wrangler.multipole_to_local( + traversal.level_start_target_or_target_parent_box_nrs, + traversal.target_or_target_parent_boxes, + traversal.sep_siblings_starts, + traversal.sep_siblings_lists, + mpole_exps) + + # local_exps represents both Gamma and Delta in [1] + + # }}} + + # {{{ "Stage 5:" evaluate sep. smaller mpoles ("list 3") at particles + + logger.debug("evaluate sep. smaller mpoles at particles ('list 3 far')") + + # (the point of aiming this stage at particles is specifically to keep its + # contribution *out* of the downward-propagating local expansions) + + potentials = wrangler.eval_multipoles( + traversal.level_start_target_box_nrs, + traversal.target_boxes, + traversal.sep_smaller_by_level, + mpole_exps) + + # these potentials are called beta in [1] + # }}} + + # {{{ "Stage 6:" form locals for separated bigger mpoles ("list 4") + + logger.debug("form locals for separated bigger mpoles ('list 4 far')") + + local_exps = local_exps + wrangler.form_locals( + traversal.level_start_target_or_target_parent_box_nrs, + traversal.target_or_target_parent_boxes, + traversal.sep_bigger_starts, + traversal.sep_bigger_lists, + src_weights) + + # }}} + + # {{{ "Stage 7:" propagate local_exps downward + + logger.debug("propagate local_exps downward") + + wrangler.refine_locals( + traversal.level_start_target_or_target_parent_box_nrs, + traversal.target_or_target_parent_boxes, + local_exps) + + # }}} + + # {{{ "Stage 8:" evaluate locals + + logger.debug("evaluate locals") + + potentials = potentials + wrangler.eval_locals( + traversal.level_start_target_box_nrs, + traversal.target_boxes, + local_exps) + + # }}} + + logger.debug("reorder potentials") + result = wrangler.reorder_potentials(potentials) + result = result + (pot_local - pot_diff) + + logger.info("fmm complete") + return result + + # {{{ expansion wrangler interface class ExpansionWranglerInterface: diff --git a/boxtree/mesh_interaction.py b/boxtree/mesh_interaction.py new file mode 100644 index 0000000000000000000000000000000000000000..1ac062847d4c382010db43a573bc52bc5ad026bf --- /dev/null +++ b/boxtree/mesh_interaction.py @@ -0,0 +1,86 @@ +import numpy as np + + +def mesh_and_target_area_intersection_box(mesh, tgt_areas): + # Returns intersection of box tgt_area and mesh + for i in range(len(mesh.groups)): + #for now, just lump all mesh groups together + elem = np.r_['1', mesh.groups[i].nodes] + #tgt_area consists of two arrays, the upper and lower bounds for each dim + tgt_lb = tgt_areas[0] + tgt_ub = tgt_areas[1] + + na = np.newaxis + intersecting = True + for i in range(mesh.dim): + lb = np.amin(elem[i, :, :], axis=1) + ub = np.amax(elem[i, :, :], axis=1) + intersecting = (intersecting & + (((tgt_lb[:, i, na] < lb) & (lb < tgt_ub[:, i, na])) | + ((tgt_lb[:, i, na] < ub) & (ub < tgt_ub[:, i, na])))) + + intersecting = np.argwhere(intersecting) + intersection = intersecting[:, 1] + intersection_starts = np.argwhere(np.diff(intersecting[:, 0]))[:, 0]+1 + intersection_starts = np.r_[0, intersection_starts, intersection.size] + #returns mesh element numbers, associated with target elem number + return intersection, intersection_starts + + +def target_and_neighbor_sources(tree, trav, queue): + #ordered like target_boxes + bss = tree.box_source_starts.get(queue) + bscc = tree.box_source_counts_cumul.get(queue) + nsbs = trav.neighbor_source_boxes_starts.get(queue) + nsbl = trav.neighbor_source_boxes_lists.get(queue) + + tgt_l1_srcs = np.array([], dtype=np.int) + tgt_l1_srcs_starts = np.array([0], dtype=np.int) + for box in range(trav.target_boxes.size): + for neighbor in range(nsbs[box], nsbs[box+1]): + l1 = nsbl[neighbor] + tgt_l1_srcs = np.append(tgt_l1_srcs, bss[l1] + np.arange(bscc[l1])) + tgt_l1_srcs_starts = np.append(tgt_l1_srcs_starts, tgt_l1_srcs.size) + return tgt_l1_srcs, tgt_l1_srcs_starts + + +def target_area_fmm_neighbors_difference(tree, trav, queue, nunit_nodes, + intersection, intersection_starts): + target_boxes = trav.target_boxes.get(queue) + bsst = tree.box_source_starts.get(queue)[target_boxes] + bscct = tree.box_source_counts_cumul.get(queue)[target_boxes] + user_source_ids = tree.user_source_ids.get(queue) + + tgt_l1_srcs, tgt_l1_srcs_starts = target_and_neighbor_sources(tree, trav, queue) + + #Preprocessing step to determine the space needed for each tgt's 'area_l1_diff' + l1_box_size = np.diff(tgt_l1_srcs_starts) + l1_size = np.zeros(tree.nsources) + area_size = nunit_nodes*np.diff(intersection_starts) + for box in range(trav.target_boxes.size): + l1_size[bsst[box] + np.arange(bscct[box])] = l1_box_size[box] + area_l1_diff_starts = np.append(0, np.cumsum(area_size - l1_size).astype(int)) + + area_l1_diff = np.zeros(area_l1_diff_starts[-1], dtype=np.int) + elem_by_node_num = np.arange(nunit_nodes) + for box in range(trav.target_boxes.size): + l1_set = tgt_l1_srcs[tgt_l1_srcs_starts[box]:tgt_l1_srcs_starts[box+1]] + l1_set = user_source_ids[l1_set] # convert to user source ordering + for tgt in bsst[box] + np.arange(bscct[box]): + #'intersection' gives us the *element nums* that are in the set + #associated with the 'tgt' point (in tree ordering), but we need the + #source nums of sources contained in each of the elements. This is easy + #to get in user ordering, since the user source order ascends in the same + #order as element ordering. + area_set = np.ravel(elem_by_node_num + nunit_nodes * intersection[ + intersection_starts[tgt]:intersection_starts[tgt+1], np.newaxis]) + + area_l1_diff[area_l1_diff_starts[tgt]:area_l1_diff_starts[tgt+1]] =\ + np.setdiff1d(area_set, l1_set, assume_unique=True) + subset = np.setdiff1d(l1_set, area_set, assume_unique=True).size + assert subset == 0, 'Area too small for l1, for tgt num: ' + str(tgt) +\ + ', box: ' + str(box) + '. ' + str(subset) + '/' + str(area_set.size)\ + + ' points not in area' + + #returns, in tree target ordering, the user source ids in the diff for the tgt + return area_l1_diff, area_l1_diff_starts diff --git a/requirements.txt b/requirements.txt index e36a43aed0596e0722e1639ab80a1b1355dd02f6..7ce323ed8b36ebe1089e673f32aec9198141f36d 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,3 +3,5 @@ git+git://github.com/pyopencl/pyopencl git+git://github.com/inducer/islpy git+git://github.com/inducer/loopy git+git://github.com/inducer/pyfmmlib +git+https://github.com/inducer/meshmode.git +git+https://github.com/inducer/sumpy.git diff --git a/test/blob-2d.msh b/test/blob-2d.msh new file mode 100644 index 0000000000000000000000000000000000000000..d9d2f7c74f9a903c7ceaaa4af3ba17a6620793fd --- /dev/null +++ b/test/blob-2d.msh @@ -0,0 +1,446 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +146 +1 -0.4725430905820004 0.1646517366169589 0 +2 -0.4767799430995929 0.1159304601168752 0 +3 -0.4767322166399771 0.06703688281240093 0 +4 -0.4702937700635583 0.01860086453797516 0 +5 -0.4538183776589282 -0.02729702900844069 0 +6 -0.4227369483017478 -0.0645819673749298 0 +7 -0.378438447476652 -0.08434169531943977 0 +8 -0.3297463531825587 -0.08755782429564593 0 +9 -0.2809922952332325 -0.08369971986917121 0 +10 -0.2323017027563748 -0.07902694381423922 0 +11 -0.1834319137022549 -0.07737353247625005 0 +12 -0.1347544746852848 -0.08171236732517156 0 +13 -0.08738693003876633 -0.09373390024415561 0 +14 -0.04169550328514776 -0.1111600381215554 0 +15 0.002930509478649007 -0.1311869920451043 0 +16 0.04743753303695913 -0.1514793872925139 0 +17 0.09323408996811799 -0.1685627310701821 0 +18 0.1414012341904861 -0.1710396922413012 0 +19 0.1732607869944924 -0.1361155792046744 0 +20 0.1831467818070197 -0.08838035357878156 0 +21 0.1831714242748896 -0.03950825402407249 0 +22 0.1771778625795523 0.009016005361491941 0 +23 0.166382378736668 0.05670673201188711 0 +24 0.1509706129502917 0.1031089155841279 0 +25 0.130397032983188 0.147451038094475 0 +26 0.103223173029645 0.1880361965098842 0 +27 0.06689182739503687 0.2204692212154625 0 +28 0.02197249056993197 0.23937220789002 0 +29 -0.02607919267527611 0.2483400416739637 0 +30 -0.07461773801025579 0.2544025145406534 0 +31 -0.1228489382886651 0.26248503488488 0 +32 -0.1703906535398324 0.2739771483445165 0 +33 -0.217649103296364 0.2866029236958204 0 +34 -0.2651766562306084 0.2981559866619584 0 +35 -0.3134270431714214 0.3060437840417993 0 +36 -0.3622390050029438 0.3059949681103042 0 +37 -0.4084081447638059 0.2909021121876902 0 +38 -0.4431589393839811 0.2570920965607918 0 +39 -0.462998838028793 0.2125712775436452 0 +40 -0.2933233764464115 0.1138267430790992 0 +41 -0.02847814951822331 0.06951342048909848 0 +42 -0.1507175526648517 0.1305814351685241 0 +43 0.05671910460881946 -0.0196559247139545 0 +44 -0.208826701719056 0.03246116157417323 0 +45 -0.3556479917501454 0.2018958929113452 0 +46 -0.3549866101923704 0.01732350767456248 0 +47 -0.1074529364327065 0.01319266311729878 0 +48 -0.2267260446319173 0.1934357802301528 0 +49 -0.0610745157004064 0.1600077687266293 0 +50 -0.3883524208181273 0.111183209576144 0 +51 0.05889041499476744 0.0776338687461553 0 +52 -0.03286645100321878 -0.02369398759946384 0 +53 0.02221409155934051 0.1462759232954491 0 +54 0.1017024117139075 -0.08755561321413979 0 +55 -0.2818538361061502 -0.003939447723841805 0 +56 -0.2930198526593978 0.2348511079107199 0 +57 -0.157597301662713 0.1988376537868406 0 +58 -0.09688546668750414 0.08681191013630563 0 +59 -0.2142058556152331 0.110169343185449 0 +60 0.1156514499997027 0.02848427806311571 0 +61 -0.4128986346885033 0.05443782465991698 0 +62 -0.4105147942411019 0.167280714103395 0 +63 0.01444774030469731 -0.05474763972018154 0 +64 -0.1634807173903056 -0.01695751457175493 0 +65 0.02472101487767853 0.02682173602829646 0 +66 -0.2982134796793256 0.1589279384008412 0 +67 -0.1550565803827874 0.0684071791682978 0 +68 -0.2570833955478836 0.04937180456916568 0 +69 -0.3352936721186352 0.07765851422549841 0 +70 -0.3305251100090268 -0.02912650765304481 0 +71 0.1232712319791713 -0.03998024696762692 0 +72 -0.1076089234619594 0.2134617374872078 0 +73 -0.01516041919738775 0.1973989129999587 0 +74 -0.3385876211475933 0.2528280561940174 0 +75 -0.2298161207061463 -0.02377906621915205 0 +76 0.07896208292736451 0.1298042755458195 0 +77 -0.4003774279055767 -0.01332706007358847 0 +78 -0.08435996183459216 -0.03781976387024838 0 +79 -0.3928732020052436 0.2331762588758936 0 +80 -0.05414129333359385 0.02275864500842628 0 +81 -0.3357226788150336 0.1314877635332329 0 +82 -0.02350537513687817 0.1234217083348019 0 +83 -0.1983295584050745 0.2336241273464102 0 +84 0.01304169938750516 0.09946768053857757 0 +85 -0.1977627582096242 0.1488792465547518 0 +86 -0.2499853987658275 0.1446668696974264 0 +87 0.04461051284959683 -0.1014584649369309 0 +88 -0.3086956668329468 0.03722270104251224 0 +89 -0.1055760958322025 0.1255049994839639 0 +90 0.1130605683101332 0.07003350494559785 0 +91 -0.4378396886896197 0.09157841975878761 0 +92 0.1316517540122663 -0.1218607719810022 0 +93 -0.05636558772024425 0.213898693008732 0 +94 0.0569057946782392 0.1670282746351865 0 +95 -0.251903797242819 0.244749207547547 0 +96 -0.04259864331817997 -0.06864874578755542 0 +97 0.07341201162712796 0.02794716069702885 0 +98 -0.1264008156441787 0.1664595192975996 0 +99 0.01109293989880011 -0.0129221658506799 0 +100 -0.1598541394438701 0.02896921228600524 0 +101 -0.2637389593339061 -0.04467477554652227 0 +102 0.1440120687448451 -0.07491671079570661 0 +103 0.0694145640978272 -0.05811944733703039 0 +104 -0.06343932078857029 0.1043590166527328 0 +105 -0.1271358308635849 -0.0395834938425037 0 +106 -0.4376053423594887 0.1328175270655886 0 +107 -0.3683443825353412 0.1605091809389714 0 +108 -0.3699274481945085 -0.04333102157057589 0 +109 0.08372726031616951 -0.1261833936989538 0 +110 -0.4297533132363489 0.01689241587734225 0 +111 -0.3733449446096971 0.06072851240081606 0 +112 -0.2385189574430036 0.01359841853098009 0 +113 -0.2945221907608757 0.2701665757127336 0 +114 0.02456850236467112 0.1911285890934324 0 +115 -0.4233094992106001 0.2071281720893554 0 +116 -0.1195246723542033 0.0567883060803262 0 +117 -0.0123133767320886 0.02402921316156392 0 +118 0.106456142433149 0.1056063205832351 0 +119 -0.2722512494827684 0.1964713698426759 0 +120 -0.3764171768348171 0.2685581126800726 0 +121 -0.07454773947258309 0.0594726504849552 0 +122 -0.1525098712549517 0.2376580258697706 0 +123 -0.1770334590520359 0.1008255546664961 0 +124 -0.3147863335002259 0.2020469892040668 0 +125 0.1413610205651505 0.0005033077372428307 0 +126 0.09390072603782196 -0.01260533241814538 0 +127 -0.201453729815424 0.06870113952150303 0 +128 -0.2973713107729749 -0.04979965501764522 0 +129 -0.2514828412664952 0.09296096822316524 0 +130 -0.1940130779682059 -0.04195144039536036 0 +131 -0.1318637160510444 0.09692995937266657 0 +132 -0.2909730970379276 0.07265145874631751 0 +133 -0.1916138009097272 0.1860347080318873 0 +134 -0.3942721861264993 0.02721104010780984 0 +135 -0.004461076794076924 -0.09344037612226547 0 +136 -0.164671330921184 0.1652013944062097 0 +137 -0.3901379739484865 0.1939980437837921 0 +138 -0.0107330904392276 0.1584094356536605 0 +139 -0.3190153057851235 0.005370063335047026 0 +140 -0.3335300684063315 0.1699261292740086 0 +141 -0.4472878422202458 0.04745912661657253 0 +142 -0.1964707115921482 -0.002085237742636636 0 +143 -0.0716416676821785 -0.005570624691587089 0 +144 0.04333753573999888 0.1182635314722547 0 +145 -0.1594972907009092 -0.05220175449839765 0 +146 0.005791678686025703 0.06169763089370611 0 +$EndNodes +$Elements +291 +1 15 2 0 1 1 +2 1 2 0 1 1 2 +3 1 2 0 1 2 3 +4 1 2 0 1 3 4 +5 1 2 0 1 4 5 +6 1 2 0 1 5 6 +7 1 2 0 1 6 7 +8 1 2 0 1 7 8 +9 1 2 0 1 8 9 +10 1 2 0 1 9 10 +11 1 2 0 1 10 11 +12 1 2 0 1 11 12 +13 1 2 0 1 12 13 +14 1 2 0 1 13 14 +15 1 2 0 1 14 15 +16 1 2 0 1 15 16 +17 1 2 0 1 16 17 +18 1 2 0 1 17 18 +19 1 2 0 1 18 19 +20 1 2 0 1 19 20 +21 1 2 0 1 20 21 +22 1 2 0 1 21 22 +23 1 2 0 1 22 23 +24 1 2 0 1 23 24 +25 1 2 0 1 24 25 +26 1 2 0 1 25 26 +27 1 2 0 1 26 27 +28 1 2 0 1 27 28 +29 1 2 0 1 28 29 +30 1 2 0 1 29 30 +31 1 2 0 1 30 31 +32 1 2 0 1 31 32 +33 1 2 0 1 32 33 +34 1 2 0 1 33 34 +35 1 2 0 1 34 35 +36 1 2 0 1 35 36 +37 1 2 0 1 36 37 +38 1 2 0 1 37 38 +39 1 2 0 1 38 39 +40 1 2 0 1 39 1 +41 2 2 0 1 52 63 96 +42 2 2 0 1 46 70 108 +43 2 2 0 1 19 92 20 +44 2 2 0 1 46 108 77 +45 2 2 0 1 49 98 72 +46 2 2 0 1 20 92 102 +47 2 2 0 1 23 60 90 +48 2 2 0 1 23 125 60 +49 2 2 0 1 1 115 62 +50 2 2 0 1 25 76 26 +51 2 2 0 1 49 89 98 +52 2 2 0 1 13 105 78 +53 2 2 0 1 48 119 95 +54 2 2 0 1 12 105 13 +55 2 2 0 1 31 72 122 +56 2 2 0 1 48 95 83 +57 2 2 0 1 30 72 31 +58 2 2 0 1 51 90 97 +59 2 2 0 1 60 97 90 +60 2 2 0 1 40 81 66 +61 2 2 0 1 50 106 62 +62 2 2 0 1 26 76 94 +63 2 2 0 1 63 135 96 +64 2 2 0 1 54 103 71 +65 2 2 0 1 68 129 127 +66 2 2 0 1 44 68 127 +67 2 2 0 1 1 39 115 +68 2 2 0 1 45 79 74 +69 2 2 0 1 40 66 86 +70 2 2 0 1 57 122 72 +71 2 2 0 1 50 61 91 +72 2 2 0 1 47 105 64 +73 2 2 0 1 40 69 81 +74 2 2 0 1 52 99 63 +75 2 2 0 1 43 63 99 +76 2 2 0 1 35 74 36 +77 2 2 0 1 50 111 61 +78 2 2 0 1 67 100 127 +79 2 2 0 1 47 116 121 +80 2 2 0 1 1 62 106 +81 2 2 0 1 50 81 69 +82 2 2 0 1 28 73 29 +83 2 2 0 1 37 79 38 +84 2 2 0 1 127 100 44 +85 2 2 0 1 47 121 80 +86 2 2 0 1 43 103 63 +87 2 2 0 1 71 103 126 +88 2 2 0 1 50 62 107 +89 2 2 0 1 15 135 87 +90 2 2 0 1 47 64 100 +91 2 2 0 1 51 65 146 +92 2 2 0 1 8 108 70 +93 2 2 0 1 51 76 118 +94 2 2 0 1 5 77 6 +95 2 2 0 1 46 139 70 +96 2 2 0 1 55 88 68 +97 2 2 0 1 43 65 97 +98 2 2 0 1 63 87 135 +99 2 2 0 1 51 97 65 +100 2 2 0 1 43 99 65 +101 2 2 0 1 21 71 125 +102 2 2 0 1 22 125 23 +103 2 2 0 1 59 127 129 +104 2 2 0 1 46 69 88 +105 2 2 0 1 15 87 16 +106 2 2 0 1 41 117 80 +107 2 2 0 1 51 118 90 +108 2 2 0 1 6 77 108 +109 2 2 0 1 50 91 106 +110 2 2 0 1 49 72 93 +111 2 2 0 1 37 120 79 +112 2 2 0 1 49 93 73 +113 2 2 0 1 45 74 124 +114 2 2 0 1 56 124 74 +115 2 2 0 1 21 102 71 +116 2 2 0 1 47 78 105 +117 2 2 0 1 41 146 117 +118 2 2 0 1 66 119 86 +119 2 2 0 1 7 108 8 +120 2 2 0 1 35 113 74 +121 2 2 0 1 53 84 82 +122 2 2 0 1 57 83 122 +123 2 2 0 1 50 69 111 +124 2 2 0 1 55 128 70 +125 2 2 0 1 10 75 130 +126 2 2 0 1 51 146 84 +127 2 2 0 1 56 95 119 +128 2 2 0 1 41 82 84 +129 2 2 0 1 57 133 83 +130 2 2 0 1 30 93 72 +131 2 2 0 1 55 68 112 +132 2 2 0 1 44 112 68 +133 2 2 0 1 76 144 94 +134 2 2 0 1 53 94 144 +135 2 2 0 1 29 73 93 +136 2 2 0 1 48 86 119 +137 2 2 0 1 17 92 18 +138 2 2 0 1 38 79 115 +139 2 2 0 1 32 83 33 +140 2 2 0 1 23 90 24 +141 2 2 0 1 21 125 22 +142 2 2 0 1 28 114 73 +143 2 2 0 1 59 86 85 +144 2 2 0 1 2 91 3 +145 2 2 0 1 54 71 102 +146 2 2 0 1 10 101 75 +147 2 2 0 1 57 72 98 +148 2 2 0 1 56 74 113 +149 2 2 0 1 6 108 7 +150 2 2 0 1 46 111 69 +151 2 2 0 1 8 70 128 +152 2 2 0 1 40 132 69 +153 2 2 0 1 48 85 86 +154 2 2 0 1 55 70 139 +155 2 2 0 1 41 104 82 +156 2 2 0 1 49 73 138 +157 2 2 0 1 53 138 114 +158 2 2 0 1 5 110 77 +159 2 2 0 1 55 75 101 +160 2 2 0 1 47 143 78 +161 2 2 0 1 13 78 96 +162 2 2 0 1 36 74 120 +163 2 2 0 1 51 144 76 +164 2 2 0 1 46 77 134 +165 2 2 0 1 52 96 78 +166 2 2 0 1 55 112 75 +167 2 2 0 1 3 91 141 +168 2 2 0 1 73 114 138 +169 2 2 0 1 54 102 92 +170 2 2 0 1 52 78 143 +171 2 2 0 1 44 100 142 +172 2 2 0 1 45 137 79 +173 2 2 0 1 50 107 81 +174 2 2 0 1 25 118 76 +175 2 2 0 1 41 121 104 +176 2 2 0 1 20 102 21 +177 2 2 0 1 10 130 11 +178 2 2 0 1 52 143 80 +179 2 2 0 1 32 122 83 +180 2 2 0 1 38 115 39 +181 2 2 0 1 49 82 104 +182 2 2 0 1 63 103 87 +183 2 2 0 1 74 79 120 +184 2 2 0 1 58 121 116 +185 2 2 0 1 18 92 19 +186 2 2 0 1 33 83 95 +187 2 2 0 1 48 83 133 +188 2 2 0 1 17 109 92 +189 2 2 0 1 42 123 85 +190 2 2 0 1 59 85 123 +191 2 2 0 1 47 80 143 +192 2 2 0 1 43 126 103 +193 2 2 0 1 61 141 91 +194 2 2 0 1 33 95 34 +195 2 2 0 1 59 129 86 +196 2 2 0 1 49 138 82 +197 2 2 0 1 29 93 30 +198 2 2 0 1 66 81 140 +199 2 2 0 1 49 104 89 +200 2 2 0 1 52 80 117 +201 2 2 0 1 41 80 121 +202 2 2 0 1 40 86 129 +203 2 2 0 1 53 114 94 +204 2 2 0 1 13 96 14 +205 2 2 0 1 65 117 146 +206 2 2 0 1 51 84 144 +207 2 2 0 1 26 94 27 +208 2 2 0 1 55 101 128 +209 2 2 0 1 9 101 10 +210 2 2 0 1 4 110 5 +211 2 2 0 1 53 144 84 +212 2 2 0 1 34 95 113 +213 2 2 0 1 64 142 100 +214 2 2 0 1 42 98 89 +215 2 2 0 1 16 87 109 +216 2 2 0 1 16 109 17 +217 2 2 0 1 54 87 103 +218 2 2 0 1 56 113 95 +219 2 2 0 1 53 82 138 +220 2 2 0 1 1 106 2 +221 2 2 0 1 54 109 87 +222 2 2 0 1 2 106 91 +223 2 2 0 1 43 97 126 +224 2 2 0 1 58 89 104 +225 2 2 0 1 27 114 28 +226 2 2 0 1 14 96 135 +227 2 2 0 1 55 139 88 +228 2 2 0 1 8 128 9 +229 2 2 0 1 58 131 89 +230 2 2 0 1 42 89 131 +231 2 2 0 1 68 88 132 +232 2 2 0 1 60 126 97 +233 2 2 0 1 58 104 121 +234 2 2 0 1 75 112 142 +235 2 2 0 1 69 132 88 +236 2 2 0 1 41 84 146 +237 2 2 0 1 54 92 109 +238 2 2 0 1 44 142 112 +239 2 2 0 1 27 94 114 +240 2 2 0 1 48 133 85 +241 2 2 0 1 42 85 136 +242 2 2 0 1 34 113 35 +243 2 2 0 1 24 118 25 +244 2 2 0 1 36 120 37 +245 2 2 0 1 24 90 118 +246 2 2 0 1 31 122 32 +247 2 2 0 1 46 88 139 +248 2 2 0 1 67 131 116 +249 2 2 0 1 14 135 15 +250 2 2 0 1 11 145 12 +251 2 2 0 1 71 126 125 +252 2 2 0 1 60 125 126 +253 2 2 0 1 58 116 131 +254 2 2 0 1 3 141 4 +255 2 2 0 1 56 119 124 +256 2 2 0 1 66 124 119 +257 2 2 0 1 67 116 100 +258 2 2 0 1 52 117 99 +259 2 2 0 1 47 100 116 +260 2 2 0 1 65 99 117 +261 2 2 0 1 45 107 137 +262 2 2 0 1 62 137 107 +263 2 2 0 1 67 123 131 +264 2 2 0 1 42 131 123 +265 2 2 0 1 68 132 129 +266 2 2 0 1 40 129 132 +267 2 2 0 1 46 134 111 +268 2 2 0 1 61 111 134 +269 2 2 0 1 9 128 101 +270 2 2 0 1 57 98 136 +271 2 2 0 1 42 136 98 +272 2 2 0 1 12 145 105 +273 2 2 0 1 45 140 107 +274 2 2 0 1 81 107 140 +275 2 2 0 1 67 127 123 +276 2 2 0 1 77 110 134 +277 2 2 0 1 61 134 110 +278 2 2 0 1 64 105 145 +279 2 2 0 1 4 141 110 +280 2 2 0 1 61 110 141 +281 2 2 0 1 79 137 115 +282 2 2 0 1 62 115 137 +283 2 2 0 1 59 123 127 +284 2 2 0 1 45 124 140 +285 2 2 0 1 66 140 124 +286 2 2 0 1 64 130 142 +287 2 2 0 1 75 142 130 +288 2 2 0 1 57 136 133 +289 2 2 0 1 85 133 136 +290 2 2 0 1 64 145 130 +291 2 2 0 1 11 130 145 +$EndElements diff --git a/test/test_mesh_and_target_area_intersection.pckl b/test/test_mesh_and_target_area_intersection.pckl new file mode 100644 index 0000000000000000000000000000000000000000..97d62cd80bbd37f0e194eb35da0fbc0b66e9afe0 Binary files /dev/null and b/test/test_mesh_and_target_area_intersection.pckl differ diff --git a/test/test_mesh_interaction.py b/test/test_mesh_interaction.py new file mode 100644 index 0000000000000000000000000000000000000000..7a056a3384f128eec30513b74174d31db65458c4 --- /dev/null +++ b/test/test_mesh_interaction.py @@ -0,0 +1,309 @@ +import numpy as np +import numpy.linalg as la +import pyopencl as cl +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) +from sumpy.kernel import one_kernel_2d +from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansion +from sumpy.expansion.local import VolumeTaylorLocalExpansion + +import pytest + +import logging +logger = logging.getLogger(__name__) + +try: + import faulthandler +except ImportError: + pass +else: + faulthandler.enable() + + +# See comment on SumpyExpansionWranglerLocalQuad below +from test_fmm import ConstantOneExpansionWrangler + + +class LocalQuadratureConstantOneExpansionWrangler(ConstantOneExpansionWrangler): + def __init__(self, tree, queue): + self.tree = tree + self.queue = queue + + def local_quadrature(self, intersection, intersection_starts): + return 2*np.diff(intersection_starts)[self.tree.sorted_target_ids] + + def pot_diff_fmm_local(self, area_l1_diff, area_l1_diff_starts, weights): + result = np.zeros(area_l1_diff_starts.size-1) + for i in range(area_l1_diff_starts.size-1): + remove = area_l1_diff[area_l1_diff_starts[i]:area_l1_diff_starts[i+1]] + result[i] = np.sum(weights[remove]) + return result[self.tree.sorted_target_ids] + + +# Modified code container that redirects to our modified local quad wrangler +from sumpy.fmm import SumpyExpansionWranglerCodeContainer + + +class SumpyExpansionWranglerLQCodeContainer(SumpyExpansionWranglerCodeContainer): + def get_wrangler(self, queue, tree, dtype, fmm_level_to_order, + source_extra_kwargs={}, + kernel_extra_kwargs=None): + return SumpyExpansionWranglerLocalQuad(self, queue, tree, dtype, + fmm_level_to_order, source_extra_kwargs, kernel_extra_kwargs) + + +# Local quad FMM driver expects two new funcs: one for local quadrature, and one for +# pointwise evaluation of sources in the diff set. Info provided to the FMM driver +# and therefore available to the wrangler are CSR-like lists (for each target point) +# of: 1) elements in the mesh/area intersection and 2) source ids in the diff set +from sumpy.fmm import SumpyExpansionWrangler + + +class SumpyExpansionWranglerLocalQuad(SumpyExpansionWrangler): + def local_quadrature(self, intersection, intersection_starts): + # drive_fmm expects potentials as device object arrays, convert before return + sorted_target_ids = self.tree.sorted_target_ids.get(self.queue) + result = 2*np.diff(intersection_starts)[sorted_target_ids] + result = cl.array.to_device(self.queue, result) + from pytools.obj_array import make_obj_array + result = make_obj_array([result]) + return result + + def pot_diff_fmm_local(self, area_l1_diff, area_l1_diff_starts, weights): + result = np.zeros(area_l1_diff_starts.size-1) + for i in range(area_l1_diff_starts.size-1): + remove = area_l1_diff[area_l1_diff_starts[i]:area_l1_diff_starts[i+1]] + result[i] = np.sum(weights.get(self.queue)[remove]) + # drive_fmm expects potentials as device object arrays, convert before return + sorted_target_ids = self.tree.sorted_target_ids.get(self.queue) + result = cl.array.to_device(self.queue, result[sorted_target_ids]) + from pytools.obj_array import make_obj_array + result = make_obj_array([result]) + return result + + +# Define a test area that could be used: a region over which local quadrature is +# necessary. This example has the quad area proportional to the local mesh size, a +# decent choice for locally refined meshes to attempt to bound the number of mesh +# elements intersecting the quad area (likely need more QBX expansion terms to work) +def local_quadrature_area(mesh, discr, user_source_ids, nodes, k): + for i in range(len(mesh.groups)): + # for now, just lump all mesh groups together + elem = np.r_['1', mesh.groups[i].nodes] + + nunit_nodes = discr.groups[0].nunit_nodes + # nodes = discr.nodes().get(queue) + tgt_lb = np.zeros((discr.nnodes, mesh.dim)) + tgt_ub = tgt_lb.copy() + for i in range(mesh.dim): + lb = np.amin(elem[i, :, :], axis=1) + ub = np.amax(elem[i, :, :], axis=1) + h = np.repeat(ub-lb, nunit_nodes) + + tgt_lb[:, i] = nodes[i, :]-k*h + tgt_ub[:, i] = nodes[i, :]+k*h + + # reorder from user to tree ordering, what mesh_interaction expects + # is user_source_ids correct or is inverse of sorted_target_ids better? + tgt_lb = tgt_lb[user_source_ids, :] + tgt_ub = tgt_ub[user_source_ids, :] + + return [tgt_lb, tgt_ub] + + +def test_mesh_and_target_area_intersection(ctx_getter): + discr_order = 3 + + ctx = ctx_getter() + queue = cl.CommandQueue(ctx) + + from meshmode.mesh.io import read_gmsh + mesh = read_gmsh("blob-2d.msh", force_ambient_dim=2) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory) + discr = Discretization(ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(discr_order)) + + user_source_ids = np.arange(discr.nnodes) + dnodes = discr.nodes().get(queue) + tgt_areas = local_quadrature_area(mesh, discr, user_source_ids, dnodes, k=1.7) + from boxtree.mesh_interaction import mesh_and_target_area_intersection_box + intersection, intersection_starts = mesh_and_target_area_intersection_box(mesh, + tgt_areas) + + import pickle + f = open('test_mesh_and_target_area_intersection.pckl', 'rb') + test_intersection, test_intersection_starts = pickle.load(f) + f.close() + + assert(all(test_intersection == intersection)) + assert(all(test_intersection_starts == intersection_starts)) + + +def test_constant_one_fmm(ctx_getter): + discr_order = 3 + + ctx = ctx_getter() + queue = cl.CommandQueue(ctx) + + from meshmode.mesh.io import read_gmsh + mesh = read_gmsh("blob-2d.msh", force_ambient_dim=2) + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory) + discr = Discretization( + ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(discr_order)) + + from boxtree import TreeBuilder + tb = TreeBuilder(ctx) + + from pytools.obj_array import make_obj_array + sources = make_obj_array([ # noqa + discr.nodes()[i, ].with_queue(queue).copy() + for i in range(2)]) + nsources = len(sources[0]) + nunit_nodes = discr.groups[0].nunit_nodes + + tree, _ = tb(queue, sources, max_particles_in_box=nunit_nodes) + + from boxtree.traversal import FMMTraversalBuilder + tbuild = FMMTraversalBuilder(ctx) + trav, _ = tbuild(queue, tree, debug=True) + if trav.sep_close_smaller_starts is not None: + trav = trav.merge_close_lists(queue) + + weights = discr.quad_weights(queue).get(queue) + weights_sum = np.sum(weights) + + host_trav = trav.get(queue=queue) + host_tree = host_trav.tree + + wrangler = LocalQuadratureConstantOneExpansionWrangler(host_tree, queue) + + user_source_ids = host_tree.user_source_ids + dnodes = discr.nodes().get(queue) + tgt_areas = local_quadrature_area(mesh, discr, user_source_ids, dnodes, k=1.7) + from boxtree.mesh_interaction import mesh_and_target_area_intersection_box + intersection, intersection_starts = mesh_and_target_area_intersection_box( + mesh, tgt_areas) + + from boxtree.mesh_interaction import target_area_fmm_neighbors_difference + area_l1_diff, area_l1_diff_starts =\ + target_area_fmm_neighbors_difference(tree, trav, queue, nunit_nodes, + intersection, intersection_starts) + + from boxtree.fmm import drive_fmm_local_quad + pot = drive_fmm_local_quad(intersection, intersection_starts, area_l1_diff, + area_l1_diff_starts, tree, host_trav, wrangler, + weights) + + rel_err = la.norm((pot - weights_sum) / nsources) + good = rel_err < 1e-8 + + assert good + + +@pytest.mark.parametrize("knl, local_expn_class, mpole_expn_class", [ + (one_kernel_2d, VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion), + ]) +def test_sumpy_fmm(ctx_getter, knl, local_expn_class, mpole_expn_class): + logging.basicConfig(level=logging.INFO) + + ctx = ctx_getter() + queue = cl.CommandQueue(ctx) + + from meshmode.mesh.io import read_gmsh + mesh = read_gmsh("blob-2d.msh", force_ambient_dim=2) + + discr_order = 3 + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory) + discr = Discretization( + ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(discr_order)) + + from boxtree import TreeBuilder + tb = TreeBuilder(ctx) + + from pytools.obj_array import make_obj_array + sources = make_obj_array([ # noqa + discr.nodes()[i, ].with_queue(queue).copy() + for i in range(2)]) + nunit_nodes = discr.groups[0].nunit_nodes + + tree, _ = tb(queue, sources, max_particles_in_box=nunit_nodes, debug=True) + + from boxtree.traversal import FMMTraversalBuilder + tbuild = FMMTraversalBuilder(ctx) + trav, _ = tbuild(queue, tree, debug=True) + + weights = discr.quad_weights(queue).get(queue) + weights = cl.array.to_device(queue, weights) + + logger.info("computing direct (reference) result") + + extra_kwargs = {} + dtype = np.float64 + order_values = [1, 2, 3] + + # Mesh-interaction variables that don't change if source locs don't + # We need to only calc these once, each new FMM won't need these updated + tgt_areas = local_quadrature_area(mesh, discr, tree.user_source_ids.get(queue), + discr.nodes().get(queue), k=1.7) + + from boxtree.mesh_interaction import mesh_and_target_area_intersection_box + intersection, intersection_starts = mesh_and_target_area_intersection_box( + mesh, tgt_areas) + + from boxtree.mesh_interaction import target_area_fmm_neighbors_difference + area_l1_diff, area_l1_diff_starts =\ + target_area_fmm_neighbors_difference(tree, trav, queue, nunit_nodes, + intersection, intersection_starts) + + from functools import partial + for order in order_values: + out_kernels = [knl] + + wcc = SumpyExpansionWranglerLQCodeContainer( + ctx, + partial(mpole_expn_class, knl), + partial(local_expn_class, knl), + out_kernels) + wrangler = wcc.get_wrangler(queue, tree, dtype, + fmm_level_to_order=lambda lev: order, + kernel_extra_kwargs=extra_kwargs) + + from boxtree.fmm import drive_fmm_local_quad + pot, = drive_fmm_local_quad(intersection, intersection_starts, area_l1_diff, + area_l1_diff_starts, tree, trav, wrangler, + weights) + + from sumpy import P2P + p2p = P2P(ctx, out_kernels, exclude_self=False) + evt, (ref_pot,) = p2p(queue, sources, sources, (weights,), + **extra_kwargs) + + pot = pot.get() + ref_pot = ref_pot.get() + + rel_err = la.norm(pot - ref_pot) / la.norm(ref_pot) + logger.info("order %d -> relative l2 error: %g" % (order, rel_err)) + + assert rel_err < 1e-8 + + +# You can test individual routines by typing +# $ python test_fmm.py 'test_routine(cl.create_some_context)' + +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