test_grudge.py

# }}}


# {{{ functions

def test_external_call(actx_factory):
    actx = actx_factory()

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

    from grudge.function_registry import \
        base_function_registry, register_external_function

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

    dims = 2

    mesh = mgen.generate_regular_rect_mesh(
            a=(0,) * dims, b=(1,) * dims, nelements_per_axis=(4,) * dims)
    dcoll = DiscretizationCollection(actx, mesh, order=1)

    ones = dcoll.discr_from_dd(dof_desc.DD_VOLUME).zeros(actx) + 1.0

    result = 3*ones + freg["ext_double_fct"](actx, ones)

    assert actx.to_numpy(flatten(result) == 5).all()


@pytest.mark.parametrize("array_type", ["scalar", "vector"])
def test_function_array(actx_factory, array_type):
    """Test if functions distribute properly over object arrays."""

    actx = actx_factory()

    dim = 2
    mesh = mgen.generate_regular_rect_mesh(
            a=(-0.5,)*dim, b=(0.5,)*dim,
            nelements_per_axis=(8,)*dim, order=4)
    dcoll = DiscretizationCollection(actx, mesh, order=4)
    nodes = op.nodes(dcoll)

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

    norm = op.norm(dcoll, x, 2)
    assert isinstance(norm, float)

# }}}


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

    actx = actx_factory()

    dim = 2
    mesh = mgen.generate_regular_rect_mesh(
            a=(-0.5,)*dim, b=(0.5,)*dim,
            nelements_per_axis=(8,)*dim, order=1)
    dcoll = DiscretizationCollection(actx, mesh, order=4)

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

    # {{ scalar

    norm = op.norm(dcoll, w[0], p)

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

    # }}}

    # {{{ vector

    norm = op.norm(dcoll, w, p)

    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_empty_boundary(actx_factory):
    # https://github.com/inducer/grudge/issues/54

    from meshmode.mesh import BTAG_NONE

    actx = actx_factory()

    dim = 2
    mesh = mgen.generate_regular_rect_mesh(
            a=(-0.5,)*dim, b=(0.5,)*dim,
            nelements_per_axis=(8,)*dim, order=4)
    dcoll = DiscretizationCollection(actx, mesh, order=4)
    normal = op.normal(dcoll, BTAG_NONE)
    from meshmode.dof_array import DOFArray
    for component in normal:
        assert isinstance(component, DOFArray)
        assert len(component) == len(dcoll.discr_from_dd(BTAG_NONE).groups)


# {{{ Deprecated symbolic tests

def test_op_collector_order_determinism():
    from grudge import sym

    class TestOperator(sym.Operator):

        def __init__(self):
            sym.Operator.__init__(self, dof_desc.DD_VOLUME, dof_desc.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]


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

    actx = actx_factory()

    dim = 2
    mesh = mgen.generate_regular_rect_mesh(
            a=(-0.5,)*dim, b=(0.5,)*dim,
            nelements_per_axis=(8,)*dim, order=4)
    dcoll = DiscretizationCollection(actx, mesh, order=4)

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


def test_operator_compiler_overwrite(actx_factory):
    """Tests that the same expression in ``eval_code`` and ``discr_code``
    does not confuse the OperatorCompiler in grudge/symbolic/compiler.py.
    """
    from grudge import sym, bind

    actx = actx_factory()

    ambient_dim = 2
    target_order = 4

    from meshmode.mesh.generation import generate_regular_rect_mesh
    mesh = generate_regular_rect_mesh(
            a=(-0.5,)*ambient_dim, b=(0.5,)*ambient_dim,
            n=(8,)*ambient_dim, order=1)
    dcoll = DiscretizationCollection(actx, mesh, order=target_order)

    # {{{ test

    sym_u = sym.nodes(ambient_dim)
    sym_div_u = sum(d(u) for d, u in zip(sym.nabla(ambient_dim), sym_u))

    div_u = bind(dcoll, sym_div_u)(actx)
    error = bind(dcoll, sym.norm(2, sym.var("x")))(actx, x=div_u - dcoll.dim)
    logger.info("error: %.5e", error)

    # }}}


@pytest.mark.parametrize("ambient_dim", [
    2,
    # FIXME, cf. https://github.com/inducer/grudge/pull/78/
    pytest.param(3, marks=pytest.mark.xfail)
    ])
def test_incorrect_assignment_aggregation(actx_factory, ambient_dim):
    """Tests that the greedy assignemnt aggregation code works on a non-trivial
    expression (on which it didn't work at the time of writing).
    """
    from grudge import sym, bind

    actx = actx_factory()

    target_order = 4

    from meshmode.mesh.generation import generate_regular_rect_mesh
    mesh = generate_regular_rect_mesh(
            a=(-0.5,)*ambient_dim, b=(0.5,)*ambient_dim,
            n=(8,)*ambient_dim, order=1)
    dcoll = DiscretizationCollection(actx, mesh, order=target_order)

    # {{{ test with a relative norm

    from grudge.dof_desc import DD_VOLUME
    dd = DD_VOLUME
    sym_x = sym.make_sym_array("y", ambient_dim, dd=dd)
    sym_y = sym.make_sym_array("y", ambient_dim, dd=dd)

    sym_norm_y = sym.norm(2, sym_y, dd=dd)
    sym_norm_d = sym.norm(2, sym_x - sym_y, dd=dd)
    sym_op = sym_norm_d / sym_norm_y
    logger.info("%s", sym.pretty(sym_op))

    # FIXME: this shouldn't raise a RuntimeError
    with pytest.raises(RuntimeError):
        bind(dcoll, sym_op)(actx, x=1.0, y=dcoll.discr_from_dd(dd).nodes())

    # }}}

    # {{{ test with repeated mass inverses

    sym_minv_y = sym.cse(sym.InverseMassOperator()(sym_y), "minv_y")

    sym_u = make_obj_array([0.5 * sym.Ones(dd), 0.0, 0.0])[:ambient_dim]
    sym_div_u = sum(d(u) for d, u in zip(sym.nabla(ambient_dim), sym_u))

    sym_op = sym.MassOperator(dd)(sym_u) \
            + sym.MassOperator(dd)(sym_minv_y * sym_div_u)
    logger.info("%s", sym.pretty(sym_op))

    # FIXME: this shouldn't raise a RuntimeError either
    bind(dcoll, sym_op)(actx, y=dcoll.discr_from_dd(dd).nodes())

    # }}}

# }}}


# 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