test_grudge.py

# }}}


# {{{ function symbol

def test_external_call(actx_factory):
    actx = actx_factory()

    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(actx, 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(actx, double=double)
    assert actx.to_numpy(flatten(result) == 5).all()


@pytest.mark.parametrize("array_type", ["scalar", "vector"])
def test_function_symbol_array(actx_factory, array_type):
    """Test if `FunctionSymbol` distributed properly over object arrays."""

    actx = actx_factory()

    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(actx, mesh, order=4)
    volume_discr = discr.discr_from_dd(sym.DD_VOLUME)

    if array_type == "scalar":
        sym_x = sym.var("x")
        x = thaw(actx, actx.np.cos(volume_discr.nodes()[0]))
    elif array_type == "vector":
        sym_x = sym.make_sym_array("x", dim)
        x = thaw(actx, volume_discr.nodes())
    else:
        raise ValueError("unknown array type")

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

# }}}


@pytest.mark.parametrize("p", [2, np.inf])
def test_norm_obj_array(actx_factory, p):
    """Test :func:`grudge.symbolic.operators.norm` for object arrays."""

    actx = actx_factory()

    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=1)
    discr = DGDiscretizationWithBoundaries(actx, mesh, order=4)

    w = make_obj_array([1.0, 2.0, 3.0])[:dim]

    # {{ scalar

    sym_w = sym.var("w")
    norm = bind(discr, sym.norm(p, sym_w))(actx, w=w[0])

    norm_exact = w[0]
    logger.info("norm: %.5e %.5e", norm, norm_exact)
    assert abs(norm - norm_exact) < 1.0e-14

    # }}}

    # {{{ vector

    sym_w = sym.make_sym_array("w", dim)
    norm = bind(discr, sym.norm(p, sym_w))(actx, w=w)

    norm_exact = np.sqrt(np.sum(w**2)) if p == 2 else np.max(w)
    logger.info("norm: %.5e %.5e", norm, norm_exact)
    assert abs(norm - norm_exact) < 1.0e-14

    # }}}


def test_map_if(actx_factory):
    """Test :meth:`grudge.symbolic.execution.ExecutionMapper.map_if` handling
    of scalar conditions.
    """

    actx = actx_factory()

    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(actx, mesh, order=4)

    sym_if = sym.If(sym.Comparison(2.0, "<", 1.0e-14), 1.0, 2.0)
    bind(discr, sym_if)(actx)


# 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