test.py

import numpy as np
import numpy.linalg as la  # noqa: F401
import pyopencl as cl  # noqa: F401
import pyopencl.array  # noqa
import pyopencl.tools  # noqa
import pyopencl.clrandom  # noqa
import loopy as lp  # noqa

import sys
import logging

import pytest
from pyopencl.tools import (  # noqa
        pytest_generate_tests_for_pyopencl
        as pytest_generate_tests)

import utilities as u


@pytest.mark.parametrize("lam_pointwise_str,lam_roe_str,lam_expected_str", [
    ("1 2 3 4 5,2 4 6 8 10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
    ("1 2 3 4 5,-2 -4 -6 -8 -10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"),
    ("1 2 3 4 5,-2 -4 -6 -8 -10", "3 6 9 12 15", "3.3 6.6 9.9 13.2 16.5"),
    ("1 2 3 4 5,2 4 6 8 10", "-3 -6 -9 -12 -15", "3.3 6.6 9.9 13.2 16.5"),
    ("3 2 9 4 5,2 6 6 12 10", "-1 -4 -3 -8 -15", "3.3 6.6 9.9 13.2 16.5")
    ])
def test_lax_wavespeeds(
        ctx_factory, lam_pointwise_str, lam_roe_str, lam_expected_str):
    prg = u.get_weno_program_with_root_kernel("lax_wavespeeds")
    queue = u.get_queue(ctx_factory)

    nvars = 5

    lam_pointwise = u.expand_to_6(u.transposed_array_from_string(lam_pointwise_str))
    lam_roe = u.array_from_string(lam_roe_str)
    lam_dev = u.empty_array_on_device(queue, nvars)

    prg(queue, nvars=nvars, lambda_pointwise=lam_pointwise,
            lambda_roe=lam_roe, wavespeeds=lam_dev)

    lam_expected = u.array_from_string(lam_expected_str)
    u.compare_arrays(lam_dev.get(), lam_expected)


@pytest.mark.parametrize("states_str,direction,lam_expected_str", [
    ("2 4 4 4 20,1 1 1 1 5.5", "x",
        "2 2 2 3.49666295 0.503337045,1 1 1 2.49666295 -0.496662955"),
    ("2 4 4 4 20,1 1 1 1 5.5", "y",
        "2 2 2 3.49666295 0.503337045,1 1 1 2.49666295 -0.496662955"),
    ("2 4 4 4 20,1 1 1 1 5.5", "z",
        "2 2 2 3.49666295 0.503337045,1 1 1 2.49666295 -0.496662955"),
    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "x",
        "-1 -1 -1 0.496662955 -2.49666295,-2 -2 -2 -0.503337045 -3.49666295"),
    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "y",
        "-1 -1 -1 0.496662955 -2.49666295,-2 -2 -2 -0.503337045 -3.49666295"),
    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "z",
        "-1 -1 -1 0.496662955 -2.49666295,-2 -2 -2 -0.503337045 -3.49666295"),
    ("2 4 8 12 64,1 1 2 3 11", "x",
        "2 2 2 3.49666295 0.503337045,1 1 1 2.49666295 -0.496662955"),
    ("2 4 8 12 64,1 1 2 3 11", "y",
        "4 4 4 5.49666295 2.50333705,2 2 2 3.49666295 0.503337045"),
    ("2 4 8 12 64,1 1 2 3 11", "z",
        "6 6 6 7.49666295 4.50333705,3 3 3 4.49666295 1.503337045"),
    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "x",
        "-1 -1 -1 0.496662955 -2.49666295,-2 -2 -2 -0.503337045 -3.49666295"),
    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "y",
        "-2 -2 -2 -0.503337045 -3.49666295,-4 -4 -4 -2.50333705 -5.49666295"),
    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "z",
        "-3 -3 -3 -1.50333705 -4.49666295,-6 -6 -6 -4.50333705 -7.49666295")
    ])
def test_pointwise_eigenvalues_ideal_gas(
        ctx_factory, states_str, direction, lam_expected_str):
    prg = u.get_weno_program_with_root_kernel("pointwise_eigenvalues")
    queue = u.get_queue(ctx_factory)

    nvars = 5
    dirs = {"x": 1, "y": 2, "z": 3}

    states = u.expand_to_6(u.transposed_array_from_string(states_str))
    lam_dev = u.empty_array_on_device(queue, nvars, 6)

    prg(queue, nvars=nvars, d=dirs[direction],
            states=states, lambda_pointwise=lam_dev)

    lam_expected = u.expand_to_6(u.transposed_array_from_string(lam_expected_str))
    u.compare_arrays(lam_dev.get(), lam_expected)


@pytest.mark.parametrize("states_str,fluxes_str,direction", [
    ("2 4 4 4 20,1 1 1 1 5.5", "4 11.2 8 8 46.4,1 2.6 1 1 7.1", "x"),
    ("2 4 4 4 20,1 1 1 1 5.5", "4 8 11.2 8 46.4,1 1 2.6 1 7.1", "y"),
    ("2 4 4 4 20,1 1 1 1 5.5", "4 8 8 11.2 46.4,1 1 1 2.6 7.1", "z"),
    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "-1 2.6 1 1 -7.1,-4 11.2 8 8 -46.4", "x"),
    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "-1 1 2.6 1 -7.1,-4 8 11.2 8 -46.4", "y"),
    ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "-1 1 1 2.6 -7.1,-4 8 8 11.2 -46.4", "z"),
    ("2 4 8 12 64,1 1 2 3 11", "4 11.2 16 24 134.4,1 2.6 2 3 12.6", "x"),
    ("2 4 8 12 64,1 1 2 3 11", "8 16 35.2 48 268.8,2 2 5.6 6 25.2", "y"),
    ("2 4 8 12 64,1 1 2 3 11", "12 24 48 75.2 403.2,3 3 6 10.6 37.8", "z"),
    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "-1 2.6 2 3 -12.6,-4 11.2 16 24 -134.4", "x"),
    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "-2 2 5.6 6 -25.2,-8 16 35.2 48 -268.8", "y"),
    ("1 -1 -2 -3 11,2 -4 -8 -12 64", "-3 3 6 10.6 -37.8,-12 24 48 75.2 -403.2", "z")
    ])
def test_roe_uniform_grid_ideal_gas(ctx_factory, states_str, fluxes_str, direction):
    def identity_matrix(n):
        return np.identity(n).astype(np.float32).copy(order="F")

    def check_roe_identity(states, r, r_inv):
        d_state = states[:, 1] - states[:, 0]
        u.compare_arrays(r@(r_inv@d_state), d_state)

    def check_roe_property(states, fluxes, r, r_inv, lam):
        d_state = states[:, 1] - states[:, 0]
        d_flux = fluxes[:, 1] - fluxes[:, 0]

        temp = r_inv@d_state
        temp = np.multiply(lam, temp)
        u.compare_arrays(r@temp, d_flux)

    prg = u.get_weno_program_with_root_kernel("roe_eigensystem")
    queue = u.get_queue(ctx_factory)

    nvars = 5
    ndim = 3
    dirs = {"x": 1, "y": 2, "z": 3}

    states = u.transposed_array_from_string(states_str)
    fluxes = u.transposed_array_from_string(fluxes_str)
    metrics_frozen = identity_matrix(ndim)

    r_dev = u.empty_array_on_device(queue, nvars, nvars)
    r_inv_dev = u.empty_array_on_device(queue, nvars, nvars)
    lam_dev = u.empty_array_on_device(queue, nvars)

    prg(queue, nvars=nvars, ndim=ndim, d=dirs[direction],
            states=states, metrics_frozen=metrics_frozen,
            R=r_dev, R_inv=r_inv_dev, lambda_roe=lam_dev)

    r = r_dev.get()
    r_inv = r_inv_dev.get()
    lam = lam_dev.get()

    print(lam)

    check_roe_identity(states, r, r_inv)
    check_roe_property(states, fluxes, r, r_inv, lam)


def test_matvec(ctx_factory):
    prg = u.get_weno_program_with_root_kernel("mult_mat_vec")
    queue = u.get_queue(ctx_factory)

    a = u.random_array_on_device(queue, 10, 10)
    b = u.random_array_on_device(queue, 10)

    c = u.empty_array_on_device(queue, 10)

    prg(queue, alpha=1.0, a=a, b=b, c=c)

    u.compare_arrays(a.get()@b.get(), c.get())


@pytest.mark.slow
def test_compute_flux_derivatives(ctx_factory):
    prg = u.get_weno_program()

    queue = u.get_queue(ctx_factory)
    prg = prg.copy(target=lp.PyOpenCLTarget(queue.device))

    lp.auto_test_vs_ref(prg, ctx_factory(), warmup_rounds=1,
            parameters=dict(ndim=3, nvars=5, nx=16, ny=16, nz=16))


@pytest.mark.slow
def test_compute_flux_derivatives_gpu(ctx_factory, write_code=False):
    prg = u.get_weno_program()
    prg = u.transform_weno_for_gpu(prg)

    queue = u.get_queue(ctx_factory)
    prg = prg.copy(target=lp.PyOpenCLTarget(queue.device))
    prg = lp.set_options(prg, no_numpy=True)

    if write_code:
        u.write_target_device_code(prg)

    lp.auto_test_vs_ref(prg, ctx_factory(), warmup_rounds=1,
            parameters=dict(ndim=3, nvars=5, nx=16, ny=16, nz=16))


# This lets you run 'python test.py test_case(cl._csc)' without pytest.
if __name__ == "__main__":
    if len(sys.argv) > 1:
        logging.basicConfig(level="INFO")
        exec(sys.argv[1])
    else:
        pytest.main([__file__])