From 784f2b0cb086ca552544e80d3643f7033dfa687f Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sat, 19 May 2018 18:57:53 -0500 Subject: [PATCH 01/16] connection: add method to compute chained resampling matrix. Just like in the case of the DirectDiscretizationConnection, this is going to be very, very slow and should only be used when absolutely necessary. --- .../discretization/connection/__init__.py | 19 +++ meshmode/mesh/refinement/no_adjacency.py | 2 +- test/test_chained.py | 160 ++++++++++++++++++ 3 files changed, 180 insertions(+), 1 deletion(-) create mode 100644 test/test_chained.py diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 30e886ed..223ac9c5 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -202,6 +202,9 @@ class DiscretizationConnection(object): self.is_surjective = is_surjective + def full_resample_matrix(self): + raise NotImplementedError() + def __call__(self, queue, vec): raise NotImplementedError() @@ -236,6 +239,22 @@ class ChainedDiscretizationConnection(DiscretizationConnection): self.connections = connections + @memoize_method + def full_resample_matrix(self, queue): + if not self.connections: + result = cl.array.to_device(queue, np.eye(self.to_discr.nnodes)) + return result + + acc = self.connections[0].full_resample_matrix(queue).get(queue) + print(acc.shape) + for cnx in self.connections[1:]: + resampler = cnx.full_resample_matrix(queue).get(queue) + print(resampler.shape) + acc = resampler.dot(acc) + print(acc.shape) + + return cl.array.to_device(queue, acc) + def __call__(self, queue, vec): for cnx in self.connections: vec = cnx(queue, vec) diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index 4413117d..7ab02e70 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -134,7 +134,7 @@ class RefinerWithoutAdjacency(object): # }}} def refine_uniformly(self): - flags = np.ones(self._last_mesh.nelements, dtype=bool) + flags = np.ones(self._current_mesh.nelements, dtype=bool) self.refine(flags) # {{{ refinement top-level diff --git a/test/test_chained.py b/test/test_chained.py new file mode 100644 index 00000000..16b2fdaa --- /dev/null +++ b/test/test_chained.py @@ -0,0 +1,160 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2018 Alexandru Fikl" + +__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 os + +import numpy as np +import numpy.linalg as la +import matplotlib.pyplot as pt + +import pyopencl as cl +import pyopencl.array # noqa +import pyopencl.clmath # noqa + +import pytest +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl + as pytest_generate_tests) + +import logging +logger = logging.getLogger(__name__) + +def create_chained_connection(queue, ndim, + nelements=30, + mesh_order=4, + discr_order=4, + visualize=False): + ctx = queue.context + connections = [None, None] + + # construct base mesh + if ndim == 2: + from functools import partial + from meshmode.mesh.generation import make_curve_mesh, ellipse + mesh = make_curve_mesh(partial(ellipse, 2), + np.linspace(0.0, 1.0, nelements + 1), + order=mesh_order) + elif ndim == 3: + from meshmode.mesh.generation import generate_torus + mesh = generate_torus(10.0, 5.0, order=mesh_order) + else: + raise ValueError("Unsupported dimension: {}".format(ndim)) + + # create base discretization + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + discr = Discretization(ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(discr_order)) + + # connection 1: refine + from meshmode.mesh.refinement import RefinerWithoutAdjacency + from meshmode.discretization.connection import make_refinement_connection + refiner = RefinerWithoutAdjacency(mesh) + + if ndim == 2: + flags = np.round(np.random.rand(refiner.get_current_mesh().nelements)) + refiner.refine(flags) + else: + refiner.refine_uniformly() + connections[0] = make_refinement_connection(refiner, discr, + InterpolatoryQuadratureSimplexGroupFactory(discr_order)) + + if visualize and ndim == 2: + from_nodes = connections[0].from_discr.nodes().get(queue) + to_nodes = connections[0].to_discr.nodes().get(queue) + + pt.plot(to_nodes[0], to_nodes[1], 'k') + pt.plot(from_nodes[0], from_nodes[1], 'o', mfc='none') + pt.plot(to_nodes[0], to_nodes[1], '.') + filename = os.path.join(os.path.dirname(__file__), + "test_chained_0_{}d.png".format(ndim)) + pt.tight_layout() + pt.savefig(filename, dpi=300) + pt.clf() + + # connection 2: refine / restrict to face + if ndim == 2: + flags = np.round(np.random.rand(refiner.get_current_mesh().nelements)) + refiner.refine(flags) + connections[1] = make_refinement_connection(refiner, + connections[0].to_discr, + InterpolatoryQuadratureSimplexGroupFactory(discr_order)) + else: + from meshmode.discretization.connection import make_face_restriction + from meshmode.discretization.connection import FACE_RESTR_INTERIOR + connections[1] = make_face_restriction(connections[0].to_discr, + InterpolatoryQuadratureSimplexGroupFactory(discr_order), + FACE_RESTR_INTERIOR, per_face_groups=True) + + if visualize and ndim == 2: + from_nodes = connections[1].from_discr.nodes().get(queue) + to_nodes = connections[1].to_discr.nodes().get(queue) + + pt.plot(to_nodes[0], to_nodes[1], 'k') + pt.plot(from_nodes[0], from_nodes[1], 'o', mfc='none') + pt.plot(to_nodes[0], to_nodes[1], '.') + filename = os.path.join(os.path.dirname(__file__), + "test_chained_1_{}d.png".format(ndim)) + pt.tight_layout() + pt.savefig(filename, dpi=300) + pt.clf() + + from meshmode.discretization.connection import ChainedDiscretizationConnection + chained = ChainedDiscretizationConnection(connections) + + return chained + + +@pytest.mark.parametrize("ndim", [2, 3]) +def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + chained = create_chained_connection(queue, ndim, visualize=visualize) + connections = chained.connections + + x = connections[0].from_discr.nodes().with_queue(queue) + def f(x): + from six.moves import reduce + return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) + + fx = f(x) + resample_mat = chained.full_resample_matrix(queue).get(queue) + f1 = np.dot(resample_mat, fx.get(queue)) + f2 = chained(queue, fx).get(queue) + f3 = connections[1](queue, connections[0](queue, fx)).get(queue) + + assert np.allclose(f1, f2) + assert np.allclose(f2, f3) + + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) + +# vim: fdm=marker -- GitLab From 809aef7d002def7389d54ae3d41afd6e9850279d Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sat, 19 May 2018 19:08:16 -0500 Subject: [PATCH 02/16] tests: make test run a bit faster --- .../discretization/connection/__init__.py | 1 + test/test_chained.py | 31 +++++++------------ 2 files changed, 12 insertions(+), 20 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 223ac9c5..1001c985 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -245,6 +245,7 @@ class ChainedDiscretizationConnection(DiscretizationConnection): result = cl.array.to_device(queue, np.eye(self.to_discr.nnodes)) return result + print(len(self.connections)) acc = self.connections[0].full_resample_matrix(queue).get(queue) print(acc.shape) for cnx in self.connections[1:]: diff --git a/test/test_chained.py b/test/test_chained.py index 16b2fdaa..8f675f43 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -40,9 +40,9 @@ import logging logger = logging.getLogger(__name__) def create_chained_connection(queue, ndim, - nelements=30, - mesh_order=4, - discr_order=4, + nelements=42, + mesh_order=2, + discr_order=2, visualize=False): ctx = queue.context connections = [None, None] @@ -72,11 +72,9 @@ def create_chained_connection(queue, ndim, from meshmode.discretization.connection import make_refinement_connection refiner = RefinerWithoutAdjacency(mesh) - if ndim == 2: - flags = np.round(np.random.rand(refiner.get_current_mesh().nelements)) - refiner.refine(flags) - else: - refiner.refine_uniformly() + threshold = 0.2 if ndim == 3 else 0.5 + flags = np.random.rand(refiner.get_current_mesh().nelements) < threshold + refiner.refine(flags) connections[0] = make_refinement_connection(refiner, discr, InterpolatoryQuadratureSimplexGroupFactory(discr_order)) @@ -94,18 +92,11 @@ def create_chained_connection(queue, ndim, pt.clf() # connection 2: refine / restrict to face - if ndim == 2: - flags = np.round(np.random.rand(refiner.get_current_mesh().nelements)) - refiner.refine(flags) - connections[1] = make_refinement_connection(refiner, - connections[0].to_discr, - InterpolatoryQuadratureSimplexGroupFactory(discr_order)) - else: - from meshmode.discretization.connection import make_face_restriction - from meshmode.discretization.connection import FACE_RESTR_INTERIOR - connections[1] = make_face_restriction(connections[0].to_discr, - InterpolatoryQuadratureSimplexGroupFactory(discr_order), - FACE_RESTR_INTERIOR, per_face_groups=True) + flags = np.random.rand(refiner.get_current_mesh().nelements) < threshold + refiner.refine(flags) + connections[1] = make_refinement_connection(refiner, + connections[0].to_discr, + InterpolatoryQuadratureSimplexGroupFactory(discr_order)) if visualize and ndim == 2: from_nodes = connections[1].from_discr.nodes().get(queue) -- GitLab From 07f75a8ec8cbb49321731723d7558ba08a41fbfa Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 21 May 2018 21:03:13 -0500 Subject: [PATCH 03/16] connection: add method to collapse chained connection to a direct connection --- .../discretization/connection/__init__.py | 151 ++++++++++++++++- test/test_chained.py | 153 +++++++++++++++++- 2 files changed, 289 insertions(+), 15 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 1001c985..bcae44b1 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -1,6 +1,9 @@ from __future__ import division, print_function, absolute_import -__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" +__copyright__ = """ +Copyright (C) 2014 Andreas Kloeckner +Copyright (C) 2018 Alexandru Fikl +""" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -29,8 +32,11 @@ import pyopencl as cl import pyopencl.array # noqa from loopy.version import MOST_RECENT_LANGUAGE_VERSION +from pytools import Record from pytools import memoize_method, memoize_in +import modepy as mp + from meshmode.discretization.connection.same_mesh import \ make_same_mesh_connection from meshmode.discretization.connection.face import ( @@ -81,6 +87,10 @@ Implementation details """ +class _ConnectionGroupData(Record): + pass + + # {{{ interpolation batch class InterpolationBatch(object): @@ -205,6 +215,9 @@ class DiscretizationConnection(object): def full_resample_matrix(self): raise NotImplementedError() + def direct_connection(self): + raise NotImplementedError() + def __call__(self, queue, vec): raise NotImplementedError() @@ -239,23 +252,145 @@ class ChainedDiscretizationConnection(DiscretizationConnection): self.connections = connections - @memoize_method def full_resample_matrix(self, queue): if not self.connections: result = cl.array.to_device(queue, np.eye(self.to_discr.nnodes)) return result - print(len(self.connections)) acc = self.connections[0].full_resample_matrix(queue).get(queue) - print(acc.shape) for cnx in self.connections[1:]: resampler = cnx.full_resample_matrix(queue).get(queue) - print(resampler.shape) acc = resampler.dot(acc) - print(acc.shape) return cl.array.to_device(queue, acc) + def _element_to_batch(self, queue, conn): + nelements = np.sum([g.nelements for g in conn.from_discr.groups]) + el_map = [[] for _ in range(nelements)] + + for i, igrp in enumerate(conn.groups): + for j, ibatch in enumerate(igrp.batches): + for p, k in enumerate(ibatch.from_element_indices.get(queue)): + el_map[k].append((i, j, p)) + + return el_map + + def _merge_groups(self, from_discr, from_groups, to_conn): + def iterbatches(groups): + for igrp, group in enumerate(groups): + for ibatch, batch in enumerate(group.batches): + yield (igrp, ibatch), (group, batch) + + def find_batch(nodes, gtb): + for igrp, batches in enumerate(gtb): + for ibatch, batch in enumerate(batches): + if np.allclose(nodes, batch.result_unit_nodes): + return (igrp, ibatch) + return (-1, -1) + + nfrom_groups = len(from_groups) + nto_groups = len(to_conn.groups) + + grp_to_grp = {} + grp_to_batch = [[] for i in range(nfrom_groups * nto_groups)] + for (igrp, ibatch), (fgrp, fbatch) in iterbatches(from_groups): + for (jgrp, jbatch), (tgrp, tbatch) in iterbatches(to_conn.groups): + # compute result_unit_nodes + ffgrp = from_discr.groups[fbatch.from_group_index] + from_matrix = mp.resampling_matrix( + ffgrp.basis(), + fbatch.result_unit_nodes, + ffgrp.unit_nodes) + result_unit_nodes = from_matrix.dot(ffgrp.unit_nodes.T) + + tfgrp = to_conn.from_discr.groups[tbatch.from_group_index] + to_matrix = mp.resampling_matrix( + tfgrp.basis(), + tbatch.result_unit_nodes, + tfgrp.unit_nodes) + result_unit_nodes = to_matrix.dot(result_unit_nodes).T + + # find new (group, batch) + (mgroup, mbatch) = find_batch(result_unit_nodes, grp_to_batch) + if mgroup < 0: + mgroup = nto_groups * igrp + jgrp + + grp_to_batch[mgroup].append(_ConnectionGroupData( + from_group_index=fbatch.from_group_index, + from_batch_index=(igrp, ibatch), + to_batch_index=(jgrp, jbatch), + result_unit_nodes=result_unit_nodes, + to_element_face=tbatch.to_element_face)) + mbatch = len(grp_to_batch[mgroup]) - 1 + + grp_to_grp[(igrp, ibatch, jgrp, jbatch)] = (mgroup, mbatch) + + return grp_to_grp, grp_to_batch + + def _construct_batches(self, queue, from_bins, to_bins, grp_to_batch): + def to_device(x): + return cl.array.to_device(queue, np.asarray(x)) + + for ibatch, (from_bin, to_bin) in enumerate(zip(from_bins, to_bins)): + info = grp_to_batch[ibatch] + + yield InterpolationBatch( + from_group_index=info.from_group_index, + from_element_indices=to_device(from_bin), + to_element_indices=to_device(to_bin), + result_unit_nodes=info.result_unit_nodes, + to_element_face=info.to_element_face) + + def direct_connection(self, queue): + def iterbatches(groups): + for igrp, grp in enumerate(groups): + for ibatch, batch in enumerate(grp.batches): + yield (igrp, ibatch), (grp, batch) + + if not self.connections: + return make_same_mesh_connection(self.to_discr, self.from_discr) + + connections = [] + for conn in self.connections: + connections.append(conn.direct_connection()) + + groups = connections[0].groups + for to_conn in connections[1:]: + # create an inverse mapping for the elements in to_conn.from_discr + el_to_batch = self._element_to_batch(queue, to_conn) + # create new group structure to go from groups -> to_conn.to_discr + grp_to_grp, grp_to_batch = self._merge_groups(self.from_discr, + groups, + to_conn) + + # distribute the indices to new groups and batches + from_bins = [[[] for _ in g] for g in grp_to_batch] + to_bins = [[[] for _ in g] for g in grp_to_batch] + + for (igrp, ibatch), (_, batch) in iterbatches(groups): + for ifrom, imid in zip(batch.from_element_indices.get(queue), + batch.to_element_indices.get(queue)): + for jgrp, jbatch, iel in el_to_batch[imid]: + mgrp, mbatch = grp_to_grp[(igrp, ibatch, jgrp, jbatch)] + ito = to_conn.groups[jgrp].batches[jbatch] \ + .to_element_indices[iel].get(queue) + + from_bins[mgrp][mbatch].append(ifrom) + to_bins[mgrp][mbatch].append(ito) + + # construct new groups + groups = [] + for igrp, (from_bin, to_bin) in enumerate(zip(from_bins, to_bins)): + groups.append(DiscretizationConnectionElementGroup( + list(self._construct_batches(queue, from_bin, to_bin, + grp_to_batch[igrp])))) + + return DirectDiscretizationConnection( + from_discr=self.from_discr, + to_discr=self.to_discr, + groups=groups, + is_surjective=self.is_surjective) + def __call__(self, queue, vec): for cnx in self.connections: vec = cnx(queue, vec) @@ -296,7 +431,6 @@ class DirectDiscretizationConnection(DiscretizationConnection): @memoize_method def _resample_matrix(self, to_group_index, ibatch_index): - import modepy as mp ibatch = self.groups[to_group_index].batches[ibatch_index] from_grp = self.from_discr.groups[ibatch.from_group_index] @@ -418,6 +552,9 @@ class DirectDiscretizationConnection(DiscretizationConnection): return result + def direct_connection(self): + return self + def __call__(self, queue, vec): @memoize_in(self, "resample_by_mat_knl") def mat_knl(): diff --git a/test/test_chained.py b/test/test_chained.py index 8f675f43..73a53cdb 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -24,7 +24,6 @@ THE SOFTWARE. import os import numpy as np -import numpy.linalg as la import matplotlib.pyplot as pt import pyopencl as cl @@ -39,10 +38,11 @@ from pyopencl.tools import ( # noqa import logging logger = logging.getLogger(__name__) + def create_chained_connection(queue, ndim, nelements=42, - mesh_order=2, - discr_order=2, + mesh_order=5, + discr_order=5, visualize=False): ctx = queue.context connections = [None, None] @@ -86,7 +86,7 @@ def create_chained_connection(queue, ndim, pt.plot(from_nodes[0], from_nodes[1], 'o', mfc='none') pt.plot(to_nodes[0], to_nodes[1], '.') filename = os.path.join(os.path.dirname(__file__), - "test_chained_0_{}d.png".format(ndim)) + "test_chained_nodes_0.png".format(ndim)) pt.tight_layout() pt.savefig(filename, dpi=300) pt.clf() @@ -106,7 +106,20 @@ def create_chained_connection(queue, ndim, pt.plot(from_nodes[0], from_nodes[1], 'o', mfc='none') pt.plot(to_nodes[0], to_nodes[1], '.') filename = os.path.join(os.path.dirname(__file__), - "test_chained_1_{}d.png".format(ndim)) + "test_chained_nodes_1.png".format(ndim)) + pt.tight_layout() + pt.savefig(filename, dpi=300) + pt.clf() + + if visualize and ndim == 2: + from_nodes = connections[0].from_discr.nodes().get(queue) + to_nodes = connections[1].to_discr.nodes().get(queue) + + pt.plot(to_nodes[0], to_nodes[1], 'k') + pt.plot(from_nodes[0], from_nodes[1], 'o', mfc='none') + pt.plot(to_nodes[0], to_nodes[1], '.') + filename = os.path.join(os.path.dirname(__file__), + "test_chained_nodes_2.png".format(ndim)) pt.tight_layout() pt.savefig(filename, dpi=300) pt.clf() @@ -118,20 +131,83 @@ def create_chained_connection(queue, ndim, @pytest.mark.parametrize("ndim", [2, 3]) -def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): +def test_chained_batch_map(ctx_factory, ndim, visualize=False): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + chained = create_chained_connection(queue, ndim, visualize=visualize) + conn = chained.connections[1] + el_to_batch = chained._element_to_batch(queue, conn) + + for igrp, grp in enumerate(conn.groups): + for ibatch, batch in enumerate(grp.batches): + for p, k in enumerate(batch.from_element_indices.get(queue)): + assert (igrp, ibatch, p) in el_to_batch[k] + + for i, batches in enumerate(el_to_batch): + for igrp, ibatch, p in batches: + batch = conn.groups[igrp].batches[ibatch] + assert i == batch.from_element_indices[p] + + +@pytest.mark.parametrize("ndim", [2, 3]) +def test_chained_merge_groups(ctx_factory, ndim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) chained = create_chained_connection(queue, ndim, visualize=visualize) connections = chained.connections + grp_to_grp, grp_to_batch = chained._merge_groups(chained.from_discr, + connections[0].groups, + connections[1]) + + if visualize: + import matplotlib.pyplot as pt + + pt.figure(figsize=(10, 8)) + for k, v in grp_to_grp.items(): + igrp, ibatch, jgrp, jbatch = k + mgroup, mbatch = v + from_group_index = connections[0].groups[igrp] \ + .batches[ibatch].from_group_index + + from_unit_nodes = connections[0].from_discr \ + .groups[from_group_index].unit_nodes + to_unit_nodes = grp_to_batch[mgroup][mbatch].result_unit_nodes + + print(k) + print(v) + + if ndim == 2: + pt.plot(from_unit_nodes, 'o') + pt.plot(to_unit_nodes, '^') + else: + pt.plot(from_unit_nodes[0], from_unit_nodes[1], 'o') + pt.plot(to_unit_nodes[0], to_unit_nodes[1], '^') + + pt.savefig('test_grp_to_grp_{}d_{:05d}_{:05d}.png' + .format(ndim, mgroup, mbatch), dpi=300) + pt.clf() + + +@pytest.mark.parametrize("ndim", [2, 3]) +def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + chained = create_chained_connection(queue, ndim, + mesh_order=2, discr_order=2, + visualize=visualize) + connections = chained.connections - x = connections[0].from_discr.nodes().with_queue(queue) def f(x): from six.moves import reduce return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) - fx = f(x) resample_mat = chained.full_resample_matrix(queue).get(queue) + + x = connections[0].from_discr.nodes().with_queue(queue) + fx = f(x) f1 = np.dot(resample_mat, fx.get(queue)) f2 = chained(queue, fx).get(queue) f3 = connections[1](queue, connections[0](queue, fx)).get(queue) @@ -140,6 +216,67 @@ def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): assert np.allclose(f2, f3) +@pytest.mark.parametrize("ndim", [2, 3]) +def test_chained_to_direct(ctx_factory, ndim, visualize=False): + import time + + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + chained = create_chained_connection(queue, ndim, visualize=visualize) + connections = chained.connections + direct = chained.direct_connection(queue) + + to_element_indices = np.full(direct.to_discr.mesh.nelements, 0, + dtype=np.int) + for grp in direct.groups: + for batch in grp.batches: + for i in batch.to_element_indices.get(queue): + to_element_indices[i] += 1 + assert np.min(to_element_indices) > 0 + + def f(x): + from six.moves import reduce + return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) + + x = connections[0].from_discr.nodes().with_queue(queue) + fx = f(x) + + t_start = time.time() + f1 = direct(queue, fx).get(queue) + t_end = time.time() + if visualize: + print('[TIME] Direct: {:.5e}'.format(t_end - t_start)) + + t_start = time.time() + f2 = chained(queue, fx).get(queue) + t_end = time.time() + if visualize: + print('[TIME] Chained: {:.5e}'.format(t_end - t_start)) + + t_start = time.time() + f3 = connections[1](queue, connections[0](queue, fx)).get(queue) + t_end = time.time() + if visualize: + print('[TIME] Chained: {:.5e}'.format(t_end - t_start)) + + if visualize and ndim == 2: + import matplotlib.pyplot as pt + + pt.figure(figsize=(10, 8)) + pt.plot(f1, label='Direct') + pt.plot(f2, label='Chained') + pt.ylim([np.min(f2) - 0.1, np.max(f2) + 0.1]) + pt.legend() + filename = os.path.join(os.path.dirname(__file__), + "test_chained_to_direct_fn.png".format(ndim)) + pt.savefig(filename, dpi=300) + pt.clf() + + assert np.allclose(f1, f2) + assert np.allclose(f2, f3) + + if __name__ == "__main__": import sys if len(sys.argv) > 1: -- GitLab From b2ba0bedbd9c045533b44a7955cbb976076786e3 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 21 May 2018 21:22:19 -0500 Subject: [PATCH 04/16] remove stray matplotlib import --- test/test_chained.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_chained.py b/test/test_chained.py index 73a53cdb..1ceabd87 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -24,7 +24,6 @@ THE SOFTWARE. import os import numpy as np -import matplotlib.pyplot as pt import pyopencl as cl import pyopencl.array # noqa -- GitLab From 2057537ff58144c8a06501a1b0a501fe8edf1396 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 21 May 2018 21:57:36 -0500 Subject: [PATCH 05/16] clean up matplotlib usage --- test/test_chained.py | 23 ++++++++++------------- 1 file changed, 10 insertions(+), 13 deletions(-) diff --git a/test/test_chained.py b/test/test_chained.py index 1ceabd87..bf3e20ea 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -21,7 +21,6 @@ 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 os import numpy as np @@ -78,16 +77,16 @@ def create_chained_connection(queue, ndim, InterpolatoryQuadratureSimplexGroupFactory(discr_order)) if visualize and ndim == 2: + import matplotlib.pyplot as pt + from_nodes = connections[0].from_discr.nodes().get(queue) to_nodes = connections[0].to_discr.nodes().get(queue) pt.plot(to_nodes[0], to_nodes[1], 'k') pt.plot(from_nodes[0], from_nodes[1], 'o', mfc='none') pt.plot(to_nodes[0], to_nodes[1], '.') - filename = os.path.join(os.path.dirname(__file__), - "test_chained_nodes_0.png".format(ndim)) pt.tight_layout() - pt.savefig(filename, dpi=300) + pt.savefig("test_chained_nodes_0.png", dpi=300) pt.clf() # connection 2: refine / restrict to face @@ -98,29 +97,29 @@ def create_chained_connection(queue, ndim, InterpolatoryQuadratureSimplexGroupFactory(discr_order)) if visualize and ndim == 2: + import matplotlib.pyplot as pt + from_nodes = connections[1].from_discr.nodes().get(queue) to_nodes = connections[1].to_discr.nodes().get(queue) pt.plot(to_nodes[0], to_nodes[1], 'k') pt.plot(from_nodes[0], from_nodes[1], 'o', mfc='none') pt.plot(to_nodes[0], to_nodes[1], '.') - filename = os.path.join(os.path.dirname(__file__), - "test_chained_nodes_1.png".format(ndim)) pt.tight_layout() - pt.savefig(filename, dpi=300) + pt.savefig("test_chained_nodes_1.png", dpi=300) pt.clf() if visualize and ndim == 2: + import matplotlib.pyplot as pt + from_nodes = connections[0].from_discr.nodes().get(queue) to_nodes = connections[1].to_discr.nodes().get(queue) pt.plot(to_nodes[0], to_nodes[1], 'k') pt.plot(from_nodes[0], from_nodes[1], 'o', mfc='none') pt.plot(to_nodes[0], to_nodes[1], '.') - filename = os.path.join(os.path.dirname(__file__), - "test_chained_nodes_2.png".format(ndim)) pt.tight_layout() - pt.savefig(filename, dpi=300) + pt.savefig("test_chained_nodes_2.png", dpi=300) pt.clf() from meshmode.discretization.connection import ChainedDiscretizationConnection @@ -267,9 +266,7 @@ def test_chained_to_direct(ctx_factory, ndim, visualize=False): pt.plot(f2, label='Chained') pt.ylim([np.min(f2) - 0.1, np.max(f2) + 0.1]) pt.legend() - filename = os.path.join(os.path.dirname(__file__), - "test_chained_to_direct_fn.png".format(ndim)) - pt.savefig(filename, dpi=300) + pt.savefig('test_chained_to_direct.png', dpi=300) pt.clf() assert np.allclose(f1, f2) -- GitLab From 9dde0cfc57629485733bcc1a5433d302bcc89c7f Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 28 May 2018 10:58:47 -0500 Subject: [PATCH 06/16] connection: move chained to direct squishing to a separate file --- .../discretization/connection/__init__.py | 171 +-------------- meshmode/discretization/connection/chained.py | 197 ++++++++++++++++++ 2 files changed, 207 insertions(+), 161 deletions(-) create mode 100644 meshmode/discretization/connection/chained.py diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index bcae44b1..aecf20d9 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -1,9 +1,6 @@ from __future__ import division, print_function, absolute_import -__copyright__ = """ -Copyright (C) 2014 Andreas Kloeckner -Copyright (C) 2018 Alexandru Fikl -""" +__copyright__ = "Copyright (C) 2014 Andreas Kloeckner" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -32,11 +29,8 @@ import pyopencl as cl import pyopencl.array # noqa from loopy.version import MOST_RECENT_LANGUAGE_VERSION -from pytools import Record from pytools import memoize_method, memoize_in -import modepy as mp - from meshmode.discretization.connection.same_mesh import \ make_same_mesh_connection from meshmode.discretization.connection.face import ( @@ -46,7 +40,8 @@ from meshmode.discretization.connection.opposite_face import \ make_opposite_face_connection from meshmode.discretization.connection.refinement import \ make_refinement_connection - +from meshmode.discretization.connection.chained import \ + make_direct_connection, make_full_resample_matrix import logging logger = logging.getLogger(__name__) @@ -59,7 +54,9 @@ __all__ = [ "make_face_restriction", "make_face_to_all_faces_embedding", "make_opposite_face_connection", - "make_refinement_connection" + "make_refinement_connection", + "make_direct_connection", + "make_full_resample_matrix" ] __doc__ = """ @@ -78,6 +75,9 @@ __doc__ = """ .. autofunction:: make_refinement_connection +.. autofunction:: make_direct_connection +.. autofunction:: make_full_resample_matrix + Implementation details ^^^^^^^^^^^^^^^^^^^^^^ @@ -87,10 +87,6 @@ Implementation details """ -class _ConnectionGroupData(Record): - pass - - # {{{ interpolation batch class InterpolationBatch(object): @@ -212,12 +208,6 @@ class DiscretizationConnection(object): self.is_surjective = is_surjective - def full_resample_matrix(self): - raise NotImplementedError() - - def direct_connection(self): - raise NotImplementedError() - def __call__(self, queue, vec): raise NotImplementedError() @@ -252,145 +242,6 @@ class ChainedDiscretizationConnection(DiscretizationConnection): self.connections = connections - def full_resample_matrix(self, queue): - if not self.connections: - result = cl.array.to_device(queue, np.eye(self.to_discr.nnodes)) - return result - - acc = self.connections[0].full_resample_matrix(queue).get(queue) - for cnx in self.connections[1:]: - resampler = cnx.full_resample_matrix(queue).get(queue) - acc = resampler.dot(acc) - - return cl.array.to_device(queue, acc) - - def _element_to_batch(self, queue, conn): - nelements = np.sum([g.nelements for g in conn.from_discr.groups]) - el_map = [[] for _ in range(nelements)] - - for i, igrp in enumerate(conn.groups): - for j, ibatch in enumerate(igrp.batches): - for p, k in enumerate(ibatch.from_element_indices.get(queue)): - el_map[k].append((i, j, p)) - - return el_map - - def _merge_groups(self, from_discr, from_groups, to_conn): - def iterbatches(groups): - for igrp, group in enumerate(groups): - for ibatch, batch in enumerate(group.batches): - yield (igrp, ibatch), (group, batch) - - def find_batch(nodes, gtb): - for igrp, batches in enumerate(gtb): - for ibatch, batch in enumerate(batches): - if np.allclose(nodes, batch.result_unit_nodes): - return (igrp, ibatch) - return (-1, -1) - - nfrom_groups = len(from_groups) - nto_groups = len(to_conn.groups) - - grp_to_grp = {} - grp_to_batch = [[] for i in range(nfrom_groups * nto_groups)] - for (igrp, ibatch), (fgrp, fbatch) in iterbatches(from_groups): - for (jgrp, jbatch), (tgrp, tbatch) in iterbatches(to_conn.groups): - # compute result_unit_nodes - ffgrp = from_discr.groups[fbatch.from_group_index] - from_matrix = mp.resampling_matrix( - ffgrp.basis(), - fbatch.result_unit_nodes, - ffgrp.unit_nodes) - result_unit_nodes = from_matrix.dot(ffgrp.unit_nodes.T) - - tfgrp = to_conn.from_discr.groups[tbatch.from_group_index] - to_matrix = mp.resampling_matrix( - tfgrp.basis(), - tbatch.result_unit_nodes, - tfgrp.unit_nodes) - result_unit_nodes = to_matrix.dot(result_unit_nodes).T - - # find new (group, batch) - (mgroup, mbatch) = find_batch(result_unit_nodes, grp_to_batch) - if mgroup < 0: - mgroup = nto_groups * igrp + jgrp - - grp_to_batch[mgroup].append(_ConnectionGroupData( - from_group_index=fbatch.from_group_index, - from_batch_index=(igrp, ibatch), - to_batch_index=(jgrp, jbatch), - result_unit_nodes=result_unit_nodes, - to_element_face=tbatch.to_element_face)) - mbatch = len(grp_to_batch[mgroup]) - 1 - - grp_to_grp[(igrp, ibatch, jgrp, jbatch)] = (mgroup, mbatch) - - return grp_to_grp, grp_to_batch - - def _construct_batches(self, queue, from_bins, to_bins, grp_to_batch): - def to_device(x): - return cl.array.to_device(queue, np.asarray(x)) - - for ibatch, (from_bin, to_bin) in enumerate(zip(from_bins, to_bins)): - info = grp_to_batch[ibatch] - - yield InterpolationBatch( - from_group_index=info.from_group_index, - from_element_indices=to_device(from_bin), - to_element_indices=to_device(to_bin), - result_unit_nodes=info.result_unit_nodes, - to_element_face=info.to_element_face) - - def direct_connection(self, queue): - def iterbatches(groups): - for igrp, grp in enumerate(groups): - for ibatch, batch in enumerate(grp.batches): - yield (igrp, ibatch), (grp, batch) - - if not self.connections: - return make_same_mesh_connection(self.to_discr, self.from_discr) - - connections = [] - for conn in self.connections: - connections.append(conn.direct_connection()) - - groups = connections[0].groups - for to_conn in connections[1:]: - # create an inverse mapping for the elements in to_conn.from_discr - el_to_batch = self._element_to_batch(queue, to_conn) - # create new group structure to go from groups -> to_conn.to_discr - grp_to_grp, grp_to_batch = self._merge_groups(self.from_discr, - groups, - to_conn) - - # distribute the indices to new groups and batches - from_bins = [[[] for _ in g] for g in grp_to_batch] - to_bins = [[[] for _ in g] for g in grp_to_batch] - - for (igrp, ibatch), (_, batch) in iterbatches(groups): - for ifrom, imid in zip(batch.from_element_indices.get(queue), - batch.to_element_indices.get(queue)): - for jgrp, jbatch, iel in el_to_batch[imid]: - mgrp, mbatch = grp_to_grp[(igrp, ibatch, jgrp, jbatch)] - ito = to_conn.groups[jgrp].batches[jbatch] \ - .to_element_indices[iel].get(queue) - - from_bins[mgrp][mbatch].append(ifrom) - to_bins[mgrp][mbatch].append(ito) - - # construct new groups - groups = [] - for igrp, (from_bin, to_bin) in enumerate(zip(from_bins, to_bins)): - groups.append(DiscretizationConnectionElementGroup( - list(self._construct_batches(queue, from_bin, to_bin, - grp_to_batch[igrp])))) - - return DirectDiscretizationConnection( - from_discr=self.from_discr, - to_discr=self.to_discr, - groups=groups, - is_surjective=self.is_surjective) - def __call__(self, queue, vec): for cnx in self.connections: vec = cnx(queue, vec) @@ -431,6 +282,7 @@ class DirectDiscretizationConnection(DiscretizationConnection): @memoize_method def _resample_matrix(self, to_group_index, ibatch_index): + import modepy as mp ibatch = self.groups[to_group_index].batches[ibatch_index] from_grp = self.from_discr.groups[ibatch.from_group_index] @@ -552,9 +404,6 @@ class DirectDiscretizationConnection(DiscretizationConnection): return result - def direct_connection(self): - return self - def __call__(self, queue, vec): @memoize_in(self, "resample_by_mat_knl") def mat_knl(): diff --git a/meshmode/discretization/connection/chained.py b/meshmode/discretization/connection/chained.py new file mode 100644 index 00000000..02588c1d --- /dev/null +++ b/meshmode/discretization/connection/chained.py @@ -0,0 +1,197 @@ +from __future__ import division, print_function, absolute_import + +__copyright__ = """Copyright (C) 2018 Alexandru Fikl""" + +__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 pyopencl as cl +import pyopencl.array # noqa + +from pytools import Record + +import modepy as mp + +from meshmode.discretization.connection import ( + DirectDiscretizationConnection, + make_same_mesh_connection) + + +class _ConnectionBatchData(Record): + pass + + +def _iterbatches(groups): + for igrp, grp in enumerate(groups): + for ibatch, batch in enumerate(grp.batches): + yield (igrp, ibatch), (grp, batch) + + +def _build_element_lookup_table(self, queue, conn): + nelements = np.sum([g.nelements for g in conn.from_discr.groups]) + el_table = [[] for _ in range(nelements)] + + for (igrp, ibatch), (_, batch) in _iterbatches(conn.groups): + for i, k in enumerate(batch.from_element_indices.get(queue)): + el_table[k].append((i, igrp, ibatch)) + + return el_table + + +def _build_new_group_table(self, from_discr, from_groups, to_conn): + def find_batch(nodes, gtb): + for igrp, batches in enumerate(gtb): + for ibatch, batch in enumerate(batches): + if np.allclose(nodes, batch.result_unit_nodes): + return (igrp, ibatch) + return (-1, -1) + + nfrom_groups = len(from_groups) + nto_groups = len(to_conn.groups) + + grp_to_grp = {} + batch_info = [[] for i in range(nfrom_groups * nto_groups)] + for (igrp, ibatch), (fgrp, fbatch) in _iterbatches(from_groups): + for (jgrp, jbatch), (tgrp, tbatch) in _iterbatches(to_conn.groups): + # compute result_unit_nodes + ffgrp = from_discr.groups[fbatch.from_group_index] + from_matrix = mp.resampling_matrix( + ffgrp.basis(), + fbatch.result_unit_nodes, + ffgrp.unit_nodes) + result_unit_nodes = from_matrix.dot(ffgrp.unit_nodes.T) + + tfgrp = to_conn.from_discr.groups[tbatch.from_group_index] + to_matrix = mp.resampling_matrix( + tfgrp.basis(), + tbatch.result_unit_nodes, + tfgrp.unit_nodes) + result_unit_nodes = to_matrix.dot(result_unit_nodes).T + + # find new (group, batch) + (igrp_new, ibatch_new) = find_batch(result_unit_nodes, batch_info) + if igrp_new < 0: + igrp_new = nto_groups * igrp + jgrp + ibatch_new = len(batch_info[igrp_new]) + 1 + + batch_info[igrp_new].append(_ConnectionBatchData( + from_group_index=fbatch.from_group_index, + result_unit_nodes=result_unit_nodes, + to_element_face=tbatch.to_element_face)) + + grp_to_grp[(igrp, ibatch, jgrp, jbatch)] = (igrp_new, ibatch_new) + + return grp_to_grp, batch_info + + +def _build_batches(self, queue, from_bins, to_bins, batch): + def to_device(x): + return cl.array.to_device(queue, np.asarray(x)) + + for ibatch, (from_bin, to_bin) in enumerate(zip(from_bins, to_bins)): + yield InterpolationBatch( + from_group_index=batch[ibatch].from_group_index, + from_element_indices=to_device(from_bin), + to_element_indices=to_device(to_bin), + result_unit_nodes=batch[ibatch].result_unit_nodes, + to_element_face=batch[ibatch].to_element_face) + + +def make_full_resample_matrix(queue, connection): + """ + """ + + if hasattr(connection, "full_resample_matrix") + return connection.full_resample_matrix(queue) + + if not hasattr(connection, 'connections'): + raise TypeError('connection is not chained') + + if not connection.connections: + result = np.eye(connection.to_discr.nnodes) + return cl.array.to_device(queue, result) + + acc = make_full_resample_matrix(queue, connection.connections[0]).get(queue) + for conn in connection.connections[1:]: + resampler = make_full_resample_matrix(queue, conn).get(queue) + acc = resampler.dot(acc) + + return cl.array.to_device(queue, acc) + + +def make_direct_connection(queue, connection): + if isinstance(connection, DirectDiscretizationConnection): + return connection + + if not hasattr(connection, 'connections') or not connection.connections: + return make_same_mesh_connection(connection.to_discr, + connection.from_discr) + + # recursively build direct connections + connections = connection.connections + direct_connections = [] + for conn in connections: + direct_connections.append(make_direct_connection(queue, conn)) + + # merge all the direct connections + groups = direct_connections[0].groups + for to_conn in direct_connections[1:]: + el_to_batch = _build_element_lookup_table(queue, to_conn) + grp_to_grp, batch_info = _build_new_group_table(connection.from_discr, + groups, to_conn) + + # distribute the indices to new groups and batches + from_bins = [[[] for _ in g] for g in batch_info] + to_bins = [[[] for _ in g] for g in batch_info] + + to_element_indices = {} + for (igrp, ibatch), (_, tbatch) in _iterbatches(to_conn.groups): + to_element_indices[(igrp, ibatch)] = \ + tbatch.to_element_indices.get(queue) + + # NOTE: notation used: + # * `ixxx` is an index in from_conn + # * `jxxx` is an index in to_conn + # * `ito` is the same as `jfrom` (on the same discr) + for (igrp, ibatch), (_, from_batch) in _iterbatches(groups): + for ifrom, ito in zip(from_batch.from_element_indices.get(queue), + from_batch.to_element_indices.get(queue)): + for j, jgrp, jbatch in el_to_batch[ito]: + igrp_new, ibatch_new = \ + grp_to_grp[(igrp, ibatch, jgrp, jbatch)] + jto = to_element_indices[(jgrp, jbatch)][j] + + from_bins[igrp_new][ibatch_new].append(ifrom) + to_bins[igrp_new][ibatch_new].append(jto) + + # build new groups + groups = [] + for igrp, (from_bin, to_bin) in enumerate(zip(from_bins, to_bins)): + groups.append(DiscretizationConnectionElementGroup( + list(_build_batches(queue, from_bin, to_bin, + batch_info[igrp])))) + + return DirectDiscretizationConnection( + from_discr=connection.from_discr, + to_discr=connection.to_discr, + groups=groups, + is_surjective=connection.is_surjective) + -- GitLab From 72b54c8da9ac28c01834cebf18d46fad8aec35b5 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 28 May 2018 11:43:57 -0500 Subject: [PATCH 07/16] connection: add some documentation --- meshmode/discretization/connection/chained.py | 50 ++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) diff --git a/meshmode/discretization/connection/chained.py b/meshmode/discretization/connection/chained.py index 02588c1d..b2908fe2 100644 --- a/meshmode/discretization/connection/chained.py +++ b/meshmode/discretization/connection/chained.py @@ -67,6 +67,9 @@ def _build_new_group_table(self, from_discr, from_groups, to_conn): nfrom_groups = len(from_groups) nto_groups = len(to_conn.groups) + # construct a table from (old groups) -> (new groups) + # NOTE: we try to reduce the number of new groups and batches by matching + # the `result_unit_nodes` and only adding a new batch if necessary grp_to_grp = {} batch_info = [[] for i in range(nfrom_groups * nto_groups)] for (igrp, ibatch), (fgrp, fbatch) in _iterbatches(from_groups): @@ -116,7 +119,22 @@ def _build_batches(self, queue, from_bins, to_bins, batch): def make_full_resample_matrix(queue, connection): - """ + """Build a dense matrix representing the discretization connection. + + This is based on + :func:`~meshmode.discretization.connection.DirectDiscretizationConnection.full_resample_matrix`. + If a chained connection is given, the matrix is constructed recursively + for each connection and multiplied left to right. + + .. warning:: + + This method will be very slow, both in terms of speed and memory + usage, and should only be used for testing or if absolutely necessary. + + :arg queue: a :class:`pyopencl.CommandQueue`. + :arg connection: a :class:`~meshmode.discretization.connection.DiscretizationConnection`. + :return: a :class:`pyopencl.array.Array` of shape + `(connection.from_discr.nnodes, connection.to_discr.nnodes)`. """ if hasattr(connection, "full_resample_matrix") @@ -138,6 +156,36 @@ def make_full_resample_matrix(queue, connection): def make_direct_connection(queue, connection): + """Collapse a connection into a direct connection. + + If the given connection is already a + :class:`~meshmode.discretization.connection.DirectDiscretizationConnection` + nothing is done. However, if the connection is a + :class:`~meshmode.discretization.connection.ChainedDiscretizationConnection`, + a new direct connection is constructed that transports from + :attr:`connection.from_discr` to :attr:`connection.to_discr`. + + The new direct connection will have a number of groups and batches that + is, at worse, the product of all the connections in the chain. For + example, if we consider a connection between a discretization and a + two-level refinement, both levels will have :math:`n` groups and + :math:`m + 1` batches per group, where :math:`m` is the number of + subdivisions of an element (exact number depends on implementation + details in :method:`meshmode.discretizationc.connection.make_refinement_connection`). + However, a direct connection from level :math:`0` to level :math:`2` + will have at worst :math:`n^2` groups and each group will have + :math:`(m + 1)^2` batches. + + .. warning:: + + If a large number of connections is chained, the number of groups and + batches can become very large. + + :arg queue: a :class:`pyopencl.CommandQueue`. + :arg connection: a :class:`~meshmode.discretization.connection.DiscretizationConnection`. + :return: a :class:`~meshmode.discretization.connection.DirectDiscretizationConnection`. + """ + if isinstance(connection, DirectDiscretizationConnection): return connection -- GitLab From 62d6aac40c1bc1b3d9bce2cf63809fb4abd0cede Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 28 May 2018 13:38:01 -0500 Subject: [PATCH 08/16] connection: make sure tests pass again --- .../discretization/connection/__init__.py | 12 +- meshmode/discretization/connection/chained.py | 130 ++++++++++-------- test/test_chained.py | 49 ++++--- 3 files changed, 110 insertions(+), 81 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index aecf20d9..bd4f1292 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -35,13 +35,14 @@ from meshmode.discretization.connection.same_mesh import \ make_same_mesh_connection from meshmode.discretization.connection.face import ( FACE_RESTR_INTERIOR, FACE_RESTR_ALL, - make_face_restriction, make_face_to_all_faces_embedding) + make_face_restriction, + make_face_to_all_faces_embedding) from meshmode.discretization.connection.opposite_face import \ make_opposite_face_connection from meshmode.discretization.connection.refinement import \ make_refinement_connection from meshmode.discretization.connection.chained import \ - make_direct_connection, make_full_resample_matrix + flatten_chained_connection import logging logger = logging.getLogger(__name__) @@ -55,12 +56,12 @@ __all__ = [ "make_face_to_all_faces_embedding", "make_opposite_face_connection", "make_refinement_connection", - "make_direct_connection", - "make_full_resample_matrix" + "flatten_chained_connection" ] __doc__ = """ .. autoclass:: DiscretizationConnection +.. autoclass:: ChainedDiscretizationConnection .. autoclass:: DirectDiscretizationConnection .. autofunction:: make_same_mesh_connection @@ -75,8 +76,7 @@ __doc__ = """ .. autofunction:: make_refinement_connection -.. autofunction:: make_direct_connection -.. autofunction:: make_full_resample_matrix +.. autofunction:: flatten_chained_connection Implementation details ^^^^^^^^^^^^^^^^^^^^^^ diff --git a/meshmode/discretization/connection/chained.py b/meshmode/discretization/connection/chained.py index b2908fe2..ef926158 100644 --- a/meshmode/discretization/connection/chained.py +++ b/meshmode/discretization/connection/chained.py @@ -30,10 +30,8 @@ from pytools import Record import modepy as mp -from meshmode.discretization.connection import ( - DirectDiscretizationConnection, - make_same_mesh_connection) +# {{{ flatten chained connection class _ConnectionBatchData(Record): pass @@ -45,7 +43,7 @@ def _iterbatches(groups): yield (igrp, ibatch), (grp, batch) -def _build_element_lookup_table(self, queue, conn): +def _build_element_lookup_table(queue, conn): nelements = np.sum([g.nelements for g in conn.from_discr.groups]) el_table = [[] for _ in range(nelements)] @@ -56,7 +54,7 @@ def _build_element_lookup_table(self, queue, conn): return el_table -def _build_new_group_table(self, from_discr, from_groups, to_conn): +def _build_new_group_table(from_conn, to_conn): def find_batch(nodes, gtb): for igrp, batches in enumerate(gtb): for ibatch, batch in enumerate(batches): @@ -64,7 +62,7 @@ def _build_new_group_table(self, from_discr, from_groups, to_conn): return (igrp, ibatch) return (-1, -1) - nfrom_groups = len(from_groups) + nfrom_groups = len(from_conn.groups) nto_groups = len(to_conn.groups) # construct a table from (old groups) -> (new groups) @@ -72,10 +70,10 @@ def _build_new_group_table(self, from_discr, from_groups, to_conn): # the `result_unit_nodes` and only adding a new batch if necessary grp_to_grp = {} batch_info = [[] for i in range(nfrom_groups * nto_groups)] - for (igrp, ibatch), (fgrp, fbatch) in _iterbatches(from_groups): + for (igrp, ibatch), (fgrp, fbatch) in _iterbatches(from_conn.groups): for (jgrp, jbatch), (tgrp, tbatch) in _iterbatches(to_conn.groups): # compute result_unit_nodes - ffgrp = from_discr.groups[fbatch.from_group_index] + ffgrp = from_conn.from_discr.groups[fbatch.from_group_index] from_matrix = mp.resampling_matrix( ffgrp.basis(), fbatch.result_unit_nodes, @@ -93,7 +91,7 @@ def _build_new_group_table(self, from_discr, from_groups, to_conn): (igrp_new, ibatch_new) = find_batch(result_unit_nodes, batch_info) if igrp_new < 0: igrp_new = nto_groups * igrp + jgrp - ibatch_new = len(batch_info[igrp_new]) + 1 + ibatch_new = len(batch_info[igrp_new]) batch_info[igrp_new].append(_ConnectionBatchData( from_group_index=fbatch.from_group_index, @@ -105,7 +103,10 @@ def _build_new_group_table(self, from_discr, from_groups, to_conn): return grp_to_grp, batch_info -def _build_batches(self, queue, from_bins, to_bins, batch): +def _build_batches(queue, from_bins, to_bins, batch): + from meshmode.discretization.connection import \ + InterpolationBatch + def to_device(x): return cl.array.to_device(queue, np.asarray(x)) @@ -118,44 +119,7 @@ def _build_batches(self, queue, from_bins, to_bins, batch): to_element_face=batch[ibatch].to_element_face) -def make_full_resample_matrix(queue, connection): - """Build a dense matrix representing the discretization connection. - - This is based on - :func:`~meshmode.discretization.connection.DirectDiscretizationConnection.full_resample_matrix`. - If a chained connection is given, the matrix is constructed recursively - for each connection and multiplied left to right. - - .. warning:: - - This method will be very slow, both in terms of speed and memory - usage, and should only be used for testing or if absolutely necessary. - - :arg queue: a :class:`pyopencl.CommandQueue`. - :arg connection: a :class:`~meshmode.discretization.connection.DiscretizationConnection`. - :return: a :class:`pyopencl.array.Array` of shape - `(connection.from_discr.nnodes, connection.to_discr.nnodes)`. - """ - - if hasattr(connection, "full_resample_matrix") - return connection.full_resample_matrix(queue) - - if not hasattr(connection, 'connections'): - raise TypeError('connection is not chained') - - if not connection.connections: - result = np.eye(connection.to_discr.nnodes) - return cl.array.to_device(queue, result) - - acc = make_full_resample_matrix(queue, connection.connections[0]).get(queue) - for conn in connection.connections[1:]: - resampler = make_full_resample_matrix(queue, conn).get(queue) - acc = resampler.dot(acc) - - return cl.array.to_device(queue, acc) - - -def make_direct_connection(queue, connection): +def flatten_chained_connection(queue, connection): """Collapse a connection into a direct connection. If the given connection is already a @@ -185,6 +149,10 @@ def make_direct_connection(queue, connection): :arg connection: a :class:`~meshmode.discretization.connection.DiscretizationConnection`. :return: a :class:`~meshmode.discretization.connection.DirectDiscretizationConnection`. """ + from meshmode.discretization.connection import ( + DirectDiscretizationConnection, + DiscretizationConnectionElementGroup, + make_same_mesh_connection) if isinstance(connection, DirectDiscretizationConnection): return connection @@ -197,14 +165,13 @@ def make_direct_connection(queue, connection): connections = connection.connections direct_connections = [] for conn in connections: - direct_connections.append(make_direct_connection(queue, conn)) + direct_connections.append(flatten_chained_connection(queue, conn)) # merge all the direct connections - groups = direct_connections[0].groups + from_conn = direct_connections[0] for to_conn in direct_connections[1:]: el_to_batch = _build_element_lookup_table(queue, to_conn) - grp_to_grp, batch_info = _build_new_group_table(connection.from_discr, - groups, to_conn) + grp_to_grp, batch_info = _build_new_group_table(from_conn, to_conn) # distribute the indices to new groups and batches from_bins = [[[] for _ in g] for g in batch_info] @@ -219,7 +186,7 @@ def make_direct_connection(queue, connection): # * `ixxx` is an index in from_conn # * `jxxx` is an index in to_conn # * `ito` is the same as `jfrom` (on the same discr) - for (igrp, ibatch), (_, from_batch) in _iterbatches(groups): + for (igrp, ibatch), (_, from_batch) in _iterbatches(from_conn.groups): for ifrom, ito in zip(from_batch.from_element_indices.get(queue), from_batch.to_element_indices.get(queue)): for j, jgrp, jbatch in el_to_batch[ito]: @@ -237,9 +204,56 @@ def make_direct_connection(queue, connection): list(_build_batches(queue, from_bin, to_bin, batch_info[igrp])))) - return DirectDiscretizationConnection( - from_discr=connection.from_discr, - to_discr=connection.to_discr, - groups=groups, - is_surjective=connection.is_surjective) + from_conn = DirectDiscretizationConnection( + from_discr=from_conn.from_discr, + to_discr=to_conn.to_discr, + groups=groups, + is_surjective=connection.is_surjective) + + return from_conn + +# }}} + + +# {{{ build chained resample matrix + +def make_full_resample_matrix(queue, connection): + """Build a dense matrix representing the discretization connection. + + This is based on + :func:`~meshmode.discretization.connection.DirectDiscretizationConnection.full_resample_matrix`. + If a chained connection is given, the matrix is constructed recursively + for each connection and multiplied left to right. + + .. warning:: + + This method will be very slow, both in terms of speed and memory + usage, and should only be used for testing or if absolutely necessary. + + :arg queue: a :class:`pyopencl.CommandQueue`. + :arg connection: a :class:`~meshmode.discretization.connection.DiscretizationConnection`. + :return: a :class:`pyopencl.array.Array` of shape + `(connection.from_discr.nnodes, connection.to_discr.nnodes)`. + """ + + if hasattr(connection, "full_resample_matrix"): + return connection.full_resample_matrix(queue) + + if not hasattr(connection, 'connections'): + raise TypeError('connection is not chained') + + if not connection.connections: + result = np.eye(connection.to_discr.nnodes) + return cl.array.to_device(queue, result) + + acc = make_full_resample_matrix(queue, connection.connections[0]).get(queue) + for conn in connection.connections[1:]: + resampler = make_full_resample_matrix(queue, conn).get(queue) + acc = resampler.dot(acc) + + return cl.array.to_device(queue, acc) + +# }}} + +# vim: foldmethod=marker diff --git a/test/test_chained.py b/test/test_chained.py index bf3e20ea..4043c481 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -129,41 +129,53 @@ def create_chained_connection(queue, ndim, @pytest.mark.parametrize("ndim", [2, 3]) -def test_chained_batch_map(ctx_factory, ndim, visualize=False): +def test_chained_batch_table(ctx_factory, ndim, visualize=False): + from meshmode.discretization.connection.chained import \ + _build_element_lookup_table + ctx = ctx_factory() queue = cl.CommandQueue(ctx) chained = create_chained_connection(queue, ndim, visualize=visualize) conn = chained.connections[1] - el_to_batch = chained._element_to_batch(queue, conn) + el_to_batch = _build_element_lookup_table(queue, conn) for igrp, grp in enumerate(conn.groups): for ibatch, batch in enumerate(grp.batches): - for p, k in enumerate(batch.from_element_indices.get(queue)): - assert (igrp, ibatch, p) in el_to_batch[k] + for i, k in enumerate(batch.from_element_indices.get(queue)): + assert (i, igrp, ibatch) in el_to_batch[k] - for i, batches in enumerate(el_to_batch): - for igrp, ibatch, p in batches: + for k, batches in enumerate(el_to_batch): + for i, igrp, ibatch in batches: batch = conn.groups[igrp].batches[ibatch] - assert i == batch.from_element_indices[p] + assert k == batch.from_element_indices[i] @pytest.mark.parametrize("ndim", [2, 3]) -def test_chained_merge_groups(ctx_factory, ndim, visualize=False): +def test_chained_new_group_table(ctx_factory, ndim, visualize=False): + from meshmode.discretization.connection.chained import \ + _build_new_group_table + ctx = ctx_factory() queue = cl.CommandQueue(ctx) - chained = create_chained_connection(queue, ndim, visualize=visualize) + chained = create_chained_connection(queue, ndim, + nelements=8, + mesh_order=2, + discr_order=2, + visualize=visualize) connections = chained.connections - grp_to_grp, grp_to_batch = chained._merge_groups(chained.from_discr, - connections[0].groups, - connections[1]) + grp_to_grp, grp_info = _build_new_group_table(connections[0], + connections[1]) if visualize: import matplotlib.pyplot as pt pt.figure(figsize=(10, 8)) for k, v in grp_to_grp.items(): + print(k) + print(v) + igrp, ibatch, jgrp, jbatch = k mgroup, mbatch = v from_group_index = connections[0].groups[igrp] \ @@ -171,10 +183,8 @@ def test_chained_merge_groups(ctx_factory, ndim, visualize=False): from_unit_nodes = connections[0].from_discr \ .groups[from_group_index].unit_nodes - to_unit_nodes = grp_to_batch[mgroup][mbatch].result_unit_nodes + to_unit_nodes = grp_info[mgroup][mbatch].result_unit_nodes - print(k) - print(v) if ndim == 2: pt.plot(from_unit_nodes, 'o') @@ -190,6 +200,9 @@ def test_chained_merge_groups(ctx_factory, ndim, visualize=False): @pytest.mark.parametrize("ndim", [2, 3]) def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): + from meshmode.discretization.connection.chained import \ + make_full_resample_matrix + ctx = ctx_factory() queue = cl.CommandQueue(ctx) @@ -202,7 +215,7 @@ def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): from six.moves import reduce return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) - resample_mat = chained.full_resample_matrix(queue).get(queue) + resample_mat = make_full_resample_matrix(queue, chained).get(queue) x = connections[0].from_discr.nodes().with_queue(queue) fx = f(x) @@ -217,13 +230,15 @@ def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): @pytest.mark.parametrize("ndim", [2, 3]) def test_chained_to_direct(ctx_factory, ndim, visualize=False): import time + from meshmode.discretization.connection.chained import \ + flatten_chained_connection ctx = ctx_factory() queue = cl.CommandQueue(ctx) chained = create_chained_connection(queue, ndim, visualize=visualize) connections = chained.connections - direct = chained.direct_connection(queue) + direct = flatten_chained_connection(queue, chained) to_element_indices = np.full(direct.to_discr.mesh.nelements, 0, dtype=np.int) -- GitLab From bc3193837ae2fb029c0f9d34b4283cdee5236a2a Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 28 May 2018 14:23:51 -0500 Subject: [PATCH 09/16] tests: extend chained connection tests a bit --- test/test_chained.py | 170 +++++++++++++++++++++++-------------------- 1 file changed, 90 insertions(+), 80 deletions(-) diff --git a/test/test_chained.py b/test/test_chained.py index 4043c481..32f09612 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -36,14 +36,11 @@ from pyopencl.tools import ( # noqa import logging logger = logging.getLogger(__name__) - -def create_chained_connection(queue, ndim, - nelements=42, - mesh_order=5, - discr_order=5, - visualize=False): +def create_discretization(queue, ndim, + nelements=42, + mesh_order=5, + discr_order=5): ctx = queue.context - connections = [None, None] # construct base mesh if ndim == 2: @@ -65,67 +62,39 @@ def create_chained_connection(queue, ndim, discr = Discretization(ctx, mesh, InterpolatoryQuadratureSimplexGroupFactory(discr_order)) - # connection 1: refine + return discr + + +def create_refined_connection(queue, discr, threshold=0.3): from meshmode.mesh.refinement import RefinerWithoutAdjacency from meshmode.discretization.connection import make_refinement_connection - refiner = RefinerWithoutAdjacency(mesh) + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory - threshold = 0.2 if ndim == 3 else 0.5 - flags = np.random.rand(refiner.get_current_mesh().nelements) < threshold + flags = np.random.rand(discr.mesh.nelements) < threshold + refiner = RefinerWithoutAdjacency(discr.mesh) refiner.refine(flags) - connections[0] = make_refinement_connection(refiner, discr, - InterpolatoryQuadratureSimplexGroupFactory(discr_order)) - - if visualize and ndim == 2: - import matplotlib.pyplot as pt - - from_nodes = connections[0].from_discr.nodes().get(queue) - to_nodes = connections[0].to_discr.nodes().get(queue) - pt.plot(to_nodes[0], to_nodes[1], 'k') - pt.plot(from_nodes[0], from_nodes[1], 'o', mfc='none') - pt.plot(to_nodes[0], to_nodes[1], '.') - pt.tight_layout() - pt.savefig("test_chained_nodes_0.png", dpi=300) - pt.clf() - - # connection 2: refine / restrict to face - flags = np.random.rand(refiner.get_current_mesh().nelements) < threshold - refiner.refine(flags) - connections[1] = make_refinement_connection(refiner, - connections[0].to_discr, + discr_order = discr.groups[0].order + connection = make_refinement_connection(refiner, discr, InterpolatoryQuadratureSimplexGroupFactory(discr_order)) - if visualize and ndim == 2: - import matplotlib.pyplot as pt - - from_nodes = connections[1].from_discr.nodes().get(queue) - to_nodes = connections[1].to_discr.nodes().get(queue) - - pt.plot(to_nodes[0], to_nodes[1], 'k') - pt.plot(from_nodes[0], from_nodes[1], 'o', mfc='none') - pt.plot(to_nodes[0], to_nodes[1], '.') - pt.tight_layout() - pt.savefig("test_chained_nodes_1.png", dpi=300) - pt.clf() - - if visualize and ndim == 2: - import matplotlib.pyplot as pt + return connection - from_nodes = connections[0].from_discr.nodes().get(queue) - to_nodes = connections[1].to_discr.nodes().get(queue) - pt.plot(to_nodes[0], to_nodes[1], 'k') - pt.plot(from_nodes[0], from_nodes[1], 'o', mfc='none') - pt.plot(to_nodes[0], to_nodes[1], '.') - pt.tight_layout() - pt.savefig("test_chained_nodes_2.png", dpi=300) - pt.clf() +def create_face_connection(queue, discr): + from meshmode.discretization.connection import FACE_RESTR_ALL + from meshmode.discretization.connection import make_face_restriction + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory - from meshmode.discretization.connection import ChainedDiscretizationConnection - chained = ChainedDiscretizationConnection(connections) + discr_order = discr.groups[0].order + connection = make_face_restriction(discr, + InterpolatoryQuadratureSimplexGroupFactory(discr_order), + FACE_RESTR_ALL, + per_face_groups=True) - return chained + return connection @pytest.mark.parametrize("ndim", [2, 3]) @@ -136,10 +105,18 @@ def test_chained_batch_table(ctx_factory, ndim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - chained = create_chained_connection(queue, ndim, visualize=visualize) + discr = create_discretization(queue, ndim) + connections = [] + conn = create_refined_connection(queue, discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + + from meshmode.discretization.connection import ChainedDiscretizationConnection + chained = ChainedDiscretizationConnection(connections) + conn = chained.connections[1] el_to_batch = _build_element_lookup_table(queue, conn) - for igrp, grp in enumerate(conn.groups): for ibatch, batch in enumerate(grp.batches): for i, k in enumerate(batch.from_element_indices.get(queue)): @@ -159,15 +136,18 @@ def test_chained_new_group_table(ctx_factory, ndim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - chained = create_chained_connection(queue, ndim, - nelements=8, - mesh_order=2, - discr_order=2, - visualize=visualize) - connections = chained.connections + discr = create_discretization(queue, ndim, + nelements=8, + mesh_order=2, + discr_order=2) + connections = [] + conn = create_refined_connection(queue, discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + grp_to_grp, grp_info = _build_new_group_table(connections[0], connections[1]) - if visualize: import matplotlib.pyplot as pt @@ -206,10 +186,17 @@ def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - chained = create_chained_connection(queue, ndim, - mesh_order=2, discr_order=2, - visualize=visualize) - connections = chained.connections + discr = create_discretization(queue, ndim, + mesh_order=2, + discr_order=2) + connections = [] + conn = create_refined_connection(queue, discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + + from meshmode.discretization.connection import ChainedDiscretizationConnection + chained = ChainedDiscretizationConnection(connections) def f(x): from six.moves import reduce @@ -228,7 +215,7 @@ def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): @pytest.mark.parametrize("ndim", [2, 3]) -def test_chained_to_direct(ctx_factory, ndim, visualize=False): +def test_chained_to_direct(ctx_factory, ndim, chain_type=1, visualize=False): import time from meshmode.discretization.connection.chained import \ flatten_chained_connection @@ -236,17 +223,40 @@ def test_chained_to_direct(ctx_factory, ndim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - chained = create_chained_connection(queue, ndim, visualize=visualize) - connections = chained.connections + discr = create_discretization(queue, ndim) + connections = [] + if chain_type == 1: + conn = create_refined_connection(queue, discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + elif chain_type == 2: + conn = create_refined_connection(queue, discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + elif chain_type == 3 and ndim == 3: + conn = create_refined_connection(queue, discr, threshold=np.inf) + connections.append(conn) + conn = create_face_connection(queue, conn.to_discr) + connections.append(conn) + else: + raise ValueError('unknown test case') + + from meshmode.discretization.connection import ChainedDiscretizationConnection + chained = ChainedDiscretizationConnection(connections) direct = flatten_chained_connection(queue, chained) - to_element_indices = np.full(direct.to_discr.mesh.nelements, 0, - dtype=np.int) - for grp in direct.groups: - for batch in grp.batches: - for i in batch.to_element_indices.get(queue): - to_element_indices[i] += 1 - assert np.min(to_element_indices) > 0 + if chain_type < 3: + to_element_indices = np.full(direct.to_discr.mesh.nelements, 0, + dtype=np.int) + for grp in direct.groups: + for batch in grp.batches: + for i in batch.to_element_indices.get(queue): + to_element_indices[i] += 1 + assert np.min(to_element_indices) > 0 def f(x): from six.moves import reduce -- GitLab From e2f22fc819b836568104ca883a4e9674dc7e2578 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 28 May 2018 14:30:29 -0500 Subject: [PATCH 10/16] connection: fix docs a bit --- meshmode/discretization/connection/chained.py | 20 +++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/meshmode/discretization/connection/chained.py b/meshmode/discretization/connection/chained.py index ef926158..e557c58c 100644 --- a/meshmode/discretization/connection/chained.py +++ b/meshmode/discretization/connection/chained.py @@ -132,10 +132,11 @@ def flatten_chained_connection(queue, connection): The new direct connection will have a number of groups and batches that is, at worse, the product of all the connections in the chain. For example, if we consider a connection between a discretization and a - two-level refinement, both levels will have :math:`n` groups and + two-level refinement, both levels will have :math:`n` groups and :math:`m + 1` batches per group, where :math:`m` is the number of subdivisions of an element (exact number depends on implementation - details in :method:`meshmode.discretizationc.connection.make_refinement_connection`). + details in + :func:`~meshmode.discretizationc.connection.make_refinement_connection`). However, a direct connection from level :math:`0` to level :math:`2` will have at worst :math:`n^2` groups and each group will have :math:`(m + 1)^2` batches. @@ -145,9 +146,11 @@ def flatten_chained_connection(queue, connection): If a large number of connections is chained, the number of groups and batches can become very large. - :arg queue: a :class:`pyopencl.CommandQueue`. - :arg connection: a :class:`~meshmode.discretization.connection.DiscretizationConnection`. - :return: a :class:`~meshmode.discretization.connection.DirectDiscretizationConnection`. + :arg queue: An instance of :class:`pyopencl.CommandQueue`. + :arg connection: An instance of + :class:`~meshmode.discretization.connection.DiscretizationConnection`. + :return: An instance of + :class:`~meshmode.discretization.connection.DirectDiscretizationConnection`. """ from meshmode.discretization.connection import ( DirectDiscretizationConnection, @@ -220,18 +223,19 @@ def flatten_chained_connection(queue, connection): def make_full_resample_matrix(queue, connection): """Build a dense matrix representing the discretization connection. - This is based on + This is based on :func:`~meshmode.discretization.connection.DirectDiscretizationConnection.full_resample_matrix`. If a chained connection is given, the matrix is constructed recursively for each connection and multiplied left to right. .. warning:: - This method will be very slow, both in terms of speed and memory + This method will be very slow, both in terms of speed and memory usage, and should only be used for testing or if absolutely necessary. :arg queue: a :class:`pyopencl.CommandQueue`. - :arg connection: a :class:`~meshmode.discretization.connection.DiscretizationConnection`. + :arg connection: a + :class:`~meshmode.discretization.connection.DiscretizationConnection`. :return: a :class:`pyopencl.array.Array` of shape `(connection.from_discr.nnodes, connection.to_discr.nnodes)`. """ -- GitLab From 4d874d82562754d54f6eba731e0e14e50d5d7195 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 28 May 2018 14:31:15 -0500 Subject: [PATCH 11/16] flake8 fixes --- test/test_chained.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_chained.py b/test/test_chained.py index 32f09612..dbad0960 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -36,6 +36,7 @@ from pyopencl.tools import ( # noqa import logging logger = logging.getLogger(__name__) + def create_discretization(queue, ndim, nelements=42, mesh_order=5, @@ -165,7 +166,6 @@ def test_chained_new_group_table(ctx_factory, ndim, visualize=False): .groups[from_group_index].unit_nodes to_unit_nodes = grp_info[mgroup][mbatch].result_unit_nodes - if ndim == 2: pt.plot(from_unit_nodes, 'o') pt.plot(to_unit_nodes, '^') -- GitLab From 3ae6b9184dddc6bfc4b14094d4267103ffae0cec Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Mon, 28 May 2018 14:37:50 -0500 Subject: [PATCH 12/16] tests: enable all chained tests --- test/test_chained.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/test/test_chained.py b/test/test_chained.py index dbad0960..450a76ec 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -214,8 +214,9 @@ def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): assert np.allclose(f2, f3) -@pytest.mark.parametrize("ndim", [2, 3]) -def test_chained_to_direct(ctx_factory, ndim, chain_type=1, visualize=False): +@pytest.mark.parametrize(("ndim", "chain_type"), [ + (2, 1), (2, 2), (3, 1), (3, 3)]) +def test_chained_to_direct(ctx_factory, ndim, chain_type, visualize=False): import time from meshmode.discretization.connection.chained import \ flatten_chained_connection @@ -277,12 +278,6 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type=1, visualize=False): if visualize: print('[TIME] Chained: {:.5e}'.format(t_end - t_start)) - t_start = time.time() - f3 = connections[1](queue, connections[0](queue, fx)).get(queue) - t_end = time.time() - if visualize: - print('[TIME] Chained: {:.5e}'.format(t_end - t_start)) - if visualize and ndim == 2: import matplotlib.pyplot as pt @@ -295,7 +290,6 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type=1, visualize=False): pt.clf() assert np.allclose(f1, f2) - assert np.allclose(f2, f3) if __name__ == "__main__": -- GitLab From 0c1edde6f3f08ca3cb2551f70f3b613c46fce2a9 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Thu, 31 May 2018 10:18:30 -0500 Subject: [PATCH 13/16] style fix --- meshmode/discretization/connection/chained.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/meshmode/discretization/connection/chained.py b/meshmode/discretization/connection/chained.py index e557c58c..d89ccf3f 100644 --- a/meshmode/discretization/connection/chained.py +++ b/meshmode/discretization/connection/chained.py @@ -98,7 +98,7 @@ def _build_new_group_table(from_conn, to_conn): result_unit_nodes=result_unit_nodes, to_element_face=tbatch.to_element_face)) - grp_to_grp[(igrp, ibatch, jgrp, jbatch)] = (igrp_new, ibatch_new) + grp_to_grp[igrp, ibatch, jgrp, jbatch] = (igrp_new, ibatch_new) return grp_to_grp, batch_info @@ -182,7 +182,7 @@ def flatten_chained_connection(queue, connection): to_element_indices = {} for (igrp, ibatch), (_, tbatch) in _iterbatches(to_conn.groups): - to_element_indices[(igrp, ibatch)] = \ + to_element_indices[igrp, ibatch] = \ tbatch.to_element_indices.get(queue) # NOTE: notation used: @@ -194,8 +194,8 @@ def flatten_chained_connection(queue, connection): from_batch.to_element_indices.get(queue)): for j, jgrp, jbatch in el_to_batch[ito]: igrp_new, ibatch_new = \ - grp_to_grp[(igrp, ibatch, jgrp, jbatch)] - jto = to_element_indices[(jgrp, jbatch)][j] + grp_to_grp[igrp, ibatch, jgrp, jbatch] + jto = to_element_indices[jgrp, jbatch][j] from_bins[igrp_new][ibatch_new].append(ifrom) to_bins[igrp_new][ibatch_new].append(jto) -- GitLab From 89d2580a0c00a6baa5849df744c04303129585b6 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Thu, 31 May 2018 10:20:34 -0500 Subject: [PATCH 14/16] added a fixme note to improve slow python loops --- meshmode/discretization/connection/chained.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/meshmode/discretization/connection/chained.py b/meshmode/discretization/connection/chained.py index d89ccf3f..13a70e00 100644 --- a/meshmode/discretization/connection/chained.py +++ b/meshmode/discretization/connection/chained.py @@ -47,6 +47,7 @@ def _build_element_lookup_table(queue, conn): nelements = np.sum([g.nelements for g in conn.from_discr.groups]) el_table = [[] for _ in range(nelements)] + # FIXME: O(n) in number of elements loop in Python should be improved for (igrp, ibatch), (_, batch) in _iterbatches(conn.groups): for i, k in enumerate(batch.from_element_indices.get(queue)): el_table[k].append((i, igrp, ibatch)) @@ -189,6 +190,7 @@ def flatten_chained_connection(queue, connection): # * `ixxx` is an index in from_conn # * `jxxx` is an index in to_conn # * `ito` is the same as `jfrom` (on the same discr) + # FIXME: O(n) in number of elements loop in Python should be improved for (igrp, ibatch), (_, from_batch) in _iterbatches(from_conn.groups): for ifrom, ito in zip(from_batch.from_element_indices.get(queue), from_batch.to_element_indices.get(queue)): -- GitLab From eb0ace27fde304e7fb9b5680c5df3391d8d69fa6 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Thu, 31 May 2018 10:46:31 -0500 Subject: [PATCH 15/16] move chainedconnection test from test_meshmode to test_chained --- meshmode/discretization/connection/chained.py | 4 +- test/test_chained.py | 44 ++++++++++++-- test/test_meshmode.py | 59 ------------------- 3 files changed, 42 insertions(+), 65 deletions(-) diff --git a/meshmode/discretization/connection/chained.py b/meshmode/discretization/connection/chained.py index 13a70e00..13a098d3 100644 --- a/meshmode/discretization/connection/chained.py +++ b/meshmode/discretization/connection/chained.py @@ -158,10 +158,10 @@ def flatten_chained_connection(queue, connection): DiscretizationConnectionElementGroup, make_same_mesh_connection) - if isinstance(connection, DirectDiscretizationConnection): + if not hasattr(connection, 'connections'): return connection - if not hasattr(connection, 'connections') or not connection.connections: + if not connection.connections: return make_same_mesh_connection(connection.to_discr, connection.from_discr) diff --git a/test/test_chained.py b/test/test_chained.py index 450a76ec..112e1c6e 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -43,7 +43,7 @@ def create_discretization(queue, ndim, discr_order=5): ctx = queue.context - # construct base mesh + # construct mesh if ndim == 2: from functools import partial from meshmode.mesh.generation import make_curve_mesh, ellipse @@ -56,7 +56,7 @@ def create_discretization(queue, ndim, else: raise ValueError("Unsupported dimension: {}".format(ndim)) - # create base discretization + # create discretization from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ InterpolatoryQuadratureSimplexGroupFactory @@ -98,6 +98,7 @@ def create_face_connection(queue, discr): return connection +@pytest.mark.skip(reason='implementation detail') @pytest.mark.parametrize("ndim", [2, 3]) def test_chained_batch_table(ctx_factory, ndim, visualize=False): from meshmode.discretization.connection.chained import \ @@ -129,6 +130,7 @@ def test_chained_batch_table(ctx_factory, ndim, visualize=False): assert k == batch.from_element_indices[i] +@pytest.mark.skip(reason='implementation detail') @pytest.mark.parametrize("ndim", [2, 3]) def test_chained_new_group_table(ctx_factory, ndim, visualize=False): from meshmode.discretization.connection.chained import \ @@ -178,6 +180,38 @@ def test_chained_new_group_table(ctx_factory, ndim, visualize=False): pt.clf() +@pytest.mark.parametrize("ndim", [2, 3]) +def test_chained_connection(ctx_factory, ndim, visualize=False): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + discr = create_discretization(queue, ndim, + nelements=10, + mesh_order=5, + discr_order=5) + connections = [] + conn = create_refined_connection(queue, discr, threshold=np.inf) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr, threshold=np.inf) + connections.append(conn) + + from meshmode.discretization.connection import \ + ChainedDiscretizationConnection + chained = ChainedDiscretizationConnection(connections) + + def f(x): + from six.moves import reduce + return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) + + x = connections[0].from_discr.nodes().with_queue(queue) + fx = f(x) + f1 = chained(queue, fx).get(queue) + f2 = connections[1](queue, connections[0](queue, fx)).get(queue) + + assert np.allclose(f1, f2) + + +@pytest.mark.slowtest @pytest.mark.parametrize("ndim", [2, 3]) def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): from meshmode.discretization.connection.chained import \ @@ -195,7 +229,8 @@ def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): conn = create_refined_connection(queue, conn.to_discr) connections.append(conn) - from meshmode.discretization.connection import ChainedDiscretizationConnection + from meshmode.discretization.connection import \ + ChainedDiscretizationConnection chained = ChainedDiscretizationConnection(connections) def f(x): @@ -246,7 +281,8 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type, visualize=False): else: raise ValueError('unknown test case') - from meshmode.discretization.connection import ChainedDiscretizationConnection + from meshmode.discretization.connection import \ + ChainedDiscretizationConnection chained = ChainedDiscretizationConnection(connections) direct = flatten_chained_connection(queue, chained) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 7173fc29..29dc2ed1 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1039,65 +1039,6 @@ def test_quad_multi_element(): plt.show() -# {{{ ChainedDiscretizationConnection - -def test_ChainedDiscretizationConnection(ctx_getter): # noqa - mesh_order = 5 - order = 5 - npanels = 10 - group_factory = InterpolatoryQuadratureSimplexGroupFactory - - def refine_flags(mesh): - return np.ones(mesh.nelements) - - cl_ctx = ctx_getter() - queue = cl.CommandQueue(cl_ctx) - - from functools import partial - from meshmode.discretization import Discretization - from meshmode.discretization.connection import make_refinement_connection - from meshmode.mesh.generation import make_curve_mesh, ellipse - - mesh = make_curve_mesh( - partial(ellipse, 1), np.linspace(0, 1, npanels + 1), - order=mesh_order) - - discr = Discretization(cl_ctx, mesh, group_factory(order)) - - connections = [] - - from meshmode.mesh.refinement import Refiner - refiner = Refiner(mesh) - - def refine_discr(discr): - mesh = discr.mesh - flags = refine_flags(mesh) - refiner.refine(flags) - connections.append( - make_refinement_connection(refiner, discr, group_factory(order))) - return connections[-1].to_discr - - discr = refine_discr(discr) - refine_discr(discr) - - from meshmode.discretization.connection import ( - ChainedDiscretizationConnection) - - chained_conn = ChainedDiscretizationConnection(connections) - - def f(x): - from six.moves import reduce - return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) - - x = connections[0].from_discr.nodes().with_queue(queue) - - assert np.allclose( - chained_conn(queue, f(x)).get(queue), - connections[1](queue, connections[0](queue, f(x))).get(queue)) - -# }}} - - if __name__ == "__main__": import sys if len(sys.argv) > 1: -- GitLab From e16ddcfaa4ed6331627686595aeb84d6f2d760dc Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Thu, 31 May 2018 16:19:23 -0500 Subject: [PATCH 16/16] chained: vectorize building some of the element tables --- meshmode/discretization/connection/chained.py | 53 +++++++++---------- test/test_chained.py | 26 +++++---- 2 files changed, 39 insertions(+), 40 deletions(-) diff --git a/meshmode/discretization/connection/chained.py b/meshmode/discretization/connection/chained.py index 13a098d3..221ce2d3 100644 --- a/meshmode/discretization/connection/chained.py +++ b/meshmode/discretization/connection/chained.py @@ -44,13 +44,12 @@ def _iterbatches(groups): def _build_element_lookup_table(queue, conn): - nelements = np.sum([g.nelements for g in conn.from_discr.groups]) - el_table = [[] for _ in range(nelements)] + el_table = [np.full(g.nelements, -1, dtype=np.int) + for g in conn.to_discr.groups] - # FIXME: O(n) in number of elements loop in Python should be improved - for (igrp, ibatch), (_, batch) in _iterbatches(conn.groups): - for i, k in enumerate(batch.from_element_indices.get(queue)): - el_table[k].append((i, igrp, ibatch)) + for (igrp, _), (_, batch) in _iterbatches(conn.groups): + el_table[igrp][batch.to_element_indices.get(queue)] = \ + batch.from_element_indices.get(queue) return el_table @@ -174,33 +173,29 @@ def flatten_chained_connection(queue, connection): # merge all the direct connections from_conn = direct_connections[0] for to_conn in direct_connections[1:]: - el_to_batch = _build_element_lookup_table(queue, to_conn) + el_table = _build_element_lookup_table(queue, from_conn) grp_to_grp, batch_info = _build_new_group_table(from_conn, to_conn) # distribute the indices to new groups and batches - from_bins = [[[] for _ in g] for g in batch_info] - to_bins = [[[] for _ in g] for g in batch_info] - - to_element_indices = {} - for (igrp, ibatch), (_, tbatch) in _iterbatches(to_conn.groups): - to_element_indices[igrp, ibatch] = \ - tbatch.to_element_indices.get(queue) - - # NOTE: notation used: - # * `ixxx` is an index in from_conn - # * `jxxx` is an index in to_conn - # * `ito` is the same as `jfrom` (on the same discr) - # FIXME: O(n) in number of elements loop in Python should be improved + from_bins = [[np.empty(0, dtype=np.int) for _ in g] for g in batch_info] + to_bins = [[np.empty(0, dtype=np.int) for _ in g] for g in batch_info] + for (igrp, ibatch), (_, from_batch) in _iterbatches(from_conn.groups): - for ifrom, ito in zip(from_batch.from_element_indices.get(queue), - from_batch.to_element_indices.get(queue)): - for j, jgrp, jbatch in el_to_batch[ito]: - igrp_new, ibatch_new = \ - grp_to_grp[igrp, ibatch, jgrp, jbatch] - jto = to_element_indices[jgrp, jbatch][j] - - from_bins[igrp_new][ibatch_new].append(ifrom) - to_bins[igrp_new][ibatch_new].append(jto) + from_to_element_indices = from_batch.to_element_indices.get(queue) + + for (jgrp, jbatch), (_, to_batch) in _iterbatches(to_conn.groups): + igrp_new, ibatch_new = grp_to_grp[igrp, ibatch, jgrp, jbatch] + + jfrom = to_batch.from_element_indices.get(queue) + jto = to_batch.to_element_indices.get(queue) + + mask = np.isin(jfrom, from_to_element_indices) + from_bins[igrp_new][ibatch_new] = \ + np.hstack([from_bins[igrp_new][ibatch_new], + el_table[igrp][jfrom[mask]]]) + to_bins[igrp_new][ibatch_new] = \ + np.hstack([to_bins[igrp_new][ibatch_new], + jto[mask]]) # build new groups groups = [] diff --git a/test/test_chained.py b/test/test_chained.py index 112e1c6e..c0f5a09e 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -117,17 +117,15 @@ def test_chained_batch_table(ctx_factory, ndim, visualize=False): from meshmode.discretization.connection import ChainedDiscretizationConnection chained = ChainedDiscretizationConnection(connections) - conn = chained.connections[1] - el_to_batch = _build_element_lookup_table(queue, conn) + conn = chained.connections[0] + el_table = _build_element_lookup_table(queue, conn) for igrp, grp in enumerate(conn.groups): for ibatch, batch in enumerate(grp.batches): - for i, k in enumerate(batch.from_element_indices.get(queue)): - assert (i, igrp, ibatch) in el_to_batch[k] + ifrom = batch.from_element_indices.get(queue) + jfrom = el_table[igrp][batch.to_element_indices.get(queue)] - for k, batches in enumerate(el_to_batch): - for i, igrp, ibatch in batches: - batch = conn.groups[igrp].batches[ibatch] - assert k == batch.from_element_indices[i] + assert np.all(ifrom == jfrom) + assert np.min(el_table[igrp]) >= 0 @pytest.mark.skip(reason='implementation detail') @@ -211,7 +209,7 @@ def test_chained_connection(ctx_factory, ndim, visualize=False): assert np.allclose(f1, f2) -@pytest.mark.slowtest +@pytest.mark.skip(reason='slow test') @pytest.mark.parametrize("ndim", [2, 3]) def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): from meshmode.discretization.connection.chained import \ @@ -251,7 +249,8 @@ def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): @pytest.mark.parametrize(("ndim", "chain_type"), [ (2, 1), (2, 2), (3, 1), (3, 3)]) -def test_chained_to_direct(ctx_factory, ndim, chain_type, visualize=False): +def test_chained_to_direct(ctx_factory, ndim, chain_type, + nelements=128, visualize=False): import time from meshmode.discretization.connection.chained import \ flatten_chained_connection @@ -259,7 +258,7 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - discr = create_discretization(queue, ndim) + discr = create_discretization(queue, ndim, nelements=nelements) connections = [] if chain_type == 1: conn = create_refined_connection(queue, discr) @@ -284,7 +283,12 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type, visualize=False): from meshmode.discretization.connection import \ ChainedDiscretizationConnection chained = ChainedDiscretizationConnection(connections) + + t_start = time.time() direct = flatten_chained_connection(queue, chained) + t_end = time.time() + if visualize: + print('[TIME] Flatten: {:.5e}'.format(t_end - t_start)) if chain_type < 3: to_element_indices = np.full(direct.to_discr.mesh.nelements, 0, -- GitLab