From bef920fc674e7fc8fbeb03cd29e510835d848a17 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 5 May 2019 20:51:56 -0500 Subject: [PATCH 1/3] separate connections into multiple files --- .../discretization/connection/__init__.py | 614 +----------------- meshmode/discretization/connection/chained.py | 45 +- meshmode/discretization/connection/direct.py | 412 ++++++++++++ .../discretization/connection/projection.py | 223 +++++++ meshmode/mesh/__init__.py | 1 - 5 files changed, 689 insertions(+), 606 deletions(-) create mode 100644 meshmode/discretization/connection/direct.py create mode 100644 meshmode/discretization/connection/projection.py diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index cc0417e6..c59511b7 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -25,14 +25,19 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -from six.moves import range, zip - import numpy as np import pyopencl as cl import pyopencl.array # noqa -from loopy.version import MOST_RECENT_LANGUAGE_VERSION -from pytools import memoize_method, memoize_in +from meshmode.discretization.connection.direct import ( + InterpolationBatch, + DiscretizationConnectionElementGroup, + DiscretizationConnection, + DirectDiscretizationConnection) +from meshmode.discretization.connection.chained import \ + ChainedDiscretizationConnection +from meshmode.discretization.connection.projection import \ + L2ProjectionInverseDiscretizationConnection from meshmode.discretization.connection.same_mesh import \ make_same_mesh_connection @@ -106,605 +111,6 @@ Implementation details .. autoclass:: DiscretizationConnectionElementGroup """ - -# {{{ interpolation batch - -class InterpolationBatch(object): - """One interpolation batch captures how a batch of elements *within* an - element group should be an interpolated. Note that while it's possible that - an interpolation batch takes care of interpolating an entire element group - from source to target, that's not *necessarily* the case. Consider the case - of extracting boundary values of a discretization. For, say, a triangle, at - least three different interpolation batches are needed to cover boundary - edges that fall onto each of the three edges of the unit triangle. - - .. attribute:: from_group_index - - An integer indicating from which element group in the *from* discretization - the data should be interpolated. - - .. attribute:: from_element_indices - - ``element_id_t [nelements]``. (a :class:`pyopencl.array.Array`) - This contains the (group-local) element index (relative to - :attr:`from_group_index` from which this "*to*" element's data will be - interpolated. - - .. attribute:: to_element_indices - - ``element_id_t [nelements]``. (a :class:`pyopencl.array.Array`) - This contains the (group-local) element index to which this "*to*" - element's data will be interpolated. - - .. attribute:: result_unit_nodes - - A :class:`numpy.ndarray` of shape - ``(from_group.dim,to_group.nunit_nodes)`` - storing the coordinates of the nodes (in unit coordinates - of the *from* reference element) from which the node - locations of this element should be interpolated. - - .. autoattribute:: nelements - - .. attribute:: to_element_face - - *int* or *None*. (a :class:`pyopencl.array.Array` if existent) If this - interpolation batch targets interpolation *to* a face, then this number - captures the face number (on all elements referenced by - :attr:`from_element_indices` to which this batch interpolates. (Since - there is a fixed set of "from" unit nodes per batch, one batch will - always go to a single face index.) - """ - - def __init__(self, from_group_index, from_element_indices, - to_element_indices, result_unit_nodes, to_element_face): - self.from_group_index = from_group_index - self.from_element_indices = from_element_indices - self.to_element_indices = to_element_indices - self.result_unit_nodes = result_unit_nodes - self.to_element_face = to_element_face - - @property - def nelements(self): - return len(self.from_element_indices) - -# }}} - - -# {{{ connection element group - -class DiscretizationConnectionElementGroup(object): - """ - .. attribute:: batches - - A list of :class:`InterpolationBatch` instances. - """ - def __init__(self, batches): - self.batches = batches - -# }}} - - -# {{{ connection classes - -class DiscretizationConnection(object): - """Abstract interface for transporting a DOF vector from one - :class:`meshmode.discretization.Discretization` to another. - Possible applications include: - - * upsampling/downsampling on the same mesh - * restricition to the boundary - * interpolation to a refined/coarsened mesh - * interpolation onto opposing faces - - .. attribute:: from_discr - - .. attribute:: to_discr - - .. attribute:: is_surjective - - A :class:`bool` indicating whether every output degree - of freedom is set by the connection. - - .. automethod:: __call__ - """ - def __init__(self, from_discr, to_discr, is_surjective): - if from_discr.cl_context != to_discr.cl_context: - raise ValueError("from_discr and to_discr must live in the " - "same OpenCL context") - - self.cl_context = from_discr.cl_context - - if from_discr.mesh.vertex_id_dtype != to_discr.mesh.vertex_id_dtype: - raise ValueError("from_discr and to_discr must agree on the " - "vertex_id_dtype") - - if from_discr.mesh.element_id_dtype != to_discr.mesh.element_id_dtype: - raise ValueError("from_discr and to_discr must agree on the " - "element_id_dtype") - - self.from_discr = from_discr - self.to_discr = to_discr - - self.is_surjective = is_surjective - - def __call__(self, queue, vec): - raise NotImplementedError() - - -class ChainedDiscretizationConnection(DiscretizationConnection): - """Aggregates multiple :class:`DiscretizationConnection` instances - into a single one. - - .. attribute:: connections - """ - - def __init__(self, connections, from_discr=None): - if connections: - if from_discr is not None: - assert from_discr is connections[0].from_discr - else: - from_discr = connections[0].from_discr - is_surjective = all(cnx.is_surjective for cnx in connections) - to_discr = connections[-1].to_discr - else: - if from_discr is None: - raise ValueError("connections may not be empty if from_discr " - "is not specified") - - to_discr = from_discr - - # It's an identity - is_surjective = True - - super(ChainedDiscretizationConnection, self).__init__( - from_discr, to_discr, is_surjective=is_surjective) - - self.connections = connections - - def __call__(self, queue, vec): - for cnx in self.connections: - vec = cnx(queue, vec) - - return vec - - -class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): - """Creates an inverse :class:`DiscretizationConnection` from an existing - connection to allow transporting from the original connection's - *to_discr* to *from_discr*. - - .. attribute:: from_discr - .. attribute:: to_discr - .. attribute:: is_surjective - - .. attribute:: conn - .. automethod:: __call__ - - """ - - def __new__(cls, connections, is_surjective=False): - if isinstance(connections, DirectDiscretizationConnection): - return DiscretizationConnection.__new__(cls) - elif isinstance(connections, ChainedDiscretizationConnection): - return cls(connections.connections, is_surjective=is_surjective) - else: - conns = [] - for cnx in reversed(connections): - conns.append(cls(cnx, is_surjective=is_surjective)) - - return ChainedDiscretizationConnection(conns) - - def __init__(self, conn, is_surjective=False): - if conn.from_discr.dim != conn.to_discr.dim: - raise RuntimeError("cannot transport from face to element") - - if not all(g.is_orthogonal_basis() for g in conn.to_discr.groups): - raise RuntimeError("`to_discr` must have an orthogonal basis") - - self.conn = conn - super(L2ProjectionInverseDiscretizationConnection, self).__init__( - from_discr=self.conn.to_discr, - to_discr=self.conn.from_discr, - is_surjective=is_surjective) - - @memoize_method - def _batch_weights(self): - """Computes scaled quadrature weights for each interpolation batch in - :attr:`conn`. The quadrature weights can be used to integrate over - child elements in the domain of the parent element, by a change of - variables. - - :return: a dictionary with keys ``(group_id, batch_id)``. - """ - - from pymbolic.geometric_algebra import MultiVector - from functools import reduce - from operator import xor - - def det(v): - nnodes = v[0].shape[0] - det_v = np.empty(nnodes) - - for i in range(nnodes): - outer_product = reduce(xor, [MultiVector(x[i, :].T) for x in v]) - det_v[i] = abs((outer_product.I | outer_product).as_scalar()) - - return det_v - - weights = {} - jac = np.empty(self.to_discr.dim, dtype=np.object) - - for igrp, grp in enumerate(self.to_discr.groups): - for ibatch, batch in enumerate(self.conn.groups[igrp].batches): - for iaxis in range(grp.dim): - mat = grp.diff_matrices()[iaxis] - jac[iaxis] = mat.dot(batch.result_unit_nodes.T) - - weights[igrp, ibatch] = det(jac) * grp.weights - - return weights - - def __call__(self, queue, vec): - @memoize_in(self, "conn_projection_knl") - def kproj(): - import loopy as lp - knl = lp.make_kernel([ - "{[k]: 0 <= k < nelements}", - "{[j]: 0 <= j < n_from_nodes}" - ], - """ - for k - <> element_dot = \ - sum(j, vec[from_element_indices[k], j] * \ - basis[j] * weights[j]) - - result[to_element_indices[k], ibasis] = \ - result[to_element_indices[k], ibasis] + element_dot - end - """, - [ - lp.GlobalArg("vec", None, - shape=("n_from_elements", "n_from_nodes")), - lp.GlobalArg("result", None, - shape=("n_to_elements", "n_to_nodes")), - lp.GlobalArg("basis", None, - shape="n_from_nodes"), - lp.GlobalArg("weights", None, - shape="n_from_nodes"), - lp.ValueArg("n_from_elements", np.int32), - lp.ValueArg("n_to_elements", np.int32), - lp.ValueArg("n_to_nodes", np.int32), - lp.ValueArg("ibasis", np.int32), - '...' - ], - name="conn_projection_knl", - lang_version=MOST_RECENT_LANGUAGE_VERSION) - - return knl - - @memoize_in(self, "conn_evaluation_knl") - def keval(): - import loopy as lp - knl = lp.make_kernel([ - "{[k]: 0 <= k < nelements}", - "{[j]: 0 <= j < n_to_nodes}" - ], - """ - result[k, j] = result[k, j] + \ - coefficients[k, ibasis] * basis[j] - """, - [ - lp.GlobalArg("coefficients", None, - shape=("nelements", "n_to_nodes")), - lp.ValueArg("ibasis", np.int32), - '...' - ], - name="conn_evaluate_knl", - default_offset=lp.auto, - lang_version=MOST_RECENT_LANGUAGE_VERSION) - - return knl - - if not isinstance(vec, cl.array.Array): - raise TypeError("non-array passed to discretization connection") - - if vec.shape != (self.from_discr.nnodes,): - raise ValueError("invalid shape of incoming resampling data") - - # compute weights on each refinement of the reference element - weights = self._batch_weights() - - # perform dot product (on reference element) to get basis coefficients - c = self.to_discr.zeros(queue, dtype=vec.dtype) - for igrp, (tgrp, cgrp) in enumerate( - zip(self.to_discr.groups, self.conn.groups)): - for ibatch, batch in enumerate(cgrp.batches): - sgrp = self.from_discr.groups[batch.from_group_index] - - for ibasis, basis_fn in enumerate(sgrp.basis()): - basis = basis_fn(batch.result_unit_nodes).flatten() - - # NOTE: batch.*_element_indices are reversed here because - # they are from the original forward connection, but - # we are going in reverse here. a bit confusing, but - # saves on recreating the connection groups and batches. - kproj()(queue, - ibasis=ibasis, - vec=sgrp.view(vec), - basis=basis, - weights=weights[igrp, ibatch], - result=tgrp.view(c), - from_element_indices=batch.to_element_indices, - to_element_indices=batch.from_element_indices) - - # evaluate at unit_nodes to get the vector on to_discr - result = self.to_discr.zeros(queue, dtype=vec.dtype) - for igrp, grp in enumerate(self.to_discr.groups): - for ibasis, basis_fn in enumerate(grp.basis()): - basis = basis_fn(grp.unit_nodes).flatten() - - keval()(queue, - ibasis=ibasis, - result=grp.view(result), - basis=basis, - coefficients=grp.view(c)) - - return result - - -class DirectDiscretizationConnection(DiscretizationConnection): - """A concrete :class:`DiscretizationConnection` supported by interpolation - data. - - .. attribute:: from_discr - - .. attribute:: to_discr - - .. attribute:: groups - - a list of :class:`DiscretizationConnectionElementGroup` - instances, with a one-to-one correspondence to the groups in - :attr:`to_discr`. - - .. attribute:: is_surjective - - A :class:`bool` indicating whether every output degree - of freedom is set by the connection. - - .. automethod:: __call__ - - .. automethod:: full_resample_matrix - - """ - - def __init__(self, from_discr, to_discr, groups, is_surjective): - super(DirectDiscretizationConnection, self).__init__( - from_discr, to_discr, is_surjective) - - self.groups = groups - - @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] - - nfrom_unit_nodes = from_grp.unit_nodes.shape[1] - if np.array_equal(from_grp.unit_nodes, ibatch.result_unit_nodes): - # Nodes are exactly identical? We can 'interpolate' even when there - # isn't a basis. - - result = np.eye(nfrom_unit_nodes) - - else: - if len(from_grp.basis()) != nfrom_unit_nodes: - from meshmode.discretization import NoninterpolatoryElementGroupError - raise NoninterpolatoryElementGroupError( - "%s does not support interpolation because it is not " - "unisolvent (its unit node count does not match its " - "number of basis functions). Using connections requires " - "the ability to interpolate." % type(from_grp).__name__) - - result = mp.resampling_matrix( - from_grp.basis(), - ibatch.result_unit_nodes, from_grp.unit_nodes) - - with cl.CommandQueue(self.cl_context) as queue: - return cl.array.to_device(queue, result).with_queue(None) - - @memoize_method - def _resample_point_pick_indices(self, to_group_index, ibatch_index, - tol_multiplier=None): - """If :meth:`_resample_matrix` *R* is a row subset of a permutation matrix *P*, - return the index subset I so that, loosely, ``x[I] == R @ x``. - - Will return *None* if no such index array exists, or a - :class:`pyopencl.array.Array` containing the index subset. - """ - - with cl.CommandQueue(self.cl_context) as queue: - mat = self._resample_matrix(to_group_index, ibatch_index).get( - queue=queue) - - nrows, ncols = mat.shape - result = np.zeros(nrows, dtype=self.to_discr.mesh.element_id_dtype) - - if tol_multiplier is None: - tol_multiplier = 50 - - tol = np.finfo(mat.dtype).eps * tol_multiplier - - for irow in range(nrows): - one_indices, = np.where(np.abs(mat[irow] - 1) < tol) - zero_indices, = np.where(np.abs(mat[irow]) < tol) - - if len(one_indices) != 1: - return None - if len(zero_indices) != ncols - 1: - return None - - one_index, = one_indices - result[irow] = one_index - - with cl.CommandQueue(self.cl_context) as queue: - return cl.array.to_device(queue, result).with_queue(None) - - def full_resample_matrix(self, queue): - """Build a dense matrix representing this discretization connection. - - .. warning:: - - On average, this will be exceedingly expensive (:math:`O(N^2)` in - the number *N* of discretization points) in terms of memory usage - and thus not what you'd typically want. - """ - - @memoize_in(self, "oversample_mat_knl") - def knl(): - import loopy as lp - knl = lp.make_kernel( - """{[k,i,j]: - 0<=k element_dot = \ + sum(j, vec[from_element_indices[k], j] * \ + basis[j] * weights[j]) + + result[to_element_indices[k], ibasis] = \ + result[to_element_indices[k], ibasis] + element_dot + end + """, + [ + lp.GlobalArg("vec", None, + shape=("n_from_elements", "n_from_nodes")), + lp.GlobalArg("result", None, + shape=("n_to_elements", "n_to_nodes")), + lp.GlobalArg("basis", None, + shape="n_from_nodes"), + lp.GlobalArg("weights", None, + shape="n_from_nodes"), + lp.ValueArg("n_from_elements", np.int32), + lp.ValueArg("n_to_elements", np.int32), + lp.ValueArg("n_to_nodes", np.int32), + lp.ValueArg("ibasis", np.int32), + '...' + ], + name="conn_projection_knl", + lang_version=MOST_RECENT_LANGUAGE_VERSION) + + return knl + + @memoize_in(self, "conn_evaluation_knl") + def keval(): + import loopy as lp + knl = lp.make_kernel([ + "{[k]: 0 <= k < nelements}", + "{[j]: 0 <= j < n_to_nodes}" + ], + """ + result[k, j] = result[k, j] + \ + coefficients[k, ibasis] * basis[j] + """, + [ + lp.GlobalArg("coefficients", None, + shape=("nelements", "n_to_nodes")), + lp.ValueArg("ibasis", np.int32), + '...' + ], + name="conn_evaluate_knl", + default_offset=lp.auto, + lang_version=MOST_RECENT_LANGUAGE_VERSION) + + return knl + + if not isinstance(vec, cl.array.Array): + raise TypeError("non-array passed to discretization connection") + + if vec.shape != (self.from_discr.nnodes,): + raise ValueError("invalid shape of incoming resampling data") + + # compute weights on each refinement of the reference element + weights = self._batch_weights() + + # perform dot product (on reference element) to get basis coefficients + c = self.to_discr.zeros(queue, dtype=vec.dtype) + for igrp, (tgrp, cgrp) in enumerate( + zip(self.to_discr.groups, self.conn.groups)): + for ibatch, batch in enumerate(cgrp.batches): + sgrp = self.from_discr.groups[batch.from_group_index] + + for ibasis, basis_fn in enumerate(sgrp.basis()): + basis = basis_fn(batch.result_unit_nodes).flatten() + + # NOTE: batch.*_element_indices are reversed here because + # they are from the original forward connection, but + # we are going in reverse here. a bit confusing, but + # saves on recreating the connection groups and batches. + kproj()(queue, + ibasis=ibasis, + vec=sgrp.view(vec), + basis=basis, + weights=weights[igrp, ibatch], + result=tgrp.view(c), + from_element_indices=batch.to_element_indices, + to_element_indices=batch.from_element_indices) + + # evaluate at unit_nodes to get the vector on to_discr + result = self.to_discr.zeros(queue, dtype=vec.dtype) + for igrp, grp in enumerate(self.to_discr.groups): + for ibasis, basis_fn in enumerate(grp.basis()): + basis = basis_fn(grp.unit_nodes).flatten() + + keval()(queue, + ibasis=ibasis, + result=grp.view(result), + basis=basis, + coefficients=grp.view(c)) + + return result + + +# vim: foldmethod=marker diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 2c799478..6d6191aa 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -640,7 +640,6 @@ class Mesh(Record): .. automethod:: __eq__ .. automethod:: __ne__ - .. automethos:: adjacency_list """ face_id_dtype = np.int8 -- GitLab From a7101d0e21f8e670a6aaa43d258c05d0127453e6 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 5 May 2019 20:55:48 -0500 Subject: [PATCH 2/3] flake8 issues --- meshmode/discretization/connection/__init__.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index c59511b7..5942c7e2 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -25,7 +25,6 @@ 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 @@ -58,6 +57,10 @@ logger = logging.getLogger(__name__) __all__ = [ "DiscretizationConnection", + "DirectDiscretizationConnection", + "ChainedDiscretizationConnection", + "L2ProjectionInverseDiscretizationConnection", + "make_same_mesh_connection", "FACE_RESTR_INTERIOR", "FACE_RESTR_ALL", "make_face_restriction", @@ -66,6 +69,9 @@ __all__ = [ "make_partition_connection", "make_refinement_connection", "flatten_chained_connection", + + "InterpolationBatch", + "DiscretizationConnectionElementGroup", ] __doc__ = """ @@ -111,6 +117,7 @@ Implementation details .. autoclass:: DiscretizationConnectionElementGroup """ + # {{{ check connection def check_connection(connection): -- GitLab From bd46136a7fb515c592d8899f9c37c0e892170691 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 5 May 2019 20:59:58 -0500 Subject: [PATCH 3/3] import directly from file with base classes --- meshmode/discretization/connection/chained.py | 2 +- meshmode/discretization/connection/face.py | 5 +++-- meshmode/discretization/connection/opposite_face.py | 2 +- meshmode/discretization/connection/refinement.py | 2 +- meshmode/discretization/connection/same_mesh.py | 5 +++-- 5 files changed, 9 insertions(+), 7 deletions(-) diff --git a/meshmode/discretization/connection/chained.py b/meshmode/discretization/connection/chained.py index 33457d5b..bd1da1f4 100644 --- a/meshmode/discretization/connection/chained.py +++ b/meshmode/discretization/connection/chained.py @@ -147,7 +147,7 @@ def _build_new_group_table(from_conn, to_conn): def _build_batches(queue, from_bins, to_bins, batch): - from meshmode.discretization.connection import \ + from meshmode.discretization.connection.direct import \ InterpolationBatch def to_device(x): diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py index 6016a146..a0d4e71e 100644 --- a/meshmode/discretization/connection/face.py +++ b/meshmode/discretization/connection/face.py @@ -65,8 +65,9 @@ class _ConnectionBatchData(Record): def _build_boundary_connection(queue, vol_discr, bdry_discr, connection_data, per_face_groups): - from meshmode.discretization.connection import ( - InterpolationBatch, DiscretizationConnectionElementGroup, + from meshmode.discretization.connection.direct import ( + InterpolationBatch, + DiscretizationConnectionElementGroup, DirectDiscretizationConnection) ibdry_grp = 0 diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py index ebb51526..46db90b7 100644 --- a/meshmode/discretization/connection/opposite_face.py +++ b/meshmode/discretization/connection/opposite_face.py @@ -246,7 +246,7 @@ def _make_cross_face_batches(queue, tgt_bdry_discr, src_bdry_discr, - template_unit_nodes.reshape(dim, 1, -1)), axis=2), axis=0) - from meshmode.discretization.connection import InterpolationBatch + from meshmode.discretization.connection.direct import InterpolationBatch yield InterpolationBatch( from_group_index=i_src_grp, from_element_indices=to_dev(src_bdry_element_indices[close_els]), diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index fb06dc5d..9a2d60d8 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -72,7 +72,7 @@ def _build_interpolation_batches_for_group( batches because their unit nodes are mapped from the same part of the reference element. """ - from meshmode.discretization.connection import InterpolationBatch + from meshmode.discretization.connection.direct import InterpolationBatch num_children = len(record.tesselation.children) \ if record.tesselation else 0 diff --git a/meshmode/discretization/connection/same_mesh.py b/meshmode/discretization/connection/same_mesh.py index a470e929..e7b781cd 100644 --- a/meshmode/discretization/connection/same_mesh.py +++ b/meshmode/discretization/connection/same_mesh.py @@ -30,8 +30,9 @@ import pyopencl.array # noqa # {{{ same-mesh constructor def make_same_mesh_connection(to_discr, from_discr): - from meshmode.discretization.connection import ( - InterpolationBatch, DiscretizationConnectionElementGroup, + from meshmode.discretization.connection.direct import ( + InterpolationBatch, + DiscretizationConnectionElementGroup, DirectDiscretizationConnection) if from_discr.mesh is not to_discr.mesh: -- GitLab