diff --git a/WENO.F90 b/WENO.F90 index e8618ef65764d2f64d310a023855a0651cfe8652..995d3b7d86fcb177684c4e97071db6d5b6b249ec 100644 --- a/WENO.F90 +++ b/WENO.F90 @@ -350,7 +350,7 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam c = sqrt((1.4 - 1.0)*(H - 0.5*q)) b1 = (1.4 - 1.0)/(c**2) - b2 = 1.0 + b1*q**2 - b1*H + b2 = 1.0 + b1*q - b1*H u_tilde(1) = 0.0 do i=1,ndim @@ -370,11 +370,11 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam alpha = rho/(2.0*c) beta = 1.0/(2.0*alpha) - lambda_roe(1) = u(1) - lambda_roe(2) = u(1) - lambda_roe(3) = u(1) - lambda_roe(4) = u(1) + c - lambda_roe(5) = u(1) - c + lambda_roe(1) = u_tilde(1)*metric_norm(ik) + lambda_roe(2) = u_tilde(1)*metric_norm(ik) + lambda_roe(3) = u_tilde(1)*metric_norm(ik) + lambda_roe(4) = u_tilde(1)*metric_norm(ik) + c*metric_norm(ik) + lambda_roe(5) = u_tilde(1)*metric_norm(ik) - c*metric_norm(ik) R(1,1) = 1.0 R(2,1) = u(1) diff --git a/comparison_fixtures.py b/comparison_fixtures.py index e21ec0c2b9d40bcc8c45dca96a9b2eddf5c9e193..04c743236fc597b1ec48ed24d5974360e6da6c7e 100644 --- a/comparison_fixtures.py +++ b/comparison_fixtures.py @@ -1,5 +1,20 @@ +import numpy as np from pytest import approx def arrays(a, b): assert a == approx(b) + + +def roe_identity(states, R, Rinv): + dState = states[:,1] - states[:,0] + arrays(R@(Rinv@dState), dState) + + +def roe_property(states, fluxes, R, Rinv, lam): + dState = states[:,1] - states[:,0] + dFlux = fluxes[:,1] - fluxes[:,0] + + temp = Rinv@dState + temp = np.multiply(lam, temp) + arrays(R@temp, dFlux) diff --git a/kernel_fixtures.py b/kernel_fixtures.py index 1ff71ac5ab413959f8f27e4a1509e1638a2a66e6..7f3dff465c30ab8fc9101ab092ef4c134bcfa951 100644 --- a/kernel_fixtures.py +++ b/kernel_fixtures.py @@ -15,6 +15,19 @@ def with_root_kernel(prg, root_name): return new_prg +def roe_eigensystem(queue, prg, params, states, metrics_frozen): + R_dev = setup.empty_array_on_device(queue, params.mat_bounds()) + Rinv_dev = setup.empty_array_on_device(queue, params.mat_bounds()) + lam_dev = setup.empty_array_on_device(queue, params.vec_bounds()) + + prg = with_root_kernel(prg, "roe_eigensystem") + prg(queue, nvars=params.nvars, ndim=params.ndim, d=params.d, + states=states, metrics_frozen=metrics_frozen, + R=R_dev, R_inv=Rinv_dev, lambda_roe=lam_dev) + + return R_dev.get(), Rinv_dev.get(), lam_dev.get() + + def mult_mat_vec(queue, prg, alpha, a, b): c_dev = setup.empty_array_on_device(queue, b.shape) diff --git a/setup_fixtures.py b/setup_fixtures.py index bcecb9eb354d4abaa81c7e5a929b6048cd0b3756..6f1debcb3e2bc1b8acb16f45d58bd3a16771a25a 100644 --- a/setup_fixtures.py +++ b/setup_fixtures.py @@ -5,6 +5,19 @@ import pyopencl.array # noqa import device_fixtures as device +class RoeParams: + def __init__(self, nvars, ndim, d): + self.nvars = nvars + self.ndim = ndim + self.d = d + + def mat_bounds(self): + return self.nvars, self.nvars + + def vec_bounds(self): + return self.nvars + + class FluxDerivativeParams: def __init__(self, nvars, ndim, nx, ny, nz): self.nvars = nvars @@ -40,6 +53,11 @@ class FluxDerivativeArrays: self.metric_jacobians = metric_jacobians +def roe_params(nvars, ndim, direction): + dirs = {"x" : 1, "y" : 2, "z" : 3} + return RoeParams(nvars, ndim, dirs[direction]) + + def flux_derivative_params(nvars, ndim, n): return FluxDerivativeParams(nvars, ndim, n, n, n) @@ -48,6 +66,10 @@ def empty_array_on_device(queue, shape): return cl.array.empty(queue, shape, dtype=np.float32, order="F") +def identity(n): + return np.identity(n).astype(np.float32).copy(order="F") + + def random_array(*shape): return np.random.random_sample(shape).astype(np.float32).copy(order="F") @@ -74,3 +96,43 @@ def random_flux_derivative_arrays_on_device(ctx_factory, params): metric_jacobians = random_array_on_device(queue, *params.jacobian_bounds()) return FluxDerivativeArrays(states, fluxes, metrics, metric_jacobians) + + +def arrays_from_string(string_arrays): + return split_map_to_list(string_arrays, array_from_string, ":") + + +def array_from_string(string_array): + if ";" not in string_array: + if "," not in string_array: + array = array_from_string_1d(string_array) + else: + array = array_from_string_2d(string_array) + else: + array = array_from_string_3d(string_array) + return array.copy(order="F") + + +def array_from_string_3d(string_array): + if string_array[0] == ";": + return array_from_string_1d(string_array[1:]).reshape((-1, 1, 1)) + else: + return np.array(split_map_to_list(string_array, array_from_string_2d, ";")) + + +def array_from_string_2d(string_array): + if string_array[0] == ",": + return array_from_string_1d(string_array[1:]).reshape((-1, 1)) + else: + return np.array(split_map_to_list(string_array, array_from_string_1d, ",")) + + +def array_from_string_1d(string_array): + if string_array[0] == "i": + return np.array(split_map_to_list(string_array[1:], int, " ")) + else: + return np.array(split_map_to_list(string_array, float, " "), dtype=np.float32) + + +def split_map_to_list(string, map_func, splitter): + return list(map(map_func, string.split(splitter))) diff --git a/test.py b/test.py index 146e052c994dd5b9f24a078562debb3ee2198bce..924171f561ad96e1eec9984bf6f46c7d8cae589a 100644 --- a/test.py +++ b/test.py @@ -2,6 +2,7 @@ import sys import logging import pytest +import pyopencl as cl from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) @@ -13,6 +14,33 @@ import setup_fixtures as setup import kernel_fixtures as kernel import comparison_fixtures as compare +@pytest.mark.xfail +@pytest.mark.parametrize("states_str,fluxes_str,direction", [ + ("2 1,4 1,4 1,4 1,20 5.5", "4 1,11.2 2.6,8 1,8 1,46.4 7.1", "x"), + ("2 1,4 1,4 1,4 1,20 5.5", "4 1,8 1,11.2 2.6,8 1,46.4 7.1", "y"), + ("2 1,4 1,4 1,4 1,20 5.5", "4 1,8 1,8 1,11.2 2.6,46.4 7.1", "z"), + ("1 2,-1 -4,-1 -4,-1 -4,5.5 20", "-1 -4,2.6 11.2,1 8,1 8,-7.1 -46.4", "x"), + ("1 2,-1 -4,-1 -4,-1 -4,5.5 20", "-1 -4,1 8,2.6 11.2,1 8,-7.1 -46.4", "y"), + ("1 2,-1 -4,-1 -4,-1 -4,5.5 20", "-1 -4,1 8,1 8,2.6 11.2,-7.1 -46.4", "z"), + ("2 1,4 1,8 2,12 3,64 11", "4 1,11.2 2.6,16 2,24 3,134.4 12.6", "x"), + ("2 1,4 1,8 2,12 3,64 11", "8 2,16 2,35.2 5.6,48 6,268.8 25.2", "y"), + ("2 1,4 1,8 2,12 3,64 11", "12 3,24 3,48 6,75.2 10.6,403.2 37.8", "z") + ]) +def test_roe_uniform_grid(ctx_factory, states_str, fluxes_str, direction): + queue = device.get_queue(ctx_factory) + prg = program.get_weno() + + params = setup.roe_params(nvars=5, ndim=3, direction=direction) + states = setup.array_from_string(states_str) + metrics_frozen = setup.identity(params.ndim) + R, Rinv, lam = kernel.roe_eigensystem(queue, prg, params, states, metrics_frozen) + + compare.roe_identity(states, R, Rinv) + + fluxes = setup.array_from_string(fluxes_str) + compare.roe_property(states, fluxes, R, Rinv, lam) + + def test_matvec(ctx_factory): queue = device.get_queue(ctx_factory) prg = program.get_weno() @@ -25,6 +53,7 @@ def test_matvec(ctx_factory): compare.arrays(a@b, c) +@pytest.mark.slow def test_compute_flux_derivatives(ctx_factory): queue = device.get_queue(ctx_factory) prg = program.get_weno() @@ -36,6 +65,7 @@ def test_compute_flux_derivatives(ctx_factory): kernel.compute_flux_derivatives(queue, prg, params, arrays) +@pytest.mark.slow def test_compute_flux_derivatives_gpu(ctx_factory): queue = device.get_queue(ctx_factory) prg = program.get_weno()