diff --git a/README.md b/README.md index 21420b4f8f81203d4ca1bf03c2886999a88f526e..72b75ea8a39af4a7d4bcceb533bc1be541e4f0c9 100644 --- a/README.md +++ b/README.md @@ -6,15 +6,13 @@ Standalone Fortran WENO code with Loopy transformations. * The WENO scheme is taken from [1] * The Roe matrix implementation follows the internal PlasCom2 documentation -* Boundary conditions are not really handled correctly at all +* Boundary conditions are not handled by the main flux derivative kernel * A declaration of "real" is not intended to suggest single precision * This implementation does not include the severe shock correction proposed by [1] * Current assumption for metrics is that metrics(i,j) = (\xi_i)_(x_j) -* It's not immediately obvious, but there might be some wasted work in the - frozen generalized flux computation -* Right now the specific gas constant gamma is implemented as a magic number - (1.4) -- this should be generalized +* Right now an ideal gas formulation is assumed and the specific gas constant + gamma is implemented as a magic number (1.4) ## Fortran dialect wishlist @@ -25,25 +23,6 @@ Standalone Fortran WENO code with Loopy transformations. - (e.g. 2:0 instead of 2:0:-1) - Scalars with intent(out) -Note: Tag 'nice' references the code at a point where several of these items are used - -## Testing plan - -Plan to test the following routines directly: - -- `compute_flux_derivatives` -- `convert_to_generalized` -- `convert_to_generalized_frozen` -- `convert_from_generalized` -- `roe_eigensystem` -- `pointwise_eigenvalues` -- `lax_wavespeeds` -- `split_characteristic_fluxes` -- `weno_flux` -- `consistent_part` -- `dissipation_part_pos` -- `dissipation_part_neg` - ## References [1] Nonomura, T., Terakado, D., Abe, Y., & Fujii, K. (2015). diff --git a/WENO.F90 b/WENO.F90 index 5509694bdf0340a3d92dcd03de17e0ada875209c..fa0d0a4559a84bcb465e291d9e0a4af591659735 100644 --- a/WENO.F90 +++ b/WENO.F90 @@ -52,9 +52,8 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, & do k=1,nz do j=1,ny do i=0,nx - call convert_to_generalized_frozen(nvars, ndim, & + call convert_to_generalized_frozen(nvars, ndim, d, & 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, & @@ -104,9 +103,8 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, & do k=1,nz do i=1,nx do j=0,ny - call convert_to_generalized_frozen(nvars, ndim, & + call convert_to_generalized_frozen(nvars, ndim, d, & 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, & @@ -157,9 +155,8 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, & do j=1,ny do i=1,nx do k=0,nz - call convert_to_generalized_frozen(nvars, ndim, & + call convert_to_generalized_frozen(nvars, ndim, d, & 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, & @@ -426,19 +423,46 @@ subroutine lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) end do end subroutine +subroutine compute_ideal_gas_fluxes(nvars, ndim, states, fluxes) + integer, intent(in) :: nvars, ndim + real*8, intent(in) :: states(nvars, -2:3) + real*8, intent(out) :: fluxes(nvars, ndim, -2:3) + + integer v, d, k + real*8 velocity(ndim) + real*8 kinetic_energy, pressure, gm + + gm = 1.4d0 + + do k=-2,3 + kinetic_energy = 0.0d0 + do v=1,ndim + velocity(v) = states(v+1,k)/states(1,k) + kinetic_energy = kinetic_energy + 0.5d0*states(1,k)*velocity(v)*velocity(v) + end do + pressure = (states(ndim+2,k) - kinetic_energy)*(gm - 1.0d0) + + do d=1,ndim + do v=1,nvars + fluxes(v,d,k) = states(v,k)*velocity(d) + end do + fluxes(d+1,d,k) = fluxes(d+1,d,k) + pressure + fluxes(ndim+2,d,k) = fluxes(ndim+2,d,k) + pressure*velocity(d) + end do + end do +end subroutine + subroutine convert_to_generalized_frozen( & - nvars, ndim, & + nvars, ndim, d, & states, & - fluxes, & metrics, & metric_jacobians, & metrics_frozen, & combined_frozen_metrics, & generalized_states_frozen, & generalized_fluxes_frozen) - integer, intent(in) :: nvars, ndim + integer, intent(in) :: nvars, ndim, d real*8, intent(in) :: states(nvars, -2:3) - real*8, intent(in) :: fluxes(nvars, -2:3) real*8, intent(in) :: metrics(ndim, ndim, 2) real*8, intent(in) :: metric_jacobians(2) real*8, intent(out) :: metrics_frozen(ndim, ndim) @@ -447,10 +471,13 @@ subroutine convert_to_generalized_frozen( & real*8, intent(out) :: generalized_fluxes_frozen(nvars, -2:3) real*8 jacobian_frozen + real*8 fluxes(nvars, ndim, -2:3) integer c, k, v integer i, j + call compute_ideal_gas_fluxes(nvars, ndim, states, fluxes) + do j=1,ndim do i=1,ndim metrics_frozen(i,j) = 0.0d0 @@ -476,10 +503,11 @@ subroutine convert_to_generalized_frozen( & do k=-2,3 do v=1,nvars - ! FIXME: this is a bug - ! ... we actually need to recompute the fluxes using frozen metrics - ! ... but this works when the metric matrix is identity - generalized_fluxes_frozen(v,k) = fluxes(v,k) + generalized_fluxes_frozen(v,k) = 0.0d0 + do i=1,ndim + generalized_fluxes_frozen(v,k) = generalized_fluxes_frozen(v,k) & + + metrics_frozen(d,i)*fluxes(v,i,k) + end do end do end do end subroutine @@ -625,12 +653,10 @@ subroutine weno_weights_pos(nvars, characteristic_fluxes, combined_frozen_metric call oscillation_pos(nvars, characteristic_fluxes, IS) do i=1,nvars - sum_alpha(1) = 0.0d0 do j=1,3 alpha(j) = C(j)/(IS(i, j) + eps)**p - sum_alpha(1) = sum_alpha(1) + alpha(j) end do - !call sum_array(3, alpha, sum_alpha) + call sum_array(3, alpha, sum_alpha) do j=1,3 w(i,j) = alpha(j)/sum_alpha(1) end do diff --git a/test/conftest.py b/test/conftest.py index 994bf9111094552282711469499ff1a6cd161602..db422088e07a742450b81da0cd73b76cab7c71ae 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -196,7 +196,7 @@ 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 + halo = int((nxhalo - nx)/2) states = np.empty( (ref.gas.nvars, nxhalo, nyhalo, nzhalo), dtype=np.float64, order="F") @@ -286,6 +286,17 @@ def grid_states_generator(): state_right = ref.gas.conserved(rho_right, vel_right, pressure_right) set_states_double_hump(state_base, state_left, state_right) + elif states_str == "sine": + def sin_func(i, n): + return np.sin(2*np.pi*(i - halo)/n) + + vel = np.array([1.0, 2.0, 3.0]) + pressure = 3.2 + for i in range(nxhalo): + for j in range(nyhalo): + for k in range(nzhalo): + rho = sin_func(i, nx)*sin_func(j, ny)*sin_func(k, nz) + 1.1 + states[:, i, j, k] = ref.gas.conserved(rho, vel, pressure) else: raise ValueError("States string {} not supported".format(states_str)) return states @@ -298,7 +309,10 @@ 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)) + nx, ny, nz = parse_grid_sizes(grid_sizes_str) + nxhalo, nyhalo, nzhalo = halo_sizes(nx, ny, nz) + halo = int((nxhalo - nx)/2) + if (metric_str == "uniform"): metrics = np.full( (ndim, ndim, nxhalo, nyhalo, nzhalo), @@ -311,6 +325,9 @@ def grid_metrics_generator(): jacobians = np.full( (nxhalo, nyhalo, nzhalo), 1.0, dtype=np.float64, order="F") + elif (metric_str == "wavy3d"): + metrics = wavy3d.grid_metrics(nx, ny, nz, halo) + jacobians = wavy3d.grid_jacobians(nx, ny, nz, halo) else: raise ValueError("Metric string {} not supported".format(metric_str)) return metrics, jacobians @@ -334,7 +351,9 @@ def grid_metrics_generator(): ("double_hump", "25x25x25", "x", "uniform"), ("double_hump", "25x25x25", "y", "uniform"), ("double_hump", "25x25x25", "z", "uniform"), - # FIXME: add remaining cases + ("sine", "10x10x10", "x", "wavy3d"), + ("sine", "10x10x10", "y", "wavy3d"), + ("sine", "10x10x10", "z", "wavy3d") ]) def grid_data(request, grid_states_generator, grid_metrics_generator): states_str = request.param[0] diff --git a/test/data_for_test.py b/test/data_for_test.py index 2cf9c2fea54857bec0cbc5dcac083c512960b620..dc4b6836b564c353473f6c5ca8589ba99898c8e0 100644 --- a/test/data_for_test.py +++ b/test/data_for_test.py @@ -1,4 +1,4 @@ -__copyright__ = "Copyright (C) 2019-2020 Timothy A. Smith" +__copyright__ = "Copyright (C) 2019 Timothy A. Smith" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -70,13 +70,15 @@ class WindowResults: self.states = data.states self.state_pair = pair.state_pair + self.fluxes = data.fluxes 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.combined_frozen_metrics = pair.combined_frozen_metrics + self.combined_frozen_metric = self.combined_frozen_metrics[self.dir_internal] self.generalized_states_frozen = self.states/self.frozen_jacobian diff --git a/test/test_frozen_metrics.py b/test/test_frozen_metrics.py new file mode 100644 index 0000000000000000000000000000000000000000..cf15153d0ecfae4971ececf881e33eaf24171b14 --- /dev/null +++ b/test/test_frozen_metrics.py @@ -0,0 +1,94 @@ +__copyright__ = "Copyright (C) 2019 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. +""" + +import numpy as np # noqa: F401 +import numpy.linalg as la # noqa: F401 +import pyopencl as cl # noqa: F401 +import pyopencl.array # noqa +import pyopencl.tools # noqa +import pyopencl.clrandom # noqa +import loopy as lp # noqa + +import sys +import logging + +import pytest + +import utilities as u + + +def test_compute_ideal_gas_fluxes(queue, window_results): + data = window_results + + prg = u.get_weno_program_with_root_kernel("compute_ideal_gas_fluxes") + + fluxes_dev = u.empty_array_on_device(queue, data.nvars, data.ndim, 6) + + prg(queue, + nvars=data.nvars, + ndim=data.ndim, + states=data.states, + fluxes=fluxes_dev) + + u.compare_arrays(fluxes_dev.get(), np.moveaxis(data.fluxes, 0, -1)) + + +def test_frozen_metrics_generalized_flux(queue, window_results): + data = window_results + + prg = u.get_weno_program_with_root_kernel("convert_to_generalized_frozen") + + pair_metrics = np.moveaxis(data.metrics[2:4], 0, -1).copy(order="F") + pair_jacobians = data.jacobians[2:4].copy() + + frozen_metrics_dev = u.empty_array_on_device(queue, data.ndim, data.ndim) + combined_frozen_metrics_dev = u.empty_array_on_device(queue, data.ndim) + generalized_states_frozen_dev = u.empty_array_on_device(queue, data.nvars, 6) + generalized_fluxes_frozen_dev = u.empty_array_on_device(queue, data.nvars, 6) + + prg(queue, + nvars=data.nvars, + ndim=data.ndim, + d=data.direction, + states=data.states, + metrics=pair_metrics, + metric_jacobians=pair_jacobians, + metrics_frozen=frozen_metrics_dev, + combined_frozen_metrics=combined_frozen_metrics_dev, + generalized_states_frozen=generalized_states_frozen_dev, + generalized_fluxes_frozen=generalized_fluxes_frozen_dev) + + u.compare_arrays(frozen_metrics_dev.get(), data.frozen_metrics) + u.compare_arrays(combined_frozen_metrics_dev.get(), data.combined_frozen_metrics) + u.compare_arrays(generalized_states_frozen_dev.get(), + data.generalized_states_frozen) + u.compare_arrays(generalized_fluxes_frozen_dev.get(), + data.generalized_fluxes_frozen) + + +# This lets you run 'python test.py test_case(cl._csc)' without pytest. +if __name__ == "__main__": + if len(sys.argv) > 1: + logging.basicConfig(level="INFO") + exec(sys.argv[1]) + else: + pytest.main([__file__]) diff --git a/test/wavy_metrics_3d.py b/test/wavy_metrics_3d.py index 93011a17d6ad8aadef288888c0e19340e10bee97..77a9c293043c738e046422f0f1fb81c458bbd96f 100644 --- a/test/wavy_metrics_3d.py +++ b/test/wavy_metrics_3d.py @@ -1,4 +1,4 @@ -__copyright__ = "Copyright (C) 2020 Timothy A. Smith" +__copyright__ = "Copyright (C) 2019 Timothy A. Smith" __license__ = """ Permission is hereby granted, free of charge, to any person obtaining a copy @@ -149,3 +149,26 @@ def jacobians(dim): def metrics(dim): return weighted_metrics(dim)*jacobians(dim)[:, None, None] + + +def grid(func, nx, ny, nz, halo): + return np.array([[[func(i, j, k) + for i in range(nx + 2*halo)] + for j in range(ny + 2*halo)] + for k in range(nz + 2*halo)], + dtype=np.float64, + order="F") + + +def grid_weighted_metrics(nx, ny, nz, halo): + return grid(weighted_metric_mat, nx, ny, nz, halo) + + +def grid_jacobians(nx, ny, nz, halo): + return np.reciprocal(grid(inverse_jacobian, nx, ny, nz, halo)) + + +def grid_metrics(nx, ny, nz, halo): + return ((np.moveaxis(grid_weighted_metrics(nx, ny, nz, halo), + [3, 4, 0, 1, 2], [0, 1, 2, 3, 4]) + * grid_jacobians(nx, ny, nz, halo)[None, None, :, :, :]).copy(order="F"))