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 2f12e9c0541307476073380e778dbac254695cac..9e62317061554bc5de818673987fcdb11b256121 100644 --- a/test/data_for_test.py +++ b/test/data_for_test.py @@ -31,7 +31,7 @@ class WindowData: ndim = ref.gas.ndim dirs = ref.gas.dirs - def __init__(self, states_str, direction): + def __init__(self, states_str, direction, window_metrics): self.direction = self.dirs[direction] self.dir_internal = self.direction-1 @@ -47,30 +47,74 @@ class WindowData: self.fluxes = ref.pointwise_fluxes( self.states)[:,:,self.dir_internal].T.copy(order="F") - 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) + self.metrics, self.jacobians = window_metrics + + +@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 @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") + ("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): - return WindowData(*request.param) +def window_data(request, metrics_generator): + states_str = request.param[0] + dir_str = request.param[1] + metric_str = request.param[2] + + return WindowData(states_str, dir_str, metrics_generator(metric_str)) + + +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] + + # FIXME: these should be generalized fluxes + # FIXME: make a clear distinction between fluxes in physical and + # generalized coordinates + self.flux_pair = data.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/test_eigensystem.py b/test/test_eigensystem.py index e8f3a9cf30863d8d5126960d584a733acfeb7d20..c077e6822beaf2c035da81eb92c11d161a2c915c 100644 --- a/test/test_eigensystem.py +++ b/test/test_eigensystem.py @@ -35,30 +35,7 @@ import pytest import utilities as u import weno_reference_implementation as ref -from data_for_test import window_data - - -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] - - # FIXME: these should be generalized fluxes - # FIXME: make a clear distinction between fluxes in physical and - # generalized coordinates - self.flux_pair = data.fluxes[:,2:4] - - # FIXME: is this the right computation? - self.frozen_metrics = np.mean(data.metrics[2:4], axis=0) - - -@pytest.fixture(scope="session") -def pair_data(window_data): - return PairData(window_data) +from data_for_test import pair_data, window_data, metrics_generator def test_roe_uniform_grid_ideal_gas(queue, pair_data): diff --git a/test/test_flux_window_ops.py b/test/test_flux_window_ops.py index b48e98a508ac3bd79cbf9d7cc554b248d59e8f2b..0f4742a4822cc5e977cb7cc64fca43d05d0834be 100644 --- a/test/test_flux_window_ops.py +++ b/test/test_flux_window_ops.py @@ -35,11 +35,13 @@ import pytest import utilities as u import weno_reference_implementation as ref -from data_for_test import window_data +from data_for_test import PairData, window_data, metrics_generator class WindowResults: def __init__(self, queue, data): + pair = PairData(data) + self.nvars = data.nvars self.ndim = data.ndim self.direction = data.direction @@ -47,17 +49,14 @@ class WindowResults: self.states = data.states self.fluxes = data.fluxes - # FIXME: duplicate of code in PairData - self.state_pair = data.states[:,2:4] + self.state_pair = pair.state_pair self.metrics = data.metrics self.jacobians = data.jacobians - # FIXME: should be computed directly from the metrics and jacobians - # FIXME: partial duplicate of code in PairData - 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 + self.frozen_metrics = pair.frozen_metrics + self.frozen_jacobian = pair.frozen_jacobian + self.combined_frozen_metric = pair.combined_frozen_metrics[self.dir_internal] self.lam_pointwise = ref.lambda_pointwise( self.states, self.metrics, self.dir_internal) @@ -74,9 +73,9 @@ class WindowResults: 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.oscillation_pos, self.combined_frozen_metric) self.weno_weights_neg = ref.weno_weights( - self.oscillation_neg, self.combined_frozen_metrics) + self.oscillation_neg, self.combined_frozen_metric) self.consistent = ref.consistent_part(self.fluxes) self.dissipation_pos = ref.dissipation_part( @@ -137,7 +136,7 @@ def test_weno_flux(queue, window_results): generalized_fluxes=data.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) @@ -169,7 +168,7 @@ def test_dissipation_part_pos(queue, window_results): 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) @@ -186,7 +185,7 @@ def test_dissipation_part_neg(queue, window_results): 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) @@ -203,7 +202,7 @@ def test_weno_weights_pos(queue, window_results): 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) @@ -222,7 +221,7 @@ def test_weno_weights_neg(queue, window_results): 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) diff --git a/test/weno_reference_implementation.py b/test/weno_reference_implementation.py index 179402a66011c2b967ae51d7aad6c9f1675396f7..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])