diff --git a/WENO.F90 b/WENO.F90 index 5a7a6f44e3d9594df952764cf52d7019192c7875..80780a4e6dbb90344c0074ebe0c8f64abb3f8505 100644 --- a/WENO.F90 +++ b/WENO.F90 @@ -61,7 +61,8 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, & generalized_states_frozen, & generalized_fluxes_frozen) - call pointwise_eigenvalues(nvars, 1, states(:, i-2:i+3, j, k), lambda_pointwise) + 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) @@ -111,7 +112,8 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, & generalized_states_frozen, & generalized_fluxes_frozen) - call pointwise_eigenvalues(nvars, 2, states(:, i, j-2:j+3, k), lambda_pointwise) + 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) @@ -162,7 +164,8 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, & generalized_states_frozen, & generalized_fluxes_frozen) - call pointwise_eigenvalues(nvars, 3, states(:, i, j, k-2:k+3), lambda_pointwise) + 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) @@ -199,32 +202,37 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, & !$loopy end tagged: flux_z_diff end subroutine -subroutine pointwise_eigenvalues(nvars, d, states, lambda_pointwise) +subroutine pointwise_eigenvalues(nvars, ndim, d, states, metrics, lambda_pointwise) implicit none integer, intent(in) :: nvars + integer, intent(in) :: ndim integer, intent(in) :: d real*8, intent(in) :: states(nvars, -2:3) + real*8, intent(in) :: metrics(ndim, ndim, -2:3) real*8, intent(out) :: lambda_pointwise(nvars, -2:3) - real*8 u(3), c, p, rho, ke + real*8 u(3), c, p, rho, ke, metric_norm integer v, k, i do k=-2,3 rho = states(1,k) ke = 0.0d0 + metric_norm = 0.0d0 do i=1,3 u(i) = states(i+1,k)/rho ke = ke + u(i)**2 + metric_norm = metric_norm + metrics(d,i,k)**2 end do + metric_norm = sqrt(metric_norm) p = (1.4d0 - 1)*(states(nvars,k) - 0.5d0*rho*ke) c = sqrt(1.4d0*p/rho) do v=1,nvars-2 lambda_pointwise(v,k) = u(d) end do - lambda_pointwise(nvars-1,k) = u(d) + c - lambda_pointwise(nvars,k) = u(d) - c + lambda_pointwise(nvars-1,k) = u(d) + c*metric_norm + lambda_pointwise(nvars,k) = u(d) - c*metric_norm end do end subroutine diff --git a/test/data_for_test.py b/test/data_for_test.py index 483d50b65f311b8cabea7497b1b1626590b62032..99feb032eb6658eac88027473c7cada9928018c6 100644 --- a/test/data_for_test.py +++ b/test/data_for_test.py @@ -1,4 +1,4 @@ -__copyright__ = "Copyright (C) 2019 Timothy A. Smith" +__copyright__ = "Copyright (C) 2019-2020 Timothy A. Smith" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -24,6 +24,7 @@ import numpy as np import pytest import utilities as u import weno_reference_implementation as ref +import wavy_metrics_3d as wavy3d class WindowData: @@ -66,12 +67,15 @@ def states_generator(): def metrics_generator(): ndim = ref.gas.ndim - def window_metrics(metric_str): + 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 @@ -80,6 +84,9 @@ def metrics_generator(): @pytest.fixture(scope="session", params=[ + ("sine", "x", "wavy3d"), + ("sine", "y", "wavy3d"), + ("sine", "z", "wavy3d"), ("sine", "x", "uniform"), ("sine", "y", "uniform"), ("sine", "z", "uniform"), @@ -106,7 +113,7 @@ def window_data(request, states_generator, metrics_generator): return WindowData(states_generator(states_str), dir_str, - metrics_generator(metric_str)) + metrics_generator(metric_str, dir_str)) class PairData: diff --git a/test/test_flux_window_ops.py b/test/test_flux_window_ops.py index 4bdb27ec925db6dda7cb451d36c812c528b3a4f4..8c22c1370cd8a2c15c9d9117430eae9565e2b9c2 100644 --- a/test/test_flux_window_ops.py +++ b/test/test_flux_window_ops.py @@ -58,6 +58,8 @@ class WindowResults: 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") @@ -101,10 +103,12 @@ def test_pointwise_eigenvalues(queue, window_results): prg = u.get_weno_program_with_root_kernel("pointwise_eigenvalues") + metric_arg = np.moveaxis(data.metrics, 0, -1).copy(order="F") 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) + prg(queue, nvars=data.nvars, ndim=data.ndim, d=data.direction, + states=data.states, metrics=metric_arg, + lambda_pointwise=lam_dev) u.compare_arrays(lam_dev.get(), data.lam_pointwise) @@ -118,7 +122,7 @@ def test_flux_splitting(queue, window_results): fluxes_neg_dev = u.empty_array_on_device(queue, data.nvars, 6) prg(queue, nvars=data.nvars, - generalized_states_frozen=data.states, + generalized_states_frozen=data.generalized_states_frozen, generalized_fluxes_frozen=data.generalized_fluxes_frozen, R_inv=data.R_inv, wavespeeds=data.wavespeeds, @@ -199,6 +203,8 @@ def test_dissipation_part_neg(queue, window_results): @pytest.mark.slow def test_weno_weights_pos(queue, window_results): data = window_results + if not np.array_equal(data.metrics[0], data.metrics[1]): + pytest.xfail("expect failure with non-uniform metrics") prg = u.get_weno_program_with_root_kernel("weno_weights_pos") @@ -218,6 +224,8 @@ def test_weno_weights_pos(queue, window_results): @pytest.mark.slow def test_weno_weights_neg(queue, window_results): data = window_results + if not np.array_equal(data.metrics[0], data.metrics[1]): + pytest.xfail("expect failure with non-uniform metrics") prg = u.get_weno_program_with_root_kernel("weno_weights_neg") diff --git a/test/wavy_metrics_3d.py b/test/wavy_metrics_3d.py new file mode 100644 index 0000000000000000000000000000000000000000..adfdf0cb059bc5a24f3cc413ab5f01d62dafc458 --- /dev/null +++ b/test/wavy_metrics_3d.py @@ -0,0 +1,144 @@ +__copyright__ = "Copyright (C) 2020 Timothy A. Smith" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +# The wavy grid and curvilinear metric formulation used here was taken from +# T. Nonomura, et al., Comput. Fluids 107 (2015) 242-255 + +import numpy as np +from numpy import pi + +# grid parameters +Lx = 4 +Ly = 4 +Lz = 4 + +Lprod = Lx*Ly*Lz + +# grid indices +# built-in assumption is j goes from 1 to j_max, and similarly with k, l +j_max = 21 +k_max = 21 +l_max = 21 + +dx = Lx/(j_max - 1) +dy = Ly/(k_max - 1) +dz = Lz/(l_max - 1) + + +def cos_x(xi): + return np.cos((4*pi*dx*xi)/Lx) + +def cos_y(eta): + return np.cos((4*pi*dy*eta)/Ly) + +def cos_z(zeta): + return np.cos((4*pi*dz*zeta)/Lz) + +def sin_x(xi): + return np.sin((4*pi*dx*xi)/Lx) + +def sin_y(eta): + return np.sin((4*pi*dy*eta)/Ly) + +def sin_z(zeta): + return np.sin((4*pi*dz*zeta)/Lz) + + +def blue_x(xi, eta, zeta): + return 16*(pi**2)*Lx*dz*dy*cos_y(eta)*cos_z(zeta)*sin_x(xi)**2 + +def blue_y(xi, eta, zeta): + return 16*(pi**2)*Ly*dx*dz*cos_z(zeta)*cos_x(xi)*sin_y(eta)**2 + +def blue_z(xi, eta, zeta): + return 16*(pi**2)*Lz*dx*dy*cos_y(eta)*cos_x(xi)*sin_z(zeta)**2 + + +def coords_from_indices(j, k, l): + return (j-1, k-1, l-1) + + +def weighted_metric_mat(j, k, l): + """ returns weighted metric matrix at a single grid point with indices j,k,l""" + xi, eta, zeta = coords_from_indices(j, k, l) + arr = np.zeros([3,3]) + + arr[0,0] = dy*dz*(Lprod - blue_x(xi, eta, zeta))/Lprod + arr[1,0] = 4*pi*dx*dy*dz*cos_x(xi)*(-Lz*sin_z(zeta) + + 4*pi*dz*cos_z(zeta)*sin_y(eta)*sin_x(xi))/(Lx*Lz) + arr[2,0] = 4*pi*dx*dy*dz*cos_x(xi)*(-Ly*sin_y(eta) + + 4*pi*dy*cos_y(eta)*sin_z(zeta)*sin_x(xi))/(Lx*Ly) + + arr[0,1] = 4*pi*dx*dy*dz*cos_y(eta)*(-Lz*sin_z(zeta) + + 4*pi*dz*cos_z(zeta)*sin_x(xi)*sin_y(eta))/(Ly*Lz) + arr[1,1] = dx*dz*(Lprod - blue_y(xi, eta, zeta))/Lprod + arr[2,1] = 4*pi*dx*dy*dz*cos_y(eta)*(-Lx*sin_x(xi) + + 4*pi*dx*cos_x(xi)*sin_z(zeta)*sin_y(eta))/(Ly*Lx) + + arr[0,2] = 4*pi*dx*dy*dz*cos_z(zeta)*(-Ly*sin_y(eta) + + 4*pi*dy*cos_y(eta)*sin_x(xi)*sin_z(zeta))/(Lz*Ly) + arr[1,2] = 4*pi*dx*dy*dz*cos_z(zeta)*(-Lx*sin_x(xi) + + 4*pi*dx*cos_x(xi)*sin_y(eta)*sin_z(zeta))/(Lz*Lx) + arr[2,2] = dx*dy*(Lprod - blue_z(xi, eta, zeta))/Lprod + + return arr + + +def inverse_jacobian(j, k, l): + """ returns inverse jacobian at a single grid point with indices j,k,l""" + xi, eta, zeta = coords_from_indices(j, k, l) + return dx*dy*dz*( + Lprod + - blue_x(xi, eta, zeta) + - blue_y(xi, eta, zeta) + - blue_z(xi, eta, zeta) + +16*(pi**3)*dx*dy*dz*sin_x(2*xi)*sin_y(2*eta)*sin_z(2*zeta) + )/Lprod + + +def array_along_dim(func, dim): + if dim == "x": + return np.array([func(j+1,1,1) for j in range(6)], + dtype=np.float64, + order="F") + elif dim == "y": + return np.array([func(1,k+1,1) for k in range(6)], + dtype=np.float64, + order="F") + elif dim == "z": + return np.array([func(1,1,l+1) for l in range(6)], + dtype=np.float64, + order="F") + else: + raise ValueError("Dimension {} not supported".format(dim)) + + +def weighted_metrics(dim): + return array_along_dim(weighted_metric_mat, dim) + + +def jacobians(dim): + return np.reciprocal(array_along_dim(inverse_jacobian, dim)) + + +def metrics(dim): + return weighted_metrics(dim)*jacobians(dim)[:,None,None]