Merge branch 'roe-tests' into 'WENO-testing'

Roe tests

See merge request !10

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)

comparison_fixtures.py

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

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)
 

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

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