diff --git a/WENO.F90 b/WENO.F90 index f14987f2d0f58773540698295b695b4d4851edb8..5509694bdf0340a3d92dcd03de17e0ada875209c 100644 --- a/WENO.F90 +++ b/WENO.F90 @@ -47,159 +47,163 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, & ! for the next three loops, note carefully that the inner loop indices have a different range - !$loopy begin tagged: flux_x_compute - do k=1,nz - do j=1,ny - do i=0,nx - call convert_to_generalized_frozen(nvars, ndim, & - states(:, i-2:i+3, j, k), & - fluxes(:, i-2:i+3, j, k), & - metrics(:, :, i:i+1, j, k), & - metric_jacobians(i:i+1, j, k), & - metrics_frozen, & - combined_frozen_metrics, & - generalized_states_frozen, & - generalized_fluxes_frozen) - - call pointwise_eigenvalues(nvars, ndim, 1, states(:, i-2:i+3, j, k), & - metrics(:, :, i-2:i+3, j, k), lambda_pointwise) - call roe_eigensystem(nvars, ndim, 1, states(:, i:i+1, j, k), & - metrics_frozen, R, R_inv, lambda_roe) - - call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) - - call split_characteristic_fluxes(nvars, & - generalized_states_frozen, & - generalized_fluxes_frozen, & - R_inv, & - wavespeeds, & - characteristic_fluxes_pos, & - characteristic_fluxes_neg) - - call weno_flux(nvars, & - fluxes(:, i-2:i+3, j, k), & - characteristic_fluxes_pos, characteristic_fluxes_neg, & - combined_frozen_metrics(1), R, weno_flux_tmp(:, i, j, k)) + if (d == 1) then + !$loopy begin tagged: flux_x_compute + do k=1,nz + do j=1,ny + do i=0,nx + call convert_to_generalized_frozen(nvars, ndim, & + states(:, i-2:i+3, j, k), & + fluxes(:, i-2:i+3, j, k), & + metrics(:, :, i:i+1, j, k), & + metric_jacobians(i:i+1, j, k), & + metrics_frozen, & + combined_frozen_metrics, & + generalized_states_frozen, & + generalized_fluxes_frozen) + + call pointwise_eigenvalues(nvars, ndim, 1, states(:, i-2:i+3, j, k), & + metrics(:, :, i-2:i+3, j, k), lambda_pointwise) + call roe_eigensystem(nvars, ndim, 1, states(:, i:i+1, j, k), & + metrics_frozen, R, R_inv, lambda_roe) + + call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) + + call split_characteristic_fluxes(nvars, & + generalized_states_frozen, & + generalized_fluxes_frozen, & + R_inv, & + wavespeeds, & + characteristic_fluxes_pos, & + characteristic_fluxes_neg) + + call weno_flux(nvars, & + fluxes(:, i-2:i+3, j, k), & + characteristic_fluxes_pos, characteristic_fluxes_neg, & + combined_frozen_metrics(1), R, weno_flux_tmp(:, i, j, k)) + end do + end do end do - end do - end do - !$loopy end tagged: flux_x_compute - - !$loopy begin tagged: flux_x_diff - do k=1,nz - do j=1,ny - do i=1,nx - do v=1,nvars - flux_derivatives(v, i, j, k) = & - (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i-1, j, k)) - end do + !$loopy end tagged: flux_x_compute + + !$loopy begin tagged: flux_x_diff + do k=1,nz + do j=1,ny + do i=1,nx + do v=1,nvars + flux_derivatives(v, i, j, k) = & + (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i-1, j, k)) + end do + end do + end do end do - end do - end do - !$loopy end tagged: flux_x_diff - - !$loopy begin tagged: flux_y_compute - do k=1,nz - do i=1,nx - do j=0,ny - call convert_to_generalized_frozen(nvars, ndim, & - states(:, i, j-2:j+3, k), & - fluxes(:, i, j-2:j+3, k), & - metrics(:, :, i, j:j+1, k), & - metric_jacobians(i, j:j+1, k), & - metrics_frozen, & - combined_frozen_metrics, & - generalized_states_frozen, & - generalized_fluxes_frozen) - - call pointwise_eigenvalues(nvars, ndim, 2, states(:, i, j-2:j+3, k), & - metrics(:, :, i, j-2:j+3, k), lambda_pointwise) - call roe_eigensystem(nvars, ndim, 2, states(:, i, j:j+1, k), & - metrics_frozen, R, R_inv, lambda_roe) - - call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) - - call split_characteristic_fluxes(nvars, & - generalized_states_frozen, & - generalized_fluxes_frozen, & - R_inv, & - wavespeeds, & - characteristic_fluxes_pos, & - characteristic_fluxes_neg) - - call weno_flux(nvars, & - fluxes(:, i, j-2:j+3, k), & - characteristic_fluxes_pos, characteristic_fluxes_neg, & - combined_frozen_metrics(2), R, weno_flux_tmp(:, i, j, k)) - + !$loopy end tagged: flux_x_diff + + else if (d == 2) then + !$loopy begin tagged: flux_y_compute + do k=1,nz + do i=1,nx + do j=0,ny + call convert_to_generalized_frozen(nvars, ndim, & + states(:, i, j-2:j+3, k), & + fluxes(:, i, j-2:j+3, k), & + metrics(:, :, i, j:j+1, k), & + metric_jacobians(i, j:j+1, k), & + metrics_frozen, & + combined_frozen_metrics, & + generalized_states_frozen, & + generalized_fluxes_frozen) + + call pointwise_eigenvalues(nvars, ndim, 2, states(:, i, j-2:j+3, k), & + metrics(:, :, i, j-2:j+3, k), lambda_pointwise) + call roe_eigensystem(nvars, ndim, 2, states(:, i, j:j+1, k), & + metrics_frozen, R, R_inv, lambda_roe) + + call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) + + call split_characteristic_fluxes(nvars, & + generalized_states_frozen, & + generalized_fluxes_frozen, & + R_inv, & + wavespeeds, & + characteristic_fluxes_pos, & + characteristic_fluxes_neg) + + call weno_flux(nvars, & + fluxes(:, i, j-2:j+3, k), & + characteristic_fluxes_pos, characteristic_fluxes_neg, & + combined_frozen_metrics(2), R, weno_flux_tmp(:, i, j, k)) + + end do + end do end do - end do - end do - !$loopy end tagged: flux_y_compute - - !$loopy begin tagged: flux_y_diff - do k=1,nz - do j=1,ny - do i=1,nx - do v=1,nvars - flux_derivatives(v, i, j, k) = & - (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j-1, k)) - end do + !$loopy end tagged: flux_y_compute + + !$loopy begin tagged: flux_y_diff + do k=1,nz + do j=1,ny + do i=1,nx + do v=1,nvars + flux_derivatives(v, i, j, k) = & + (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j-1, k)) + end do + end do + end do end do - end do - end do - !$loopy end tagged: flux_y_diff - - !$loopy begin tagged: flux_z_compute - do j=1,ny - do i=1,nx - do k=0,nz - call convert_to_generalized_frozen(nvars, ndim, & - states(:, i, j, k-2:k+3), & - fluxes(:, i, j, k-2:k+3), & - metrics(:, :, i, j, k:k+1), & - metric_jacobians(i, j, k:k+1), & - metrics_frozen, & - combined_frozen_metrics, & - generalized_states_frozen, & - generalized_fluxes_frozen) - - call pointwise_eigenvalues(nvars, ndim, 3, states(:, i, j, k-2:k+3), & - metrics(:, :, i, j, k-2:k+3), lambda_pointwise) - call roe_eigensystem(nvars, ndim, 3, states(:, i, j, k:k+1), & - metrics_frozen, R, R_inv, lambda_roe) - - call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) - - call split_characteristic_fluxes(nvars, & - generalized_states_frozen, & - generalized_fluxes_frozen, & - R_inv, & - wavespeeds, & - characteristic_fluxes_pos, & - characteristic_fluxes_neg) - - call weno_flux(nvars, & - fluxes(:, i, j, k-2:k+3), & - characteristic_fluxes_pos, characteristic_fluxes_neg, & - combined_frozen_metrics(3), R, weno_flux_tmp(:, i, j, k)) + !$loopy end tagged: flux_y_diff + + else + !$loopy begin tagged: flux_z_compute + do j=1,ny + do i=1,nx + do k=0,nz + call convert_to_generalized_frozen(nvars, ndim, & + states(:, i, j, k-2:k+3), & + fluxes(:, i, j, k-2:k+3), & + metrics(:, :, i, j, k:k+1), & + metric_jacobians(i, j, k:k+1), & + metrics_frozen, & + combined_frozen_metrics, & + generalized_states_frozen, & + generalized_fluxes_frozen) + + call pointwise_eigenvalues(nvars, ndim, 3, states(:, i, j, k-2:k+3), & + metrics(:, :, i, j, k-2:k+3), lambda_pointwise) + call roe_eigensystem(nvars, ndim, 3, states(:, i, j, k:k+1), & + metrics_frozen, R, R_inv, lambda_roe) + + call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) + + call split_characteristic_fluxes(nvars, & + generalized_states_frozen, & + generalized_fluxes_frozen, & + R_inv, & + wavespeeds, & + characteristic_fluxes_pos, & + characteristic_fluxes_neg) + + call weno_flux(nvars, & + fluxes(:, i, j, k-2:k+3), & + characteristic_fluxes_pos, characteristic_fluxes_neg, & + combined_frozen_metrics(3), R, weno_flux_tmp(:, i, j, k)) + end do + end do end do - end do - end do - !$loopy end tagged: flux_z_compute - - !$loopy begin tagged: flux_z_diff - do k=1,nz - do j=1,ny - do i=1,nx - do v=1,nvars - flux_derivatives(v, i, j, k) = & - (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j, k-1)) - end do + !$loopy end tagged: flux_z_compute + + !$loopy begin tagged: flux_z_diff + do k=1,nz + do j=1,ny + do i=1,nx + do v=1,nvars + flux_derivatives(v, i, j, k) = & + (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j, k-1)) + end do + end do + end do end do - end do - end do - !$loopy end tagged: flux_z_diff + !$loopy end tagged: flux_z_diff + end if end subroutine subroutine pointwise_eigenvalues(nvars, ndim, d, states, metrics, lambda_pointwise) diff --git a/test/conftest.py b/test/conftest.py index 67f8e77d0ea2bbee068affa1be83330ff052724c..994bf9111094552282711469499ff1a6cd161602 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -30,7 +30,14 @@ from pyopencl.tools import ( clear_first_arg_caches ) +import numpy as np import utilities as u +from data_for_test import ( + PairData, WindowData, WindowResults, + GridData, GridResults, + halo_sizes, parse_grid_sizes) +import weno_reference_implementation as ref +import wavy_metrics_3d as wavy3d # {{{ mark for slow tests @@ -104,4 +111,246 @@ def queue(ctx_factory): return u.get_queue(ctx_factory) +@pytest.fixture(scope="session") +def pair_data(window_data): + return PairData(window_data) + + +@pytest.fixture(scope="session") +def window_states_generator(): + def window_states(states_str): + if states_str == "sine": + rho = np.sin(0.01*np.arange(6)+1.55) + vel = np.repeat(np.array([1.0, 2.0, 3.0])[None, :], 6, axis=0) + pressure = np.repeat(3.2, 6) + return u.transposed_array([ref.gas.conserved(r, v, p) + for r, v, p + in zip(rho, vel, pressure)]) + else: + state_pair = u.transposed_array_from_string(states_str) + return u.expand_to_n(state_pair, 6) + + return window_states + + +@pytest.fixture(scope="session") +def window_metrics_generator(): + ndim = ref.gas.ndim + + def window_metrics(metric_str, dim): + if (metric_str == "uniform"): + metrics = np.array([np.identity(ndim) for i in range(6)], + dtype=np.float64, + order="F") + jacobians = np.repeat(1.0, 6) + elif (metric_str == "wavy3d"): + metrics = wavy3d.metrics(dim) + jacobians = wavy3d.jacobians(dim) + else: + raise ValueError("Metric string {} not supported".format(metric_str)) + return metrics, jacobians + + return window_metrics + + +@pytest.fixture(scope="session", params=[ + ("sine", "x", "wavy3d"), + ("sine", "y", "wavy3d"), + ("sine", "z", "wavy3d"), + ("sine", "x", "uniform"), + ("sine", "y", "uniform"), + ("sine", "z", "uniform"), + ("1 1 1 1 5.5,1 1 1 1 5.5", "x", "uniform"), + ("1 1 1 1 5.5,1 1 1 1 5.5", "y", "uniform"), + ("1 1 1 1 5.5,1 1 1 1 5.5", "z", "uniform"), + ("2 4 4 4 20,1 1 1 1 5.5", "x", "uniform"), + ("2 4 4 4 20,1 1 1 1 5.5", "y", "uniform"), + ("2 4 4 4 20,1 1 1 1 5.5", "z", "uniform"), + ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "x", "uniform"), + ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "y", "uniform"), + ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "z", "uniform"), + ("2 4 8 12 64,1 1 2 3 11", "x", "uniform"), + ("2 8 12 4 64,1 2 3 1 11", "y", "uniform"), + ("2 12 4 8 64,1 3 1 2 11", "z", "uniform"), + ("1 -1 -2 -3 11,2 -4 -8 -12 64", "x", "uniform"), + ("1 -2 -3 -1 11,2 -8 -12 -4 64", "y", "uniform"), + ("1 -3 -1 -2 11,2 -12 -4 -8 64", "z", "uniform") + ]) +def window_data(request, window_states_generator, window_metrics_generator): + states_str = request.param[0] + dir_str = request.param[1] + metric_str = request.param[2] + + return WindowData(window_states_generator(states_str), + dir_str, + window_metrics_generator(metric_str, dir_str)) + + +@pytest.fixture(scope="session") +def window_results(queue, window_data): + return WindowResults(queue, window_data) + + +@pytest.fixture(scope="session") +def grid_states_generator(): + def grid_states(states_str, grid_sizes_str, dir_str): + nx, ny, nz = parse_grid_sizes(grid_sizes_str) + nxhalo, nyhalo, nzhalo = halo_sizes(nx, ny, nz) + halo = (nxhalo - nx)/2 + states = np.empty( + (ref.gas.nvars, nxhalo, nyhalo, nzhalo), + dtype=np.float64, order="F") + + def is_right(i, j, k): + right_x = (dir_str == "x" and i >= nxhalo/2) + right_y = (dir_str == "y" and j >= nyhalo/2) + right_z = (dir_str == "z" and k >= nzhalo/2) + return right_x or right_y or right_z + + def set_states_jump(state_left, state_right): + for i in range(nxhalo): + for j in range(nyhalo): + for k in range(nzhalo): + if is_right(i, j, k): + states[:, i, j, k] = state_right + else: + states[:, i, j, k] = state_left + + def is_left_hump(i, j, k): + hump_x = (dir_str == "x" and i >= nx/5 + halo and i <= 2*nx/5 + halo) + hump_y = (dir_str == "y" and j >= ny/5 + halo and j <= 2*ny/5 + halo) + hump_z = (dir_str == "z" and k >= nz/5 + halo and k <= 2*nz/5 + halo) + return hump_x or hump_y or hump_z + + def is_right_hump(i, j, k): + hump_x = (dir_str == "x" and i >= 3*nx/5 + halo and i <= 4*nx/5 + halo) + hump_y = (dir_str == "y" and j >= 3*ny/5 + halo and j <= 4*ny/5 + halo) + hump_z = (dir_str == "z" and k >= 3*nz/5 + halo and k <= 4*nz/5 + halo) + return hump_x or hump_y or hump_z + + def set_states_double_hump(state_base, state_left, state_right): + for i in range(nxhalo): + for j in range(nyhalo): + for k in range(nzhalo): + if is_left_hump(i, j, k): + states[:, i, j, k] = state_left + elif is_right_hump(i, j, k): + states[:, i, j, k] = state_right + else: + states[:, i, j, k] = state_base + + if states_str == "flat": + rho = 1.0 + vel = np.array([1.0, 1.0, 1.0]) + pressure = 1.6 + state = ref.gas.conserved(rho, vel, pressure) + set_states_jump(state, state) + elif states_str == "single_jump": + rho_left = 2.0 + vel_left = np.array([2.0, 2.0, 2.0]) + pressure_left = 3.2 + state_left = ref.gas.conserved(rho_left, vel_left, pressure_left) + + rho_right = 1.0 + vel_right = np.array([1.0, 1.0, 1.0]) + pressure_right = 1.6 + state_right = ref.gas.conserved(rho_right, vel_right, pressure_right) + + set_states_jump(state_left, state_right) + elif states_str == "single_jump_varying": + rho_left = 2.0 + vel_left = np.array([2.0, 4.0, 6.0]) + pressure_left = 3.2 + state_left = ref.gas.conserved(rho_left, vel_left, pressure_left) + + rho_right = 1.0 + vel_right = np.array([1.0, 2.0, 3.0]) + pressure_right = 1.6 + state_right = ref.gas.conserved(rho_right, vel_right, pressure_right) + + set_states_jump(state_left, state_right) + elif states_str == "double_hump": + rho_base = 1.0 + vel_base = np.array([1.0, 1.0, 1.0]) + pressure_base = 1.6 + state_base = ref.gas.conserved(rho_base, vel_base, pressure_base) + + rho_left = 2.0 + vel_left = np.array([2.0, 2.0, 2.0]) + pressure_left = 3.2 + state_left = ref.gas.conserved(rho_left, vel_left, pressure_left) + + rho_right = 3.0 + vel_right = np.array([3.0, 3.0, 3.0]) + pressure_right = 4.8 + state_right = ref.gas.conserved(rho_right, vel_right, pressure_right) + + set_states_double_hump(state_base, state_left, state_right) + else: + raise ValueError("States string {} not supported".format(states_str)) + return states + + return grid_states + + +@pytest.fixture(scope="session") +def grid_metrics_generator(): + ndim = ref.gas.ndim + + def grid_metrics(metric_str, grid_sizes_str, dim): + nxhalo, nyhalo, nzhalo = halo_sizes(*parse_grid_sizes(grid_sizes_str)) + if (metric_str == "uniform"): + metrics = np.full( + (ndim, ndim, nxhalo, nyhalo, nzhalo), + 0.0, dtype=np.float64, order="F") + for i in range(nxhalo): + for j in range(nyhalo): + for k in range(nzhalo): + for m in range(ndim): + metrics[m, m, i, j, k] = 1.0 + jacobians = np.full( + (nxhalo, nyhalo, nzhalo), + 1.0, dtype=np.float64, order="F") + else: + raise ValueError("Metric string {} not supported".format(metric_str)) + return metrics, jacobians + + return grid_metrics + + +@pytest.fixture(scope="session", params=[ + ("flat", "10x10x10", "x", "uniform"), + ("flat", "10x10x10", "y", "uniform"), + ("flat", "10x10x10", "z", "uniform"), + ("single_jump", "10x10x10", "x", "uniform"), + ("single_jump", "10x10x10", "y", "uniform"), + ("single_jump", "10x10x10", "z", "uniform"), + ("single_jump", "12x14x16", "x", "uniform"), + ("single_jump", "12x14x16", "y", "uniform"), + ("single_jump", "12x14x16", "z", "uniform"), + ("single_jump_varying", "10x10x10", "x", "uniform"), + ("single_jump_varying", "10x10x10", "y", "uniform"), + ("single_jump_varying", "10x10x10", "z", "uniform"), + ("double_hump", "25x25x25", "x", "uniform"), + ("double_hump", "25x25x25", "y", "uniform"), + ("double_hump", "25x25x25", "z", "uniform"), + # FIXME: add remaining cases + ]) +def grid_data(request, grid_states_generator, grid_metrics_generator): + states_str = request.param[0] + grid_sizes_str = request.param[1] + dir_str = request.param[2] + metric_str = request.param[3] + + return GridData(grid_states_generator(states_str, grid_sizes_str, dir_str), + grid_sizes_str, + dir_str, + grid_metrics_generator(metric_str, grid_sizes_str, dir_str)) + + +@pytest.fixture(scope="session") +def grid_results(queue, grid_data): + return GridResults(queue, grid_data) + + # vim: foldmethod=marker diff --git a/test/data_for_test.py b/test/data_for_test.py index 990c5eef0042c6e001cf94015eb7331610f5adb7..2cf9c2fea54857bec0cbc5dcac083c512960b620 100644 --- a/test/data_for_test.py +++ b/test/data_for_test.py @@ -21,10 +21,23 @@ THE SOFTWARE. """ import numpy as np -import pytest -import utilities as u import weno_reference_implementation as ref -import wavy_metrics_3d as wavy3d + + +class PairData: + def __init__(self, data): + self.nvars = data.nvars + self.ndim = data.ndim + self.direction = data.direction + self.dir_internal = data.dir_internal + + self.state_pair = data.states[:, 2:4] + self.generalized_flux_pair = data.generalized_fluxes[:, 2:4] + + weighted_metric_slice = data.metrics[2:4]/data.jacobians[2:4, None, None] + self.frozen_metrics = np.mean(weighted_metric_slice, axis=0) + self.frozen_jacobian = np.mean(data.jacobians[2:4]) + self.combined_frozen_metrics = np.sum(self.frozen_metrics**2, axis=1) class WindowData: @@ -32,11 +45,11 @@ class WindowData: ndim = ref.gas.ndim dirs = ref.gas.dirs - def __init__(self, states, direction, window_metrics): + def __init__(self, window_states, direction, window_metrics): self.direction = self.dirs[direction] self.dir_internal = self.direction-1 - self.states = states + self.states = window_states self.metrics, self.jacobians = window_metrics self.fluxes = ref.pointwise_fluxes(self.states) @@ -46,95 +59,164 @@ class WindowData: )[:, :, self.dir_internal].T.copy(order="F") -@pytest.fixture(scope="session") -def states_generator(): - def states(states_str): - if states_str == "sine": - rho = np.sin(0.01*np.arange(6)+1.55) - vel = np.repeat(np.array([1.0, 2.0, 3.0])[None, :], 6, axis=0) - pressure = np.repeat(3.2, 6) - return u.transposed_array([ref.gas.conserved(r, v, p) - for r, v, p - in zip(rho, vel, pressure)]) - else: - state_pair = u.transposed_array_from_string(states_str) - return u.expand_to_n(state_pair, 6) +class WindowResults: + def __init__(self, queue, data): + pair = PairData(data) + + self.nvars = data.nvars + self.ndim = data.ndim + self.direction = data.direction + self.dir_internal = data.dir_internal + + self.states = data.states + self.state_pair = pair.state_pair + + self.metrics = data.metrics + self.jacobians = data.jacobians + + self.frozen_metrics = pair.frozen_metrics + self.frozen_jacobian = pair.frozen_jacobian + self.combined_frozen_metric = pair.combined_frozen_metrics[self.dir_internal] + + self.generalized_states_frozen = self.states/self.frozen_jacobian + + self.generalized_fluxes = data.generalized_fluxes + self.generalized_fluxes_frozen = np.array([np.dot(flux, self.frozen_metrics) + for flux in data.fluxes])[:, :, self.dir_internal].T.copy(order="F") + + self.lam_pointwise = ref.lambda_pointwise( + self.states, self.metrics, self.dir_internal) + self.R, self.R_inv, self.lam_roe = ref.roe_eigensystem( + queue, self.state_pair, self.frozen_metrics, self.direction) + self.wavespeeds = ref.wavespeeds(self.lam_pointwise, self.lam_roe) + + self.char_fluxes_pos, self.char_fluxes_neg = ref.split_char_fluxes( + self.states, self.wavespeeds, + self.frozen_metrics[self.dir_internal], self.frozen_jacobian, + self.R_inv) - return states + self.oscillation_pos = ref.oscillation(self.char_fluxes_pos) + self.oscillation_neg = ref.oscillation(self.char_fluxes_neg[:, ::-1]) + self.weno_weights_pos = ref.weno_weights( + self.oscillation_pos, self.combined_frozen_metric) + self.weno_weights_neg = ref.weno_weights( + self.oscillation_neg, self.combined_frozen_metric) -@pytest.fixture(scope="session") -def metrics_generator(): + self.consistent = ref.consistent_part(self.generalized_fluxes) + self.dissipation_pos = ref.dissipation_part( + self.R, self.char_fluxes_pos, self.weno_weights_pos, 1) + self.dissipation_neg = ref.dissipation_part( + self.R, self.char_fluxes_neg, self.weno_weights_neg, -1) + + self.weno_flux = ref.weno_flux( + self.consistent, self.dissipation_pos, self.dissipation_neg) + + +halo = 3 + + +def parse_grid_sizes(grid_sizes_str): + return map(int, grid_sizes_str.split("x")) + + +def halo_sizes(*grid_sizes): + return map(lambda n: n + 2*halo, grid_sizes) + + +class GridData: ndim = ref.gas.ndim + nvars = ref.gas.nvars + dirs = ref.gas.dirs - def window_metrics(metric_str, dim): - if (metric_str == "uniform"): - metrics = np.array([np.identity(ndim) for i in range(6)], - dtype=np.float64, - order="F") - jacobians = np.repeat(1.0, 6) - elif (metric_str == "wavy3d"): - metrics = wavy3d.metrics(dim) - jacobians = wavy3d.jacobians(dim) - else: - raise ValueError("Metric string {} not supported".format(metric_str)) - return metrics, jacobians - - return window_metrics - - -@pytest.fixture(scope="session", params=[ - ("sine", "x", "wavy3d"), - ("sine", "y", "wavy3d"), - ("sine", "z", "wavy3d"), - ("sine", "x", "uniform"), - ("sine", "y", "uniform"), - ("sine", "z", "uniform"), - ("1 1 1 1 5.5,1 1 1 1 5.5", "x", "uniform"), - ("1 1 1 1 5.5,1 1 1 1 5.5", "y", "uniform"), - ("1 1 1 1 5.5,1 1 1 1 5.5", "z", "uniform"), - ("2 4 4 4 20,1 1 1 1 5.5", "x", "uniform"), - ("2 4 4 4 20,1 1 1 1 5.5", "y", "uniform"), - ("2 4 4 4 20,1 1 1 1 5.5", "z", "uniform"), - ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "x", "uniform"), - ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "y", "uniform"), - ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "z", "uniform"), - ("2 4 8 12 64,1 1 2 3 11", "x", "uniform"), - ("2 8 12 4 64,1 2 3 1 11", "y", "uniform"), - ("2 12 4 8 64,1 3 1 2 11", "z", "uniform"), - ("1 -1 -2 -3 11,2 -4 -8 -12 64", "x", "uniform"), - ("1 -2 -3 -1 11,2 -8 -12 -4 64", "y", "uniform"), - ("1 -3 -1 -2 11,2 -12 -4 -8 64", "z", "uniform") - ]) -def window_data(request, states_generator, metrics_generator): - states_str = request.param[0] - dir_str = request.param[1] - metric_str = request.param[2] - - return WindowData(states_generator(states_str), - dir_str, - metrics_generator(metric_str, dir_str)) + def __init__(self, grid_states, grid_sizes_str, dir_str, grid_metrics): + self.dir_str = dir_str + self.direction = self.dirs[self.dir_str] + self.dir_internal = self.direction-1 + self.set_grid_sizes(grid_sizes_str) + self.states = grid_states + self.metrics, self.jacobians = grid_metrics -class PairData: - def __init__(self, data): + self.generalized_fluxes = np.empty( + (self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), + dtype=np.float64, order="F") + for i in range(self.nxhalo): + for j in range(self.nyhalo): + for k in range(self.nzhalo): + flux = ref.gas.flux(self.states[:, i, j, k]) + self.generalized_fluxes[:, i, j, k] = np.dot( + flux, + self.metrics[:, :, i, j, k]/self.jacobians[i, j, k] + )[:, self.dir_internal] + + def set_grid_sizes(self, grid_sizes_str): + self.nx, self.ny, self.nz = parse_grid_sizes(grid_sizes_str) + + self.nxhalo, self.nyhalo, self.nzhalo = halo_sizes(self.nx, self.ny, self.nz) + + self.output_dims = (self.nvars, self.nx, self.ny, self.nz) + + +class GridResults: + def __init__(self, queue, data): + self.queue = queue + self.nx = data.nx + self.ny = data.ny + self.nz = data.nz self.nvars = data.nvars self.ndim = data.ndim + self.output_dims = data.output_dims + self.dir_str = data.dir_str self.direction = data.direction self.dir_internal = data.dir_internal - self.state_pair = data.states[:, 2:4] - self.generalized_flux_pair = data.generalized_fluxes[:, 2:4] - - weighted_metric_slice = data.metrics[2:4]/data.jacobians[2:4, None, None] - self.frozen_metrics = np.mean(weighted_metric_slice, axis=0) - self.frozen_jacobian = np.mean(data.jacobians[2:4]) - self.combined_frozen_metrics = np.sum(self.frozen_metrics**2, axis=1) - - -@pytest.fixture(scope="session") -def pair_data(window_data): - return PairData(window_data) + self.states = data.states + self.generalized_fluxes = data.generalized_fluxes + self.metrics = data.metrics + self.jacobians = data.jacobians + + self.generalized_flux_derivatives = np.full(self.output_dims, 0.0, + dtype=np.float64, order="F") + for i in range(self.nx): + for j in range(self.ny): + for k in range(self.nz): + flux_pos = self.local_flux( + *self.window_variables(i, j, k)) + neg_loc = ( + i - (1 if 0 == self.dir_internal else 0), + j - (1 if 1 == self.dir_internal else 0), + k - (1 if 2 == self.dir_internal else 0)) + flux_neg = self.local_flux( + *self.window_variables(*neg_loc)) + self.generalized_flux_derivatives[:, i, j, k] = ( + flux_pos - flux_neg) + + def window_variables(self, i, j, k): + def lo(n, d): + return n + halo - (2 if d == self.dir_internal else 0) + + def hi(n, d): + return n + halo + (4 if d == self.dir_internal else 1) + + xlo = lo(i, 0) + xhi = hi(i, 0) + ylo = lo(j, 1) + yhi = hi(j, 1) + zlo = lo(k, 2) + zhi = hi(k, 2) + local_states = self.states[:, xlo:xhi, ylo:yhi, zlo:zhi] + local_metrics = self.metrics[:, :, xlo:xhi, ylo:yhi, zlo:zhi] + local_jacobians = self.jacobians[xlo:xhi, ylo:yhi, zlo:zhi] + return (np.squeeze(local_states), + np.moveaxis(np.squeeze(local_metrics), -1, 0), + np.squeeze(local_jacobians)) + + def local_flux(self, state_window, metric_window, jacobian_window): + local_data = WindowData( + state_window, self.dir_str, + (metric_window, jacobian_window)) + return WindowResults(self.queue, local_data).weno_flux # vim: foldmethod=marker diff --git a/test/test_eigensystem.py b/test/test_eigensystem.py index b1a8448505f5f64a6264ca6db8f67413ddfb50e8..e96c05a48e03d2e0acaee1987279d29eea77e0fd 100644 --- a/test/test_eigensystem.py +++ b/test/test_eigensystem.py @@ -34,12 +34,9 @@ import logging import pytest import utilities as u -import weno_reference_implementation as ref # noqa -from data_for_test import ( # noqa - pair_data, window_data, states_generator, metrics_generator) -def test_roe_uniform_grid_ideal_gas(queue, pair_data): # noqa +def test_roe_uniform_grid_ideal_gas(queue, pair_data): data = pair_data def check_roe_identity(states, R, R_inv): diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index e80cf3d59ad0d8aa21df52a172076f9fd6af6452..b52c59b8d864706a940dbaa3a07d508b7fc767f5 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -20,7 +20,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -import numpy as np +import numpy as np # noqa: F401 import numpy.linalg as la # noqa: F401 import pyopencl as cl # noqa: F401 import pyopencl.array # noqa @@ -34,133 +34,15 @@ import logging import pytest import utilities as u -import weno_reference_implementation as ref -# FIXME: is there a better way to divide responsibilities with these fixture classes? -class FluxDataVector: - nvars = 5 - ndim = 3 - dirs = {"x": 1, "y": 2, "z": 3} - halo = 3 - - def __init__(self, nx, ny, nz, states_str, direction): - self.direction = self.dirs[direction] - self.dir_internal = self.direction-1 - - self.nx = nx - self.ny = ny - self.nz = nz - - self.nxhalo = self.nx + 2*self.halo - self.nyhalo = self.ny + 2*self.halo - self.nzhalo = self.nz + 2*self.halo - - self.flux_dims = (self.nvars, self.nx, self.ny, self.nz) - - self.metrics = np.stack( - [np.stack( - [np.stack( - [np.identity(self.ndim) for i in range(self.nxhalo)], - axis=-1) for j in range(self.nyhalo)], - axis=-1) for k in range(self.nzhalo)], - axis=-1).copy(order="F") - self.jacobians = np.ones((self.nxhalo, self.nyhalo, self.nzhalo), - order="F") - - # FIXME: Move array_from_string stuff outside WindowResults - # -- just pass an array & have external utilities that generate - # Riemann, sine wave, etc. initial conditions - # FIXME: Consider handling row swapping outside as well? - # FIXME: Do we even need to swap rows? - self.state_pair = self.swap_array_rows( - u.transposed_array_from_string(states_str), self.dir_internal) - # NOTE: dimensions are nvars x nxhalo x nyhalo x nzhalo - self.states = self.fill_from_pair() - - # NOTE: dimensions are nvars x nxhalo x nyhalo x nzhalo - # FIXME: these should be generalized fluxes - # FIXME: make a clear distinction between fluxes in physical and - # generalized coordinates - npoints = self.nxhalo*self.nyhalo*self.nzhalo - flat_states = self.states.reshape((self.nvars, npoints)) - self.fluxes = ref.pointwise_fluxes( - flat_states)[:, :, self.dir_internal].T.reshape( - (self.nvars, self.nxhalo, self.nyhalo, self.nzhalo) - ).copy(order="F") - - # FIXME: use reference implementation - # NOTE: dimensions are nvars x nx x ny x nz - self.flux_derivatives = np.zeros((self.nvars, self.nx, self.ny, - self.nz), order="F") - - def swap_array_rows(self, arr, d): - p = self.permutation(d) - arr[p, :] = arr[[1, 2, 3], :] - return arr - - def permutation(self, d): - return [(d-1+i) % 3 + 1 for i in range(3)] - - def fill_from_pair(self): - d = self.dir_internal - nx_arr = np.array([self.nxhalo, self.nyhalo, self.nzhalo]) - result = u.expand_to_n(self.state_pair, nx_arr[d]) - - for i in range(d): - result = self.add_dimension(result, nx_arr[i]) - result = np.swapaxes(result, -2, -1) - for i in range(d+1, self.ndim): - result = self.add_dimension(result, nx_arr[i]) - - return result.copy(order="F") - - def add_dimension(self, arr, n): - return np.stack([arr for i in range(n)], axis=-1) - - -vector_data = {} - -vector_data["Case flat:x"] = FluxDataVector( - nx=6, ny=2, nz=2, - states_str="1 1 1 1 5.5,1 1 1 1 5.5", - direction="x") -vector_data["Case flat:y"] = FluxDataVector( - nx=2, ny=6, nz=2, - states_str="1 1 1 1 5.5,1 1 1 1 5.5", - direction="y") -vector_data["Case flat:z"] = FluxDataVector( - nx=2, ny=2, nz=6, - states_str="1 1 1 1 5.5,1 1 1 1 5.5", - direction="z") - -vector_data["Case a:x"] = FluxDataVector( - nx=6, ny=2, nz=2, - states_str="2 4 4 4 20,1 1 1 1 5.5", - direction="x") -vector_data["Case a:y"] = FluxDataVector( - nx=2, ny=6, nz=2, - states_str="2 4 4 4 20,1 1 1 1 5.5", - direction="y") -vector_data["Case a:z"] = FluxDataVector( - nx=2, ny=2, nz=6, - states_str="2 4 4 4 20,1 1 1 1 5.5", - direction="z") - - -@pytest.fixture(scope="session", params=[ - "Case flat:x", "Case flat:y", "Case flat:z"]) -def cfd_test_data_fixture(request): - return vector_data[request.param] - - -def test_compute_flux_derivatives_uniform_grid(queue, cfd_test_data_fixture): - data = cfd_test_data_fixture +def test_compute_flux_derivatives(queue, grid_results): + data = grid_results prg = u.get_weno_program() prg = prg.copy(target=lp.PyOpenCLTarget(queue.device)) - flux_derivatives_dev = u.empty_array_on_device(queue, *data.flux_dims) + flux_derivatives_dev = u.empty_array_on_device(queue, *data.output_dims) prg(queue, nvars=data.nvars, @@ -170,12 +52,12 @@ def test_compute_flux_derivatives_uniform_grid(queue, cfd_test_data_fixture): nz=data.nz, d=data.direction, states=data.states, - fluxes=data.fluxes, + fluxes=data.generalized_fluxes, metrics=data.metrics, metric_jacobians=data.jacobians, flux_derivatives=flux_derivatives_dev) - u.compare_arrays(flux_derivatives_dev.get(), data.flux_derivatives) + u.compare_arrays(flux_derivatives_dev.get(), data.generalized_flux_derivatives) # This lets you run 'python test.py test_case(cl._csc)' without pytest. diff --git a/test/test_flux_window_ops.py b/test/test_flux_window_ops.py index dbd47415f89c710d751e0ba5324864d906e0d8f4..367bb77ba716cc0aeacb6a3283b1b080df3ee498 100644 --- a/test/test_flux_window_ops.py +++ b/test/test_flux_window_ops.py @@ -34,68 +34,6 @@ import logging import pytest import utilities as u -import weno_reference_implementation as ref -from data_for_test import ( # noqa - PairData, window_data, states_generator, metrics_generator) - - -class WindowResults: - def __init__(self, queue, data): - pair = PairData(data) - - self.nvars = data.nvars - self.ndim = data.ndim - self.direction = data.direction - self.dir_internal = data.dir_internal - - self.states = data.states - self.state_pair = pair.state_pair - - self.metrics = data.metrics - self.jacobians = data.jacobians - - self.frozen_metrics = pair.frozen_metrics - self.frozen_jacobian = pair.frozen_jacobian - self.combined_frozen_metric = pair.combined_frozen_metrics[self.dir_internal] - - self.generalized_states_frozen = self.states/self.frozen_jacobian - - self.generalized_fluxes = data.generalized_fluxes - self.generalized_fluxes_frozen = np.array([np.dot(flux, self.frozen_metrics) - for flux in data.fluxes])[:, :, self.dir_internal].T.copy(order="F") - - self.lam_pointwise = ref.lambda_pointwise( - self.states, self.metrics, self.dir_internal) - self.R, self.R_inv, self.lam_roe = ref.roe_eigensystem( - queue, self.state_pair, self.frozen_metrics, self.direction) - self.wavespeeds = ref.wavespeeds(self.lam_pointwise, self.lam_roe) - - self.char_fluxes_pos, self.char_fluxes_neg = ref.split_char_fluxes( - self.states, self.wavespeeds, - self.frozen_metrics[self.dir_internal], self.frozen_jacobian, - self.R_inv) - - self.oscillation_pos = ref.oscillation(self.char_fluxes_pos) - self.oscillation_neg = ref.oscillation(self.char_fluxes_neg[:, ::-1]) - - self.weno_weights_pos = ref.weno_weights( - self.oscillation_pos, self.combined_frozen_metric) - self.weno_weights_neg = ref.weno_weights( - self.oscillation_neg, self.combined_frozen_metric) - - self.consistent = ref.consistent_part(self.generalized_fluxes) - self.dissipation_pos = ref.dissipation_part( - self.R, self.char_fluxes_pos, self.weno_weights_pos, 1) - self.dissipation_neg = ref.dissipation_part( - self.R, self.char_fluxes_neg, self.weno_weights_neg, -1) - - self.weno_flux = ref.weno_flux( - self.consistent, self.dissipation_pos, self.dissipation_neg) - - -@pytest.fixture(scope="session") -def window_results(queue, window_data): # noqa - return WindowResults(queue, window_data) def test_pointwise_eigenvalues(queue, window_results): diff --git a/test/utilities.py b/test/utilities.py index 3f1fa0714c28f6c4051b07167842f638759cfa64..793a40721cedebdfc042d67dc73c86c22cccbdfd 100644 --- a/test/utilities.py +++ b/test/utilities.py @@ -35,7 +35,7 @@ import os.path # {{{ arrays def compare_arrays(a, b): - assert a == approx(b, rel=1e-12, abs=1e-14) + assert a == approx(b, rel=1e-12, abs=1e-13) def transposed_array(a): diff --git a/test/weno_reference_implementation.py b/test/weno_reference_implementation.py index 6215ba28150f164714ec442b5a2ed51cbe4497c8..30f1d80784e8d29e6c98d00ab63ac843da70637c 100644 --- a/test/weno_reference_implementation.py +++ b/test/weno_reference_implementation.py @@ -52,7 +52,8 @@ def roe_eigensystem(queue, state_pair, frozen_metrics, direction): lam_dev = u.empty_array_on_device(queue, nvars) prg(queue, nvars=nvars, ndim=ndim, d=direction, - states=state_pair, metrics_frozen=frozen_metrics, + states=state_pair.copy(order='F'), + metrics_frozen=frozen_metrics.copy(order='F'), R=R_dev, R_inv=R_inv_dev, lambda_roe=lam_dev) return R_dev.get(), R_inv_dev.get(), lam_dev.get()