test_grudge.py

    from grudge.discretization import PointsDiscretization
    pdiscr = PointsDiscretization(cl.array.to_device(queue, cp.points))

    bind(pdiscr, sym.nodes(dim)**2)(queue)

# }}}


# {{{ operator collector determinism

def test_op_collector_order_determinism():
    class TestOperator(sym.Operator):

        def __init__(self):
            sym.Operator.__init__(self, sym.DD_VOLUME, sym.DD_VOLUME)

        mapper_method = "map_test_operator"

    from grudge.symbolic.mappers import BoundOperatorCollector

    class TestBoundOperatorCollector(BoundOperatorCollector):

        def map_test_operator(self, expr):
            return self.map_operator(expr)

    v0 = sym.var("v0")
    ob0 = sym.OperatorBinding(TestOperator(), v0)

    v1 = sym.var("v1")
    ob1 = sym.OperatorBinding(TestOperator(), v1)

    # The output order isn't significant, but it should always be the same.
    assert list(TestBoundOperatorCollector(TestOperator)(ob0 + ob1)) == [ob0, ob1]

# }}}


# {{{ bessel

def test_bessel(ctx_factory):
    cl_ctx = ctx_factory()
    queue = cl.CommandQueue(cl_ctx)

    dims = 2

    from meshmode.mesh.generation import generate_regular_rect_mesh
    mesh = generate_regular_rect_mesh(
            a=(0.1,)*dims,
            b=(1.0,)*dims,
            n=(8,)*dims)

    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=3)

    nodes = sym.nodes(dims)
    r = sym.cse(sym.sqrt(nodes[0]**2 + nodes[1]**2))

    # https://dlmf.nist.gov/10.6.1
    n = 3
    bessel_zero = (
            sym.bessel_j(n+1, r)
            + sym.bessel_j(n-1, r)
            - 2*n/r * sym.bessel_j(n, r))

    z = bind(discr, sym.norm(2, bessel_zero))(queue)

    assert z < 1e-15

# }}}


# {{{ function symbol

def test_external_call(ctx_factory):
    cl_ctx = ctx_factory()
    queue = cl.CommandQueue(cl_ctx)

    def double(queue, x):
        return 2 * x

    from meshmode.mesh.generation import generate_regular_rect_mesh

    dims = 2

    mesh = generate_regular_rect_mesh(a=(0,) * dims, b=(1,) * dims, n=(4,) * dims)
    discr = DGDiscretizationWithBoundaries(cl_ctx, mesh, order=1)

    ones = sym.Ones(sym.DD_VOLUME)
    op = (
            ones * 3
            + sym.FunctionSymbol("double")(ones))

    from grudge.function_registry import (
            base_function_registry, register_external_function)

    freg = register_external_function(
            base_function_registry,
            "double",
            implementation=double,
            dd=sym.DD_VOLUME)

    bound_op = bind(discr, op, function_registry=freg)

    result = bound_op(queue, double=double)
    assert (result == 5).get().all()


@pytest.mark.parametrize("array_type", ["scalar", "vector"])
def test_function_symbol_array(ctx_factory, array_type):
    ctx = ctx_factory()
    queue = cl.CommandQueue(ctx)

    from meshmode.mesh.generation import generate_regular_rect_mesh
    dim = 2
    mesh = generate_regular_rect_mesh(
            a=(-0.5,)*dim, b=(0.5,)*dim,
            n=(8,)*dim, order=4)
    discr = DGDiscretizationWithBoundaries(ctx, mesh, order=4)
    nnodes = discr.discr_from_dd(sym.DD_VOLUME).nnodes

    import pyopencl.clrandom        # noqa: F401
    if array_type == "scalar":
        sym_x = sym.var("x")
        x = cl.clrandom.rand(queue, nnodes, dtype=np.float)
    elif array_type == "vector":
        sym_x = sym.make_sym_array("x", dim)
        x = make_obj_array([
            cl.clrandom.rand(queue, nnodes, dtype=np.float)
            for _ in range(dim)
            ])
    else:
        raise ValueError("unknown array type")

    norm = bind(discr, sym.norm(2, sym_x))(queue, x=x)
    assert isinstance(norm, float)

# }}}


# 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