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

examples/maxwell/testing.py → test/test_trace_pair.py

-"""Minimal example of a grudge driver."""
-
-from __future__ import division, print_function
-
-__copyright__ = "Copyright (C) 2015 Andreas Kloeckner"
+__copyright__ = "Copyright (C) 2021 University of Illinois Board of Trustees"
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -25,61 +21,53 @@ THE SOFTWARE.
 """
 
 
-import numpy as np
-import pyopencl as cl
-import sumpy.point_calculus as spy
-from grudge import bind
+import logging
 
+import numpy as np
 
-def sumpy_test_3D(order, cl_ctx, queue):
-
-    epsilon = 1
-    mu = 1
+import meshmode.mesh.generation as mgen
+from arraycontext import pytest_generate_tests_for_array_contexts
+from meshmode.dof_array import DOFArray
 
-    kx, ky, kz = factors = [n*np.pi for n in (1, 2, 2)]
-    k = np.sqrt(sum(f**2 for f in factors))
+from grudge.array_context import PytestPyOpenCLArrayContextFactory
+from grudge.discretization import make_discretization_collection
+from grudge.trace_pair import TracePair
 
-    patch = spy.CalculusPatch((0.4, 0.3, 0.4))
 
-    from grudge.discretization import PointsDiscretization
-    pdiscr = PointsDiscretization(cl.array.to_device(queue, patch.points))
+logger = logging.getLogger(__name__)
+pytest_generate_tests = pytest_generate_tests_for_array_contexts(
+        [PytestPyOpenCLArrayContextFactory])
 
-    from grudge.models.em import get_rectangular_cavity_mode
-    sym_mode = get_rectangular_cavity_mode(1, (1, 2, 2))
-    fields = bind(pdiscr, sym_mode)(queue, t=0, epsilon=epsilon, mu=mu)
 
-    for i in range(len(fields)):
-        if isinstance(fields[i], (int, float)):
-            fields[i] = np.zeros(patch.points.shape[1])
-        else:
-            fields[i] = fields[i].get()
+def test_trace_pair(actx_factory):
+    """Simple smoke test for :class:`grudge.trace_pair.TracePair`."""
+    actx = actx_factory()
+    dim = 3
+    order = 1
+    n = 4
 
-    e = np.array([fields[0], fields[1], fields[2]])
-    h = np.array([fields[3], fields[4], fields[5]])
-    print(frequency_domain_maxwell(patch, e, h, k))
-    return
+    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 frequency_domain_maxwell(cpatch, e, h, k):
-    mu = 1
-    epsilon = 1
-    c = 1/np.sqrt(mu*epsilon)
-    omega = k*c
-    b = mu*h
-    d = epsilon*e
-    # https://en.wikipedia.org/w/index.php?title=Maxwell%27s_equations&oldid=798940325#Macroscopic_formulation
-    # assumed time dependence exp(-1j*omega*t)
-    # Agrees with Jackson, Third Ed., (8.16)
-    resid_faraday = np.vstack(cpatch.curl(e)) - 1j * omega * b
-    resid_ampere = np.vstack(cpatch.curl(h)) + 1j * omega * d
-    resid_div_e = cpatch.div(e)
-    resid_div_h = cpatch.div(h)
-    return (resid_faraday, resid_ampere,  resid_div_e,    resid_div_h)
+    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))
 
-if __name__ == "__main__":
-    cl_ctx = cl.create_some_context()
-    queue = cl.CommandQueue(cl_ctx)
+    from grudge.dof_desc import DD_VOLUME_ALL
+    interior = rand()
+    exterior = rand()
+    tpair = TracePair(DD_VOLUME_ALL, interior=interior, exterior=exterior)
 
-    print("Doing Sumpy Test")
-    sumpy_test_3D(4, cl_ctx, queue)
+    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