diff --git a/WENO.F90 b/WENO.F90 index d798cdbe7f82a973d88be3506e1a6a9055187929..5a7a6f44e3d9594df952764cf52d7019192c7875 100644 --- a/WENO.F90 +++ b/WENO.F90 @@ -509,14 +509,14 @@ subroutine split_characteristic_fluxes(nvars, & end subroutine subroutine weno_flux(nvars, generalized_fluxes, characteristic_fluxes_pos, & - characteristic_fluxes_neg, combined_frozen_metrics, R, flux) + characteristic_fluxes_neg, combined_frozen_metric, R, flux) implicit none integer, intent(in) :: nvars real*8, intent(in) :: generalized_fluxes(nvars, -2:3) real*8, intent(in) :: characteristic_fluxes_pos(nvars, -2:3) real*8, intent(in) :: characteristic_fluxes_neg(nvars, -2:3) - real*8, intent(in) :: combined_frozen_metrics + real*8, intent(in) :: combined_frozen_metric real*8, intent(in) :: R(nvars, nvars) real*8, intent(out) :: flux(nvars) @@ -527,9 +527,9 @@ subroutine weno_flux(nvars, generalized_fluxes, characteristic_fluxes_pos, & call consistent_part(nvars, generalized_fluxes, consistent) call dissipation_part_pos(nvars, characteristic_fluxes_pos, & - combined_frozen_metrics, R, dissipation_pos) + combined_frozen_metric, R, dissipation_pos) call dissipation_part_neg(nvars, characteristic_fluxes_neg, & - combined_frozen_metrics, R, dissipation_neg) + combined_frozen_metric, R, dissipation_neg) do v=1,nvars flux(v) = consistent(v) + dissipation_pos(v) + dissipation_neg(v) @@ -554,36 +554,36 @@ subroutine consistent_part(nvars, generalized_fluxes, consistent) end subroutine subroutine dissipation_part_pos(nvars, characteristic_fluxes, & - combined_frozen_metrics, R, dissipation_pos) + combined_frozen_metric, R, dissipation_pos) implicit none integer, intent(in) :: nvars real*8, intent(in) :: characteristic_fluxes(nvars, -2:3) - real*8, intent(in) :: combined_frozen_metrics + real*8, intent(in) :: combined_frozen_metric real*8, intent(in) :: R(nvars, nvars) real*8, intent(out) :: dissipation_pos(nvars) real*8 combined_fluxes(nvars) - call weno_combination_pos(nvars, characteristic_fluxes, combined_frozen_metrics, combined_fluxes) + call weno_combination_pos(nvars, characteristic_fluxes, combined_frozen_metric, combined_fluxes) call mult_mat_vec(nvars, nvars, -1.0d0/60, R, combined_fluxes, dissipation_pos) end subroutine subroutine weno_combination_pos(nvars, characteristic_fluxes, & - combined_frozen_metrics, combined_fluxes) + combined_frozen_metric, combined_fluxes) implicit none integer, intent(in) :: nvars real*8, intent(in) :: characteristic_fluxes(nvars, -2:3) - real*8, intent(in) :: combined_frozen_metrics + real*8, intent(in) :: combined_frozen_metric real*8, intent(out) :: combined_fluxes(nvars) real*8 w(nvars, 3) real*8 flux_differences(nvars, 3) integer v - call weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metrics, w) + call weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric, w) call flux_differences_pos(nvars, characteristic_fluxes, flux_differences) do v=1,nvars @@ -593,12 +593,12 @@ subroutine weno_combination_pos(nvars, characteristic_fluxes, & end do end subroutine -subroutine weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metrics, w) +subroutine weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric, w) implicit none integer, intent(in) :: nvars real*8, intent(in) :: characteristic_fluxes(nvars, -2:3) - real*8, intent(in) :: combined_frozen_metrics + real*8, intent(in) :: combined_frozen_metric real*8, intent(out) :: w(nvars, 3) real*8 :: C(3) = (/0.1d0, 0.6d0, 0.3d0/) @@ -607,7 +607,7 @@ subroutine weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric integer p, i, j real*8 sum_alpha(1) - eps = 1.0d-6!*combined_frozen_metrics + eps = 1.0d-6!*combined_frozen_metric p = 2 call oscillation_pos(nvars, characteristic_fluxes, IS) @@ -671,36 +671,36 @@ subroutine flux_differences_pos(nvars, characteristic_fluxes, flux_differences) end subroutine subroutine dissipation_part_neg(nvars, characteristic_fluxes, & - combined_frozen_metrics, R, dissipation_neg) + combined_frozen_metric, R, dissipation_neg) implicit none integer, intent(in) :: nvars real*8, intent(in) :: characteristic_fluxes(nvars, -2:3) - real*8, intent(in) :: combined_frozen_metrics + real*8, intent(in) :: combined_frozen_metric real*8, intent(in) :: R(nvars, nvars) real*8, intent(out) :: dissipation_neg(nvars) real*8 combined_fluxes(nvars) - call weno_combination_neg(nvars, characteristic_fluxes, combined_frozen_metrics, combined_fluxes) + call weno_combination_neg(nvars, characteristic_fluxes, combined_frozen_metric, combined_fluxes) call mult_mat_vec(nvars, nvars, 1.0d0/60, R, combined_fluxes, dissipation_neg) end subroutine subroutine weno_combination_neg(nvars, characteristic_fluxes, & - combined_frozen_metrics, combined_fluxes) + combined_frozen_metric, combined_fluxes) implicit none integer, intent(in) :: nvars real*8, intent(in) :: characteristic_fluxes(nvars, -2:3) - real*8, intent(in) :: combined_frozen_metrics + real*8, intent(in) :: combined_frozen_metric real*8, intent(out) :: combined_fluxes(nvars) real*8 w(nvars, 3) real*8 flux_differences(nvars, 3) integer v - call weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metrics, w) + call weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric, w) call flux_differences_neg(nvars, characteristic_fluxes, flux_differences) do v=1,nvars @@ -710,12 +710,12 @@ subroutine weno_combination_neg(nvars, characteristic_fluxes, & end do end subroutine -subroutine weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metrics, w) +subroutine weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric, w) implicit none integer, intent(in) :: nvars real*8, intent(in) :: characteristic_fluxes(nvars, -2:3) - real*8, intent(in) :: combined_frozen_metrics + real*8, intent(in) :: combined_frozen_metric real*8, intent(out) :: w(nvars, 3) real*8 :: C(3) = (/0.1d0, 0.6d0, 0.3d0/) @@ -724,7 +724,7 @@ subroutine weno_weights_neg(nvars, characteristic_fluxes, combined_frozen_metric integer p, i, j real*8 sum_alpha(1) - eps = 1.0d-6!*combined_frozen_metrics + eps = 1.0d-6!*combined_frozen_metric p = 2 call oscillation_neg(nvars, characteristic_fluxes, IS) diff --git a/test/data_for_test.py b/test/data_for_test.py index 3bba73debd927eb3fd65b504cbf0f5bf494ecfe8..483d50b65f311b8cabea7497b1b1626590b62032 100644 --- a/test/data_for_test.py +++ b/test/data_for_test.py @@ -26,226 +26,108 @@ import utilities as u import weno_reference_implementation as ref -# {{{ FluxDataSingle - -class FluxDataSingle: - # FIXME: can we set some of these constants from ref.gas? - # -- if all nvars references come from there, it's relatively easy to - # introduce a new gas with more (e.g. scalar) variables - nvars = 5 - ndim = 3 - dirs = {"x": 1, "y": 2, "z": 3} - - def __init__(self, queue, states_str, direction): +class WindowData: + nvars = ref.gas.nvars + ndim = ref.gas.ndim + dirs = ref.gas.dirs + def __init__(self, states, direction, window_metrics): self.direction = self.dirs[direction] self.dir_internal = self.direction-1 - self.metrics = np.array([np.identity(self.ndim) for i in range(6)], - dtype=np.float64, order="F") - self.jacobians = np.repeat(1.0, 6) - - # FIXME: should be computed directly from the metrics and jacobians - self.frozen_metrics = np.mean(self.metrics[2:4], axis=0) - self.frozen_jacobian = np.mean(self.jacobians[2:4], axis=0) - self.combined_frozen_metrics = 1.0 - - # FIXME: Move array_from_string stuff outside FluxDataSingle - # -- 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) - self.states = u.expand_to_n(self.state_pair, 6) - - # FIXME: these should be generalized fluxes - # FIXME: make a clear distinction between fluxes in physical and - # generalized coordinates - self.flux_pair = ref.pointwise_fluxes( - self.state_pair)[:,:,self.dir_internal].T.copy(order="F") - self.fluxes = ref.pointwise_fluxes( - self.states)[:,:,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_metrics) - self.weno_weights_neg = ref.weno_weights( - self.oscillation_neg, self.combined_frozen_metrics) - - self.consistent = ref.consistent_part(self.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) - - 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+i) % 3 + 1 for i in range(3)] - -# }}} - -# {{{ FluxDataVector - -# 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.states = states + self.metrics, self.jacobians = window_metrics + + self.fluxes = ref.pointwise_fluxes(self.states) + self.generalized_fluxes = np.array([np.dot(flux, metrics/jacobian) + for flux, metrics, jacobian + in zip(self.fluxes, self.metrics, self.jacobians)] + )[:,:,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) + + return states + + +@pytest.fixture(scope="session") +def metrics_generator(): + ndim = ref.gas.ndim + + def window_metrics(metric_str): + 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) + else: + raise ValueError("Metric string {} not supported".format(metric_str)) + return metrics, jacobians + + return window_metrics - 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") - - state_pair = self.swap_array_rows( - u.transposed_array_from_string(states_str), self.direction) - # FIXME: Move array_from_string stuff outside FluxDataSingle - # -- 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) - -# }}} @pytest.fixture(scope="session", params=[ - ("1 1 1 1 5.5,1 1 1 1 5.5", "x"), - ("1 1 1 1 5.5,1 1 1 1 5.5", "y"), - ("1 1 1 1 5.5,1 1 1 1 5.5", "z"), - ("2 4 4 4 20,1 1 1 1 5.5", "x"), - ("2 4 4 4 20,1 1 1 1 5.5", "y"), - ("2 4 4 4 20,1 1 1 1 5.5", "z"), - ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "x"), - ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "y"), - ("1 -1 -1 -1 5.5,2 -4 -4 -4 20", "z"), - ("2 4 8 12 64,1 1 2 3 11", "x"), - ("2 8 12 4 64,1 2 3 1 11", "y"), - ("2 12 4 8 64,1 3 1 2 11", "z"), - ("1 -1 -2 -3 11,2 -4 -8 -12 64", "x"), - ("1 -2 -3 -1 11,2 -8 -12 -4 64", "y"), - ("1 -3 -1 -2 11,2 -12 -4 -8 64", "z") + ("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 flux_test_data_fixture(request, queue): - return FluxDataSingle(queue, *request.param) - - -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") +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)) -@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] + +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) + + +@pytest.fixture(scope="session") +def pair_data(window_data): + return PairData(window_data) # vim: foldmethod=marker diff --git a/test/ideal_gas.py b/test/ideal_gas.py index a062acb14260e027962b5d6bbb96b200233f8764..85aa4b737a985ef7c0b83d1cec7fa9c9ca59d41a 100644 --- a/test/ideal_gas.py +++ b/test/ideal_gas.py @@ -29,8 +29,9 @@ import pyopencl.clrandom # noqa import loopy as lp # noqa -ndim = 3 -nvars = 3 +dirs = {"x": 1, "y": 2, "z": 3} +ndim = len(dirs) +nvars = ndim + 2 gamma = 1.4 @@ -46,13 +47,12 @@ def energy(state): return state[ndim+1] -def kinetic_energy(state): - vel = velocity(state) - return 0.5*rho(state)*np.dot(vel, vel) +def kinetic_energy(rho, vel): + return 0.5*rho*np.dot(vel, vel) def internal_energy(state): - return energy(state) - kinetic_energy(state) + return energy(state) - kinetic_energy(rho(state), velocity(state)) def pressure(state): @@ -74,3 +74,14 @@ def flux(state): result[ndim+1,i] += vel[i]*p return result + + +def conserved(rho, vel, p): + def energy(rho, vel, p): + return p/(gamma - 1) + kinetic_energy(rho, vel) + + state = np.empty(nvars) + state[0] = rho + state[1:nvars-1] = rho*vel + state[nvars-1] = energy(rho, vel, p) + return state diff --git a/test/test_eigensystem.py b/test/test_eigensystem.py index bc7df0a0e44a89b1cf3451576be91ed1c4ecf23a..8c1ba9369290ae16a5245a886c6c64af49874b01 100644 --- a/test/test_eigensystem.py +++ b/test/test_eigensystem.py @@ -34,26 +34,13 @@ import logging import pytest import utilities as u -from data_for_test import ( # noqa: F401 - flux_test_data_fixture - ) +import weno_reference_implementation as ref +from data_for_test import ( + pair_data, window_data, states_generator, metrics_generator) -def test_pointwise_eigenvalues_ideal_gas(queue, flux_test_data_fixture): - data = flux_test_data_fixture - - prg = u.get_weno_program_with_root_kernel("pointwise_eigenvalues") - - lam_dev = u.empty_array_on_device(queue, data.nvars, 6) - - prg(queue, nvars=data.nvars, d=data.direction, - states=data.states, lambda_pointwise=lam_dev) - - u.compare_arrays(lam_dev.get(), data.lam_pointwise) - - -def test_roe_uniform_grid_ideal_gas(queue, flux_test_data_fixture): - data = flux_test_data_fixture +def test_roe_uniform_grid_ideal_gas(queue, pair_data): + data = pair_data def check_roe_identity(states, R, R_inv): d_state = states[:, 1] - states[:, 0] @@ -82,7 +69,7 @@ def test_roe_uniform_grid_ideal_gas(queue, flux_test_data_fixture): lam = lam_dev.get() check_roe_identity(data.state_pair, R, R_inv) - check_roe_property(data.state_pair, data.flux_pair, R, R_inv, lam) + check_roe_property(data.state_pair, data.generalized_flux_pair, R, R_inv, lam) @pytest.mark.parametrize("lam_pointwise_str,lam_roe_str,lam_expected_str", [ diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 33b84303b3af166d13feb0970890d9b7c5e2c63f..2c4ec9ca439f3559c9f7670ec14b48c0828a0d70 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -34,9 +34,126 @@ import logging import pytest import utilities as u -from data_for_test import ( # noqa: F401 - cfd_test_data_fixture - ) +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") + + state_pair = self.swap_array_rows( + u.transposed_array_from_string(states_str), self.direction) + # 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): diff --git a/test/test_weno_flux.py b/test/test_flux_window_ops.py similarity index 55% rename from test/test_weno_flux.py rename to test/test_flux_window_ops.py index 40f9b2f5636e0ba833b7c257a4bcf7f84b6b3cb4..4bdb27ec925db6dda7cb451d36c812c528b3a4f4 100644 --- a/test/test_weno_flux.py +++ b/test/test_flux_window_ops.py @@ -34,45 +34,137 @@ import logging import pytest import utilities as u -from data_for_test import ( # noqa: F401 - flux_test_data_fixture - ) +import weno_reference_implementation as ref +from data_for_test import ( + PairData, window_data, states_generator, metrics_generator) -def test_weno_flux_uniform_grid(queue, flux_test_data_fixture): - data = flux_test_data_fixture +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_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): + return WindowResults(queue, window_data) + + +def test_pointwise_eigenvalues(queue, window_results): + data = window_results + + prg = u.get_weno_program_with_root_kernel("pointwise_eigenvalues") + + lam_dev = u.empty_array_on_device(queue, data.nvars, 6) + + prg(queue, nvars=data.nvars, d=data.direction, + states=data.states, lambda_pointwise=lam_dev) + + u.compare_arrays(lam_dev.get(), data.lam_pointwise) + + +def test_flux_splitting(queue, window_results): + data = window_results + + prg = u.get_weno_program_with_root_kernel("split_characteristic_fluxes") + + fluxes_pos_dev = u.empty_array_on_device(queue, data.nvars, 6) + fluxes_neg_dev = u.empty_array_on_device(queue, data.nvars, 6) + + prg(queue, nvars=data.nvars, + generalized_states_frozen=data.states, + generalized_fluxes_frozen=data.generalized_fluxes_frozen, + R_inv=data.R_inv, + wavespeeds=data.wavespeeds, + characteristic_fluxes_pos=fluxes_pos_dev, + characteristic_fluxes_neg=fluxes_neg_dev) + + u.compare_arrays(fluxes_pos_dev.get(), data.char_fluxes_pos) + u.compare_arrays(fluxes_neg_dev.get(), data.char_fluxes_neg) + + +def test_weno_flux(queue, window_results): + data = window_results prg = u.get_weno_program_with_root_kernel("weno_flux") flux_dev = u.empty_array_on_device(queue, data.nvars) prg(queue, nvars=data.nvars, - generalized_fluxes=data.fluxes, + generalized_fluxes=data.generalized_fluxes, characteristic_fluxes_pos=data.char_fluxes_pos, characteristic_fluxes_neg=data.char_fluxes_neg, - combined_frozen_metrics=data.combined_frozen_metrics, + combined_frozen_metric=data.combined_frozen_metric, R=data.R, flux=flux_dev) u.compare_arrays(flux_dev.get(), data.weno_flux) -def test_consistent_part_uniform_grid(queue, flux_test_data_fixture): - data = flux_test_data_fixture +@pytest.mark.slow +def test_consistent_part(queue, window_results): + data = window_results prg = u.get_weno_program_with_root_kernel("consistent_part") consistent_dev = u.empty_array_on_device(queue, data.nvars) prg(queue, nvars=data.nvars, - generalized_fluxes=data.fluxes, + generalized_fluxes=data.generalized_fluxes, consistent=consistent_dev) u.compare_arrays(consistent_dev.get(), data.consistent) -def test_dissipation_part_pos_uniform_grid(queue, flux_test_data_fixture): - data = flux_test_data_fixture +@pytest.mark.slow +def test_dissipation_part_pos(queue, window_results): + data = window_results prg = u.get_weno_program_with_root_kernel("dissipation_part_pos") @@ -80,15 +172,16 @@ def test_dissipation_part_pos_uniform_grid(queue, flux_test_data_fixture): prg(queue, nvars=data.nvars, characteristic_fluxes=data.char_fluxes_pos, - combined_frozen_metrics=data.combined_frozen_metrics, + combined_frozen_metric=data.combined_frozen_metric, R=data.R, dissipation_pos=dissipation_dev) u.compare_arrays(dissipation_dev.get(), data.dissipation_pos) -def test_dissipation_part_neg_uniform_grid(queue, flux_test_data_fixture): - data = flux_test_data_fixture +@pytest.mark.slow +def test_dissipation_part_neg(queue, window_results): + data = window_results prg = u.get_weno_program_with_root_kernel("dissipation_part_neg") @@ -96,15 +189,16 @@ def test_dissipation_part_neg_uniform_grid(queue, flux_test_data_fixture): prg(queue, nvars=data.nvars, characteristic_fluxes=data.char_fluxes_neg, - combined_frozen_metrics=data.combined_frozen_metrics, + combined_frozen_metric=data.combined_frozen_metric, R=data.R, dissipation_neg=dissipation_dev) u.compare_arrays(dissipation_dev.get(), data.dissipation_neg) -def test_weno_weights_pos_uniform_grid(queue, flux_test_data_fixture): - data = flux_test_data_fixture +@pytest.mark.slow +def test_weno_weights_pos(queue, window_results): + data = window_results prg = u.get_weno_program_with_root_kernel("weno_weights_pos") @@ -112,7 +206,7 @@ def test_weno_weights_pos_uniform_grid(queue, flux_test_data_fixture): prg(queue, nvars=data.nvars, characteristic_fluxes=data.char_fluxes_pos, - combined_frozen_metrics=data.combined_frozen_metrics, + combined_frozen_metric=data.combined_frozen_metric, w=weights_dev) sum_weights = np.sum(weights_dev.get(), axis=1) @@ -121,8 +215,9 @@ def test_weno_weights_pos_uniform_grid(queue, flux_test_data_fixture): u.compare_arrays(weights_dev.get(), data.weno_weights_pos) -def test_weno_weights_neg_uniform_grid(queue, flux_test_data_fixture): - data = flux_test_data_fixture +@pytest.mark.slow +def test_weno_weights_neg(queue, window_results): + data = window_results prg = u.get_weno_program_with_root_kernel("weno_weights_neg") @@ -130,7 +225,7 @@ def test_weno_weights_neg_uniform_grid(queue, flux_test_data_fixture): prg(queue, nvars=data.nvars, characteristic_fluxes=data.char_fluxes_neg, - combined_frozen_metrics=data.combined_frozen_metrics, + combined_frozen_metric=data.combined_frozen_metric, w=weights_dev) sum_weights = np.sum(weights_dev.get(), axis=1) @@ -139,8 +234,9 @@ def test_weno_weights_neg_uniform_grid(queue, flux_test_data_fixture): u.compare_arrays(weights_dev.get(), data.weno_weights_neg) -def test_oscillation_pos_uniform_grid(queue, flux_test_data_fixture): - data = flux_test_data_fixture +@pytest.mark.slow +def test_oscillation_pos(queue, window_results): + data = window_results prg = u.get_weno_program_with_root_kernel("oscillation_pos") @@ -153,8 +249,9 @@ def test_oscillation_pos_uniform_grid(queue, flux_test_data_fixture): u.compare_arrays(oscillation_dev.get(), data.oscillation_pos) -def test_oscillation_neg_uniform_grid(queue, flux_test_data_fixture): - data = flux_test_data_fixture +@pytest.mark.slow +def test_oscillation_neg(queue, window_results): + data = window_results prg = u.get_weno_program_with_root_kernel("oscillation_neg") @@ -167,26 +264,6 @@ def test_oscillation_neg_uniform_grid(queue, flux_test_data_fixture): u.compare_arrays(oscillation_dev.get(), data.oscillation_neg) -def test_flux_splitting_uniform_grid(queue, flux_test_data_fixture): - data = flux_test_data_fixture - - prg = u.get_weno_program_with_root_kernel("split_characteristic_fluxes") - - fluxes_pos_dev = u.empty_array_on_device(queue, data.nvars, 6) - fluxes_neg_dev = u.empty_array_on_device(queue, data.nvars, 6) - - prg(queue, nvars=data.nvars, - generalized_states_frozen=data.states, - generalized_fluxes_frozen=data.fluxes, - R_inv=data.R_inv, - wavespeeds=data.wavespeeds, - characteristic_fluxes_pos=fluxes_pos_dev, - characteristic_fluxes_neg=fluxes_neg_dev) - - u.compare_arrays(fluxes_pos_dev.get(), data.char_fluxes_pos) - u.compare_arrays(fluxes_neg_dev.get(), data.char_fluxes_neg) - - # This lets you run 'python test.py test_case(cl._csc)' without pytest. if __name__ == "__main__": if len(sys.argv) > 1: diff --git a/test/weno_reference_implementation.py b/test/weno_reference_implementation.py index bd4f7bcf7f62521eb9229a50a48ca262c6b69900..a7b4095acb3a91b16ad0d5eeb6e73f9db7366d38 100644 --- a/test/weno_reference_implementation.py +++ b/test/weno_reference_implementation.py @@ -31,6 +31,7 @@ import loopy as lp # noqa import utilities as u import ideal_gas as gas +### WENO algorithm taken from Nonomura, et al., Comput. Fluids 107 (2015) 242-255 def pointwise_fluxes(states): return np.array([gas.flux(s) for s in states.T]) @@ -124,18 +125,16 @@ def flux_differences(fluxes): return u.transposed_array([np.dot(w, f) for f in fluxes.T[indices]]) -def combination_weighting(w): - return np.array( - [20*w[:,0] - 1, -10*(w[:,0] + w[:,1]) + 5, np.ones(w.shape[0])] - ).T - - -def combine_fluxes(w, f): - cw = combination_weighting(w) - return np.sum(cw*f, axis=1) +def dissipation_part(R, char_fluxes, w, sign): + def combine_fluxes(w, f): + def combination_weighting(w): + return np.array( + [20*w[:,0] - 1, -10*(w[:,0] + w[:,1]) + 5, np.ones(w.shape[0])] + ).T + cw = combination_weighting(w) + return np.sum(cw*f, axis=1) -def dissipation_part(R, char_fluxes, w, sign): flux_diff = flux_differences(char_fluxes)[:,::sign] flux_comb = combine_fluxes(w, flux_diff)