diff --git a/WENO.F90 b/WENO.F90 index 829080e5f56556ffcc276e33a5ad8d0aa0628ab6..f829dfe057a663f93ef02f815ce7c22737072dc5 100644 --- a/WENO.F90 +++ b/WENO.F90 @@ -38,7 +38,7 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, & real metric_solve_tmp(ndim) real R(nvars, nvars), R_inv(nvars, nvars), lambda_roe(nvars) real lambda_pointwise(nvars, -2:3) - real lambda(nvars) + real wavespeeds(nvars) integer i, j, k real delta_xi, delta_eta, delta_zeta real characteristic_fluxes_pos(nvars, -2:3) @@ -89,13 +89,13 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, & 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, lambda) + call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) call split_characteristic_fluxes(nvars, & generalized_states_frozen, & generalized_fluxes_frozen(:, 1, :), & R_inv, & - lambda, & + wavespeeds, & characteristic_fluxes_pos, & characteristic_fluxes_neg) @@ -139,13 +139,13 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, & 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, lambda) + call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) call split_characteristic_fluxes(nvars, & generalized_states_frozen, & generalized_fluxes_frozen(:, 2, :), & R_inv, & - lambda, & + wavespeeds, & characteristic_fluxes_pos, & characteristic_fluxes_neg) @@ -190,13 +190,13 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, & 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, lambda) + call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) call split_characteristic_fluxes(nvars, & generalized_states_frozen, & generalized_fluxes_frozen(:, 3, :), & R_inv, & - lambda, & + wavespeeds, & characteristic_fluxes_pos, & characteristic_fluxes_neg) @@ -437,13 +437,13 @@ subroutine roe_eigensystem(nvars, ndim, d, states, metrics_frozen, R, R_inv, lam R_inv(5,5) = beta*b1 end subroutine -subroutine lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, lambda) +subroutine lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) implicit none integer, intent(in) :: nvars real, intent(in) :: lambda_pointwise(nvars, -2:3) real, intent(in) :: lambda_roe(nvars) - real, intent(out) :: lambda(nvars) + real, intent(out) :: wavespeeds(nvars) real kappa integer v @@ -452,11 +452,11 @@ subroutine lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, lambda) kappa = 1.1 do v=1,nvars - lambda(v) = abs(lambda_roe(v)) + wavespeeds(v) = abs(lambda_roe(v)) do k=-2,3 - lambda(v) = max(lambda(v), abs(lambda_pointwise(v,k))) + wavespeeds(v) = max(wavespeeds(v), abs(lambda_pointwise(v,k))) end do - lambda(v) = kappa*lambda(v) + wavespeeds(v) = kappa*wavespeeds(v) end do end subroutine @@ -615,7 +615,7 @@ subroutine split_characteristic_fluxes(nvars, & generalized_states_frozen, & generalized_fluxes_frozen, & R_inv, & - lambda, & + wavespeeds, & characteristic_fluxes_pos, & characteristic_fluxes_neg) implicit none @@ -624,7 +624,7 @@ subroutine split_characteristic_fluxes(nvars, & real, intent(in) :: generalized_states_frozen(nvars, -2:3) real, intent(in) :: generalized_fluxes_frozen(nvars, -2:3) real, intent(in) :: R_inv(nvars, nvars) - real, intent(in) :: lambda(nvars) + real, intent(in) :: wavespeeds(nvars) real, intent(out) :: characteristic_fluxes_pos(nvars, -2:3) real, intent(out) :: characteristic_fluxes_neg(nvars, -2:3) @@ -638,10 +638,10 @@ subroutine split_characteristic_fluxes(nvars, & do l=1,nvars characteristic_fluxes_pos(m,k) = characteristic_fluxes_pos(m,k) & + 0.5*R_inv(m,l) & - *(generalized_fluxes_frozen(l,k) + lambda(m)*generalized_states_frozen(l,k)) + *(generalized_fluxes_frozen(l,k) + wavespeeds(m)*generalized_states_frozen(l,k)) characteristic_fluxes_neg(m,k) = characteristic_fluxes_neg(m,k) & + 0.5*R_inv(m,l) & - *(generalized_fluxes_frozen(l,k) - lambda(m)*generalized_states_frozen(l,k)) + *(generalized_fluxes_frozen(l,k) - wavespeeds(m)*generalized_states_frozen(l,k)) end do end do end do diff --git a/test.py b/test.py index ce10501781aa85aedaacc9d6d64f06e395e9bbe5..1bfe9c86d3ee5151c84dbb7066f3495a7ca9929b 100644 --- a/test.py +++ b/test.py @@ -18,6 +18,31 @@ from pyopencl.tools import ( # noqa from utilities import * +@pytest.mark.parametrize("lam_pointwise_str,lam_roe_str,lam_expected_str", [ + ("1 2 3 4 5,2 4 6 8 10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"), + ("1 2 3 4 5,-2 -4 -6 -8 -10", "1.5 3 4.5 6 7.5", "2.2 4.4 6.6 8.8 11"), + ("1 2 3 4 5,-2 -4 -6 -8 -10", "3 6 9 12 15", "3.3 6.6 9.9 13.2 16.5"), + ("1 2 3 4 5,2 4 6 8 10", "-3 -6 -9 -12 -15", "3.3 6.6 9.9 13.2 16.5"), + ("3 2 9 4 5,2 6 6 12 10", "-1 -4 -3 -8 -15", "3.3 6.6 9.9 13.2 16.5") + ]) +def test_lax_wavespeeds( + ctx_factory, lam_pointwise_str, lam_roe_str, lam_expected_str): + prg = get_weno_program_with_root_kernel("lax_wavespeeds") + queue = get_queue(ctx_factory) + + nvars = 5 + + lam_pointwise = expand_to_6(transposed_array_from_string(lam_pointwise_str)) + lam_roe = array_from_string(lam_roe_str) + lam_dev = empty_array_on_device(queue, nvars) + + prg(queue, nvars=nvars, lambda_pointwise=lam_pointwise, + lambda_roe=lam_roe, wavespeeds=lam_dev) + + lam_expected = array_from_string(lam_expected_str) + compare_arrays(lam_dev.get(), lam_expected) + + @pytest.mark.parametrize("states_str,direction,lam_expected_str", [ ("2 4 4 4 20,1 1 1 1 5.5", "x", "2 2 2 3.49666295 0.503337045,1 1 1 2.49666295 -0.496662955"), @@ -46,9 +71,6 @@ from utilities import * ]) def test_pointwise_eigenvalues( ctx_factory, states_str, direction, lam_expected_str): - def expand_to_6(pair): - return np.repeat(pair, 3, axis=1).copy(order="F") - prg = get_weno_program_with_root_kernel("pointwise_eigenvalues") queue = get_queue(ctx_factory) diff --git a/utilities.py b/utilities.py index b825fcbe82ee7ef0d2a991edf9ed70f136b55f80..c619310ed0ff59c260744f5b08a8590f5997a7d6 100644 --- a/utilities.py +++ b/utilities.py @@ -64,6 +64,10 @@ def array_from_string(string_array): def split_map_to_list(string, map_func, splitter): return list(map(map_func, string.split(splitter))) + +def expand_to_6(pair): + return np.repeat(pair, 3, axis=1).copy(order="F") + # }}} # {{{ device