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()