test/test_metrics.py

+__copyright__ = """
+Copyright (C) 2015 Andreas Kloeckner
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import logging
+
+import numpy as np
+import pytest
+
+import meshmode.mesh.generation as mgen
+from arraycontext import pytest_generate_tests_for_array_contexts
+from meshmode.dof_array import flat_norm
+
+from grudge.array_context import (
+    PytestNumpyArrayContextFactory,
+    PytestPyOpenCLArrayContextFactory,
+    PytestPytatoPyOpenCLArrayContextFactory,
+)
+from grudge.discretization import make_discretization_collection
+
+
+logger = logging.getLogger(__name__)
+pytest_generate_tests = pytest_generate_tests_for_array_contexts(
+        [PytestPyOpenCLArrayContextFactory,
+         PytestPytatoPyOpenCLArrayContextFactory,
+         PytestNumpyArrayContextFactory])
+
+
+# {{{ inverse metric
+
+@pytest.mark.parametrize("dim", [2, 3])
+@pytest.mark.parametrize("nonaffine", [False, True])
+@pytest.mark.parametrize("use_quad", [False, True])
+def test_inverse_metric(actx_factory, dim, nonaffine, use_quad):
+    actx = actx_factory()
+
+    order = 3
+    mesh = mgen.generate_regular_rect_mesh(a=(-0.5,)*dim, b=(0.5,)*dim,
+            nelements_per_axis=(6,)*dim, order=order)
+
+    if nonaffine:
+        def m(x):
+            result = np.empty_like(x)
+            result[0] = (
+                    1.5*x[0] + np.cos(x[0])
+                    + 0.1*np.sin(10*x[1]))
+            result[1] = (
+                    0.05*np.cos(10*x[0])
+                    + 1.3*x[1] + np.sin(x[1]))
+            if len(x) == 3:
+                result[2] = x[2]
+            return result
+
+        from meshmode.mesh.processing import map_mesh
+        mesh = map_mesh(mesh, m)
+
+    from meshmode.discretization.poly_element import (
+        QuadratureSimplexGroupFactory,
+        default_simplex_group_factory,
+    )
+
+    from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD, as_dofdesc
+
+    dcoll = make_discretization_collection(
+        actx, mesh,
+        discr_tag_to_group_factory={
+            DISCR_TAG_BASE: default_simplex_group_factory(base_dim=dim, order=order),
+            DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order + 1),
+        }
+    )
+
+    from grudge.geometry import (
+        forward_metric_derivative_mat,
+        inverse_metric_derivative_mat,
+    )
+
+    dd = as_dofdesc("vol")
+    if use_quad:
+        dd = dd.with_discr_tag(DISCR_TAG_QUAD)
+
+    mat = forward_metric_derivative_mat(
+        actx, dcoll, dd,
+        _use_geoderiv_connection=actx.supports_nonscalar_broadcasting).dot(
+        inverse_metric_derivative_mat(
+            actx, dcoll, dd,
+            _use_geoderiv_connection=actx.supports_nonscalar_broadcasting))
+
+    for i in range(mesh.dim):
+        for j in range(mesh.dim):
+            tgt = 1 if i == j else 0
+
+            err = actx.to_numpy(flat_norm(mat[i, j] - tgt, ord=np.inf))
+            logger.info("error[%d, %d]: %.5e", i, j, err)
+            assert err < 1.0e-12, (i, j, err)
+
+# }}}
+
+
+# You can test individual routines by typing
+# $ python test_metrics.py 'test_routine()'
+
+if __name__ == "__main__":
+    import sys
+    if len(sys.argv) > 1:
+        exec(sys.argv[1])
+    else:
+        pytest.main([__file__])
+
+# vim: fdm=marker

test/test_modal_connections.py

+__copyright__ = "Copyright (C) 2021 University of Illinois Board of Trustees"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+from arraycontext import pytest_generate_tests_for_array_contexts
+
+from grudge.array_context import PytestPyOpenCLArrayContextFactory
+from grudge.discretization import make_discretization_collection
+
+
+pytest_generate_tests = pytest_generate_tests_for_array_contexts(
+        [PytestPyOpenCLArrayContextFactory])
+
+import pytest
+
+import meshmode.mesh.generation as mgen
+from meshmode.discretization.poly_element import (
+    # Simplex group factories
+    InterpolatoryQuadratureSimplexGroupFactory,
+    # Tensor product group factories
+    LegendreGaussLobattoTensorProductGroupFactory,
+    ModalGroupFactory,
+    PolynomialEquidistantSimplexGroupFactory,
+    PolynomialWarpAndBlend2DRestrictingGroupFactory,
+    # Quadrature-based (non-interpolatory) group factories
+    QuadratureSimplexGroupFactory,
+)
+from meshmode.dof_array import flat_norm
+
+import grudge.dof_desc as dof_desc
+
+
+@pytest.mark.parametrize("nodal_group_factory", [
+    InterpolatoryQuadratureSimplexGroupFactory,
+    PolynomialWarpAndBlend2DRestrictingGroupFactory,
+    PolynomialEquidistantSimplexGroupFactory,
+    LegendreGaussLobattoTensorProductGroupFactory,
+    ]
+)
+def test_inverse_modal_connections(actx_factory, nodal_group_factory):
+    actx = actx_factory()
+    order = 4
+
+    def f(x):
+        return 2*actx.np.sin(20*x) + 0.5*actx.np.cos(10*x)
+
+    # Make a regular rectangle mesh
+    mesh = mgen.generate_regular_rect_mesh(
+        a=(0, 0), b=(5, 3), npoints_per_axis=(10, 6), order=order,
+        group_cls=nodal_group_factory.mesh_group_class
+    )
+
+    dcoll = make_discretization_collection(
+        actx, mesh,
+        discr_tag_to_group_factory={
+            dof_desc.DISCR_TAG_BASE: nodal_group_factory(order),
+            dof_desc.DISCR_TAG_MODAL: ModalGroupFactory(order),
+        }
+    )
+
+    dd_modal = dof_desc.DD_VOLUME_ALL_MODAL
+    dd_volume = dof_desc.DD_VOLUME_ALL
+
+    x_nodal = actx.thaw(dcoll.discr_from_dd(dd_volume).nodes()[0])
+    nodal_f = f(x_nodal)
+
+    # Map nodal coefficients of f to modal coefficients
+    forward_conn = dcoll.connection_from_dds(dd_volume, dd_modal)
+    modal_f = forward_conn(nodal_f)
+    # Now map the modal coefficients back to nodal
+    backward_conn = dcoll.connection_from_dds(dd_modal, dd_volume)
+    nodal_f_2 = backward_conn(modal_f)
+
+    # This error should be small since we composed a map with
+    # its inverse
+    err = flat_norm(nodal_f - nodal_f_2)
+
+    assert err <= 1e-13
+
+
+def test_inverse_modal_connections_quadgrid(actx_factory):
+    actx = actx_factory()
+    order = 5
+
+    def f(x):
+        return 1 + 2*x + 3*x**2
+
+    # Make a regular rectangle mesh
+    mesh = mgen.generate_regular_rect_mesh(
+        a=(0, 0), b=(5, 3), npoints_per_axis=(10, 6), order=order,
+        group_cls=QuadratureSimplexGroupFactory.mesh_group_class
+    )
+
+    dcoll = make_discretization_collection(
+        actx, mesh,
+        discr_tag_to_group_factory={
+            dof_desc.DISCR_TAG_BASE:
+            PolynomialWarpAndBlend2DRestrictingGroupFactory(order),
+            dof_desc.DISCR_TAG_QUAD: QuadratureSimplexGroupFactory(2*order),
+            dof_desc.DISCR_TAG_MODAL: ModalGroupFactory(order),
+        }
+    )
+
+    # Use dof descriptors on the quadrature grid
+    dd_modal = dof_desc.DD_VOLUME_ALL_MODAL
+    dd_quad = dof_desc.DOFDesc(dof_desc.DTAG_VOLUME_ALL,
+                               dof_desc.DISCR_TAG_QUAD)
+
+    x_quad = actx.thaw(dcoll.discr_from_dd(dd_quad).nodes()[0])
+    quad_f = f(x_quad)
+
+    # Map nodal coefficients of f to modal coefficients
+    forward_conn = dcoll.connection_from_dds(dd_quad, dd_modal)
+    modal_f = forward_conn(quad_f)
+    # Now map the modal coefficients back to nodal
+    backward_conn = dcoll.connection_from_dds(dd_modal, dd_quad)
+    quad_f_2 = backward_conn(modal_f)
+
+    # This error should be small since we composed a map with
+    # its inverse
+    err = flat_norm(quad_f - quad_f_2)
+
+    assert err <= 1e-11

test/test_mpi_communication.py

+__copyright__ = """
+Copyright (C) 2017 Ellis Hoag
+Copyright (C) 2017 Andreas Kloeckner
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import logging
+import os
+import sys
+
+import numpy as np
+import pytest
+
+from meshmode.dof_array import flat_norm
+from pytools.obj_array import flat_obj_array
+
+from grudge import dof_desc, op
+from grudge.array_context import MPIPyOpenCLArrayContext, MPIPytatoArrayContext
+from grudge.discretization import make_discretization_collection
+from grudge.shortcuts import rk4_step
+
+
+logger = logging.getLogger(__name__)
+
+
+class SimpleTag:
+    pass
+
+
+# {{{ mpi test infrastructure
+
+DISTRIBUTED_ACTXS = [MPIPyOpenCLArrayContext, MPIPytatoArrayContext]
+
+
+def run_test_with_mpi(num_ranks, f, *args):
+    import pytest
+    pytest.importorskip("mpi4py")
+
+    from base64 import b64encode
+    from pickle import dumps
+
+    invocation_info = b64encode(dumps((f, args))).decode()
+    from subprocess import check_call
+
+    # NOTE: CI uses OpenMPI; -x to pass env vars. MPICH uses -env
+    check_call([
+        "mpiexec", "-np", str(num_ranks),
+        "-x", "RUN_WITHIN_MPI=1",
+        "-x", f"INVOCATION_INFO={invocation_info}",
+        sys.executable, "-m", "mpi4py", __file__])
+
+
+def run_test_with_mpi_inner():
+    from base64 import b64decode
+    from pickle import loads
+    f, (actx_class, *args) = loads(b64decode(os.environ["INVOCATION_INFO"].encode()))
+
+    import pyopencl as cl
+
+    cl_context = cl.create_some_context()
+    queue = cl.CommandQueue(cl_context)
+
+    from mpi4py import MPI
+    comm = MPI.COMM_WORLD
+
+    if actx_class is MPIPytatoArrayContext:
+        actx = actx_class(comm, queue, mpi_base_tag=15000)
+    elif actx_class is MPIPyOpenCLArrayContext:
+        actx = actx_class(comm, queue)
+    else:
+        raise ValueError("unknown actx_class")
+
+    f(actx, *args)
+
+# }}}
+
+
+# {{{ func_comparison
+
+@pytest.mark.parametrize("actx_class", DISTRIBUTED_ACTXS)
+@pytest.mark.parametrize("num_ranks", [2])
+def test_func_comparison_mpi(actx_class, num_ranks):
+    run_test_with_mpi(
+            num_ranks, _test_func_comparison_mpi_communication_entrypoint,
+            actx_class)
+
+
+def _test_func_comparison_mpi_communication_entrypoint(actx):
+    """Discretize a function, communicate it, check that it matches the
+    function discretized by the other end.
+    """
+
+    comm = actx.mpi_communicator
+
+    from meshmode.distributed import get_partition_by_pymetis, membership_list_to_map
+    from meshmode.mesh import BTAG_ALL
+    from meshmode.mesh.processing import partition_mesh
+
+    num_parts = comm.size
+
+    if comm.rank == 0:
+        from meshmode.mesh.generation import generate_regular_rect_mesh
+        mesh = generate_regular_rect_mesh(a=(-1,)*2,
+                                          b=(1,)*2,
+                                          nelements_per_axis=(2,)*2)
+
+        part_id_to_part = partition_mesh(mesh,
+                       membership_list_to_map(
+                           get_partition_by_pymetis(mesh, num_parts)))
+        parts = [part_id_to_part[i] for i in range(num_parts)]
+        local_mesh = comm.scatter(parts)
+    else:
+        local_mesh = comm.scatter(None)
+
+    dcoll = make_discretization_collection(actx, local_mesh, order=5)
+
+    x = actx.thaw(dcoll.nodes())
+    myfunc = actx.np.sin(np.dot(x, [2, 3]))
+
+    dd_int = dof_desc.as_dofdesc("int_faces")
+    dd_vol = dof_desc.as_dofdesc("vol")
+    dd_af = dof_desc.as_dofdesc("all_faces")
+
+    all_faces_func = op.project(dcoll, dd_vol, dd_af, myfunc)
+    int_faces_func = op.project(dcoll, dd_vol, dd_int, myfunc)
+    bdry_faces_func = op.project(dcoll, BTAG_ALL, dd_af,
+                                 op.project(dcoll, dd_vol, BTAG_ALL, myfunc))
+
+    def hopefully_zero():
+        return (
+            op.project(
+                dcoll, "int_faces", "all_faces",
+                dcoll.opposite_face_connection(
+                    dof_desc.BoundaryDomainTag(
+                        dof_desc.FACE_RESTR_INTERIOR, dof_desc.VTAG_ALL)
+                    )(int_faces_func)
+            )
+            + sum(op.project(dcoll, tpair.dd, "all_faces", tpair.ext)
+                  for tpair in op.cross_rank_trace_pairs(dcoll, myfunc,
+                      comm_tag=SimpleTag))
+        ) - (all_faces_func - bdry_faces_func)
+
+    hopefully_zero_result = actx.compile(hopefully_zero)()
+
+    error = actx.to_numpy(flat_norm(hopefully_zero_result, ord=np.inf))
+
+    with np.printoptions(threshold=100000000, suppress=True):
+        logger.debug(hopefully_zero)
+    logger.info("error: %.5e", error)
+
+    assert error < 1e-14
+
+
+# }}}
+
+
+# {{{ wave operator
+
+@pytest.mark.parametrize("actx_class", DISTRIBUTED_ACTXS)
+@pytest.mark.parametrize("num_ranks", [2])
+def test_mpi_wave_op(actx_class, num_ranks):
+    run_test_with_mpi(num_ranks, _test_mpi_wave_op_entrypoint, actx_class)
+
+
+def _test_mpi_wave_op_entrypoint(actx, visualize=False):
+    comm = actx.mpi_communicator
+    num_parts = comm.size
+
+    from meshmode.distributed import get_partition_by_pymetis, membership_list_to_map
+    from meshmode.mesh.processing import partition_mesh
+
+    dim = 2
+    order = 4
+
+    if comm.rank == 0:
+        from meshmode.mesh.generation import generate_regular_rect_mesh
+        mesh = generate_regular_rect_mesh(a=(-0.5,)*dim,
+                                          b=(0.5,)*dim,
+                                          nelements_per_axis=(16,)*dim)
+
+        part_id_to_part = partition_mesh(mesh,
+                       membership_list_to_map(
+                           get_partition_by_pymetis(mesh, num_parts)))
+        parts = [part_id_to_part[i] for i in range(num_parts)]
+        local_mesh = comm.scatter(parts)
+
+        del mesh
+    else:
+        local_mesh = comm.scatter(None)
+
+    dcoll = make_discretization_collection(actx, local_mesh, order=order)
+
+    def source_f(actx, dcoll, t=0):
+        source_center = np.array([0.1, 0.22, 0.33])[:dcoll.dim]
+        source_width = 0.05
+        source_omega = 3
+        nodes = actx.thaw(dcoll.nodes())
+        source_center_dist = flat_obj_array(
+            [nodes[i] - source_center[i] for i in range(dcoll.dim)]
+        )
+        return (
+            actx.np.sin(source_omega*t)
+            * actx.np.exp(
+                -np.dot(source_center_dist, source_center_dist)
+                / source_width**2
+            )
+        )
+
+    from meshmode.mesh import BTAG_ALL, BTAG_NONE
+
+    from grudge.models.wave import WeakWaveOperator
+
+    wave_op = WeakWaveOperator(
+        dcoll,
+        0.1,
+        source_f=source_f,
+        dirichlet_tag=BTAG_NONE,
+        neumann_tag=BTAG_NONE,
+        radiation_tag=BTAG_ALL,
+        flux_type="upwind",
+        comm_tag=SimpleTag,
+    )
+
+    fields = flat_obj_array(
+        dcoll.zeros(actx),
+        [dcoll.zeros(actx) for i in range(dcoll.dim)]
+    )
+
+    dt = actx.to_numpy(
+        wave_op.estimate_rk4_timestep(actx, dcoll, fields=fields))
+
+    wave_op.check_bc_coverage(local_mesh)
+
+    from logpyle import LogManager, add_general_quantities, add_run_info
+    log_filename = None
+    # NOTE: LogManager hangs when using a file on a shared directory.
+    # log_filename = "grudge_log.dat"
+    logmgr = LogManager(log_filename, "w", comm)
+    add_run_info(logmgr)
+    add_general_quantities(logmgr)
+
+    def rhs(t, w):
+        return wave_op.operator(t, w)
+
+    compiled_rhs = actx.compile(rhs)
+
+    final_t = 4
+    nsteps = int(final_t/dt)
+    logger.info("[%04d] dt %.5e nsteps %4d", comm.rank, dt, nsteps)
+
+    step = 0
+
+    from time import time
+    t_last_step = time()
+
+    if visualize:
+        from grudge.shortcuts import make_visualizer
+        vis = make_visualizer(dcoll)
+
+    logmgr.tick_before()
+    for step in range(nsteps):
+        t = step*dt
+        fields = rk4_step(fields, t=t, h=dt, f=compiled_rhs)
+        fields = actx.thaw(actx.freeze(fields))
+
+        norm = actx.to_numpy(op.norm(dcoll, fields, 2))
+        logger.info("[%04d] t = %.5e |u| = %.5e elapsed %.5e",
+                    step, t, norm, time() - t_last_step)
+
+        if visualize:
+            vis.write_parallel_vtk_file(
+                comm,
+                f"fld-wave-mpi-{type(actx).__name__}-{{rank:03d}}-{step:04d}.vtu",
+                [
+                    ("u", fields[0]),
+                    ("v", fields[1:]),
+                ]
+            )
+        assert norm < 1
+
+        t_last_step = time()
+        logmgr.tick_after()
+        logmgr.tick_before()
+
+    logmgr.tick_after()
+    logmgr.close()
+    logger.info("Rank %d exiting", comm.rank)
+
+# }}}
+
+
+if __name__ == "__main__":
+    if "RUN_WITHIN_MPI" in os.environ:
+        run_test_with_mpi_inner()
+    elif len(sys.argv) > 1:
+        exec(sys.argv[1])
+    else:
+        from pytest import main
+        main([__file__])
+
+# vim: fdm=marker

test/test_op.py

+__copyright__ = "Copyright (C) 2021 University of Illinois Board of Trustees"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+import logging
+
+import numpy as np
+import pytest
+
+import meshmode.mesh.generation as mgen
+from arraycontext import pytest_generate_tests_for_array_contexts
+from meshmode.discretization.poly_element import (
+    InterpolatoryEdgeClusteredGroupFactory,
+    QuadratureGroupFactory,
+)
+from meshmode.mesh import BTAG_ALL
+from pytools.obj_array import make_obj_array
+
+from grudge import geometry, op
+from grudge.array_context import PytestPyOpenCLArrayContextFactory
+from grudge.discretization import make_discretization_collection
+from grudge.dof_desc import (
+    DISCR_TAG_BASE,
+    DISCR_TAG_QUAD,
+    DTAG_VOLUME_ALL,
+    FACE_RESTR_ALL,
+    VTAG_ALL,
+    BoundaryDomainTag,
+    as_dofdesc,
+)
+from grudge.trace_pair import bv_trace_pair
+
+
+logger = logging.getLogger(__name__)
+pytest_generate_tests = pytest_generate_tests_for_array_contexts(
+        [PytestPyOpenCLArrayContextFactory])
+
+
+# {{{ gradient
+
+@pytest.mark.parametrize("form", ["strong", "weak", "weak-overint"])
+@pytest.mark.parametrize("dim", [1, 2, 3])
+@pytest.mark.parametrize("order", [2, 3])
+@pytest.mark.parametrize("warp_mesh", [False, True])
+@pytest.mark.parametrize(("vectorize", "nested"), [
+    (False, False),
+    (True, False),
+    (True, True)
+    ])
+def test_gradient(actx_factory, form, dim, order, vectorize, nested,
+                  warp_mesh, visualize=False):
+    actx = actx_factory()
+
+    from pytools.convergence import EOCRecorder
+    eoc_rec = EOCRecorder()
+
+    for n in [8, 12, 16] if warp_mesh else [4, 6, 8]:
+        if warp_mesh:
+            if dim == 1:
+                pytest.skip("warped mesh in 1D not implemented")
+            mesh = mgen.generate_warped_rect_mesh(
+                          dim=dim, order=order, nelements_side=n)
+        else:
+            mesh = mgen.generate_regular_rect_mesh(
+                    a=(-1,)*dim, b=(1,)*dim,
+                    nelements_per_axis=(n,)*dim)
+
+        dcoll = make_discretization_collection(
+                   actx, mesh,
+                   discr_tag_to_group_factory={
+                       DISCR_TAG_BASE: InterpolatoryEdgeClusteredGroupFactory(order),
+                       DISCR_TAG_QUAD: QuadratureGroupFactory(3 * order)
+                   })
+
+        def f(x):
+            result = 1
+            for i in range(dim-1):
+                result = result * actx.np.sin(np.pi*x[i])
+            result = result * actx.np.cos(np.pi/2*x[dim-1])
+            return result
+
+        def grad_f(x):
+            result = make_obj_array([1 for _ in range(dim)])
+            for i in range(dim-1):
+                for j in range(i):
+                    result[i] = result[i] * actx.np.sin(np.pi*x[j])
+                result[i] = result[i] * np.pi*actx.np.cos(np.pi*x[i])
+                for j in range(i+1, dim-1):
+                    result[i] = result[i] * actx.np.sin(np.pi*x[j])
+                result[i] = result[i] * actx.np.cos(np.pi/2*x[dim-1])
+            for j in range(dim-1):
+                result[dim-1] = result[dim-1] * actx.np.sin(np.pi*x[j])
+            result[dim-1] = result[dim-1] * (-np.pi/2*actx.np.sin(np.pi/2*x[dim-1]))
+            return result
+
+        def vectorize_if_requested(vec):
+            if vectorize:
+                return make_obj_array([(i+1)*vec for i in range(dim)])
+            else:
+                return vec
+
+        def get_flux(u_tpair, dcoll=dcoll):
+            dd = u_tpair.dd
+            dd_allfaces = dd.with_domain_tag(
+                BoundaryDomainTag(FACE_RESTR_ALL, VTAG_ALL)
+                )
+            normal = geometry.normal(actx, dcoll, dd)
+            u_avg = u_tpair.avg
+            if vectorize:
+                if nested:
+                    flux = make_obj_array([u_avg_i * normal for u_avg_i in u_avg])
+                else:
+                    flux = np.outer(u_avg, normal)
+            else:
+                flux = u_avg * normal
+            return op.project(dcoll, dd, dd_allfaces, flux)
+
+        x = actx.thaw(dcoll.nodes())
+        u = vectorize_if_requested(f(x))
+
+        bdry_dd_base = as_dofdesc(BTAG_ALL)
+        bdry_x = actx.thaw(dcoll.nodes(bdry_dd_base))
+        bdry_u = vectorize_if_requested(f(bdry_x))
+
+        if form == "strong":
+            grad_u = (
+                op.local_grad(dcoll, u, nested=nested)
+                # No flux terms because u doesn't have inter-el jumps
+                )
+        elif form.startswith("weak"):
+            assert form in ["weak", "weak-overint"]
+            if "overint" in form:
+                quad_discr_tag = DISCR_TAG_QUAD
+            else:
+                quad_discr_tag = DISCR_TAG_BASE
+
+            allfaces_dd_base = as_dofdesc(FACE_RESTR_ALL, quad_discr_tag)
+            vol_dd_base = as_dofdesc(DTAG_VOLUME_ALL)
+            vol_dd_quad = vol_dd_base.with_discr_tag(quad_discr_tag)
+            bdry_dd_quad = bdry_dd_base.with_discr_tag(quad_discr_tag)
+            allfaces_dd_quad = allfaces_dd_base.with_discr_tag(quad_discr_tag)
+
+            grad_u = op.inverse_mass(dcoll,
+                -op.weak_local_grad(dcoll, vol_dd_quad,
+                        op.project(dcoll, vol_dd_base, vol_dd_quad, u),
+                        nested=nested)
+                +
+                op.face_mass(dcoll,
+                    allfaces_dd_quad,
+                    sum(get_flux(
+                        op.project_tracepair(dcoll, allfaces_dd_quad, utpair))
+                        for utpair in op.interior_trace_pairs(
+                                      dcoll, u, volume_dd=vol_dd_base))
+                    + get_flux(
+                        op.project_tracepair(dcoll, bdry_dd_quad,
+                                   bv_trace_pair(dcoll, bdry_dd_base, u, bdry_u)))
+                )
+            )
+        else:
+            raise ValueError("Invalid form argument.")
+
+        if vectorize:
+            expected_grad_u = make_obj_array(
+                [(i+1)*grad_f(x) for i in range(dim)])
+            if not nested:
+                expected_grad_u = np.stack(expected_grad_u, axis=0)
+        else:
+            expected_grad_u = grad_f(x)
+
+        if visualize:
+            # the code below does not handle the vectorized case
+            assert not vectorize
+
+            from grudge.shortcuts import make_visualizer
+            vis = make_visualizer(dcoll, vis_order=order if dim == 3 else dim+3)
+
+            filename = (f"test_gradient_{form}_{dim}_{order}"
+                f"{'_vec' if vectorize else ''}{'_nested' if nested else ''}.vtu")
+            vis.write_vtk_file(filename, [
+                ("u", u),
+                ("grad_u", grad_u),
+                ("expected_grad_u", expected_grad_u),
+                ], overwrite=True)
+
+        rel_linf_err = actx.to_numpy(
+            op.norm(dcoll, grad_u - expected_grad_u, np.inf)
+            / op.norm(dcoll, expected_grad_u, np.inf))
+        eoc_rec.add_data_point(1./n, rel_linf_err)
+
+    print("L^inf error:")
+    print(eoc_rec)
+    assert (eoc_rec.order_estimate() >= order - 0.5
+                or eoc_rec.max_error() < 1e-11)
+
+# }}}
+
+
+# {{{ divergence
+
+@pytest.mark.parametrize("form", ["strong", "weak"])
+@pytest.mark.parametrize("dim", [1, 2, 3])
+@pytest.mark.parametrize("order", [2, 3])
+@pytest.mark.parametrize(("vectorize", "nested"), [
+    (False, False),
+    (True, False),
+    (True, True)
+    ])
+def test_divergence(actx_factory, form, dim, order, vectorize, nested,
+        visualize=False):
+    actx = actx_factory()
+
+    from pytools.convergence import EOCRecorder
+    eoc_rec = EOCRecorder()
+
+    for n in [4, 6, 8]:
+        mesh = mgen.generate_regular_rect_mesh(
+                a=(-1,)*dim, b=(1,)*dim,
+                nelements_per_axis=(n,)*dim)
+
+        dcoll = make_discretization_collection(actx, mesh, order=order)
+
+        def f(x, dcoll=dcoll):
+            result = make_obj_array([dcoll.zeros(actx) + (i+1) for i in range(dim)])
+            for i in range(dim-1):
+                result = result * actx.np.sin(np.pi*x[i])
+            result = result * actx.np.cos(np.pi/2*x[dim-1])
+            return result
+
+        def div_f(x, dcoll=dcoll):
+            result = dcoll.zeros(actx)
+            for i in range(dim-1):
+                deriv = dcoll.zeros(actx) + (i+1)
+                for j in range(i):
+                    deriv = deriv * actx.np.sin(np.pi*x[j])
+                deriv = deriv * np.pi*actx.np.cos(np.pi*x[i])
+                for j in range(i+1, dim-1):
+                    deriv = deriv * actx.np.sin(np.pi*x[j])
+                deriv = deriv * actx.np.cos(np.pi/2*x[dim-1])
+                result = result + deriv
+            deriv = dcoll.zeros(actx) + dim
+            for j in range(dim-1):
+                deriv = deriv * actx.np.sin(np.pi*x[j])
+            deriv = deriv * (-np.pi/2*actx.np.sin(np.pi/2*x[dim-1]))
+            result = result + deriv
+            return result
+
+        x = actx.thaw(dcoll.nodes())
+
+        if vectorize:
+            u = make_obj_array([(i+1)*f(x) for i in range(dim)])
+            if not nested:
+                u = np.stack(u, axis=0)
+        else:
+            u = f(x)
+
+        def get_flux(u_tpair, dcoll=dcoll):
+            dd = u_tpair.dd
+            dd_allfaces = dd.with_domain_tag(
+                BoundaryDomainTag(FACE_RESTR_ALL, VTAG_ALL)
+                )
+            normal = geometry.normal(actx, dcoll, dd)
+            flux = u_tpair.avg @ normal
+            return op.project(dcoll, dd, dd_allfaces, flux)
+
+        dd_allfaces = as_dofdesc(FACE_RESTR_ALL)
+
+        if form == "strong":
+            div_u = (
+                op.local_div(dcoll, u)
+                # No flux terms because u doesn't have inter-el jumps
+                )
+        elif form == "weak":
+            div_u = op.inverse_mass(dcoll,
+                -op.weak_local_div(dcoll, u)
+                +
+                op.face_mass(dcoll,
+                    dd_allfaces,
+                    # Note: no boundary flux terms here because u_ext == u_int == 0
+                    sum(get_flux(utpair)
+                        for utpair in op.interior_trace_pairs(dcoll, u))
+                )
+            )
+        else:
+            raise ValueError("Invalid form argument.")
+
+        if vectorize:
+            expected_div_u = make_obj_array([(i+1)*div_f(x) for i in range(dim)])
+        else:
+            expected_div_u = div_f(x)
+
+        if visualize:
+            from grudge.shortcuts import make_visualizer
+            vis = make_visualizer(dcoll, vis_order=order if dim == 3 else dim+3)
+
+            filename = (f"test_divergence_{form}_{dim}_{order}"
+                f"{'_vec' if vectorize else ''}{'_nested' if nested else ''}.vtu")
+            vis.write_vtk_file(filename, [
+                ("u", u),
+                ("div_u", div_u),
+                ("expected_div_u", expected_div_u),
+                ], overwrite=True)
+
+        rel_linf_err = actx.to_numpy(
+            op.norm(dcoll, div_u - expected_div_u, np.inf)
+            / op.norm(dcoll, expected_div_u, np.inf))
+        eoc_rec.add_data_point(1./n, rel_linf_err)
+
+    print("L^inf error:")
+    print(eoc_rec)
+    assert (eoc_rec.order_estimate() >= order - 0.5
+                or eoc_rec.max_error() < 1e-11)
+
+# }}}
+
+
+# You can test individual routines by typing
+# $ python test_grudge.py 'test_routine()'
+
+if __name__ == "__main__":
+    import sys
+    if len(sys.argv) > 1:
+        exec(sys.argv[1])
+    else:
+        pytest.main([__file__])
+
+# vim: fdm=marker

test/test_reductions.py

+__copyright__ = """
+Copyright (C) 2021 University of Illinois Board of Trustees
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import logging
+from dataclasses import dataclass
+
+import mesh_data
+import numpy as np
+import pytest
+
+from arraycontext import (
+    dataclass_array_container,
+    flatten,
+    pytest_generate_tests_for_array_contexts,
+    with_container_arithmetic,
+)
+from meshmode.dof_array import DOFArray
+from pytools.obj_array import make_obj_array
+
+from grudge import op
+from grudge.array_context import PytestPyOpenCLArrayContextFactory
+from grudge.discretization import make_discretization_collection
+
+
+logger = logging.getLogger(__name__)
+pytest_generate_tests = pytest_generate_tests_for_array_contexts(
+        [PytestPyOpenCLArrayContextFactory])
+
+
+@pytest.mark.parametrize(("mesh_size", "with_initial"), [
+    (4, False),
+    (4, True),
+    (0, False),
+    (0, True)
+])
+def test_nodal_reductions(actx_factory, mesh_size, with_initial):
+    actx = actx_factory()
+
+    builder = mesh_data.BoxMeshBuilder1D()
+
+    mesh = builder.get_mesh(mesh_size)
+    dcoll = make_discretization_collection(actx, mesh, order=4)
+    x = actx.thaw(dcoll.nodes())
+
+    def f(x):
+        return -actx.np.sin(10*x[0])
+
+    def g(x):
+        return actx.np.cos(2*x[0])
+
+    def h(x):
+        return -actx.np.tan(5*x[0])
+
+    fields = make_obj_array([f(x), g(x), h(x)])
+
+    f_ref = actx.to_numpy(flatten(fields[0], actx))
+    g_ref = actx.to_numpy(flatten(fields[1], actx))
+    h_ref = actx.to_numpy(flatten(fields[2], actx))
+    concat_fields = np.concatenate([f_ref, g_ref, h_ref])
+
+    for grudge_op, np_op in [(op.nodal_max, np.max),
+                             (op.nodal_min, np.min),
+                             (op.nodal_sum, np.sum)]:
+        extra_kwargs = {}
+        if with_initial:
+            if grudge_op is op.nodal_max:
+                extra_kwargs["initial"] = -100.
+            elif grudge_op is op.nodal_min:
+                extra_kwargs["initial"] = 100.
+
+        # nodal_min/nodal_max have default initial values, so they behave
+        # differently from numpy in the empty case
+        extra_np_only_kwargs = {}
+        if mesh_size == 0 and not with_initial:
+            if grudge_op is op.nodal_max:
+                extra_np_only_kwargs["initial"] = -np.inf
+            elif grudge_op is op.nodal_min:
+                extra_np_only_kwargs["initial"] = np.inf
+
+        # Componentwise reduction checks
+        assert np.isclose(
+            actx.to_numpy(grudge_op(dcoll, "vol", fields[0], **extra_kwargs)),
+            np_op(f_ref, **extra_kwargs, **extra_np_only_kwargs),
+            rtol=1e-13)
+        assert np.isclose(
+            actx.to_numpy(grudge_op(dcoll, "vol", fields[1], **extra_kwargs)),
+            np_op(g_ref, **extra_kwargs, **extra_np_only_kwargs),
+            rtol=1e-13)
+        assert np.isclose(
+            actx.to_numpy(grudge_op(dcoll, "vol", fields[2], **extra_kwargs)),
+            np_op(h_ref, **extra_kwargs, **extra_np_only_kwargs),
+            rtol=1e-13)
+
+        # Test nodal reductions work on object arrays
+        assert np.isclose(
+            actx.to_numpy(grudge_op(dcoll, "vol", fields, **extra_kwargs)),
+            np_op(concat_fields, **extra_kwargs, **extra_np_only_kwargs),
+            rtol=1e-13)
+
+
+def test_elementwise_reductions(actx_factory):
+    actx = actx_factory()
+
+    builder = mesh_data.BoxMeshBuilder1D()
+
+    nelements = 4
+    mesh = builder.get_mesh(nelements)
+    dcoll = make_discretization_collection(actx, mesh, order=4)
+    x = actx.thaw(dcoll.nodes())
+
+    def f(x):
+        return actx.np.sin(x[0])
+
+    field = f(x)
+    mins = []
+    maxs = []
+    sums = []
+    for grp_f in field:
+        min_res = np.empty(grp_f.shape)
+        max_res = np.empty(grp_f.shape)
+        sum_res = np.empty(grp_f.shape)
+        for eidx in range(mesh.nelements):
+            element_data = actx.to_numpy(grp_f[eidx])
+            min_res[eidx, :] = np.min(element_data)
+            max_res[eidx, :] = np.max(element_data)
+            sum_res[eidx, :] = np.sum(element_data)
+        mins.append(actx.from_numpy(min_res))
+        maxs.append(actx.from_numpy(max_res))
+        sums.append(actx.from_numpy(sum_res))
+
+    ref_mins = DOFArray(actx, data=tuple(mins))
+    ref_maxs = DOFArray(actx, data=tuple(maxs))
+    ref_sums = DOFArray(actx, data=tuple(sums))
+
+    elem_mins = op.elementwise_min(dcoll, field)
+    elem_maxs = op.elementwise_max(dcoll, field)
+    elem_sums = op.elementwise_sum(dcoll, field)
+
+    assert actx.to_numpy(op.norm(dcoll, elem_mins - ref_mins, np.inf)) < 1.e-15
+    assert actx.to_numpy(op.norm(dcoll, elem_maxs - ref_maxs, np.inf)) < 1.e-15
+    assert actx.to_numpy(op.norm(dcoll, elem_sums - ref_sums, np.inf)) < 1.e-15
+
+
+# {{{ Array container tests
+
+@with_container_arithmetic(bcasts_across_obj_array=False,
+                           eq_comparison=False,
+                           rel_comparison=False,
+                           _cls_has_array_context_attr=True,
+                           )
+@dataclass_array_container
+@dataclass(frozen=True)
+class MyContainer:
+    name: str
+    mass: DOFArray
+    momentum: np.ndarray
+    enthalpy: DOFArray
+
+    # NOTE: disable numpy doing any array math
+    __array_ufunc__ = None
+
+    @property
+    def array_context(self):
+        return self.mass.array_context
+
+
+def test_nodal_reductions_with_container(actx_factory):
+    actx = actx_factory()
+
+    builder = mesh_data.BoxMeshBuilder2D()
+
+    mesh = builder.get_mesh(4)
+    dcoll = make_discretization_collection(actx, mesh, order=4)
+    x = actx.thaw(dcoll.nodes())
+
+    def f(x):
+        return -actx.np.sin(10*x[0]) * actx.np.cos(2*x[1])
+
+    def g(x):
+        return actx.np.cos(2*x[0]) * actx.np.sin(10*x[1])
+
+    def h(x):
+        return -actx.np.tan(5*x[0]) * actx.np.tan(0.5*x[1])
+
+    mass = f(x) + g(x)
+    momentum = make_obj_array([f(x)/g(x), h(x)])
+    enthalpy = h(x) - g(x)
+
+    ary_container = MyContainer(name="container",
+                                mass=mass,
+                                momentum=momentum,
+                                enthalpy=enthalpy)
+
+    mass_ref = actx.to_numpy(flatten(mass, actx))
+    momentum_ref = np.concatenate([
+            actx.to_numpy(mom_i)
+            for mom_i in flatten(momentum, actx, leaf_class=DOFArray)])
+    enthalpy_ref = actx.to_numpy(flatten(enthalpy, actx))
+    concat_fields = np.concatenate([mass_ref, momentum_ref, enthalpy_ref])
+
+    for grudge_op, np_op in [(op.nodal_sum, np.sum),
+                             (op.nodal_max, np.max),
+                             (op.nodal_min, np.min)]:
+
+        assert np.isclose(actx.to_numpy(grudge_op(dcoll, "vol", ary_container)),
+                          np_op(concat_fields), rtol=1e-13)
+
+    # Check norm reduction
+    assert np.isclose(actx.to_numpy(op.norm(dcoll, ary_container, np.inf)),
+                      np.linalg.norm(concat_fields, ord=np.inf),
+                      rtol=1e-13)
+
+
+def test_elementwise_reductions_with_container(actx_factory):
+    actx = actx_factory()
+
+    builder = mesh_data.BoxMeshBuilder2D()
+
+    nelements = 4
+    mesh = builder.get_mesh(nelements)
+    dcoll = make_discretization_collection(actx, mesh, order=4)
+    x = actx.thaw(dcoll.nodes())
+
+    def f(x):
+        return actx.np.sin(x[0]) * actx.np.sin(x[1])
+
+    def g(x):
+        return actx.np.cos(x[0]) * actx.np.cos(x[1])
+
+    def h(x):
+        return actx.np.cos(x[0]) * actx.np.sin(x[1])
+
+    mass = 2*f(x) + 0.5*g(x)
+    momentum = make_obj_array([f(x)/g(x), h(x)])
+    enthalpy = 3*h(x) - g(x)
+
+    ary_container = MyContainer(name="container",
+                                mass=mass,
+                                momentum=momentum,
+                                enthalpy=enthalpy)
+
+    def _get_ref_data(field):
+        mins = []
+        maxs = []
+        sums = []
+        for grp_f in field:
+            min_res = np.empty(grp_f.shape)
+            max_res = np.empty(grp_f.shape)
+            sum_res = np.empty(grp_f.shape)
+            for eidx in range(mesh.nelements):
+                element_data = actx.to_numpy(grp_f[eidx])
+                min_res[eidx, :] = np.min(element_data)
+                max_res[eidx, :] = np.max(element_data)
+                sum_res[eidx, :] = np.sum(element_data)
+            mins.append(actx.from_numpy(min_res))
+            maxs.append(actx.from_numpy(max_res))
+            sums.append(actx.from_numpy(sum_res))
+        min_field = DOFArray(actx, data=tuple(mins))
+        max_field = DOFArray(actx, data=tuple(maxs))
+        sums_field = DOFArray(actx, data=tuple(sums))
+        return min_field, max_field, sums_field
+
+    min_mass, max_mass, sums_mass = _get_ref_data(mass)
+    min_enthalpy, max_enthalpy, sums_enthalpy = _get_ref_data(enthalpy)
+    min_mom_x, max_mom_x, sums_mom_x = _get_ref_data(momentum[0])
+    min_mom_y, max_mom_y, sums_mom_y = _get_ref_data(momentum[1])
+    min_momentum = make_obj_array([min_mom_x, min_mom_y])
+    max_momentum = make_obj_array([max_mom_x, max_mom_y])
+    sums_momentum = make_obj_array([sums_mom_x, sums_mom_y])
+
+    reference_min = MyContainer(
+        name="Reference min",
+        mass=min_mass,
+        momentum=min_momentum,
+        enthalpy=min_enthalpy
+    )
+
+    reference_max = MyContainer(
+        name="Reference max",
+        mass=max_mass,
+        momentum=max_momentum,
+        enthalpy=max_enthalpy
+    )
+
+    reference_sum = MyContainer(
+        name="Reference sums",
+        mass=sums_mass,
+        momentum=sums_momentum,
+        enthalpy=sums_enthalpy
+    )
+
+    elem_mins = op.elementwise_min(dcoll, ary_container)
+    elem_maxs = op.elementwise_max(dcoll, ary_container)
+    elem_sums = op.elementwise_sum(dcoll, ary_container)
+
+    assert actx.to_numpy(op.norm(dcoll, elem_mins - reference_min, np.inf)) < 1.e-14
+    assert actx.to_numpy(op.norm(dcoll, elem_maxs - reference_max, np.inf)) < 1.e-14
+    assert actx.to_numpy(op.norm(dcoll, elem_sums - reference_sum, np.inf)) < 1.e-14
+
+# }}}
+
+
+# You can test individual routines by typing
+# $ python test_grudge.py 'test_routine()'
+
+if __name__ == "__main__":
+    import sys
+    if len(sys.argv) > 1:
+        exec(sys.argv[1])
+    else:
+        pytest.main([__file__])

test/test_tools.py

+__copyright__ = """
+Copyright (C) 2022 University of Illinois Board of Trustees
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+from dataclasses import dataclass
+
+import numpy as np
+import numpy.linalg as la  # noqa
+
+from arraycontext import pytest_generate_tests_for_array_contexts
+
+from grudge.array_context import PytestPyOpenCLArrayContextFactory
+
+
+pytest_generate_tests = pytest_generate_tests_for_array_contexts(
+        [PytestPyOpenCLArrayContextFactory])
+
+import logging
+
+import pytest
+
+from pytools.obj_array import make_obj_array
+
+
+logger = logging.getLogger(__name__)
+
+
+# {{{ map_subarrays and rec_map_subarrays
+
+@dataclass(frozen=True, eq=True)
+class _DummyScalar:
+    val: int
+
+
+def test_map_subarrays():
+    """Test map_subarrays."""
+    from grudge.tools import map_subarrays
+
+    # Scalar
+    result = map_subarrays(
+        lambda x: np.array([x, 2*x]), (), (2,), 1)
+    assert result.dtype == int
+    assert np.all(result == np.array([1, 2]))
+
+    # Scalar, nested
+    result = map_subarrays(
+        lambda x: np.array([x, 2*x]), (), (2,), 1, return_nested=True)
+    assert result.dtype == int
+    assert np.all(result == np.array([1, 2]))  # Same as non-nested
+
+    # in_shape is whole array
+    result = map_subarrays(
+        lambda x: 2*x, (2,), (2,), np.array([1, 2]))
+    assert result.dtype == int
+    assert np.all(result == np.array([2, 4]))
+
+    # in_shape is whole array, nested
+    result = map_subarrays(
+        lambda x: 2*x, (2,), (2,), np.array([1, 2]), return_nested=True)
+    assert result.dtype == int
+    assert np.all(result == np.array([2, 4]))  # Same as non-nested
+
+    # len(out_shape) == 0
+    result = map_subarrays(
+        np.sum, (2,), (), np.array([[1, 2], [2, 4]]))
+    assert result.dtype == int
+    assert np.all(result == np.array([3, 6]))
+
+    # len(out_shape) == 0, nested output
+    result = map_subarrays(
+        np.sum, (2,), (), np.array([[1, 2], [2, 4]]), return_nested=True)
+    assert np.all(result == np.array([3, 6]))  # Same as non-nested output
+
+    # len(out_shape) == 0, non-numerical scalars
+    result = map_subarrays(
+        lambda x: _DummyScalar(x[0].val + x[1].val), (2,), (),
+        np.array([
+            [_DummyScalar(1), _DummyScalar(2)],
+            [_DummyScalar(2), _DummyScalar(4)]]))
+    assert result.dtype == object
+    assert result.shape == (2,)
+    assert result[0] == _DummyScalar(3)
+    assert result[1] == _DummyScalar(6)
+
+    # len(out_shape) != 0
+    result = map_subarrays(
+        lambda x: np.array([x, 2*x]), (), (2,), np.array([1, 2]))
+    assert result.dtype == int
+    assert np.all(result == np.array([[1, 2], [2, 4]]))
+
+    # len(out_shape) != 0, nested
+    result = map_subarrays(
+        lambda x: np.array([x, 2*x]), (), (2,), np.array([1, 2]), return_nested=True)
+    assert result.dtype == object
+    assert result.shape == (2,)
+    assert np.all(result[0] == np.array([1, 2]))
+    assert np.all(result[1] == np.array([2, 4]))
+
+    # len(out_shape) != 0, non-numerical scalars
+    result = map_subarrays(
+        lambda x: np.array([_DummyScalar(x), _DummyScalar(2*x)]), (), (2,),
+        np.array([1, 2]))
+    assert result.dtype == object
+    assert result.shape == (2, 2)
+    assert np.all(result[0] == np.array([_DummyScalar(1), _DummyScalar(2)]))
+    assert np.all(result[1] == np.array([_DummyScalar(2), _DummyScalar(4)]))
+
+    # Zero-size input array
+    result = map_subarrays(
+        lambda x: np.array([x, 2*x]), (), (2,), np.empty((2, 0)))
+    assert result.dtype == object
+    assert result.shape == (2, 0, 2)
+
+    # Zero-size input array, nested
+    result = map_subarrays(
+        lambda x: np.array([x, 2*x]), (), (2,), np.empty((2, 0)),
+        return_nested=True)
+    assert result.dtype == object
+    assert result.shape == (2, 0)
+
+
+def test_rec_map_subarrays():
+    """Test rec_map_subarrays."""
+    from grudge.tools import rec_map_subarrays
+
+    # Scalar
+    result = rec_map_subarrays(
+        lambda x: np.array([x, 2*x]), (), (2,), 1)
+    assert result.dtype == int
+    assert np.all(result == np.array([1, 2]))
+
+    # Scalar, non-numerical
+    result = rec_map_subarrays(
+        lambda x: np.array([x.val, 2*x.val]), (), (2,), _DummyScalar(1),
+        scalar_cls=_DummyScalar)
+    assert result.dtype == int
+    assert np.all(result == np.array([1, 2]))
+
+    # Array of scalars
+    result = rec_map_subarrays(
+        np.sum, (2,), (), np.array([[1, 2], [2, 4]]))
+    assert result.dtype == int
+    assert np.all(result == np.array([3, 6]))
+
+    # Array of scalars, non-numerical
+    result = rec_map_subarrays(
+        lambda x: x[0].val + x[1].val, (2,), (),
+        np.array([
+            [_DummyScalar(1), _DummyScalar(2)],
+            [_DummyScalar(2), _DummyScalar(4)]]),
+        scalar_cls=_DummyScalar)
+    assert result.dtype == int
+    assert np.all(result == np.array([3, 6]))
+
+    # Array container
+    result = rec_map_subarrays(
+        np.sum, (2,), (), make_obj_array([np.array([1, 2]), np.array([2, 4])]))
+    assert result.dtype == object
+    assert result[0] == 3
+    assert result[1] == 6
+
+    # Array container, non-numerical scalars
+    result = rec_map_subarrays(
+        lambda x: x[0].val + x[1], (2,), (),
+        make_obj_array([
+            np.array([_DummyScalar(1), 2]),
+            np.array([_DummyScalar(2), 4])]),
+        scalar_cls=_DummyScalar)
+    assert result.dtype == object
+    assert result[0] == 3
+    assert result[1] == 6
+
+# }}}
+
+
+# You can test individual routines by typing
+# $ python test_tools.py 'test_routine()'
+
+if __name__ == "__main__":
+    import sys
+    if len(sys.argv) > 1:
+        exec(sys.argv[1])
+    else:
+        pytest.main([__file__])
+
+# vim: fdm=marker

test/test_trace_pair.py

+__copyright__ = "Copyright (C) 2021 University of Illinois Board of Trustees"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+
+import logging
+
+import numpy as np
+
+import meshmode.mesh.generation as mgen
+from arraycontext import pytest_generate_tests_for_array_contexts
+from meshmode.dof_array import DOFArray
+
+from grudge.array_context import PytestPyOpenCLArrayContextFactory
+from grudge.discretization import make_discretization_collection
+from grudge.trace_pair import TracePair
+
+
+logger = logging.getLogger(__name__)
+pytest_generate_tests = pytest_generate_tests_for_array_contexts(
+        [PytestPyOpenCLArrayContextFactory])
+
+
+def test_trace_pair(actx_factory):
+    """Simple smoke test for :class:`grudge.trace_pair.TracePair`."""
+    actx = actx_factory()
+    dim = 3
+    order = 1
+    n = 4
+
+    mesh = mgen.generate_regular_rect_mesh(
+        a=(-1,)*dim, b=(1,)*dim,
+        nelements_per_axis=(n,)*dim)
+
+    dcoll = make_discretization_collection(actx, mesh, order=order)
+
+    rng = np.random.default_rng(1234)
+
+    def rand():
+        return DOFArray(
+                actx,
+                tuple(actx.from_numpy(
+                    rng.uniform(size=(grp.nelements, grp.nunit_dofs)))
+                    for grp in dcoll.discr_from_dd("vol").groups))
+
+    from grudge.dof_desc import DD_VOLUME_ALL
+    interior = rand()
+    exterior = rand()
+    tpair = TracePair(DD_VOLUME_ALL, interior=interior, exterior=exterior)
+
+    import grudge.op as op
+    assert op.norm(dcoll, tpair.avg - 0.5*(exterior + interior), np.inf) == 0
+    assert op.norm(dcoll, tpair.diff - (exterior - interior), np.inf) == 0
+    assert op.norm(dcoll, tpair.int - interior, np.inf) == 0
+    assert op.norm(dcoll, tpair.ext - exterior, np.inf) == 0