test_grudge.py

        else:
            qtag = None

        ns = [20, 25]
        for n in ns:
            mesh = mgen.generate_regular_rect_mesh(
                a=(-0.5,)*dims,
                b=(0.5,)*dims,
                nelements_per_axis=(n,)*dims,
                order=order)

            if use_quad:
                discr_tag_to_group_factory = {
                    qtag: QuadratureSimplexGroupFactory(order=4*order)
                }
            else:
                discr_tag_to_group_factory = {}

            dcoll = make_discretization_collection(
                actx, mesh, order=order,
                discr_tag_to_group_factory=discr_tag_to_group_factory
            )

            nodes = actx.thaw(dcoll.nodes())

            def zero_inflow(dtag, t=0, dcoll=dcoll):
                dd = dof_desc.as_dofdesc(dtag, qtag)
                return dcoll.discr_from_dd(dd).zeros(actx)

            adv_op = VariableCoefficientAdvectionOperator(
                dcoll,
                flat_obj_array(-1*nodes[1], nodes[0]),
                inflow_u=lambda t: zero_inflow(BTAG_ALL, t=t),
                flux_type="upwind",
                quad_tag=qtag
            )

            total_error = op.norm(
                dcoll, adv_op.operator(0, gaussian_mode(nodes)), 2
            )
            eoc_rec.add_data_point(1.0/n, actx.to_numpy(total_error))

        logger.info("\n%s", eoc_rec.pretty_print(
            abscissa_label="h",
            error_label="L2 Error"))

        return eoc_rec.order_estimate(), np.array([x[1] for x in eoc_rec.history])

    eoc, errs = conv_test("no quadrature", False)
    q_eoc, q_errs = conv_test("with quadrature", True)

    assert q_eoc > eoc
    assert (q_errs < errs).all()
    assert q_eoc > order - 0.1

# }}}


# {{{ bessel

@pytest.mark.xfail
def test_bessel(actx_factory):
    actx = actx_factory()

    dims = 2

    mesh = mgen.generate_regular_rect_mesh(
            a=(0.1,)*dims,
            b=(1.0,)*dims,
            nelements_per_axis=(8,)*dims)

    dcoll = make_discretization_collection(actx, mesh, order=3)

    nodes = actx.thaw(dcoll.nodes())
    r = actx.np.sqrt(nodes[0]**2 + nodes[1]**2)

    # FIXME: Bessel functions need to brought out of the symbolic
    # layer. Related issue: https://github.com/inducer/grudge/issues/93
    def bessel_j(actx, n, r):
        from grudge import bind, sym  # pylint: disable=no-name-in-module
        return bind(dcoll, sym.bessel_j(n, sym.var("r")))(actx, r=r)

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

    z = op.norm(dcoll, bessel_zero, 2)

    assert z < 1e-15

# }}}


# {{{ test norms

@pytest.mark.parametrize("p", [2, np.inf])
def test_norm_real(actx_factory, p):
    actx = actx_factory()

    dim = 2
    mesh = mgen.generate_regular_rect_mesh(
            a=(0,)*dim, b=(1,)*dim,
            nelements_per_axis=(8,)*dim, order=1)
    dcoll = make_discretization_collection(actx, mesh, order=4)
    nodes = actx.thaw(dcoll.nodes())

    norm = op.norm(dcoll, nodes[0], p)
    if p == 2:
        ref_norm = (1/3)**0.5
    elif p == np.inf:
        ref_norm = 1
    else:
        raise AssertionError("unsupported p")

    logger.info("norm: %.5e %.5e", norm, ref_norm)
    assert abs(norm-ref_norm) / abs(ref_norm) < 1e-13


@pytest.mark.parametrize("p", [2, np.inf])
def test_norm_complex(actx_factory, p):
    actx = actx_factory()

    dim = 2
    mesh = mgen.generate_regular_rect_mesh(
            a=(0,)*dim, b=(1,)*dim,
            nelements_per_axis=(8,)*dim, order=1)
    dcoll = make_discretization_collection(actx, mesh, order=4)
    nodes = actx.thaw(dcoll.nodes())

    norm = op.norm(dcoll, (1 + 1j)*nodes[0], p)
    if p == 2:
        ref_norm = (2/3)**0.5
    elif p == np.inf:
        ref_norm = 2**0.5
    else:
        raise AssertionError("unsupported p")

    logger.info("norm: %.5e %.5e", norm, ref_norm)
    assert abs(norm-ref_norm) / abs(ref_norm) < 1e-13


@pytest.mark.parametrize("p", [2, np.inf])
def test_norm_obj_array(actx_factory, p):
    actx = actx_factory()

    dim = 2
    mesh = mgen.generate_regular_rect_mesh(
            a=(0,)*dim, b=(1,)*dim,
            nelements_per_axis=(8,)*dim, order=1)
    dcoll = make_discretization_collection(actx, mesh, order=4)
    nodes = actx.thaw(dcoll.nodes())

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

    if p == 2:
        ref_norm = (dim/3)**0.5
    elif p == np.inf:
        ref_norm = 1
    else:
        raise AssertionError("unsupported p")

    logger.info("norm: %.5e %.5e", norm, ref_norm)
    assert abs(norm-ref_norm) / abs(ref_norm) < 1e-14

# }}}


# {{{ empty boundaries

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

# }}}


# {{{ multi-volume

def test_multiple_independent_volumes(actx_factory):
    dim = 2
    actx = actx_factory()

    mesh1 = mgen.generate_regular_rect_mesh(
            a=(-2,)*dim, b=(-1,)*dim,
            nelements_per_axis=(4,)*dim, order=4)

    mesh2 = mgen.generate_regular_rect_mesh(
            a=(1,)*dim, b=(2,)*dim,
            nelements_per_axis=(8,)*dim, order=4)

    volume_to_mesh = {
        "vol1": mesh1,
        "vol2": mesh2}

    make_discretization_collection(actx, volume_to_mesh, order=4)

# }}}


# 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