diff --git a/WENO.F90 b/WENO.F90 index 6bfad5c1d58ed39569fb0a308a426215d35a13e8..829080e5f56556ffcc276e33a5ad8d0aa0628ab6 100644 --- a/WENO.F90 +++ b/WENO.F90 @@ -303,9 +303,9 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam if (im > 3) im = im - 3 do j=1,2 - u_orig(1,j) = states(ik+1,j)/states(1,j) - u_orig(2,j) = states(il+1,j)/states(1,j) - u_orig(3,j) = states(im+1,j)/states(1,j) + do i=1,ndim + u_orig(i,j) = states(i+1,j)/states(1,j) + end do end do do j=1,2 diff --git a/test.py b/test.py index 0a27274293198443e4d3bec4a80f7094e7f75a0d..97ccfe7cc1d668e3ccfeae5e4bb4dfd5410686ca 100644 --- a/test.py +++ b/test.py @@ -18,50 +18,24 @@ from pyopencl.tools import ( # noqa from utilities import * -@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") + ("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(ctx_factory, states_str, fluxes_str, direction): - 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_bound(self): - return self.nvars - - def setup_roe_params(nvars, ndim, direction): - dirs = {"x" : 1, "y" : 2, "z" : 3} - return RoeParams(nvars, ndim, dirs[direction]) - def identity_matrix(n): return np.identity(n).astype(np.float32).copy(order="F") - def kernel_roe_eigensystem(queue, prg, params, states, metrics_frozen): - R_dev = empty_array_on_device(queue, *params.mat_bounds()) - Rinv_dev = empty_array_on_device(queue, *params.mat_bounds()) - lam_dev = empty_array_on_device(queue, params.vec_bound()) - - 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 check_roe_identity(states, R, Rinv): dState = states[:,1] - states[:,0] compare_arrays(R@(Rinv@dState), dState) @@ -74,22 +48,37 @@ def test_roe_uniform_grid(ctx_factory, states_str, fluxes_str, direction): temp = np.multiply(lam, temp) compare_arrays(R@temp, dFlux) + prg = get_weno_program_with_root_kernel("roe_eigensystem") queue = get_queue(ctx_factory) - prg = get_weno_program() - params = setup_roe_params(nvars=5, ndim=3, direction=direction) - states = array_from_string(states_str) - metrics_frozen = identity_matrix(params.ndim) - R, Rinv, lam = kernel_roe_eigensystem(queue, prg, params, states, metrics_frozen) + nvars = 5 + ndim = 3 + dirs = {"x" : 1, "y" : 2, "z" : 3} - check_roe_identity(states, R, Rinv) + states = transposed_array_from_string(states_str) + fluxes = transposed_array_from_string(fluxes_str) + metrics_frozen = identity_matrix(ndim) + + R_dev = empty_array_on_device(queue, nvars, nvars) + Rinv_dev = empty_array_on_device(queue, nvars, nvars) + lam_dev = 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=Rinv_dev, lambda_roe=lam_dev) + + R = R_dev.get() + Rinv = Rinv_dev.get() + lam = lam_dev.get() - fluxes = array_from_string(fluxes_str) + print(lam) + + check_roe_identity(states, R, Rinv) check_roe_property(states, fluxes, R, Rinv, lam) def test_matvec(ctx_factory): - prg = get_weno_program() + prg = get_weno_program_with_root_kernel("mult_mat_vec") queue = get_queue(ctx_factory) a = random_array_on_device(queue, 10, 10) @@ -97,7 +86,6 @@ def test_matvec(ctx_factory): c = empty_array_on_device(queue, 10) - prg = with_root_kernel(prg, "mult_mat_vec") prg(queue, alpha=1.0, a=a, b=b, c=c) compare_arrays(a.get()@b.get(), c.get()) diff --git a/utilities.py b/utilities.py index 22834325b4d6478490a5168143fd2937dbccfdaf..b825fcbe82ee7ef0d2a991edf9ed70f136b55f80 100644 --- a/utilities.py +++ b/utilities.py @@ -28,6 +28,10 @@ def arrays_from_string(string_arrays): return split_map_to_list(string_arrays, array_from_string, ":") +def transposed_array_from_string(string_array): + return array_from_string(string_array).transpose().copy(order="F") + + def array_from_string(string_array): def array_from_string_1d(string_array): if string_array[0] == "i": @@ -72,13 +76,7 @@ def get_queue(ctx_factory): # {{{ program / kernel -_WENO_PRG = [] - - -def get_weno_program(): - if not _WENO_PRG: - parse_weno() - return _WENO_PRG[0] +_WENO_PRG = {} def parse_weno(): @@ -88,7 +86,13 @@ def parse_weno(): infile_content = infile.read() prg = lp.parse_transformed_fortran(infile_content, filename=fn) - _WENO_PRG.append(prg) + _WENO_PRG["default"] = prg + + +def get_weno_program(): + if not _WENO_PRG: + parse_weno() + return _WENO_PRG["default"] def with_root_kernel(prg, root_name): @@ -103,6 +107,18 @@ def with_root_kernel(prg, root_name): return new_prg +def setup_weno_program_with_root_kernel(knl): + prg = get_weno_program() + prg = with_root_kernel(prg, knl) + _WENO_PRG[knl] = prg + + +def get_weno_program_with_root_kernel(knl): + if not knl in _WENO_PRG: + setup_weno_program_with_root_kernel(knl) + return _WENO_PRG[knl] + + def transform_weno_for_gpu(prg, print_kernel=False): cfd = prg["compute_flux_derivatives"]