From 09f85e059976d1681c9f6095ec193940b2c9ecdf Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sun, 25 Nov 2018 22:48:27 -0600 Subject: [PATCH 01/19] connection: add a reverse connection --- .../discretization/connection/__init__.py | 100 ++++++++++++++++++ test/test_chained.py | 66 +++++++++++- 2 files changed, 164 insertions(+), 2 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 5f136982..fab24c9d 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -251,6 +251,106 @@ class ChainedDiscretizationConnection(DiscretizationConnection): return vec +class ReversedChainedDiscretizationConnection(DiscretizationConnection): + def __init__(self, chained_connection): + self.conn = chained_connection + if len(self.conn.to_discr.groups) != 1 \ + or len(self.conn.from_discr.groups) != 1: + raise RuntimeError("multiple element groups are not supported") + + super(ReversedChainedDiscretizationConnection, self).__init__( + from_discr=self.conn.to_discr, + to_discr=self.conn.from_discr, + is_surjective=False) + + @memoize_method + def _evaluate_basis(self, queue, group, ibasis): + basis_fn = group.basis()[ibasis] + + # evaluate + basis = np.zeros(self.to_discr.nnodes) + view = group.view(basis) + for k in range(view.shape[0]): + view[k, :] = basis_fn(group.unit_nodes) + + basis = cl.array.to_device(queue, basis) + basis = self.conn(queue, basis).with_queue(None) + + return basis + + @memoize_method + def _build_element_index_map(self, queue, group): + # evaluate element indices on source + indices = np.full(self.to_discr.nnodes, -1.0) + view = group.view(indices) + for k in range(view.shape[0]): + view[k, :] = k + + # interpolate to target + indices = self.conn(queue, cl.array.to_device(queue, indices)) + indices = cl.clmath.round(indices.with_queue(queue)).astype(np.int) + + return indices.with_queue(None) + + def __call__(self, queue, vec): + @memoize_in(self, "conn_projection_kernel") + def knl(): + import loopy as lp + knl = lp.make_kernel([ + "{[k]: 0 <= k < nsources}", + "{[j]: 0 <= j < n_from_nodes}" + ], + """ + for k, j + result[indices[k, j], ibasis] = \ + result[indices[k, j], ibasis] + \ + vec[k, j] * basis[k, j] * weights[j] + end + """, + [ + lp.GlobalArg("vec", None, + shape=("nsources", "n_from_nodes")), + lp.GlobalArg("result", None, + shape=("ntargets", "n_to_nodes")), + lp.GlobalArg("indices", None, + shape=("nsources", "n_from_nodes")), + lp.GlobalArg("basis", None, + shape=("nsources", "n_from_nodes")), + lp.GlobalArg("weights", None, + shape="n_from_nodes"), + lp.ValueArg("ibasis", np.int32), + lp.ValueArg("n_to_nodes", np.int32), + lp.ValueArg("ntargets", np.int32), + '...' + ], + name="conn_projection_kernel", + 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") + + result = self.to_discr.zeros(queue, dtype=vec.dtype) + for tgrp in self.to_discr.groups: + el_indices = self._build_element_index_map(queue, tgrp) + for ibasis in range(len(tgrp.basis())): + basis = self._evaluate_basis(queue, tgrp, ibasis) + + for sgrp in self.from_discr.groups: + knl()(queue, + weights=sgrp.weights, ibasis=ibasis, + vec=sgrp.view(vec), + basis=sgrp.view(basis), + indices=sgrp.view(el_indices), + result=tgrp.view(result)) + + return result + + class DirectDiscretizationConnection(DiscretizationConnection): """A concrete :class:`DiscretizationConnection` supported by interpolation data. diff --git a/test/test_chained.py b/test/test_chained.py index c0f5a09e..f516ecb1 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -321,17 +321,79 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type, if visualize and ndim == 2: import matplotlib.pyplot as pt - pt.figure(figsize=(10, 8)) + pt.figure(figsize=(10, 8), dpi=300) pt.plot(f1, label='Direct') pt.plot(f2, label='Chained') pt.ylim([np.min(f2) - 0.1, np.max(f2) + 0.1]) pt.legend() - pt.savefig('test_chained_to_direct.png', dpi=300) + pt.savefig('test_chained_to_direct.png') pt.clf() assert np.allclose(f1, f2) +@pytest.mark.parametrize("ndim", [2, 3]) +def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + # build test connection + np.random.seed(42) + discr = create_discretization(queue, ndim, + nelements=16, + 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) + from meshmode.discretization.connection import \ + ReversedChainedDiscretizationConnection + reverse = ReversedChainedDiscretizationConnection(chained) + from meshmode.discretization.connection import \ + flatten_chained_connection + direct = flatten_chained_connection(queue, chained) + + + # create test vector + from_nodes = chained.from_discr.nodes().with_queue(queue) + to_nodes = chained.to_discr.nodes().with_queue(queue) + + from_x = 0.0 + to_x = 0.0 + for i in range(ndim): + from_x += cl.clmath.cos(from_nodes[i]) ** (i + 1.0) + to_x += cl.clmath.cos(to_nodes[i]) ** (i + 1.0) + from_x.fill(1.0) + to_x.fill(1.0) + + from_interp = reverse(queue, to_x) + print(from_interp) + print(from_x) + + if visualize and ndim == 2: + import matplotlib.pyplot as pt + + from_group = chained.from_discr.groups[0] + to_group = chained.to_discr.groups[0] + + from_t = np.linspace(0.0, 1.0, chained.from_discr.nnodes) + to_t = chained(queue, cl.array.to_device(queue, from_t)).get(queue) + + pt.plot(from_t, from_x.get(queue), '--', label="From") + pt.plot(to_t, to_x.get(queue), '--', label="To") + pt.plot(from_t, from_interp.get(queue), label="Projection") + pt.legend() + pt.savefig("test_reverse_chained_conn.png") + pt.clf() + + if __name__ == "__main__": import sys if len(sys.argv) > 1: -- GitLab From 4d2d290735b42ef768c6e9b64e8bf2a998c2e83f Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 27 Nov 2018 20:15:48 -0600 Subject: [PATCH 02/19] connection: more work on reverse connection. Seems to work in 2D --- .../discretization/connection/__init__.py | 249 +++++++++++++++--- test/test_chained.py | 46 ++-- 2 files changed, 234 insertions(+), 61 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index fab24c9d..9e984677 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -251,67 +251,188 @@ class ChainedDiscretizationConnection(DiscretizationConnection): return vec -class ReversedChainedDiscretizationConnection(DiscretizationConnection): - def __init__(self, chained_connection): - self.conn = chained_connection - if len(self.conn.to_discr.groups) != 1 \ - or len(self.conn.from_discr.groups) != 1: +class ReversedDiscretizationConnection(DiscretizationConnection): + """Creates a reverse :class:`DiscretizationConnection` to allow + interpolating from *to_discr* to *from_discr*. + """ + + def __new__(cls, connections, is_surjective=False): + if isinstance(connections, DirectDiscretizationConnection): + return DiscretizationConnection.__new__(cls) + 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 len(conn.to_discr.groups) != 1 or len(conn.from_discr.groups) != 1: raise RuntimeError("multiple element groups are not supported") - super(ReversedChainedDiscretizationConnection, self).__init__( + # TODO: add a check for element -> face connections? + + self.conn = conn + super(ReversedDiscretizationConnection, self).__init__( from_discr=self.conn.to_discr, to_discr=self.conn.from_discr, - is_surjective=False) + is_surjective=is_surjective) @memoize_method - def _evaluate_basis(self, queue, group, ibasis): + def _evaluate_target_basis(self, queue, target, igroup, ibasis): + """Evaluates the a basis functions of the target :attr:`to_discr` + on another discretization. + + :arg target: one of :attr:`to_discr` or :attr:`from_discr`, on which + to evaluate the basis functions. + :arg igroup: index of a group in :attr:`to_discr`. + :arg ibasis: index of a basis in the group. + :returns: an array of shape ``(target.nnodes,)``. The basis function + *ibasis* is evaluated on each element at the group's + reference nodes. It is then transported to the target discretization, + if not already there. + """ + + @memoize_in(self, "conn_basis_eval_kernel") + def knl(): + import loopy as lp + knl = lp.make_kernel([ + "{[i]: 0 <= i < nelements}", + "{[j]: 0 <= j < n_to_nodes}" + ], + """ + vec[i, j] = basis[j] + """, + [ + lp.GlobalArg("vec", None, + shape=("nelements", "n_to_nodes")), + lp.GlobalArg("basis", None, + shape=("n_to_nodes")), + lp.ValueArg("nelements", np.int32), + lp.ValueArg("n_to_nodes", np.int32), + '...' + ], + name="conn_basis_eval_kernel", + lang_version=MOST_RECENT_LANGUAGE_VERSION) + + return knl + + group = self.to_discr.groups[igroup] basis_fn = group.basis()[ibasis] # evaluate - basis = np.zeros(self.to_discr.nnodes) - view = group.view(basis) - for k in range(view.shape[0]): - view[k, :] = basis_fn(group.unit_nodes) + basis = self.to_discr.zeros(queue) + knl()(queue, + vec=group.view(basis), + basis=basis_fn(group.unit_nodes).flatten()) - basis = cl.array.to_device(queue, basis) - basis = self.conn(queue, basis).with_queue(None) + if target is self.from_discr: + basis = self.conn(queue, basis) - return basis + return basis.with_queue(None) @memoize_method - def _build_element_index_map(self, queue, group): + def _evaluate_target_index(self, queue, target, igroup): + """ + :arg target: one of :attr:`to_discr` or `from_discr`. + :arg igroup: index of a group in :attr:`to_discr`. + :returns: a tuple ``(indices, count)``, where *indices* is an + array of size ``(target.nnodes)`` and *count* is an array + of size ``(group.nelements)``. The *indices* array is a mapping + from target node indices to :attr:`to_discr` element indices + in the given group. The *count* array is a mapping from the + element indices in the group to the number of elements in was + refined on *target*. + """ + + @memoize_in(self, "conn_index_kernel") + def kindex(): + import loopy as lp + knl = lp.make_kernel([ + "{[i]: 0 <= i < nelements}", + "{[j]: 0 <= j < n_to_nodes}" + ], + """ + indices[i, j] = i + """, + [ + lp.GlobalArg("indices", None, + shape=("nelements", "n_to_nodes")), + lp.ValueArg("n_to_nodes", np.int32), + lp.ValueArg("nelements", np.int32), + '...' + ], + name="conn_index_kernel", + lang_version=MOST_RECENT_LANGUAGE_VERSION) + + return knl + + @memoize_in(self, "conn_count_kernel") + def kcount(): + import loopy as lp + knl = lp.make_kernel([ + "{[i]: 0 <= i < nsources}", + ], + """ + count[indices[i, 0]] = count[indices[i, 0]] + 1 + """, + [ + lp.GlobalArg("indices", None, + shape=("nsources", "n_from_nodes")), + lp.GlobalArg("count", None, + shape=("ntargets")), + lp.ValueArg("n_from_nodes", np.int32), + lp.ValueArg("ntargets", np.int32), + '...' + ], + name="conn_count_kernel", + lang_version=MOST_RECENT_LANGUAGE_VERSION) + + return knl + + from_group = self.from_discr.groups[igroup] + to_group = self.to_discr.groups[igroup] + # evaluate element indices on source - indices = np.full(self.to_discr.nnodes, -1.0) - view = group.view(indices) - for k in range(view.shape[0]): - view[k, :] = k + indices = cl.array.empty(queue, self.to_discr.nnodes, dtype=np.float) + indices.fill(-1.0) + kindex()(queue, + indices=to_group.view(indices)) - # interpolate to target - indices = self.conn(queue, cl.array.to_device(queue, indices)) - indices = cl.clmath.round(indices.with_queue(queue)).astype(np.int) + if target is self.from_discr: + indices = self.conn(queue, indices).with_queue(queue) + indices = cl.clmath.round(indices).astype(np.int) - return indices.with_queue(None) + # count how many times each element was divided during refinement + count = cl.array.zeros(queue, to_group.nelements, dtype=np.int) + kcount()(queue, + indices=from_group.view(indices), + count=count) + + return indices.with_queue(None), count.with_queue(None) def __call__(self, queue, vec): @memoize_in(self, "conn_projection_kernel") - def knl(): + def kproj(): import loopy as lp knl = lp.make_kernel([ "{[k]: 0 <= k < nsources}", "{[j]: 0 <= j < n_from_nodes}" ], """ - for k, j - result[indices[k, j], ibasis] = \ - result[indices[k, j], ibasis] + \ - vec[k, j] * basis[k, j] * weights[j] - end + for k, j + result[indices[k, j], ibasis] = \ + result[indices[k, j], ibasis] + \ + vec[k, j] * basis[k, j] * weights[j] / count[indices[k, j]] + end """, [ lp.GlobalArg("vec", None, shape=("nsources", "n_from_nodes")), lp.GlobalArg("result", None, shape=("ntargets", "n_to_nodes")), + lp.GlobalArg("count", None, + shape=("ntargets")), lp.GlobalArg("indices", None, shape=("nsources", "n_from_nodes")), lp.GlobalArg("basis", None, @@ -328,25 +449,73 @@ class ReversedChainedDiscretizationConnection(DiscretizationConnection): return knl + @memoize_in(self, "conn_evaluate_kernel") + def keval(): + import loopy as lp + knl = lp.make_kernel([ + "{[k]: 0 <= k < nresults}", + "{[j]: 0 <= j < n_to_nodes}" + ], + """ + result[k, j] = result[k, j] + \ + coefficients[k, ibasis] * basis[k, j] + """, + [ + lp.GlobalArg("result", None, + shape=("nresults", "n_to_nodes")), + lp.GlobalArg("basis", None, + shape=("nresults", "n_to_nodes")), + lp.GlobalArg("coefficients", None, + shape=("nresults", "n_to_nodes")), + lp.ValueArg("ibasis", np.int32), + lp.ValueArg("n_to_nodes", np.int32), + lp.ValueArg("nresults", np.int32), + '...' + ], + name="conn_evaluate_kernel", + 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") - result = self.to_discr.zeros(queue, dtype=vec.dtype) - for tgrp in self.to_discr.groups: - el_indices = self._build_element_index_map(queue, tgrp) + # perform dot product to get basis coefficients + # NOTE: this may work for multiple groups (untested), but it will + # set all nodes not in the group to zero. + c = self.to_discr.zeros(queue, dtype=vec.dtype) + for igrp, tgrp in enumerate(self.to_discr.groups): + el_indices, el_count = self._evaluate_target_index(queue, + self.from_discr, igrp) + for ibasis in range(len(tgrp.basis())): - basis = self._evaluate_basis(queue, tgrp, ibasis) + basis = self._evaluate_target_basis(queue, + self.from_discr, igrp, ibasis) for sgrp in self.from_discr.groups: - knl()(queue, - weights=sgrp.weights, ibasis=ibasis, - vec=sgrp.view(vec), - basis=sgrp.view(basis), - indices=sgrp.view(el_indices), - result=tgrp.view(result)) + kproj()(queue, + vec=sgrp.view(vec), + basis=sgrp.view(basis), + indices=sgrp.view(el_indices), + count=el_count, + weights=sgrp.weights, ibasis=ibasis, + result=tgrp.view(c)) + + # 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 in range(len(grp.basis())): + basis = self._evaluate_target_basis(queue, + self.to_discr, igrp, ibasis) + assert basis.size == self.to_discr.nnodes + + keval()(queue, + result=grp.view(result), + basis=grp.view(basis), ibasis=ibasis, + coefficients=grp.view(c)) return result diff --git a/test/test_chained.py b/test/test_chained.py index f516ecb1..43cdc477 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -23,6 +23,7 @@ THE SOFTWARE. """ import numpy as np +import numpy.linalg as la import pyopencl as cl import pyopencl.array # noqa @@ -337,29 +338,33 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) + if ndim == 2: + order = 5 + threshold = 0.3 + else: + order = 8 + threshold = 0.1 + # build test connection np.random.seed(42) discr = create_discretization(queue, ndim, - nelements=16, - mesh_order=2, - discr_order=2) + mesh_order=order, + discr_order=order) connections = [] - conn = create_refined_connection(queue, discr) + conn = create_refined_connection(queue, discr, threshold=threshold) connections.append(conn) - conn = create_refined_connection(queue, conn.to_discr) + conn = create_refined_connection(queue, conn.to_discr, threshold=threshold) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr, threshold=threshold) connections.append(conn) from meshmode.discretization.connection import \ ChainedDiscretizationConnection chained = ChainedDiscretizationConnection(connections) from meshmode.discretization.connection import \ - ReversedChainedDiscretizationConnection - reverse = ReversedChainedDiscretizationConnection(chained) - from meshmode.discretization.connection import \ - flatten_chained_connection - direct = flatten_chained_connection(queue, chained) - + ReversedDiscretizationConnection + reverse = ReversedDiscretizationConnection(connections) # create test vector from_nodes = chained.from_discr.nodes().with_queue(queue) @@ -370,25 +375,24 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): for i in range(ndim): from_x += cl.clmath.cos(from_nodes[i]) ** (i + 1.0) to_x += cl.clmath.cos(to_nodes[i]) ** (i + 1.0) - from_x.fill(1.0) - to_x.fill(1.0) - from_interp = reverse(queue, to_x) - print(from_interp) - print(from_x) + + from_interp = from_interp.get(queue) + from_x = from_x.get(queue) + + print(to_x.shape, from_x.shape) + print(la.norm(from_interp - from_x)) + assert la.norm(from_interp - from_x) / la.norm(from_x) < 1.0e-12 if visualize and ndim == 2: import matplotlib.pyplot as pt - from_group = chained.from_discr.groups[0] - to_group = chained.to_discr.groups[0] - from_t = np.linspace(0.0, 1.0, chained.from_discr.nnodes) to_t = chained(queue, cl.array.to_device(queue, from_t)).get(queue) - pt.plot(from_t, from_x.get(queue), '--', label="From") + pt.plot(from_t, from_x, '--', label="From") pt.plot(to_t, to_x.get(queue), '--', label="To") - pt.plot(from_t, from_interp.get(queue), label="Projection") + pt.plot(from_t, from_interp, 'o-', label="Projection") pt.legend() pt.savefig("test_reverse_chained_conn.png") pt.clf() -- GitLab From 5704cef8569e5cf1e961ded2e7b9888f21e7982c Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Wed, 5 Dec 2018 10:49:25 -0600 Subject: [PATCH 03/19] connection: clean up reverse connection impl --- .../discretization/connection/__init__.py | 236 +++++++----------- test/test_chained.py | 38 +-- 2 files changed, 120 insertions(+), 154 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 9e984677..5cabbcdb 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -252,8 +252,17 @@ class ChainedDiscretizationConnection(DiscretizationConnection): class ReversedDiscretizationConnection(DiscretizationConnection): - """Creates a reverse :class:`DiscretizationConnection` to allow - interpolating from *to_discr* to *from_discr*. + """Creates a reverse :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): @@ -267,6 +276,7 @@ class ReversedDiscretizationConnection(DiscretizationConnection): return ChainedDiscretizationConnection(conns) def __init__(self, conn, is_surjective=False): + # TODO: this may not be strictly necessary if len(conn.to_discr.groups) != 1 or len(conn.from_discr.groups) != 1: raise RuntimeError("multiple element groups are not supported") @@ -279,20 +289,7 @@ class ReversedDiscretizationConnection(DiscretizationConnection): is_surjective=is_surjective) @memoize_method - def _evaluate_target_basis(self, queue, target, igroup, ibasis): - """Evaluates the a basis functions of the target :attr:`to_discr` - on another discretization. - - :arg target: one of :attr:`to_discr` or :attr:`from_discr`, on which - to evaluate the basis functions. - :arg igroup: index of a group in :attr:`to_discr`. - :arg ibasis: index of a basis in the group. - :returns: an array of shape ``(target.nnodes,)``. The basis function - *ibasis* is evaluated on each element at the group's - reference nodes. It is then transported to the target discretization, - if not already there. - """ - + def _evaluate_to_basis(self, queue, target, igroup, ibasis): @memoize_in(self, "conn_basis_eval_kernel") def knl(): import loopy as lp @@ -300,18 +297,8 @@ class ReversedDiscretizationConnection(DiscretizationConnection): "{[i]: 0 <= i < nelements}", "{[j]: 0 <= j < n_to_nodes}" ], - """ - vec[i, j] = basis[j] - """, - [ - lp.GlobalArg("vec", None, - shape=("nelements", "n_to_nodes")), - lp.GlobalArg("basis", None, - shape=("n_to_nodes")), - lp.ValueArg("nelements", np.int32), - lp.ValueArg("n_to_nodes", np.int32), - '...' - ], + """vec[i, j] = basis[j]""", + ['...'], name="conn_basis_eval_kernel", lang_version=MOST_RECENT_LANGUAGE_VERSION) @@ -332,128 +319,93 @@ class ReversedDiscretizationConnection(DiscretizationConnection): return basis.with_queue(None) @memoize_method - def _evaluate_target_index(self, queue, target, igroup): - """ - :arg target: one of :attr:`to_discr` or `from_discr`. - :arg igroup: index of a group in :attr:`to_discr`. - :returns: a tuple ``(indices, count)``, where *indices* is an - array of size ``(target.nnodes)`` and *count* is an array - of size ``(group.nelements)``. The *indices* array is a mapping - from target node indices to :attr:`to_discr` element indices - in the given group. The *count* array is a mapping from the - element indices in the group to the number of elements in was - refined on *target*. - """ - - @memoize_in(self, "conn_index_kernel") - def kindex(): - import loopy as lp - knl = lp.make_kernel([ - "{[i]: 0 <= i < nelements}", - "{[j]: 0 <= j < n_to_nodes}" - ], - """ - indices[i, j] = i - """, - [ - lp.GlobalArg("indices", None, - shape=("nelements", "n_to_nodes")), - lp.ValueArg("n_to_nodes", np.int32), - lp.ValueArg("nelements", np.int32), - '...' - ], - name="conn_index_kernel", - lang_version=MOST_RECENT_LANGUAGE_VERSION) - - return knl - - @memoize_in(self, "conn_count_kernel") - def kcount(): + def _from_quad_weights(self, queue): + @memoize_in(self, "conn_quad_weights_knl") + def kweights(): import loopy as lp knl = lp.make_kernel([ - "{[i]: 0 <= i < nsources}", + "{[k]: 0 <= k < nelements}", + "{[j]: 0 <= j < n_from_nodes}", ], """ - count[indices[i, 0]] = count[indices[i, 0]] + 1 + result[from_element_indices[k], j] = \ + weights[j] / scaling[to_element_indices[k]] """, [ - lp.GlobalArg("indices", None, - shape=("nsources", "n_from_nodes")), - lp.GlobalArg("count", None, - shape=("ntargets")), - lp.ValueArg("n_from_nodes", np.int32), - lp.ValueArg("ntargets", np.int32), - '...' + lp.GlobalArg("result", None, + shape=("n_from_elements", "n_from_nodes")), + lp.ValueArg("n_from_elements", np.int32), + "..." ], - name="conn_count_kernel", + name="conn_quad_weights_knl", lang_version=MOST_RECENT_LANGUAGE_VERSION) return knl - from_group = self.from_discr.groups[igroup] - to_group = self.to_discr.groups[igroup] - - # evaluate element indices on source - indices = cl.array.empty(queue, self.to_discr.nnodes, dtype=np.float) - indices.fill(-1.0) - kindex()(queue, - indices=to_group.view(indices)) - - if target is self.from_discr: - indices = self.conn(queue, indices).with_queue(queue) - indices = cl.clmath.round(indices).astype(np.int) + elranges = np.cumsum([0] + + [g.nelements for g in self.to_discr.groups]) + scaling = cl.array.zeros(queue, elranges[-1], dtype=np.int) + for igrp, grp in enumerate(self.conn.groups): + for batch in grp.batches: + scaling[elranges[igrp] + batch.from_element_indices] += 1 - # count how many times each element was divided during refinement - count = cl.array.zeros(queue, to_group.nelements, dtype=np.int) - kcount()(queue, - indices=from_group.view(indices), - count=count) + weights = self.from_discr.empty(queue, dtype=np.float) + for igrp, (tgrp, cgrp) in enumerate( + zip(self.from_discr.groups, self.conn.groups)): + for batch in cgrp.batches: + kweights()(queue, + weights=tgrp.weights, + result=tgrp.view(weights), + scaling=scaling[elranges[igrp]:elranges[igrp + 1]], + to_element_indices=batch.from_element_indices, + from_element_indices=batch.to_element_indices) - return indices.with_queue(None), count.with_queue(None) + return weights.with_queue(None) def __call__(self, queue, vec): - @memoize_in(self, "conn_projection_kernel") + @memoize_in(self, "conn_projection_knl") def kproj(): import loopy as lp knl = lp.make_kernel([ - "{[k]: 0 <= k < nsources}", + "{[k]: 0 <= k < nelements}", "{[j]: 0 <= j < n_from_nodes}" ], """ - for k, j - result[indices[k, j], ibasis] = \ - result[indices[k, j], ibasis] + \ - vec[k, j] * basis[k, j] * weights[j] / count[indices[k, j]] + for k + <> element_dot = \ + sum(j, vec[from_element_indices[k], j] * \ + basis[from_element_indices[k], j] * \ + weights[from_element_indices[k], j]) + + result[to_element_indices[k], ibasis] = \ + result[to_element_indices[k], ibasis] + element_dot end """, [ lp.GlobalArg("vec", None, - shape=("nsources", "n_from_nodes")), - lp.GlobalArg("result", None, - shape=("ntargets", "n_to_nodes")), - lp.GlobalArg("count", None, - shape=("ntargets")), - lp.GlobalArg("indices", None, - shape=("nsources", "n_from_nodes")), + shape=("n_from_elements", "n_from_nodes")), lp.GlobalArg("basis", None, - shape=("nsources", "n_from_nodes")), + shape=("n_from_elements", "n_from_nodes")), lp.GlobalArg("weights", None, - shape="n_from_nodes"), - lp.ValueArg("ibasis", np.int32), + shape=("n_from_elements", "n_from_nodes")), + lp.GlobalArg("result", None, + shape=("n_to_elements", "n_to_nodes")), + lp.ValueArg("n_from_elements", np.int32), + lp.ValueArg("n_to_elements", np.int32), lp.ValueArg("n_to_nodes", np.int32), - lp.ValueArg("ntargets", np.int32), + lp.ValueArg("ibasis", np.int32), '...' ], - name="conn_projection_kernel", + name="conn_projection_knl", lang_version=MOST_RECENT_LANGUAGE_VERSION) return knl - @memoize_in(self, "conn_evaluate_kernel") + @memoize_in(self, "conn_evaluate_knl") def keval(): import loopy as lp knl = lp.make_kernel([ - "{[k]: 0 <= k < nresults}", + "{[k]: 0 <= k < nelements}", "{[j]: 0 <= j < n_to_nodes}" ], """ @@ -461,18 +413,13 @@ class ReversedDiscretizationConnection(DiscretizationConnection): coefficients[k, ibasis] * basis[k, j] """, [ - lp.GlobalArg("result", None, - shape=("nresults", "n_to_nodes")), - lp.GlobalArg("basis", None, - shape=("nresults", "n_to_nodes")), lp.GlobalArg("coefficients", None, - shape=("nresults", "n_to_nodes")), + shape=("nelements", "n_to_nodes")), lp.ValueArg("ibasis", np.int32), - lp.ValueArg("n_to_nodes", np.int32), - lp.ValueArg("nresults", np.int32), '...' ], - name="conn_evaluate_kernel", + name="conn_evaluate_knl", + default_offset=lp.auto, lang_version=MOST_RECENT_LANGUAGE_VERSION) return knl @@ -483,38 +430,49 @@ class ReversedDiscretizationConnection(DiscretizationConnection): if vec.shape != (self.from_discr.nnodes,): raise ValueError("invalid shape of incoming resampling data") + # the code proceeds as follows: + # * we evaluate the basis of to_discr on from_discr nodes + # * in each element of from_discr, we compute the dot product + # with the basis components to get the basis coefficients on + # to_discr + # * once the cofficients are computed and scaled, we evaluate them + # at the actual nodes of to_discr. + # perform dot product to get basis coefficients - # NOTE: this may work for multiple groups (untested), but it will - # set all nodes not in the group to zero. + weights = self._from_quad_weights(queue) c = self.to_discr.zeros(queue, dtype=vec.dtype) - for igrp, tgrp in enumerate(self.to_discr.groups): - el_indices, el_count = self._evaluate_target_index(queue, - self.from_discr, igrp) - + for igrp, (tgrp, cgrp) in enumerate( + zip(self.to_discr.groups, self.conn.groups)): for ibasis in range(len(tgrp.basis())): - basis = self._evaluate_target_basis(queue, + to_basis = self._evaluate_to_basis(queue, self.from_discr, igrp, ibasis) - for sgrp in self.from_discr.groups: - kproj()(queue, + for ibatch, batch in enumerate(cgrp.batches): + sgrp = self.from_discr.groups[batch.from_group_index] + + # 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=sgrp.view(basis), - indices=sgrp.view(el_indices), - count=el_count, - weights=sgrp.weights, ibasis=ibasis, - result=tgrp.view(c)) + basis=sgrp.view(to_basis), + weights=sgrp.view(weights), + 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 in range(len(grp.basis())): - basis = self._evaluate_target_basis(queue, + basis = self._evaluate_to_basis(queue, self.to_discr, igrp, ibasis) - assert basis.size == self.to_discr.nnodes - keval()(queue, + keval()(queue, ibasis=ibasis, result=grp.view(result), - basis=grp.view(basis), ibasis=ibasis, + basis=grp.view(basis), coefficients=grp.view(c)) return result diff --git a/test/test_chained.py b/test/test_chained.py index 43cdc477..535d27f2 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -53,7 +53,7 @@ def create_discretization(queue, ndim, order=mesh_order) elif ndim == 3: from meshmode.mesh.generation import generate_torus - mesh = generate_torus(10.0, 5.0, order=mesh_order) + mesh = generate_torus(2.0, 1.0, order=mesh_order) else: raise ValueError("Unsupported dimension: {}".format(ndim)) @@ -339,14 +339,13 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): queue = cl.CommandQueue(ctx) if ndim == 2: - order = 5 + order = 6 threshold = 0.3 else: - order = 8 + order = 12 threshold = 0.1 # build test connection - np.random.seed(42) discr = create_discretization(queue, ndim, mesh_order=order, discr_order=order) @@ -373,29 +372,38 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): from_x = 0.0 to_x = 0.0 for i in range(ndim): - from_x += cl.clmath.cos(from_nodes[i]) ** (i + 1.0) - to_x += cl.clmath.cos(to_nodes[i]) ** (i + 1.0) + from_x += cl.clmath.cos(from_nodes[i] / (4.0 * np.pi)) ** (i + 1.0) + to_x += cl.clmath.cos(to_nodes[i] / (4.0 * np.pi)) ** (i + 1.0) from_interp = reverse(queue, to_x) - from_interp = from_interp.get(queue) - from_x = from_x.get(queue) - - print(to_x.shape, from_x.shape) - print(la.norm(from_interp - from_x)) - assert la.norm(from_interp - from_x) / la.norm(from_x) < 1.0e-12 - if visualize and ndim == 2: import matplotlib.pyplot as pt from_t = np.linspace(0.0, 1.0, chained.from_discr.nnodes) to_t = chained(queue, cl.array.to_device(queue, from_t)).get(queue) - pt.plot(from_t, from_x, '--', label="From") + pt.figure(figsize=(10, 8), dpi=300) + pt.plot(from_t, from_x.get(queue), '--', label="From") pt.plot(to_t, to_x.get(queue), '--', label="To") - pt.plot(from_t, from_interp, 'o-', label="Projection") + pt.plot(from_t, from_interp.get(queue), 'o-', label="Projection") pt.legend() pt.savefig("test_reverse_chained_conn.png") pt.clf() + elif visualize and ndim == 3: + from meshmode.discretization.visualization import make_visualizer + vis = make_visualizer(queue, chained.from_discr, order) + + vis.write_vtk_file("test_reverse_chained_conn.vtu", [ + ("x_interp", from_interp), + ("x", from_x) + ], overwrite=True) + + from_interp = from_interp.get(queue) + from_x = from_x.get(queue) + + print(to_x.shape, from_x.shape) + print(la.norm(from_interp - from_x) / la.norm(from_x)) + assert la.norm(from_interp - from_x) / la.norm(from_x) < 1.0e-9 if __name__ == "__main__": -- GitLab From a7381c8c874446f795431c78a35a704d2b85ab4c Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 11 Dec 2018 16:22:32 -0600 Subject: [PATCH 04/19] connection: add reverse class to docs --- meshmode/discretization/connection/__init__.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 5cabbcdb..0d7ab3cc 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 @@ -63,6 +66,7 @@ __all__ = [ __doc__ = """ .. autoclass:: DiscretizationConnection .. autoclass:: ChainedDiscretizationConnection +.. autoclass:: ReversedDiscretizationConnection .. autoclass:: DirectDiscretizationConnection .. autofunction:: make_same_mesh_connection -- GitLab From aa32233d162fc4b83b5b1d18bd49d1ae93e54e4a Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 11 Dec 2018 16:29:53 -0600 Subject: [PATCH 05/19] connection: allow taking a chained connection --- meshmode/discretization/connection/__init__.py | 2 ++ test/test_chained.py | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 0d7ab3cc..150c999c 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -272,6 +272,8 @@ class ReversedDiscretizationConnection(DiscretizationConnection): 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): diff --git a/test/test_chained.py b/test/test_chained.py index 535d27f2..c2b7ee55 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -363,7 +363,7 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): chained = ChainedDiscretizationConnection(connections) from meshmode.discretization.connection import \ ReversedDiscretizationConnection - reverse = ReversedDiscretizationConnection(connections) + reverse = ReversedDiscretizationConnection(chained) # create test vector from_nodes = chained.from_discr.nodes().with_queue(queue) -- GitLab From c0c126190f1b88ad45f868c88ccf49cd946744f7 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 11 Dec 2018 16:34:59 -0600 Subject: [PATCH 06/19] connection: do not memoize with queue --- .../discretization/connection/__init__.py | 69 ++++++++++--------- 1 file changed, 36 insertions(+), 33 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 150c999c..f88a38f5 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -295,7 +295,7 @@ class ReversedDiscretizationConnection(DiscretizationConnection): is_surjective=is_surjective) @memoize_method - def _evaluate_to_basis(self, queue, target, igroup, ibasis): + def _evaluate_to_basis(self, target, igroup, ibasis): @memoize_in(self, "conn_basis_eval_kernel") def knl(): import loopy as lp @@ -313,19 +313,19 @@ class ReversedDiscretizationConnection(DiscretizationConnection): group = self.to_discr.groups[igroup] basis_fn = group.basis()[ibasis] - # evaluate - basis = self.to_discr.zeros(queue) - knl()(queue, - vec=group.view(basis), - basis=basis_fn(group.unit_nodes).flatten()) + with cl.CommandQueue(target.cl_context) as queue: + basis = self.to_discr.zeros(queue) + knl()(queue, + vec=group.view(basis), + basis=basis_fn(group.unit_nodes).flatten()) - if target is self.from_discr: - basis = self.conn(queue, basis) + if target is self.from_discr: + basis = self.conn(queue, basis) return basis.with_queue(None) @memoize_method - def _from_quad_weights(self, queue): + def _from_quad_weights(self): @memoize_in(self, "conn_quad_weights_knl") def kweights(): import loopy as lp @@ -348,23 +348,25 @@ class ReversedDiscretizationConnection(DiscretizationConnection): return knl - elranges = np.cumsum([0] - + [g.nelements for g in self.to_discr.groups]) - scaling = cl.array.zeros(queue, elranges[-1], dtype=np.int) - for igrp, grp in enumerate(self.conn.groups): - for batch in grp.batches: - scaling[elranges[igrp] + batch.from_element_indices] += 1 - - weights = self.from_discr.empty(queue, dtype=np.float) - for igrp, (tgrp, cgrp) in enumerate( - zip(self.from_discr.groups, self.conn.groups)): - for batch in cgrp.batches: - kweights()(queue, - weights=tgrp.weights, - result=tgrp.view(weights), - scaling=scaling[elranges[igrp]:elranges[igrp + 1]], - to_element_indices=batch.from_element_indices, - from_element_indices=batch.to_element_indices) + with cl.CommandQueue(self.to_discr.cl_context) as queue: + elranges = np.cumsum([0] + + [g.nelements for g in self.to_discr.groups]) + scaling = cl.array.zeros(queue, elranges[-1], dtype=np.int) + for igrp, grp in enumerate(self.conn.groups): + for batch in grp.batches: + indices = batch.from_element_indices.with_queue(queue) + scaling[elranges[igrp] + indices] += 1 + + weights = self.from_discr.empty(queue, dtype=np.float) + for igrp, (tgrp, cgrp) in enumerate( + zip(self.from_discr.groups, self.conn.groups)): + for batch in cgrp.batches: + kweights()(queue, + weights=tgrp.weights, + result=tgrp.view(weights), + scaling=scaling[elranges[igrp]:elranges[igrp + 1]], + to_element_indices=batch.from_element_indices, + from_element_indices=batch.to_element_indices) return weights.with_queue(None) @@ -445,13 +447,13 @@ class ReversedDiscretizationConnection(DiscretizationConnection): # at the actual nodes of to_discr. # perform dot product to get basis coefficients - weights = self._from_quad_weights(queue) + weights = self._from_quad_weights() c = self.to_discr.zeros(queue, dtype=vec.dtype) for igrp, (tgrp, cgrp) in enumerate( zip(self.to_discr.groups, self.conn.groups)): for ibasis in range(len(tgrp.basis())): - to_basis = self._evaluate_to_basis(queue, - self.from_discr, igrp, ibasis) + to_basis = \ + self._evaluate_to_basis(self.from_discr, igrp, ibasis) for ibatch, batch in enumerate(cgrp.batches): sgrp = self.from_discr.groups[batch.from_group_index] @@ -460,7 +462,8 @@ class ReversedDiscretizationConnection(DiscretizationConnection): # 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, + kproj()(queue, + ibasis=ibasis, vec=sgrp.view(vec), basis=sgrp.view(to_basis), weights=sgrp.view(weights), @@ -473,10 +476,10 @@ class ReversedDiscretizationConnection(DiscretizationConnection): for igrp, grp in enumerate(self.to_discr.groups): for ibasis in range(len(grp.basis())): - basis = self._evaluate_to_basis(queue, - self.to_discr, igrp, ibasis) + basis = self._evaluate_to_basis(self.to_discr, igrp, ibasis) - keval()(queue, ibasis=ibasis, + keval()(queue, + ibasis=ibasis, result=grp.view(result), basis=grp.view(basis), coefficients=grp.view(c)) -- GitLab From 88bcb4457edae63376cfcf36a7ae04aba46eba85 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Wed, 12 Dec 2018 09:50:23 -0600 Subject: [PATCH 07/19] connection: add more checks in for reverse conn --- meshmode/discretization/connection/__init__.py | 15 ++++++++++----- test/test_chained.py | 18 ++++++++++++++++-- 2 files changed, 26 insertions(+), 7 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index f88a38f5..017af349 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -66,7 +66,7 @@ __all__ = [ __doc__ = """ .. autoclass:: DiscretizationConnection .. autoclass:: ChainedDiscretizationConnection -.. autoclass:: ReversedDiscretizationConnection +.. autoclass:: L2ProjectionInverseDiscretizationConnection .. autoclass:: DirectDiscretizationConnection .. autofunction:: make_same_mesh_connection @@ -256,7 +256,7 @@ class ChainedDiscretizationConnection(DiscretizationConnection): class ReversedDiscretizationConnection(DiscretizationConnection): - """Creates a reverse :class:`DiscretizationConnection` from an existing + """Creates an inverse :class:`DiscretizationConnection` from an existing connection to allow transporting from the original connection's *to_discr* to *from_discr*. @@ -282,11 +282,15 @@ class ReversedDiscretizationConnection(DiscretizationConnection): return ChainedDiscretizationConnection(conns) def __init__(self, conn, is_surjective=False): - # TODO: this may not be strictly necessary if len(conn.to_discr.groups) != 1 or len(conn.from_discr.groups) != 1: - raise RuntimeError("multiple element groups are not supported") + raise NotImplementedError("multiple element groups are " + "implemented in principle, but not yet tested") - # TODO: add a check for element -> face connections? + 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 and orthogonal basis") self.conn = conn super(ReversedDiscretizationConnection, self).__init__( @@ -352,6 +356,7 @@ class ReversedDiscretizationConnection(DiscretizationConnection): elranges = np.cumsum([0] + [g.nelements for g in self.to_discr.groups]) scaling = cl.array.zeros(queue, elranges[-1], dtype=np.int) + for igrp, grp in enumerate(self.conn.groups): for batch in grp.batches: indices = batch.from_element_indices.with_queue(queue) diff --git a/test/test_chained.py b/test/test_chained.py index c2b7ee55..37ed076d 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -374,7 +374,19 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): for i in range(ndim): from_x += cl.clmath.cos(from_nodes[i] / (4.0 * np.pi)) ** (i + 1.0) to_x += cl.clmath.cos(to_nodes[i] / (4.0 * np.pi)) ** (i + 1.0) + + import time + t_start = time.time() + from_interp = reverse(queue, to_x) + t_end = time.time() + if visualize: + print('Elapsed: {}'.format(t_end - t_start)) + + t_start = time.time() from_interp = reverse(queue, to_x) + t_end = time.time() + if visualize: + print('Elapsed: {}'.format(t_end - t_start)) if visualize and ndim == 2: import matplotlib.pyplot as pt @@ -401,8 +413,10 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): from_interp = from_interp.get(queue) from_x = from_x.get(queue) - print(to_x.shape, from_x.shape) - print(la.norm(from_interp - from_x) / la.norm(from_x)) + if visualize: + print(to_x.shape, from_x.shape) + print(la.norm(from_interp - from_x) / la.norm(from_x)) + assert la.norm(from_interp - from_x) / la.norm(from_x) < 1.0e-9 -- GitLab From 0f35f7e822cca7e98f4f5871cb76ed5b909dbfba Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sat, 13 Apr 2019 19:46:01 -0500 Subject: [PATCH 08/19] l2conn: finish rename --- meshmode/discretization/connection/__init__.py | 4 ++-- test/test_chained.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index c96aaa2c..ac378e6f 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -269,7 +269,7 @@ class ChainedDiscretizationConnection(DiscretizationConnection): return vec -class ReversedDiscretizationConnection(DiscretizationConnection): +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*. @@ -307,7 +307,7 @@ class ReversedDiscretizationConnection(DiscretizationConnection): raise RuntimeError("`to_discr` must have and orthogonal basis") self.conn = conn - super(ReversedDiscretizationConnection, self).__init__( + super(L2ProjectionInverseDiscretizationConnection, self).__init__( from_discr=self.conn.to_discr, to_discr=self.conn.from_discr, is_surjective=is_surjective) diff --git a/test/test_chained.py b/test/test_chained.py index 37ed076d..523aa48b 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -362,8 +362,8 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): ChainedDiscretizationConnection chained = ChainedDiscretizationConnection(connections) from meshmode.discretization.connection import \ - ReversedDiscretizationConnection - reverse = ReversedDiscretizationConnection(chained) + L2ProjectionInverseDiscretizationConnection + reverse = L2ProjectionInverseDiscretizationConnection(chained) # create test vector from_nodes = chained.from_discr.nodes().with_queue(queue) -- GitLab From 7e00b9f41dbb17023811f89f1723f7c3f442b03d Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 16 Apr 2019 16:35:51 -0500 Subject: [PATCH 09/19] l2conn: rewrite dot product --- .../discretization/connection/__init__.py | 118 +++++++----------- test/test_chained.py | 17 +-- 2 files changed, 51 insertions(+), 84 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index ac378e6f..cdf7643d 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -312,9 +312,8 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): to_discr=self.conn.from_discr, is_surjective=is_surjective) - @memoize_method - def _evaluate_to_basis(self, target, igroup, ibasis): - @memoize_in(self, "conn_basis_eval_kernel") + def _evaluate_basis_on_target(self, discr, igroup, ibasis): + @memoize_in(self, "conn_basis_evaluation_knl") def knl(): import loopy as lp knl = lp.make_kernel([ @@ -329,65 +328,49 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): return knl group = self.to_discr.groups[igroup] - basis_fn = group.basis()[ibasis] + basis = group.basis()[ibasis] - with cl.CommandQueue(target.cl_context) as queue: - basis = self.to_discr.zeros(queue) + with cl.CommandQueue(discr.cl_context) as queue: + vec = self.to_discr.zeros(queue) knl()(queue, - vec=group.view(basis), - basis=basis_fn(group.unit_nodes).flatten()) + vec=group.view(vec), + basis=basis(group.unit_nodes).flatten()) - if target is self.from_discr: - basis = self.conn(queue, basis) + if discr is self.from_discr: + vec = self.conn(queue, vec) - return basis.with_queue(None) + return vec.with_queue(None) @memoize_method - def _from_quad_weights(self): - @memoize_in(self, "conn_quad_weights_knl") - def kweights(): - import loopy as lp - knl = lp.make_kernel([ - "{[k]: 0 <= k < nelements}", - "{[j]: 0 <= j < n_from_nodes}", - ], - """ - result[from_element_indices[k], j] = \ - weights[j] / scaling[to_element_indices[k]] - """, - [ - lp.GlobalArg("result", None, - shape=("n_from_elements", "n_from_nodes")), - lp.ValueArg("n_from_elements", np.int32), - "..." - ], - name="conn_quad_weights_knl", - lang_version=MOST_RECENT_LANGUAGE_VERSION) + def _batch_weights(self): + from pymbolic.geometric_algebra import MultiVector + from functools import reduce + from operator import xor - return knl + 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()) - with cl.CommandQueue(self.to_discr.cl_context) as queue: - elranges = np.cumsum([0] - + [g.nelements for g in self.to_discr.groups]) - scaling = cl.array.zeros(queue, elranges[-1], dtype=np.int) - - for igrp, grp in enumerate(self.conn.groups): - for batch in grp.batches: - indices = batch.from_element_indices.with_queue(queue) - scaling[elranges[igrp] + indices] += 1 - - weights = self.from_discr.empty(queue, dtype=np.float) - for igrp, (tgrp, cgrp) in enumerate( - zip(self.from_discr.groups, self.conn.groups)): - for batch in cgrp.batches: - kweights()(queue, - weights=tgrp.weights, - result=tgrp.view(weights), - scaling=scaling[elranges[igrp]:elranges[igrp + 1]], - to_element_indices=batch.from_element_indices, - from_element_indices=batch.to_element_indices) - - return weights.with_queue(None) + return det_v + + to_discr = self.to_discr + J = np.empty(to_discr.dim, dtype=np.object) + weights = [np.empty(len(g.batches), dtype=np.object) + for g in self.conn.groups] + + for igrp, grp in enumerate(to_discr.groups): + for ibatch, batch in enumerate(self.conn.groups[igrp].batches): + for iaxis in range(grp.dim): + D = grp.diff_matrices()[iaxis] + J[iaxis] = D.dot(batch.result_unit_nodes.T) + + weights[igrp][ibatch] = det(J) * grp.weights + + return weights def __call__(self, queue, vec): @memoize_in(self, "conn_projection_knl") @@ -402,7 +385,7 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): <> element_dot = \ sum(j, vec[from_element_indices[k], j] * \ basis[from_element_indices[k], j] * \ - weights[from_element_indices[k], j]) + weights[j]) result[to_element_indices[k], ibasis] = \ result[to_element_indices[k], ibasis] + element_dot @@ -414,7 +397,7 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): lp.GlobalArg("basis", None, shape=("n_from_elements", "n_from_nodes")), lp.GlobalArg("weights", None, - shape=("n_from_elements", "n_from_nodes")), + shape="n_from_nodes"), lp.GlobalArg("result", None, shape=("n_to_elements", "n_to_nodes")), lp.ValueArg("n_from_elements", np.int32), @@ -428,7 +411,7 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): return knl - @memoize_in(self, "conn_evaluate_knl") + @memoize_in(self, "conn_evaluation_knl") def keval(): import loopy as lp knl = lp.make_kernel([ @@ -457,22 +440,16 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): if vec.shape != (self.from_discr.nnodes,): raise ValueError("invalid shape of incoming resampling data") - # the code proceeds as follows: - # * we evaluate the basis of to_discr on from_discr nodes - # * in each element of from_discr, we compute the dot product - # with the basis components to get the basis coefficients on - # to_discr - # * once the cofficients are computed and scaled, we evaluate them - # at the actual nodes of to_discr. + # compute weights on each refinement of the reference element + weights = self._batch_weights(); - # perform dot product to get basis coefficients - weights = self._from_quad_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 ibasis in range(len(tgrp.basis())): - to_basis = \ - self._evaluate_to_basis(self.from_discr, igrp, ibasis) + to_basis = self._evaluate_basis_on_target( + self.from_discr, igrp, ibasis) for ibatch, batch in enumerate(cgrp.batches): sgrp = self.from_discr.groups[batch.from_group_index] @@ -485,7 +462,7 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): ibasis=ibasis, vec=sgrp.view(vec), basis=sgrp.view(to_basis), - weights=sgrp.view(weights), + weights=weights[igrp][ibatch], result=tgrp.view(c), from_element_indices=batch.to_element_indices, to_element_indices=batch.from_element_indices) @@ -495,7 +472,8 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): for igrp, grp in enumerate(self.to_discr.groups): for ibasis in range(len(grp.basis())): - basis = self._evaluate_to_basis(self.to_discr, igrp, ibasis) + basis = self._evaluate_basis_on_target( + self.to_discr, igrp, ibasis) keval()(queue, ibasis=ibasis, diff --git a/test/test_chained.py b/test/test_chained.py index 523aa48b..8ff70a00 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -375,18 +375,7 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): from_x += cl.clmath.cos(from_nodes[i] / (4.0 * np.pi)) ** (i + 1.0) to_x += cl.clmath.cos(to_nodes[i] / (4.0 * np.pi)) ** (i + 1.0) - import time - t_start = time.time() from_interp = reverse(queue, to_x) - t_end = time.time() - if visualize: - print('Elapsed: {}'.format(t_end - t_start)) - - t_start = time.time() - from_interp = reverse(queue, to_x) - t_end = time.time() - if visualize: - print('Elapsed: {}'.format(t_end - t_start)) if visualize and ndim == 2: import matplotlib.pyplot as pt @@ -413,11 +402,11 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): from_interp = from_interp.get(queue) from_x = from_x.get(queue) + error = la.norm(from_interp - from_x) / la.norm(from_x) if visualize: - print(to_x.shape, from_x.shape) - print(la.norm(from_interp - from_x) / la.norm(from_x)) + print('Error: {}'.format(error)) - assert la.norm(from_interp - from_x) / la.norm(from_x) < 1.0e-9 + assert error < 1.0e-9 if __name__ == "__main__": -- GitLab From cdd0c7fa95b6b3aa796ef6af716980d3287515aa Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 16 Apr 2019 16:45:12 -0500 Subject: [PATCH 10/19] l2conn: appease flake8 --- meshmode/discretization/connection/__init__.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index cdf7643d..bfdf264c 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -358,17 +358,17 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): return det_v to_discr = self.to_discr - J = np.empty(to_discr.dim, dtype=np.object) + jac = np.empty(to_discr.dim, dtype=np.object) weights = [np.empty(len(g.batches), dtype=np.object) for g in self.conn.groups] for igrp, grp in enumerate(to_discr.groups): for ibatch, batch in enumerate(self.conn.groups[igrp].batches): for iaxis in range(grp.dim): - D = grp.diff_matrices()[iaxis] - J[iaxis] = D.dot(batch.result_unit_nodes.T) + mat = grp.diff_matrices()[iaxis] + jac[iaxis] = mat.dot(batch.result_unit_nodes.T) - weights[igrp][ibatch] = det(J) * grp.weights + weights[igrp][ibatch] = det(jac) * grp.weights return weights @@ -441,7 +441,7 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): raise ValueError("invalid shape of incoming resampling data") # compute weights on each refinement of the reference element - weights = self._batch_weights(); + weights = self._batch_weights() # perform dot product (on reference element) to get basis coefficients c = self.to_discr.zeros(queue, dtype=vec.dtype) -- GitLab From 2ee232b1ac63f2a9840efcf818b6663090e38aa6 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Tue, 16 Apr 2019 16:59:29 -0500 Subject: [PATCH 11/19] l2conn: use queue from call instead of creating many others --- .../discretization/connection/__init__.py | 21 +++++++++---------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index bfdf264c..b727c2e9 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -312,7 +312,7 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): to_discr=self.conn.from_discr, is_surjective=is_surjective) - def _evaluate_basis_on_target(self, discr, igroup, ibasis): + def _evaluate_basis_on_target(self, queue, discr, igroup, ibasis): @memoize_in(self, "conn_basis_evaluation_knl") def knl(): import loopy as lp @@ -330,16 +330,15 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): group = self.to_discr.groups[igroup] basis = group.basis()[ibasis] - with cl.CommandQueue(discr.cl_context) as queue: - vec = self.to_discr.zeros(queue) - knl()(queue, - vec=group.view(vec), - basis=basis(group.unit_nodes).flatten()) + vec = self.to_discr.zeros(queue) + knl()(queue, + vec=group.view(vec), + basis=basis(group.unit_nodes).flatten()) - if discr is self.from_discr: - vec = self.conn(queue, vec) + if discr is self.from_discr: + vec = self.conn(queue, vec) - return vec.with_queue(None) + return vec @memoize_method def _batch_weights(self): @@ -449,7 +448,7 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): zip(self.to_discr.groups, self.conn.groups)): for ibasis in range(len(tgrp.basis())): to_basis = self._evaluate_basis_on_target( - self.from_discr, igrp, ibasis) + queue, self.from_discr, igrp, ibasis) for ibatch, batch in enumerate(cgrp.batches): sgrp = self.from_discr.groups[batch.from_group_index] @@ -473,7 +472,7 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): for igrp, grp in enumerate(self.to_discr.groups): for ibasis in range(len(grp.basis())): basis = self._evaluate_basis_on_target( - self.to_discr, igrp, ibasis) + queue, self.to_discr, igrp, ibasis) keval()(queue, ibasis=ibasis, -- GitLab From f0f20818b823f18d9778935110edfafeede0180a Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 23 Apr 2019 12:07:39 -0500 Subject: [PATCH 12/19] l2conn: exception to warning --- meshmode/discretization/connection/__init__.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index b727c2e9..bd67f5e0 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -297,14 +297,15 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): def __init__(self, conn, is_surjective=False): if len(conn.to_discr.groups) != 1 or len(conn.from_discr.groups) != 1: - raise NotImplementedError("multiple element groups are " - "implemented in principle, but not yet tested") + from warnings import warn + warn("multiple element groups are implemented in principle, " + "but not yet tested") 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 and orthogonal basis") + raise RuntimeError("`to_discr` must have an orthogonal basis") self.conn = conn super(L2ProjectionInverseDiscretizationConnection, self).__init__( -- GitLab From d270d5c023db1f05cab98ab92237029fd2436af3 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 25 Apr 2019 21:53:06 -0500 Subject: [PATCH 13/19] l2conn: add some docs and cleanup --- .../discretization/connection/__init__.py | 46 +++++++++++++------ test/test_chained.py | 6 +-- 2 files changed, 35 insertions(+), 17 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index bd67f5e0..314b312d 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -313,7 +313,17 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): to_discr=self.conn.from_discr, is_surjective=is_surjective) - def _evaluate_basis_on_target(self, queue, discr, igroup, ibasis): + def _evaluate_basis_on_target(self, queue, discr, igroup, ibasis_func): + """Evaluates a basis function from the target discretization + :attr:`to_discr` on all elements of `discr`. + + :arg discr: discretization on which to evaluate the basis function. + :arg igroup: group index in the array returned by + :attr:`to_discr.groups`. + :arg ibasis_func: basis index in the array returned by + :attr:`to_discr.basis()`. + """ + @memoize_in(self, "conn_basis_evaluation_knl") def knl(): import loopy as lp @@ -322,14 +332,13 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): "{[j]: 0 <= j < n_to_nodes}" ], """vec[i, j] = basis[j]""", - ['...'], name="conn_basis_eval_kernel", lang_version=MOST_RECENT_LANGUAGE_VERSION) return knl group = self.to_discr.groups[igroup] - basis = group.basis()[ibasis] + basis = group.basis()[ibasis_func] vec = self.to_discr.zeros(queue) knl()(queue, @@ -343,6 +352,15 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): @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 + refined elements in the domain of the parent element. + + :returns: a list of lists of size `(ngroups, ngroup_batches)` + containing the quadrature weights (of size `batch.nunit_nodes`) + for a given batch in a given group of :attr:`conn`. + """ + from pymbolic.geometric_algebra import MultiVector from functools import reduce from operator import xor @@ -387,8 +405,8 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): basis[from_element_indices[k], j] * \ weights[j]) - result[to_element_indices[k], ibasis] = \ - result[to_element_indices[k], ibasis] + element_dot + result[to_element_indices[k], ibasis_func] = \ + result[to_element_indices[k], ibasis_func] + element_dot end """, [ @@ -403,7 +421,7 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): 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), + lp.ValueArg("ibasis_func", np.int32), '...' ], name="conn_projection_knl", @@ -420,12 +438,12 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): ], """ result[k, j] = result[k, j] + \ - coefficients[k, ibasis] * basis[k, j] + coefficients[k, ibasis_func] * basis[k, j] """, [ lp.GlobalArg("coefficients", None, shape=("nelements", "n_to_nodes")), - lp.ValueArg("ibasis", np.int32), + lp.ValueArg("ibasis_func", np.int32), '...' ], name="conn_evaluate_knl", @@ -447,9 +465,9 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): c = self.to_discr.zeros(queue, dtype=vec.dtype) for igrp, (tgrp, cgrp) in enumerate( zip(self.to_discr.groups, self.conn.groups)): - for ibasis in range(len(tgrp.basis())): + for ibasis_func in range(len(tgrp.basis())): to_basis = self._evaluate_basis_on_target( - queue, self.from_discr, igrp, ibasis) + queue, self.from_discr, igrp, ibasis_func) for ibatch, batch in enumerate(cgrp.batches): sgrp = self.from_discr.groups[batch.from_group_index] @@ -459,7 +477,7 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): # we are going in reverse here. a bit confusing, but # saves on recreating the connection groups and batches. kproj()(queue, - ibasis=ibasis, + ibasis_func=ibasis_func, vec=sgrp.view(vec), basis=sgrp.view(to_basis), weights=weights[igrp][ibatch], @@ -471,12 +489,12 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): result = self.to_discr.zeros(queue, dtype=vec.dtype) for igrp, grp in enumerate(self.to_discr.groups): - for ibasis in range(len(grp.basis())): + for ibasis_func in range(len(grp.basis())): basis = self._evaluate_basis_on_target( - queue, self.to_discr, igrp, ibasis) + queue, self.to_discr, igrp, ibasis_func) keval()(queue, - ibasis=ibasis, + ibasis_func=ibasis_func, result=grp.view(result), basis=grp.view(basis), coefficients=grp.view(c)) diff --git a/test/test_chained.py b/test/test_chained.py index 8ff70a00..fb186aed 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -53,7 +53,7 @@ def create_discretization(queue, ndim, order=mesh_order) elif ndim == 3: from meshmode.mesh.generation import generate_torus - mesh = generate_torus(2.0, 1.0, order=mesh_order) + mesh = generate_torus(10.0, 5.0, order=mesh_order) else: raise ValueError("Unsupported dimension: {}".format(ndim)) @@ -342,7 +342,7 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): order = 6 threshold = 0.3 else: - order = 12 + order = 6 threshold = 0.1 # build test connection @@ -406,7 +406,7 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): if visualize: print('Error: {}'.format(error)) - assert error < 1.0e-9 + assert error < 1.0e-9, error if __name__ == "__main__": -- GitLab From 46ba318760f33911976c281d23a41c4abbdcd6a1 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 27 Apr 2019 20:15:03 -0500 Subject: [PATCH 14/19] l2conn: clean up basis evaluation --- .../discretization/connection/__init__.py | 91 +++++-------------- test/test_chained.py | 8 +- 2 files changed, 28 insertions(+), 71 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 314b312d..800655c4 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -313,52 +313,13 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): to_discr=self.conn.from_discr, is_surjective=is_surjective) - def _evaluate_basis_on_target(self, queue, discr, igroup, ibasis_func): - """Evaluates a basis function from the target discretization - :attr:`to_discr` on all elements of `discr`. - - :arg discr: discretization on which to evaluate the basis function. - :arg igroup: group index in the array returned by - :attr:`to_discr.groups`. - :arg ibasis_func: basis index in the array returned by - :attr:`to_discr.basis()`. - """ - - @memoize_in(self, "conn_basis_evaluation_knl") - def knl(): - import loopy as lp - knl = lp.make_kernel([ - "{[i]: 0 <= i < nelements}", - "{[j]: 0 <= j < n_to_nodes}" - ], - """vec[i, j] = basis[j]""", - name="conn_basis_eval_kernel", - lang_version=MOST_RECENT_LANGUAGE_VERSION) - - return knl - - group = self.to_discr.groups[igroup] - basis = group.basis()[ibasis_func] - - vec = self.to_discr.zeros(queue) - knl()(queue, - vec=group.view(vec), - basis=basis(group.unit_nodes).flatten()) - - if discr is self.from_discr: - vec = self.conn(queue, vec) - - return vec - @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 refined elements in the domain of the parent element. - :returns: a list of lists of size `(ngroups, ngroup_batches)` - containing the quadrature weights (of size `batch.nunit_nodes`) - for a given batch in a given group of :attr:`conn`. + :return: a dictionary with keys ``(group_id, batch_id)`` tuples. """ from pymbolic.geometric_algebra import MultiVector @@ -375,18 +336,16 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): return det_v - to_discr = self.to_discr - jac = np.empty(to_discr.dim, dtype=np.object) - weights = [np.empty(len(g.batches), dtype=np.object) - for g in self.conn.groups] + weights = {} + jac = np.empty(self.to_discr.dim, dtype=np.object) - for igrp, grp in enumerate(to_discr.groups): + 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 + weights[igrp, ibatch] = det(jac) * grp.weights return weights @@ -402,26 +361,25 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): for k <> element_dot = \ sum(j, vec[from_element_indices[k], j] * \ - basis[from_element_indices[k], j] * \ - weights[j]) + basis[j] * weights[j]) - result[to_element_indices[k], ibasis_func] = \ - result[to_element_indices[k], ibasis_func] + element_dot + 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_elements", "n_from_nodes")), + shape="n_from_nodes"), lp.GlobalArg("weights", None, shape="n_from_nodes"), - lp.GlobalArg("result", None, - shape=("n_to_elements", "n_to_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_func", np.int32), + lp.ValueArg("ibasis", np.int32), '...' ], name="conn_projection_knl", @@ -438,12 +396,12 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): ], """ result[k, j] = result[k, j] + \ - coefficients[k, ibasis_func] * basis[k, j] + coefficients[k, ibasis] * basis[j] """, [ lp.GlobalArg("coefficients", None, shape=("nelements", "n_to_nodes")), - lp.ValueArg("ibasis_func", np.int32), + lp.ValueArg("ibasis", np.int32), '...' ], name="conn_evaluate_knl", @@ -465,22 +423,20 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): c = self.to_discr.zeros(queue, dtype=vec.dtype) for igrp, (tgrp, cgrp) in enumerate( zip(self.to_discr.groups, self.conn.groups)): - for ibasis_func in range(len(tgrp.basis())): - to_basis = self._evaluate_basis_on_target( - queue, self.from_discr, igrp, ibasis_func) - + for ibasis, basis_fn in enumerate(tgrp.basis()): for ibatch, batch in enumerate(cgrp.batches): sgrp = self.from_discr.groups[batch.from_group_index] + 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_func=ibasis_func, + ibasis=ibasis, vec=sgrp.view(vec), - basis=sgrp.view(to_basis), - weights=weights[igrp][ibatch], + basis=basis, + weights=weights[igrp, ibatch], result=tgrp.view(c), from_element_indices=batch.to_element_indices, to_element_indices=batch.from_element_indices) @@ -489,14 +445,13 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): result = self.to_discr.zeros(queue, dtype=vec.dtype) for igrp, grp in enumerate(self.to_discr.groups): - for ibasis_func in range(len(grp.basis())): - basis = self._evaluate_basis_on_target( - queue, self.to_discr, igrp, ibasis_func) + for ibasis, basis_fn in enumerate(grp.basis()): + basis = basis_fn(grp.unit_nodes).flatten() keval()(queue, - ibasis_func=ibasis_func, + ibasis=ibasis, result=grp.view(result), - basis=grp.view(basis), + basis=basis, coefficients=grp.view(c)) return result diff --git a/test/test_chained.py b/test/test_chained.py index fb186aed..2934fe2d 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -371,9 +371,11 @@ def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): from_x = 0.0 to_x = 0.0 - for i in range(ndim): - from_x += cl.clmath.cos(from_nodes[i] / (4.0 * np.pi)) ** (i + 1.0) - to_x += cl.clmath.cos(to_nodes[i] / (4.0 * np.pi)) ** (i + 1.0) + from_t = cl.array.to_device(queue, np.linspace(0.0, 1.0, from_nodes.shape[1])) + to_t = cl.array.to_device(queue, np.linspace(0.0, 1.0, to_nodes.shape[1])) + for i in range(2): + from_x += cl.clmath.cos(from_t / (4.0 * np.pi)) ** (i + 1.0) + to_x += cl.clmath.cos(to_t / (4.0 * np.pi)) ** (i + 1.0) from_interp = reverse(queue, to_x) -- GitLab From 826f83b93611db352b087f9f330e379adbd7338b Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 27 Apr 2019 21:26:33 -0500 Subject: [PATCH 15/19] l2conn: check proper convergence in tests --- .../discretization/connection/__init__.py | 11 +- test/test_chained.py | 152 +++++++++--------- 2 files changed, 78 insertions(+), 85 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 800655c4..0bb4e91e 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -296,11 +296,6 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): return ChainedDiscretizationConnection(conns) def __init__(self, conn, is_surjective=False): - if len(conn.to_discr.groups) != 1 or len(conn.from_discr.groups) != 1: - from warnings import warn - warn("multiple element groups are implemented in principle, " - "but not yet tested") - if conn.from_discr.dim != conn.to_discr.dim: raise RuntimeError("cannot transport from face to element") @@ -317,9 +312,10 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): def _batch_weights(self): """Computes scaled quadrature weights for each interpolation batch in :attr:`conn`. The quadrature weights can be used to integrate over - refined elements in the domain of the parent element. + child elements in the domain of the parent element, by a change of + variables. - :return: a dictionary with keys ``(group_id, batch_id)`` tuples. + :return: a dictionary with keys ``(group_id, batch_id)``. """ from pymbolic.geometric_algebra import MultiVector @@ -443,7 +439,6 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): # 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() diff --git a/test/test_chained.py b/test/test_chained.py index 2934fe2d..12f82f3e 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -40,29 +40,47 @@ logger = logging.getLogger(__name__) def create_discretization(queue, ndim, nelements=42, - mesh_order=5, - discr_order=5): + mesh_name=None, + order=4): ctx = queue.context # construct 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) + from meshmode.mesh.generation import make_curve_mesh, ellipse, starfish + + if mesh_name is None: + mesh_name = "ellipse" + + t = np.linspace(0.0, 1.0, nelements + 1) + if mesh_name == "ellipse": + mesh = make_curve_mesh(partial(ellipse, 2), t, order=order) + elif mesh_name == "starfish": + mesh = make_curve_mesh(starfish, t, order=order) + else: + raise ValueError('unknown mesh name: {}'.format(mesh_name)) elif ndim == 3: from meshmode.mesh.generation import generate_torus - mesh = generate_torus(10.0, 5.0, order=mesh_order) + from meshmode.mesh.generation import generate_warped_rect_mesh + + if mesh_name is None: + mesh_name = "torus" + + if mesh_name == "torus": + mesh = generate_torus(10.0, 5.0, order=order) + elif mesh_name == "warp": + mesh = generate_warped_rect_mesh(ndim, order=order, n=nelements) + else: + raise ValueError("unknown mesh name: {}".format(mesh_name)) else: - raise ValueError("Unsupported dimension: {}".format(ndim)) + raise ValueError("unsupported dimension: {}".format(ndim)) # create discretization from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ InterpolatoryQuadratureSimplexGroupFactory discr = Discretization(ctx, mesh, - InterpolatoryQuadratureSimplexGroupFactory(discr_order)) + InterpolatoryQuadratureSimplexGroupFactory(order)) return discr @@ -185,9 +203,7 @@ def test_chained_connection(ctx_factory, ndim, visualize=False): queue = cl.CommandQueue(ctx) discr = create_discretization(queue, ndim, - nelements=10, - mesh_order=5, - discr_order=5) + nelements=10) connections = [] conn = create_refined_connection(queue, discr, threshold=np.inf) connections.append(conn) @@ -219,9 +235,7 @@ def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - discr = create_discretization(queue, ndim, - mesh_order=2, - discr_order=2) + discr = create_discretization(queue, ndim) connections = [] conn = create_refined_connection(queue, discr) connections.append(conn) @@ -333,83 +347,67 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type, assert np.allclose(f1, f2) -@pytest.mark.parametrize("ndim", [2, 3]) -def test_reversed_chained_connection(ctx_factory, ndim, visualize=False): +@pytest.mark.parametrize(("ndim", "mesh_name"), [ + (2, "starfish"), + (3, "warp")]) +def test_reversed_chained_connection(ctx_factory, ndim, mesh_name, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) - if ndim == 2: - order = 6 - threshold = 0.3 - else: - order = 6 - threshold = 0.1 - - # build test connection - discr = create_discretization(queue, ndim, - mesh_order=order, - discr_order=order) - - connections = [] - conn = create_refined_connection(queue, discr, threshold=threshold) - connections.append(conn) - conn = create_refined_connection(queue, conn.to_discr, threshold=threshold) - connections.append(conn) - conn = create_refined_connection(queue, conn.to_discr, threshold=threshold) - connections.append(conn) + def run(nelements, order): + discr = create_discretization(queue, ndim, + nelements=nelements, + order=order, + mesh_name=mesh_name) - from meshmode.discretization.connection import \ - ChainedDiscretizationConnection - chained = ChainedDiscretizationConnection(connections) - from meshmode.discretization.connection import \ - L2ProjectionInverseDiscretizationConnection - reverse = L2ProjectionInverseDiscretizationConnection(chained) + threshold = 1.0 + connections = [] + conn = create_refined_connection(queue, discr, threshold=threshold) + connections.append(conn) + # conn = create_refined_connection(queue, conn.to_discr, threshold=threshold) + # connections.append(conn) + # conn = create_refined_connection(queue, conn.to_discr, threshold=threshold) + # connections.append(conn) - # create test vector - from_nodes = chained.from_discr.nodes().with_queue(queue) - to_nodes = chained.to_discr.nodes().with_queue(queue) + from meshmode.discretization.connection import \ + ChainedDiscretizationConnection + chained = ChainedDiscretizationConnection(connections) + from meshmode.discretization.connection import \ + L2ProjectionInverseDiscretizationConnection + reverse = L2ProjectionInverseDiscretizationConnection(chained) - from_x = 0.0 - to_x = 0.0 - from_t = cl.array.to_device(queue, np.linspace(0.0, 1.0, from_nodes.shape[1])) - to_t = cl.array.to_device(queue, np.linspace(0.0, 1.0, to_nodes.shape[1])) - for i in range(2): - from_x += cl.clmath.cos(from_t / (4.0 * np.pi)) ** (i + 1.0) - to_x += cl.clmath.cos(to_t / (4.0 * np.pi)) ** (i + 1.0) + # create test vector + nnodes = chained.from_discr.nnodes + from_t = cl.array.to_device(queue, np.linspace(0.0, 1.0, nnodes)) + nnodes = chained.to_discr.nnodes + to_t = cl.array.to_device(queue, np.linspace(0.0, 1.0, nnodes)) - from_interp = reverse(queue, to_x) + from_x = cl.clmath.cos(2.0 * np.pi * from_t) + to_x = cl.clmath.cos(2.0 * np.pi * to_t) - if visualize and ndim == 2: - import matplotlib.pyplot as pt + from_interp = reverse(queue, to_x) - from_t = np.linspace(0.0, 1.0, chained.from_discr.nnodes) - to_t = chained(queue, cl.array.to_device(queue, from_t)).get(queue) + from_interp = from_interp.get(queue) + from_x = from_x.get(queue) - pt.figure(figsize=(10, 8), dpi=300) - pt.plot(from_t, from_x.get(queue), '--', label="From") - pt.plot(to_t, to_x.get(queue), '--', label="To") - pt.plot(from_t, from_interp.get(queue), 'o-', label="Projection") - pt.legend() - pt.savefig("test_reverse_chained_conn.png") - pt.clf() - elif visualize and ndim == 3: - from meshmode.discretization.visualization import make_visualizer - vis = make_visualizer(queue, chained.from_discr, order) + return 1.0 / nelements, la.norm(from_interp - from_x) / la.norm(from_x) - vis.write_vtk_file("test_reverse_chained_conn.vtu", [ - ("x_interp", from_interp), - ("x", from_x) - ], overwrite=True) + from pytools.convergence import EOCRecorder + eoc = EOCRecorder() - from_interp = from_interp.get(queue) - from_x = from_x.get(queue) + order = 4 + if ndim == 2: + mesh_sizes = [32, 64, 128, 256, 512] + else: + mesh_sizes = [2, 4, 6, 8, 10, 12, 14, 16] - error = la.norm(from_interp - from_x) / la.norm(from_x) - if visualize: - print('Error: {}'.format(error)) + for n in mesh_sizes: + h, error = run(n, order) + eoc.add_data_point(h, error) - assert error < 1.0e-9, error + print(eoc) + assert eoc.order_estimate() > (order - 0.5) if __name__ == "__main__": import sys -- GitLab From 4c715a7bc38b4c3d168d9519197a0ed3ce02eedb Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 27 Apr 2019 23:28:42 -0500 Subject: [PATCH 16/19] l2conn: flake8 --- test/test_chained.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/test/test_chained.py b/test/test_chained.py index 12f82f3e..55821e1d 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -362,12 +362,17 @@ def test_reversed_chained_connection(ctx_factory, ndim, mesh_name, visualize=Fal threshold = 1.0 connections = [] - conn = create_refined_connection(queue, discr, threshold=threshold) + conn = create_refined_connection(queue, + discr, threshold=threshold) connections.append(conn) - # conn = create_refined_connection(queue, conn.to_discr, threshold=threshold) - # connections.append(conn) - # conn = create_refined_connection(queue, conn.to_discr, threshold=threshold) - # connections.append(conn) + if ndim == 2: + # NOTE: additional refinement makes the 3D meshes explode in size + conn = create_refined_connection(queue, + conn.to_discr, threshold=threshold) + connections.append(conn) + conn = create_refined_connection(queue, + conn.to_discr, threshold=threshold) + connections.append(conn) from meshmode.discretization.connection import \ ChainedDiscretizationConnection @@ -409,6 +414,7 @@ def test_reversed_chained_connection(ctx_factory, ndim, mesh_name, visualize=Fal assert eoc.order_estimate() > (order - 0.5) + if __name__ == "__main__": import sys if len(sys.argv) > 1: -- GitLab From f1fc3107ffb022547fbb125b87ff91268f563270 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 2 May 2019 20:50:32 -0500 Subject: [PATCH 17/19] l2conn: use correct basis --- meshmode/discretization/connection/__init__.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 0bb4e91e..cc0417e6 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -419,9 +419,10 @@ class L2ProjectionInverseDiscretizationConnection(DiscretizationConnection): c = self.to_discr.zeros(queue, dtype=vec.dtype) for igrp, (tgrp, cgrp) in enumerate( zip(self.to_discr.groups, self.conn.groups)): - for ibasis, basis_fn in enumerate(tgrp.basis()): - for ibatch, batch in enumerate(cgrp.batches): - sgrp = self.from_discr.groups[batch.from_group_index] + 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 -- GitLab From 5508a122d3a54d23186c9e3b3f0d1213b29f1af2 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 5 May 2019 20:09:25 -0500 Subject: [PATCH 18/19] l2conn: fix h-convergence --- test/test_chained.py | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/test/test_chained.py b/test/test_chained.py index 55821e1d..6e8eb50d 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -67,7 +67,8 @@ def create_discretization(queue, ndim, mesh_name = "torus" if mesh_name == "torus": - mesh = generate_torus(10.0, 5.0, order=order) + mesh = generate_torus(10.0, 5.0, order=order, + n_inner=nelements, n_outer=nelements) elif mesh_name == "warp": mesh = generate_warped_rect_mesh(ndim, order=order, n=nelements) else: @@ -349,7 +350,7 @@ def test_chained_to_direct(ctx_factory, ndim, chain_type, @pytest.mark.parametrize(("ndim", "mesh_name"), [ (2, "starfish"), - (3, "warp")]) + (3, "torus")]) def test_reversed_chained_connection(ctx_factory, ndim, mesh_name, visualize=False): ctx = ctx_factory() queue = cl.CommandQueue(ctx) @@ -382,13 +383,14 @@ def test_reversed_chained_connection(ctx_factory, ndim, mesh_name, visualize=Fal reverse = L2ProjectionInverseDiscretizationConnection(chained) # create test vector - nnodes = chained.from_discr.nnodes - from_t = cl.array.to_device(queue, np.linspace(0.0, 1.0, nnodes)) - nnodes = chained.to_discr.nnodes - to_t = cl.array.to_device(queue, np.linspace(0.0, 1.0, nnodes)) + from_nodes = chained.from_discr.nodes().with_queue(queue) + to_nodes = chained.to_discr.nodes().with_queue(queue) - from_x = cl.clmath.cos(2.0 * np.pi * from_t) - to_x = cl.clmath.cos(2.0 * np.pi * to_t) + from_x = 0 + to_x = 0 + for d in range(ndim): + from_x += cl.clmath.cos(from_nodes[d]) ** (d + 1) + to_x += cl.clmath.cos(to_nodes[d]) ** (d + 1) from_interp = reverse(queue, to_x) @@ -401,10 +403,7 @@ def test_reversed_chained_connection(ctx_factory, ndim, mesh_name, visualize=Fal eoc = EOCRecorder() order = 4 - if ndim == 2: - mesh_sizes = [32, 64, 128, 256, 512] - else: - mesh_sizes = [2, 4, 6, 8, 10, 12, 14, 16] + mesh_sizes = [16, 32, 48, 64, 96, 128] for n in mesh_sizes: h, error = run(n, order) -- GitLab From a5724e9a1feb6da2b003643f5c5e2bff0a7ba461 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Tue, 7 May 2019 20:27:17 +0200 Subject: [PATCH 19/19] Apply suggestion to test/test_chained.py --- 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 6e8eb50d..76c416f3 100644 --- a/test/test_chained.py +++ b/test/test_chained.py @@ -411,7 +411,7 @@ def test_reversed_chained_connection(ctx_factory, ndim, mesh_name, visualize=Fal print(eoc) - assert eoc.order_estimate() > (order - 0.5) + assert eoc.order_estimate() > (order + 1 - 0.5) if __name__ == "__main__": -- GitLab