diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index ee1efd0c3c7078c298894c3825a9787532b3b3a3..cc0417e6f78559762f506c441bb4a501acb65ff6 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 @@ -65,6 +68,7 @@ Base classes ------------ .. autoclass:: DiscretizationConnection .. autoclass:: ChainedDiscretizationConnection +.. autoclass:: L2ProjectionInverseDiscretizationConnection .. autoclass:: DirectDiscretizationConnection @@ -265,6 +269,190 @@ class ChainedDiscretizationConnection(DiscretizationConnection): 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. diff --git a/test/test_chained.py b/test/test_chained.py index c0f5a09e320ce4442f1747a35c06718f4067d1dd..76c416f37f96abde8e8b43404acf8d85d4c39173 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 @@ -39,29 +40,48 @@ 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, + n_inner=nelements, n_outer=nelements) + 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 @@ -184,9 +204,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) @@ -218,9 +236,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) @@ -321,17 +337,83 @@ 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", "mesh_name"), [ + (2, "starfish"), + (3, "torus")]) +def test_reversed_chained_connection(ctx_factory, ndim, mesh_name, visualize=False): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + def run(nelements, order): + discr = create_discretization(queue, ndim, + nelements=nelements, + order=order, + mesh_name=mesh_name) + + threshold = 1.0 + connections = [] + conn = create_refined_connection(queue, + 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 + chained = ChainedDiscretizationConnection(connections) + from meshmode.discretization.connection import \ + L2ProjectionInverseDiscretizationConnection + reverse = L2ProjectionInverseDiscretizationConnection(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 + 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) + + from_interp = from_interp.get(queue) + from_x = from_x.get(queue) + + return 1.0 / nelements, la.norm(from_interp - from_x) / la.norm(from_x) + + from pytools.convergence import EOCRecorder + eoc = EOCRecorder() + + order = 4 + mesh_sizes = [16, 32, 48, 64, 96, 128] + + for n in mesh_sizes: + h, error = run(n, order) + eoc.add_data_point(h, error) + + print(eoc) + + assert eoc.order_estimate() > (order + 1 - 0.5) + + if __name__ == "__main__": import sys if len(sys.argv) > 1: