From 1e7bfc8c65576ca7a9517936b5558b770f8af4f6 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Sat, 18 Apr 2020 10:17:25 -0600 Subject: [PATCH 01/25] start refactoring RHS tests -- for now only prints direction index --- test/test_flux_derivatives.py | 312 +++++++++++++++++++--------------- 1 file changed, 171 insertions(+), 141 deletions(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 2c4ec9c..6a8b35d 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -37,147 +37,54 @@ import utilities as u 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): - data = cfd_test_data_fixture - - prg = u.get_weno_program() - prg = prg.copy(target=lp.PyOpenCLTarget(queue.device)) - - flux_derivatives_dev = u.empty_array_on_device(queue, *data.flux_dims) - - prg(queue, - nvars=data.nvars, - ndim=data.ndim, - nx=data.nx, - ny=data.ny, - nz=data.nz, - d=data.direction, - states=data.states, - fluxes=data.fluxes, - metrics=data.metrics, - metric_jacobians=data.jacobians, - flux_derivatives=flux_derivatives_dev) - - u.compare_arrays(flux_derivatives_dev.get(), data.flux_derivatives) +class GridResults: + #nx = 10 + #ny = 10 + #nz = 10 + #halo = 3 + #ndim = 3 + #nvars = 5 + dir_map = {"xi":1, "eta":2, "zeta":3} + + def __init__(self, dir_str): + self.direction = self.dir_map[dir_str] + + #self.states = pass + #self.fluxes = pass + #self.metrics = pass + #self.jacobians = pass + #self.flux_derivatives = pass + + +@pytest.fixture(scope="session", params=["xi", "eta", "zeta"]) +def grid_results(request): + return GridResults(request.param) + + +def test_compute_flux_derivatives(queue, grid_results): + data = grid_results + + print(data.direction) + + #prg = u.get_weno_program() + #prg = prg.copy(target=lp.PyOpenCLTarget(queue.device)) + + #flux_derivatives_dev = u.empty_array_on_device(queue, *data.flux_dims) + + #prg(queue, + # nvars=data.nvars, + # ndim=data.ndim, + # nx=data.nx, + # ny=data.ny, + # nz=data.nz, + # d=data.direction, + # states=data.states, + # fluxes=data.fluxes, + # metrics=data.metrics, + # metric_jacobians=data.jacobians, + # flux_derivatives=flux_derivatives_dev) + + #u.compare_arrays(flux_derivatives_dev.get(), data.flux_derivatives) # This lets you run 'python test.py test_case(cl._csc)' without pytest. @@ -187,3 +94,126 @@ if __name__ == "__main__": exec(sys.argv[1]) else: pytest.main([__file__]) + + +### old code for reference ### + +## 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] + + -- GitLab From a95a43aa3f212e38210c7e6de57533cc51e7d712 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Sat, 18 Apr 2020 11:09:00 -0600 Subject: [PATCH 02/25] flux derivative test is working with bogus values, but correct array sizes --- WENO.F90 | 1 + test/test_flux_derivatives.py | 65 ++++++++++++++++++----------------- 2 files changed, 34 insertions(+), 32 deletions(-) diff --git a/WENO.F90 b/WENO.F90 index f14987f..6e02d65 100644 --- a/WENO.F90 +++ b/WENO.F90 @@ -195,6 +195,7 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, & do v=1,nvars flux_derivatives(v, i, j, k) = & (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j, k-1)) + flux_derivatives(v, i, j, k) = 0.0d0 end do end do end do diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 6a8b35d..697478d 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -38,22 +38,25 @@ import weno_reference_implementation as ref class GridResults: - #nx = 10 - #ny = 10 - #nz = 10 - #halo = 3 - #ndim = 3 - #nvars = 5 + nx = 10 + ny = 10 + nz = 10 + halo = 3 + nxhalo = nx + 2*halo + nyhalo = ny + 2*halo + nzhalo = nz + 2*halo + ndim = 3 + nvars = 5 dir_map = {"xi":1, "eta":2, "zeta":3} def __init__(self, dir_str): self.direction = self.dir_map[dir_str] - #self.states = pass - #self.fluxes = pass - #self.metrics = pass - #self.jacobians = pass - #self.flux_derivatives = pass + self.states = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") + self.generalized_fluxes = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") + self.metrics = np.full((self.ndim, self.ndim, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") + self.jacobians = np.full((self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") + self.generalized_flux_derivatives = np.full((self.nvars, self.nx, self.ny, self.nz), 0.0) @pytest.fixture(scope="session", params=["xi", "eta", "zeta"]) @@ -64,27 +67,25 @@ def grid_results(request): def test_compute_flux_derivatives(queue, grid_results): data = grid_results - print(data.direction) - - #prg = u.get_weno_program() - #prg = prg.copy(target=lp.PyOpenCLTarget(queue.device)) - - #flux_derivatives_dev = u.empty_array_on_device(queue, *data.flux_dims) - - #prg(queue, - # nvars=data.nvars, - # ndim=data.ndim, - # nx=data.nx, - # ny=data.ny, - # nz=data.nz, - # d=data.direction, - # states=data.states, - # fluxes=data.fluxes, - # metrics=data.metrics, - # metric_jacobians=data.jacobians, - # flux_derivatives=flux_derivatives_dev) - - #u.compare_arrays(flux_derivatives_dev.get(), data.flux_derivatives) + prg = u.get_weno_program() + prg = prg.copy(target=lp.PyOpenCLTarget(queue.device)) + + flux_derivatives_dev = u.empty_array_on_device(queue, data.nvars,data.nx,data.ny,data.nz) + + prg(queue, + nvars=data.nvars, + ndim=data.ndim, + nx=data.nx, + ny=data.ny, + nz=data.nz, + d=data.direction, + states=data.states, + fluxes=data.generalized_fluxes, + metrics=data.metrics, + metric_jacobians=data.jacobians, + flux_derivatives=flux_derivatives_dev) + + u.compare_arrays(flux_derivatives_dev.get(), data.generalized_flux_derivatives) # This lets you run 'python test.py test_case(cl._csc)' without pytest. -- GitLab From eece86686d89bc1bb6f2d190730715212814f40e Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Sat, 18 Apr 2020 11:14:08 -0600 Subject: [PATCH 03/25] create a variable to hold output dimensions --- test/test_flux_derivatives.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 697478d..16d37fb 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -48,6 +48,7 @@ class GridResults: ndim = 3 nvars = 5 dir_map = {"xi":1, "eta":2, "zeta":3} + output_dims = (nvars, nx, ny, nz) def __init__(self, dir_str): self.direction = self.dir_map[dir_str] @@ -56,7 +57,7 @@ class GridResults: self.generalized_fluxes = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") self.metrics = np.full((self.ndim, self.ndim, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") self.jacobians = np.full((self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") - self.generalized_flux_derivatives = np.full((self.nvars, self.nx, self.ny, self.nz), 0.0) + self.generalized_flux_derivatives = np.full(self.output_dims, 0.0) @pytest.fixture(scope="session", params=["xi", "eta", "zeta"]) @@ -70,7 +71,7 @@ def test_compute_flux_derivatives(queue, grid_results): prg = u.get_weno_program() prg = prg.copy(target=lp.PyOpenCLTarget(queue.device)) - flux_derivatives_dev = u.empty_array_on_device(queue, data.nvars,data.nx,data.ny,data.nz) + flux_derivatives_dev = u.empty_array_on_device(queue, *data.output_dims) prg(queue, nvars=data.nvars, -- GitLab From c34be36b327da32a9433e2f6f972b52f5aaad4f7 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Mon, 20 Apr 2020 17:47:11 -0600 Subject: [PATCH 04/25] split data and results --- test/test_flux_derivatives.py | 28 ++++++++++++++++++++++++---- 1 file changed, 24 insertions(+), 4 deletions(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 16d37fb..fd42251 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -37,7 +37,7 @@ import utilities as u import weno_reference_implementation as ref -class GridResults: +class GridData: nx = 10 ny = 10 nz = 10 @@ -54,15 +54,35 @@ class GridResults: self.direction = self.dir_map[dir_str] self.states = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") - self.generalized_fluxes = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") self.metrics = np.full((self.ndim, self.ndim, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") self.jacobians = np.full((self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") + + +class GridResults: + def __init__(self, data): + self.nx = data.nx + self.ny = data.nx + self.nz = data.nx + self.nvars = data.nvars + self.ndim = data.ndim + self.output_dims = data.output_dims + self.direction = data.direction + + self.states = data.states + self.generalized_fluxes = np.full((data.nvars, data.nxhalo, data.nyhalo, data.nzhalo), 1.0, dtype=np.float64, order="F") + self.metrics = data.metrics + self.jacobians = data.jacobians self.generalized_flux_derivatives = np.full(self.output_dims, 0.0) @pytest.fixture(scope="session", params=["xi", "eta", "zeta"]) -def grid_results(request): - return GridResults(request.param) +def grid_data(request): + return GridData(request.param) + + +@pytest.fixture(scope="session") +def grid_results(grid_data): + return GridResults(grid_data) def test_compute_flux_derivatives(queue, grid_results): -- GitLab From a3e895e75c15e9fa9471eb8fe2ed34ca4b07c0d0 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Mon, 20 Apr 2020 17:52:52 -0600 Subject: [PATCH 05/25] using reference implementation for GridData --- test/test_flux_derivatives.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index fd42251..4c042fb 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -38,6 +38,10 @@ import weno_reference_implementation as ref class GridData: + ndim = ref.gas.ndim + nvars = ref.gas.nvars + dirs = ref.gas.dirs + nx = 10 ny = 10 nz = 10 @@ -45,13 +49,11 @@ class GridData: nxhalo = nx + 2*halo nyhalo = ny + 2*halo nzhalo = nz + 2*halo - ndim = 3 - nvars = 5 - dir_map = {"xi":1, "eta":2, "zeta":3} + output_dims = (nvars, nx, ny, nz) def __init__(self, dir_str): - self.direction = self.dir_map[dir_str] + self.direction = self.dirs[dir_str] self.states = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") self.metrics = np.full((self.ndim, self.ndim, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") @@ -75,7 +77,7 @@ class GridResults: self.generalized_flux_derivatives = np.full(self.output_dims, 0.0) -@pytest.fixture(scope="session", params=["xi", "eta", "zeta"]) +@pytest.fixture(scope="session", params=["x", "y", "z"]) def grid_data(request): return GridData(request.param) -- GitLab From 600a2319794752385e3709025ba5f1e442860c1d Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Mon, 20 Apr 2020 17:54:12 -0600 Subject: [PATCH 06/25] move generalized fluxes to GridData --- test/test_flux_derivatives.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 4c042fb..cb3a608 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -59,6 +59,8 @@ class GridData: self.metrics = np.full((self.ndim, self.ndim, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") self.jacobians = np.full((self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") + self.generalized_fluxes = np.full((data.nvars, data.nxhalo, data.nyhalo, data.nzhalo), 1.0, dtype=np.float64, order="F") + class GridResults: def __init__(self, data): @@ -71,7 +73,7 @@ class GridResults: self.direction = data.direction self.states = data.states - self.generalized_fluxes = np.full((data.nvars, data.nxhalo, data.nyhalo, data.nzhalo), 1.0, dtype=np.float64, order="F") + self.generalized_fluxes = data.generalized_fluxes self.metrics = data.metrics self.jacobians = data.jacobians self.generalized_flux_derivatives = np.full(self.output_dims, 0.0) -- GitLab From acbb07d753840e88d0a555c16f1a7d63ff8ebf0d Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Mon, 20 Apr 2020 18:11:44 -0600 Subject: [PATCH 07/25] rename window-related fixtures --- test/data_for_test.py | 18 +++++++++--------- test/test_eigensystem.py | 2 +- test/test_flux_derivatives.py | 2 +- test/test_flux_window_ops.py | 2 +- 4 files changed, 12 insertions(+), 12 deletions(-) diff --git a/test/data_for_test.py b/test/data_for_test.py index 99feb03..128cded 100644 --- a/test/data_for_test.py +++ b/test/data_for_test.py @@ -32,11 +32,11 @@ class WindowData: ndim = ref.gas.ndim dirs = ref.gas.dirs - def __init__(self, states, direction, window_metrics): + def __init__(self, window_states, direction, window_metrics): self.direction = self.dirs[direction] self.dir_internal = self.direction-1 - self.states = states + self.states = window_states self.metrics, self.jacobians = window_metrics self.fluxes = ref.pointwise_fluxes(self.states) @@ -47,8 +47,8 @@ class WindowData: @pytest.fixture(scope="session") -def states_generator(): - def states(states_str): +def window_states_generator(): + def window_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) @@ -60,11 +60,11 @@ def states_generator(): state_pair = u.transposed_array_from_string(states_str) return u.expand_to_n(state_pair, 6) - return states + return window_states @pytest.fixture(scope="session") -def metrics_generator(): +def window_metrics_generator(): ndim = ref.gas.ndim def window_metrics(metric_str, dim): @@ -106,14 +106,14 @@ def metrics_generator(): ("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, states_generator, metrics_generator): +def window_data(request, window_states_generator, window_metrics_generator): states_str = request.param[0] dir_str = request.param[1] metric_str = request.param[2] - return WindowData(states_generator(states_str), + return WindowData(window_states_generator(states_str), dir_str, - metrics_generator(metric_str, dir_str)) + window_metrics_generator(metric_str, dir_str)) class PairData: diff --git a/test/test_eigensystem.py b/test/test_eigensystem.py index 8c1ba93..793337f 100644 --- a/test/test_eigensystem.py +++ b/test/test_eigensystem.py @@ -36,7 +36,7 @@ import pytest import utilities as u import weno_reference_implementation as ref from data_for_test import ( - pair_data, window_data, states_generator, metrics_generator) + pair_data, window_data, window_states_generator, window_metrics_generator) def test_roe_uniform_grid_ideal_gas(queue, pair_data): diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index cb3a608..415c3aa 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -59,7 +59,7 @@ class GridData: self.metrics = np.full((self.ndim, self.ndim, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") self.jacobians = np.full((self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") - self.generalized_fluxes = np.full((data.nvars, data.nxhalo, data.nyhalo, data.nzhalo), 1.0, dtype=np.float64, order="F") + self.generalized_fluxes = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") class GridResults: diff --git a/test/test_flux_window_ops.py b/test/test_flux_window_ops.py index 9685c38..ad228da 100644 --- a/test/test_flux_window_ops.py +++ b/test/test_flux_window_ops.py @@ -36,7 +36,7 @@ import pytest import utilities as u import weno_reference_implementation as ref from data_for_test import ( - PairData, window_data, states_generator, metrics_generator) + PairData, window_data, window_states_generator, window_metrics_generator) class WindowResults: -- GitLab From 1c9e35d39acc7417cd1538cba3850c12e5176120 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Mon, 20 Apr 2020 18:36:21 -0600 Subject: [PATCH 08/25] added grid data generators --- test/test_flux_derivatives.py | 59 +++++++++++++++++++++++++++++------ 1 file changed, 50 insertions(+), 9 deletions(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 415c3aa..412c92c 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -52,16 +52,62 @@ class GridData: output_dims = (nvars, nx, ny, nz) - def __init__(self, dir_str): + def __init__(self, grid_states, dir_str, grid_metrics): self.direction = self.dirs[dir_str] - self.states = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") - self.metrics = np.full((self.ndim, self.ndim, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") - self.jacobians = np.full((self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") + self.states = grid_states + self.metrics, self.jacobians = grid_metrics self.generalized_fluxes = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") +@pytest.fixture(scope="session") +def grid_states_generator(): + def grid_states(states_str): + if states_str == "flat": + return np.full((5, 16, 16, 16), 1.0, dtype=np.float64, order="F") + #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: + raise ValueError("States string {} not supported".format(states_str)) + + return grid_states + + +@pytest.fixture(scope="session") +def grid_metrics_generator(): + ndim = ref.gas.ndim + + def grid_metrics(metric_str, dim): + if (metric_str == "uniform"): + metrics = np.full((3, 3, 16, 16, 16), 1.0, dtype=np.float64, order="F") + jacobians = np.full((16, 16, 16), 1.0, dtype=np.float64, order="F") + else: + raise ValueError("Metric string {} not supported".format(metric_str)) + return metrics, jacobians + + return grid_metrics + + +@pytest.fixture(scope="session", params=[ + ("flat", "x", "uniform"), + ("flat", "y", "uniform"), + ("flat", "z", "uniform") + ]) +def grid_data(request, grid_states_generator, grid_metrics_generator): + states_str = request.param[0] + dir_str = request.param[1] + metric_str = request.param[2] + + return GridData(grid_states_generator(states_str), + dir_str, + grid_metrics_generator(metric_str, dir_str)) + + class GridResults: def __init__(self, data): self.nx = data.nx @@ -79,11 +125,6 @@ class GridResults: self.generalized_flux_derivatives = np.full(self.output_dims, 0.0) -@pytest.fixture(scope="session", params=["x", "y", "z"]) -def grid_data(request): - return GridData(request.param) - - @pytest.fixture(scope="session") def grid_results(grid_data): return GridResults(grid_data) -- GitLab From 95a00747f74245e758cdbbdba6522f1011963832 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Mon, 20 Apr 2020 18:42:40 -0600 Subject: [PATCH 09/25] added some fixme notes --- test/test_flux_derivatives.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 412c92c..2a7af38 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -42,6 +42,7 @@ class GridData: nvars = ref.gas.nvars dirs = ref.gas.dirs + # FIXME: Try with multiple dimensions nx = 10 ny = 10 nz = 10 @@ -58,6 +59,7 @@ class GridData: self.states = grid_states self.metrics, self.jacobians = grid_metrics + # FIXME: compute fluxes / generalized fluxes from states and metrics self.generalized_fluxes = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") @@ -65,6 +67,7 @@ class GridData: def grid_states_generator(): def grid_states(states_str): if states_str == "flat": + # FIXME: use reference gas to setup states and make overall dimensions generic return np.full((5, 16, 16, 16), 1.0, dtype=np.float64, order="F") #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) @@ -84,6 +87,7 @@ def grid_metrics_generator(): def grid_metrics(metric_str, dim): if (metric_str == "uniform"): + # FIXME: fill identity matrices correctly and make overall dimensions generic metrics = np.full((3, 3, 16, 16, 16), 1.0, dtype=np.float64, order="F") jacobians = np.full((16, 16, 16), 1.0, dtype=np.float64, order="F") else: @@ -97,6 +101,7 @@ def grid_metrics_generator(): ("flat", "x", "uniform"), ("flat", "y", "uniform"), ("flat", "z", "uniform") + # FIXME: add remaining cases ]) def grid_data(request, grid_states_generator, grid_metrics_generator): states_str = request.param[0] @@ -122,6 +127,8 @@ class GridResults: self.generalized_fluxes = data.generalized_fluxes self.metrics = data.metrics self.jacobians = data.jacobians + + # FIXME: compute flux derivatives from WENO reference implementation self.generalized_flux_derivatives = np.full(self.output_dims, 0.0) -- GitLab From ed1aeaa9ba440eb12544f27ba9f89340a57f0003 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Mon, 8 Jun 2020 20:41:02 -0600 Subject: [PATCH 10/25] allow GridData to accept a specification for grid sizes --- test/test_flux_derivatives.py | 39 ++++++++++++++++++++--------------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 2a7af38..73341b1 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -42,19 +42,9 @@ class GridData: nvars = ref.gas.nvars dirs = ref.gas.dirs - # FIXME: Try with multiple dimensions - nx = 10 - ny = 10 - nz = 10 - halo = 3 - nxhalo = nx + 2*halo - nyhalo = ny + 2*halo - nzhalo = nz + 2*halo - - output_dims = (nvars, nx, ny, nz) - - def __init__(self, grid_states, dir_str, grid_metrics): + def __init__(self, grid_states, grid_sizes_str, dir_str, grid_metrics): self.direction = self.dirs[dir_str] + self.set_grid_sizes(grid_sizes_str) self.states = grid_states self.metrics, self.jacobians = grid_metrics @@ -62,6 +52,19 @@ class GridData: # FIXME: compute fluxes / generalized fluxes from states and metrics self.generalized_fluxes = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") + def set_grid_sizes(self, grid_sizes_str): + self.nx, self.ny, self.nz = self.parse_grid_sizes(grid_sizes_str) + + halo = 3 + self.nxhalo = self.nx + 2*halo + self.nyhalo = self.ny + 2*halo + self.nzhalo = self.nz + 2*halo + + self.output_dims = (self.nvars, self.nx, self.ny, self.nz) + + def parse_grid_sizes(self, grid_sizes_str): + return map(int, grid_sizes_str.split("x")) + @pytest.fixture(scope="session") def grid_states_generator(): @@ -98,17 +101,19 @@ def grid_metrics_generator(): @pytest.fixture(scope="session", params=[ - ("flat", "x", "uniform"), - ("flat", "y", "uniform"), - ("flat", "z", "uniform") + ("flat", "10x10x10", "x", "uniform"), + ("flat", "10x10x10", "y", "uniform"), + ("flat", "10x10x10", "z", "uniform") # FIXME: add remaining cases ]) def grid_data(request, grid_states_generator, grid_metrics_generator): states_str = request.param[0] - dir_str = request.param[1] - metric_str = request.param[2] + grid_sizes_str = request.param[1] + dir_str = request.param[2] + metric_str = request.param[3] return GridData(grid_states_generator(states_str), + grid_sizes_str, dir_str, grid_metrics_generator(metric_str, dir_str)) -- GitLab From 493d0e949560cb6185264041df8050ef00d628e7 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Mon, 8 Jun 2020 20:58:38 -0600 Subject: [PATCH 11/25] make state and metric generators generic with respect to grid size --- test/test_flux_derivatives.py | 49 ++++++++++++++++++++++------------- 1 file changed, 31 insertions(+), 18 deletions(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 73341b1..a997190 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -37,6 +37,15 @@ import utilities as u import weno_reference_implementation as ref +def parse_grid_sizes(grid_sizes_str): + return map(int, grid_sizes_str.split("x")) + + +def halo_sizes(*grid_sizes): + halo = 3 + return map(lambda n: n + 2*halo, grid_sizes) + + class GridData: ndim = ref.gas.ndim nvars = ref.gas.nvars @@ -50,28 +59,27 @@ class GridData: self.metrics, self.jacobians = grid_metrics # FIXME: compute fluxes / generalized fluxes from states and metrics - self.generalized_fluxes = np.full((self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), 1.0, dtype=np.float64, order="F") + self.generalized_fluxes = np.full( + (self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), + 1.0, dtype=np.float64, order="F") def set_grid_sizes(self, grid_sizes_str): - self.nx, self.ny, self.nz = self.parse_grid_sizes(grid_sizes_str) + self.nx, self.ny, self.nz = parse_grid_sizes(grid_sizes_str) - halo = 3 - self.nxhalo = self.nx + 2*halo - self.nyhalo = self.ny + 2*halo - self.nzhalo = self.nz + 2*halo + self.nxhalo, self.nyhalo, self.nzhalo = halo_sizes(self.nx, self.ny, self.nz) self.output_dims = (self.nvars, self.nx, self.ny, self.nz) - def parse_grid_sizes(self, grid_sizes_str): - return map(int, grid_sizes_str.split("x")) - @pytest.fixture(scope="session") def grid_states_generator(): - def grid_states(states_str): + def grid_states(states_str, grid_sizes_str): + nxhalo, nyhalo, nzhalo = halo_sizes(*parse_grid_sizes(grid_sizes_str)) if states_str == "flat": - # FIXME: use reference gas to setup states and make overall dimensions generic - return np.full((5, 16, 16, 16), 1.0, dtype=np.float64, order="F") + # FIXME: use reference gas to setup states + return np.full( + (ref.gas.nvars, nxhalo, nyhalo, nzhalo), + 1.0, dtype=np.float64, order="F") #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) @@ -88,11 +96,16 @@ def grid_states_generator(): def grid_metrics_generator(): ndim = ref.gas.ndim - def grid_metrics(metric_str, dim): + def grid_metrics(metric_str, grid_sizes_str, dim): + nxhalo, nyhalo, nzhalo = halo_sizes(*parse_grid_sizes(grid_sizes_str)) if (metric_str == "uniform"): - # FIXME: fill identity matrices correctly and make overall dimensions generic - metrics = np.full((3, 3, 16, 16, 16), 1.0, dtype=np.float64, order="F") - jacobians = np.full((16, 16, 16), 1.0, dtype=np.float64, order="F") + # FIXME: fill identity matrices correctly + metrics = np.full( + (ref.gas.ndim, ref.gas.ndim, nxhalo, nyhalo, nzhalo), + 1.0, dtype=np.float64, order="F") + jacobians = np.full( + (nxhalo, nyhalo, nzhalo), + 1.0, dtype=np.float64, order="F") else: raise ValueError("Metric string {} not supported".format(metric_str)) return metrics, jacobians @@ -112,10 +125,10 @@ def grid_data(request, grid_states_generator, grid_metrics_generator): dir_str = request.param[2] metric_str = request.param[3] - return GridData(grid_states_generator(states_str), + return GridData(grid_states_generator(states_str, grid_sizes_str), grid_sizes_str, dir_str, - grid_metrics_generator(metric_str, dir_str)) + grid_metrics_generator(metric_str, grid_sizes_str, dir_str)) class GridResults: -- GitLab From 966142b0a7182fc954b1695d312c5949d3df0403 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Tue, 9 Jun 2020 21:35:37 -0600 Subject: [PATCH 12/25] clean up state and metric generators so the flat/uniform cases use real values --- test/test_flux_derivatives.py | 29 +++++++++++++++++------------ 1 file changed, 17 insertions(+), 12 deletions(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index a997190..405beaa 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -76,16 +76,17 @@ def grid_states_generator(): def grid_states(states_str, grid_sizes_str): nxhalo, nyhalo, nzhalo = halo_sizes(*parse_grid_sizes(grid_sizes_str)) if states_str == "flat": - # FIXME: use reference gas to setup states - return np.full( - (ref.gas.nvars, nxhalo, nyhalo, nzhalo), - 1.0, dtype=np.float64, order="F") - #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)]) + states = np.empty( + (ref.gas.nvars, nxhalo, nyhalo, nzhalo), + dtype=np.float64, order="F") + rho = 1.0 + vel = np.array([1.0, 1.0, 1.0]) + pressure = 1.6 + state = ref.gas.conserved(rho, vel, pressure) + for i in range(nxhalo): + for j in range(nyhalo): + for k in range(nzhalo): + states[:, i, j, k] = state else: raise ValueError("States string {} not supported".format(states_str)) @@ -99,10 +100,14 @@ def grid_metrics_generator(): def grid_metrics(metric_str, grid_sizes_str, dim): nxhalo, nyhalo, nzhalo = halo_sizes(*parse_grid_sizes(grid_sizes_str)) if (metric_str == "uniform"): - # FIXME: fill identity matrices correctly metrics = np.full( (ref.gas.ndim, ref.gas.ndim, nxhalo, nyhalo, nzhalo), - 1.0, dtype=np.float64, order="F") + 0.0, dtype=np.float64, order="F") + for i in range(nxhalo): + for j in range(nyhalo): + for k in range(nzhalo): + for m in range(ref.gas.ndim): + metrics[m, m, i, j, k] = 1.0 jacobians = np.full( (nxhalo, nyhalo, nzhalo), 1.0, dtype=np.float64, order="F") -- GitLab From 1fb4bf00647b043c6a0b83525716186e61bc7c19 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Fri, 3 Jul 2020 19:56:14 -0600 Subject: [PATCH 13/25] fix generalized_fluxes so they are actually computed from states/metrics --- test/test_flux_derivatives.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 405beaa..8278b03 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -53,15 +53,22 @@ class GridData: def __init__(self, grid_states, grid_sizes_str, dir_str, grid_metrics): self.direction = self.dirs[dir_str] + self.dir_internal = self.direction-1 self.set_grid_sizes(grid_sizes_str) self.states = grid_states self.metrics, self.jacobians = grid_metrics - # FIXME: compute fluxes / generalized fluxes from states and metrics - self.generalized_fluxes = np.full( + self.generalized_fluxes = np.empty( (self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), - 1.0, dtype=np.float64, order="F") + dtype=np.float64, order="F") + for i in range(self.nxhalo): + for j in range(self.nyhalo): + for k in range(self.nzhalo): + flux = ref.gas.flux(self.states[:, i, j, k]) + self.generalized_fluxes[:,i,j,k] = np.dot( + flux, self.metrics[:,:,i,j,k]/self.jacobians[i,j,k] + )[:,self.dir_internal] def set_grid_sizes(self, grid_sizes_str): self.nx, self.ny, self.nz = parse_grid_sizes(grid_sizes_str) @@ -89,6 +96,7 @@ def grid_states_generator(): states[:, i, j, k] = state else: raise ValueError("States string {} not supported".format(states_str)) + return states return grid_states -- GitLab From 37e396b5a7e7b0307bd83469b8bbdf187713c67e Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Wed, 8 Jul 2020 18:06:48 -0600 Subject: [PATCH 14/25] move WindowResults to data_for_test --- test/data_for_test.py | 59 +++++++++++++++++++++++++++++++++ test/test_flux_window_ops.py | 63 ++---------------------------------- 2 files changed, 61 insertions(+), 61 deletions(-) diff --git a/test/data_for_test.py b/test/data_for_test.py index 128cded..0a84394 100644 --- a/test/data_for_test.py +++ b/test/data_for_test.py @@ -137,4 +137,63 @@ def pair_data(window_data): return PairData(window_data) +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_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") + + 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) + + # vim: foldmethod=marker diff --git a/test/test_flux_window_ops.py b/test/test_flux_window_ops.py index ad228da..3f30064 100644 --- a/test/test_flux_window_ops.py +++ b/test/test_flux_window_ops.py @@ -35,67 +35,8 @@ import pytest import utilities as u import weno_reference_implementation as ref -from data_for_test import ( - PairData, window_data, window_states_generator, window_metrics_generator) - - -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_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") - - 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) +from data_for_test import (window_results, window_data, + window_states_generator, window_metrics_generator) def test_pointwise_eigenvalues(queue, window_results): -- GitLab From 4d012cb63b8b722ee937b3c859879b323bd9c06c Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Thu, 24 Sep 2020 20:18:57 -0600 Subject: [PATCH 15/25] Compute expected flux derivatives from WENO reference implementation --- test/test_flux_derivatives.py | 60 +++++++++++++++++++++++---- test/weno_reference_implementation.py | 3 +- 2 files changed, 55 insertions(+), 8 deletions(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 8278b03..48f20f8 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -35,6 +35,10 @@ import pytest import utilities as u import weno_reference_implementation as ref +from data_for_test import WindowResults, WindowData + + +halo = 3 def parse_grid_sizes(grid_sizes_str): @@ -42,7 +46,6 @@ def parse_grid_sizes(grid_sizes_str): def halo_sizes(*grid_sizes): - halo = 3 return map(lambda n: n + 2*halo, grid_sizes) @@ -52,7 +55,8 @@ class GridData: dirs = ref.gas.dirs def __init__(self, grid_states, grid_sizes_str, dir_str, grid_metrics): - self.direction = self.dirs[dir_str] + self.dir_str = dir_str + self.direction = self.dirs[self.dir_str] self.dir_internal = self.direction-1 self.set_grid_sizes(grid_sizes_str) @@ -145,27 +149,69 @@ def grid_data(request, grid_states_generator, grid_metrics_generator): class GridResults: - def __init__(self, data): + def __init__(self, queue, data): + self.queue = queue self.nx = data.nx self.ny = data.nx self.nz = data.nx self.nvars = data.nvars self.ndim = data.ndim self.output_dims = data.output_dims + self.dir_str = data.dir_str self.direction = data.direction + self.dir_internal = data.dir_internal self.states = data.states self.generalized_fluxes = data.generalized_fluxes self.metrics = data.metrics self.jacobians = data.jacobians - # FIXME: compute flux derivatives from WENO reference implementation - self.generalized_flux_derivatives = np.full(self.output_dims, 0.0) + self.generalized_flux_derivatives = np.full(self.output_dims, 0.0, + dtype=np.float64, order="F") + for i in range(self.nx): + for j in range(self.ny): + for k in range(self.nz): + flux_pos = self.local_flux( + *self.window_variables(i,j,k)) + neg_loc = ( + i - (1 if 0 is self.dir_internal else 0), + j - (1 if 1 is self.dir_internal else 0), + k - (1 if 2 is self.dir_internal else 0)) + flux_neg = self.local_flux( + *self.window_variables(*neg_loc)) + self.generalized_flux_derivatives[:,i,j,k] = flux_pos - flux_neg + + def window_variables(self, i, j, k): + def lo(n, d): + return n + halo - (2 if d == self.dir_internal else 0) + + def hi(n, d): + return n + halo + (4 if d == self.dir_internal else 1) + + xlo = lo(i, 0) + xhi = hi(i, 0) + ylo = lo(j, 1) + yhi = hi(j, 1) + zlo = lo(k, 2) + zhi = hi(k, 2) + local_states = self.states[:,xlo:xhi,ylo:yhi,zlo:zhi] + local_metrics = self.metrics[:,:,xlo:xhi,ylo:yhi,zlo:zhi] + local_jacobians = self.jacobians[xlo:xhi,ylo:yhi,zlo:zhi] + return (np.squeeze(local_states), + np.moveaxis(np.squeeze(local_metrics), -1, 0), + np.squeeze(local_jacobians)) + + def local_flux(self, state_window, metric_window, jacobian_window): + local_data = WindowData( + state_window, self.dir_str, + (metric_window, jacobian_window)) + return WindowResults(self.queue, local_data).weno_flux + @pytest.fixture(scope="session") -def grid_results(grid_data): - return GridResults(grid_data) +def grid_results(queue, grid_data): + return GridResults(queue, grid_data) def test_compute_flux_derivatives(queue, grid_results): diff --git a/test/weno_reference_implementation.py b/test/weno_reference_implementation.py index a7b4095..c971bdf 100644 --- a/test/weno_reference_implementation.py +++ b/test/weno_reference_implementation.py @@ -51,7 +51,8 @@ def roe_eigensystem(queue, state_pair, frozen_metrics, direction): lam_dev = u.empty_array_on_device(queue, nvars) prg(queue, nvars=nvars, ndim=ndim, d=direction, - states=state_pair, metrics_frozen=frozen_metrics, + states=state_pair.copy(order='F'), + metrics_frozen=frozen_metrics.copy(order='F'), R=R_dev, R_inv=R_inv_dev, lambda_roe=lam_dev) return R_dev.get(), R_inv_dev.get(), lam_dev.get() -- GitLab From a9880458306b7a6612a3a4dbede718c3b224ac64 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Thu, 1 Oct 2020 20:11:22 -0600 Subject: [PATCH 16/25] get rid of old commented code --- test/test_flux_derivatives.py | 121 ---------------------------------- 1 file changed, 121 deletions(-) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 48f20f8..7a5ff76 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -247,124 +247,3 @@ if __name__ == "__main__": pytest.main([__file__]) -### old code for reference ### - -## 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] - - -- GitLab From a00c6fe05884352c6970e19de49ea5f3c0d3b0a7 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Thu, 1 Oct 2020 20:35:20 -0600 Subject: [PATCH 17/25] fix merge issues and flake errors --- test/data_for_test.py | 4 +-- test/test_flux_derivatives.py | 31 +++++++++--------- test/test_flux_window_ops.py | 60 ----------------------------------- 3 files changed, 17 insertions(+), 78 deletions(-) diff --git a/test/data_for_test.py b/test/data_for_test.py index 3c6d131..547f81c 100644 --- a/test/data_for_test.py +++ b/test/data_for_test.py @@ -160,7 +160,7 @@ class WindowResults: 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") + 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) @@ -174,7 +174,7 @@ class WindowResults: self.R_inv) self.oscillation_pos = ref.oscillation(self.char_fluxes_pos) - self.oscillation_neg = ref.oscillation(self.char_fluxes_neg[:,::-1]) + self.oscillation_neg = ref.oscillation(self.char_fluxes_neg[:, ::-1]) self.weno_weights_pos = ref.weno_weights( self.oscillation_pos, self.combined_frozen_metric) diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 7a5ff76..dd37737 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -70,9 +70,10 @@ class GridData: for j in range(self.nyhalo): for k in range(self.nzhalo): flux = ref.gas.flux(self.states[:, i, j, k]) - self.generalized_fluxes[:,i,j,k] = np.dot( - flux, self.metrics[:,:,i,j,k]/self.jacobians[i,j,k] - )[:,self.dir_internal] + self.generalized_fluxes[:, i, j, k] = np.dot( + flux, + self.metrics[:, :, i, j, k]/self.jacobians[i, j, k] + )[:, self.dir_internal] def set_grid_sizes(self, grid_sizes_str): self.nx, self.ny, self.nz = parse_grid_sizes(grid_sizes_str) @@ -113,12 +114,12 @@ def grid_metrics_generator(): nxhalo, nyhalo, nzhalo = halo_sizes(*parse_grid_sizes(grid_sizes_str)) if (metric_str == "uniform"): metrics = np.full( - (ref.gas.ndim, ref.gas.ndim, nxhalo, nyhalo, nzhalo), + (ndim, ndim, nxhalo, nyhalo, nzhalo), 0.0, dtype=np.float64, order="F") for i in range(nxhalo): for j in range(nyhalo): for k in range(nzhalo): - for m in range(ref.gas.ndim): + for m in range(ndim): metrics[m, m, i, j, k] = 1.0 jacobians = np.full( (nxhalo, nyhalo, nzhalo), @@ -172,14 +173,15 @@ class GridResults: for j in range(self.ny): for k in range(self.nz): flux_pos = self.local_flux( - *self.window_variables(i,j,k)) + *self.window_variables(i, j, k)) neg_loc = ( - i - (1 if 0 is self.dir_internal else 0), - j - (1 if 1 is self.dir_internal else 0), - k - (1 if 2 is self.dir_internal else 0)) + i - (1 if 0 == self.dir_internal else 0), + j - (1 if 1 == self.dir_internal else 0), + k - (1 if 2 == self.dir_internal else 0)) flux_neg = self.local_flux( *self.window_variables(*neg_loc)) - self.generalized_flux_derivatives[:,i,j,k] = flux_pos - flux_neg + self.generalized_flux_derivatives[:, i, j, k] = ( + flux_pos - flux_neg) def window_variables(self, i, j, k): def lo(n, d): @@ -194,9 +196,9 @@ class GridResults: yhi = hi(j, 1) zlo = lo(k, 2) zhi = hi(k, 2) - local_states = self.states[:,xlo:xhi,ylo:yhi,zlo:zhi] - local_metrics = self.metrics[:,:,xlo:xhi,ylo:yhi,zlo:zhi] - local_jacobians = self.jacobians[xlo:xhi,ylo:yhi,zlo:zhi] + local_states = self.states[:, xlo:xhi, ylo:yhi, zlo:zhi] + local_metrics = self.metrics[:, :, xlo:xhi, ylo:yhi, zlo:zhi] + local_jacobians = self.jacobians[xlo:xhi, ylo:yhi, zlo:zhi] return (np.squeeze(local_states), np.moveaxis(np.squeeze(local_metrics), -1, 0), np.squeeze(local_jacobians)) @@ -208,7 +210,6 @@ class GridResults: return WindowResults(self.queue, local_data).weno_flux - @pytest.fixture(scope="session") def grid_results(queue, grid_data): return GridResults(queue, grid_data) @@ -245,5 +246,3 @@ if __name__ == "__main__": exec(sys.argv[1]) else: pytest.main([__file__]) - - diff --git a/test/test_flux_window_ops.py b/test/test_flux_window_ops.py index 74d71e5..96117b1 100644 --- a/test/test_flux_window_ops.py +++ b/test/test_flux_window_ops.py @@ -40,66 +40,6 @@ from data_for_test import ( # noqa window_states_generator, window_metrics_generator) -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_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") - - 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): # noqa - return WindowResults(queue, window_data) ->>>>>>> master - - def test_pointwise_eigenvalues(queue, window_results): data = window_results -- GitLab From c97f54b477b0a2772aec84de402ba75fb1faeab4 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Thu, 1 Oct 2020 20:42:27 -0600 Subject: [PATCH 18/25] fix another flake bug --- test/test_flux_window_ops.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_flux_window_ops.py b/test/test_flux_window_ops.py index 96117b1..89f3bbd 100644 --- a/test/test_flux_window_ops.py +++ b/test/test_flux_window_ops.py @@ -34,7 +34,6 @@ import logging import pytest import utilities as u -import weno_reference_implementation as ref from data_for_test import ( # noqa window_results, window_data, window_states_generator, window_metrics_generator) -- GitLab From fc09b5e8d385fcb6aa962dbf80cb599cd5cbba6a Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Thu, 1 Oct 2020 20:49:51 -0600 Subject: [PATCH 19/25] trying again to make flake happy --- test/test_flux_window_ops.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_flux_window_ops.py b/test/test_flux_window_ops.py index 89f3bbd..a79b187 100644 --- a/test/test_flux_window_ops.py +++ b/test/test_flux_window_ops.py @@ -34,8 +34,7 @@ import logging import pytest import utilities as u -from data_for_test import ( # noqa - window_results, window_data, +from data_for_test import (window_results, window_data, window_states_generator, window_metrics_generator) -- GitLab From 9af6f77a66b86fc72e84e89eb2662e49587714ac Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Thu, 1 Oct 2020 21:18:28 -0600 Subject: [PATCH 20/25] rearranging stuff so hopefully flake8 will finally be happy --- test/conftest.py | 158 +++++++++++++++++++++++++ test/data_for_test.py | 216 +++++++++++++++++++--------------- test/test_eigensystem.py | 5 +- test/test_flux_derivatives.py | 179 ---------------------------- test/test_flux_window_ops.py | 2 - 5 files changed, 279 insertions(+), 281 deletions(-) diff --git a/test/conftest.py b/test/conftest.py index 67f8e77..5fb6967 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -30,7 +30,14 @@ from pyopencl.tools import ( clear_first_arg_caches ) +import numpy as np import utilities as u +from data_for_test import ( + PairData, WindowData, WindowResults, + GridData, GridResults, + halo_sizes, parse_grid_sizes) +import weno_reference_implementation as ref +import wavy_metrics_3d as wavy3d # {{{ mark for slow tests @@ -104,4 +111,155 @@ def queue(ctx_factory): return u.get_queue(ctx_factory) +@pytest.fixture(scope="session") +def pair_data(window_data): + return PairData(window_data) + + +@pytest.fixture(scope="session") +def window_states_generator(): + def window_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 window_states + + +@pytest.fixture(scope="session") +def window_metrics_generator(): + ndim = ref.gas.ndim + + 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 + + return window_metrics + + +@pytest.fixture(scope="session", params=[ + ("sine", "x", "wavy3d"), + ("sine", "y", "wavy3d"), + ("sine", "z", "wavy3d"), + ("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 window_data(request, window_states_generator, window_metrics_generator): + states_str = request.param[0] + dir_str = request.param[1] + metric_str = request.param[2] + + return WindowData(window_states_generator(states_str), + dir_str, + window_metrics_generator(metric_str, dir_str)) + + +@pytest.fixture(scope="session") +def window_results(queue, window_data): + return WindowResults(queue, window_data) + + +@pytest.fixture(scope="session") +def grid_states_generator(): + def grid_states(states_str, grid_sizes_str): + nxhalo, nyhalo, nzhalo = halo_sizes(*parse_grid_sizes(grid_sizes_str)) + if states_str == "flat": + states = np.empty( + (ref.gas.nvars, nxhalo, nyhalo, nzhalo), + dtype=np.float64, order="F") + rho = 1.0 + vel = np.array([1.0, 1.0, 1.0]) + pressure = 1.6 + state = ref.gas.conserved(rho, vel, pressure) + for i in range(nxhalo): + for j in range(nyhalo): + for k in range(nzhalo): + states[:, i, j, k] = state + else: + raise ValueError("States string {} not supported".format(states_str)) + return states + + return grid_states + + +@pytest.fixture(scope="session") +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)) + if (metric_str == "uniform"): + metrics = np.full( + (ndim, ndim, nxhalo, nyhalo, nzhalo), + 0.0, dtype=np.float64, order="F") + for i in range(nxhalo): + for j in range(nyhalo): + for k in range(nzhalo): + for m in range(ndim): + metrics[m, m, i, j, k] = 1.0 + jacobians = np.full( + (nxhalo, nyhalo, nzhalo), + 1.0, dtype=np.float64, order="F") + else: + raise ValueError("Metric string {} not supported".format(metric_str)) + return metrics, jacobians + + return grid_metrics + + +@pytest.fixture(scope="session", params=[ + ("flat", "10x10x10", "x", "uniform"), + ("flat", "10x10x10", "y", "uniform"), + ("flat", "10x10x10", "z", "uniform") + # FIXME: add remaining cases + ]) +def grid_data(request, grid_states_generator, grid_metrics_generator): + states_str = request.param[0] + grid_sizes_str = request.param[1] + dir_str = request.param[2] + metric_str = request.param[3] + + return GridData(grid_states_generator(states_str, grid_sizes_str), + grid_sizes_str, + dir_str, + grid_metrics_generator(metric_str, grid_sizes_str, dir_str)) + + +@pytest.fixture(scope="session") +def grid_results(queue, grid_data): + return GridResults(queue, grid_data) + + # vim: foldmethod=marker diff --git a/test/data_for_test.py b/test/data_for_test.py index 547f81c..7ef126a 100644 --- a/test/data_for_test.py +++ b/test/data_for_test.py @@ -21,10 +21,24 @@ THE SOFTWARE. """ import numpy as np -import pytest import utilities as u import weno_reference_implementation as ref -import wavy_metrics_3d as wavy3d + + +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) class WindowData: @@ -46,97 +60,6 @@ class WindowData: )[:, :, self.dir_internal].T.copy(order="F") -@pytest.fixture(scope="session") -def window_states_generator(): - def window_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 window_states - - -@pytest.fixture(scope="session") -def window_metrics_generator(): - ndim = ref.gas.ndim - - 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 - - return window_metrics - - -@pytest.fixture(scope="session", params=[ - ("sine", "x", "wavy3d"), - ("sine", "y", "wavy3d"), - ("sine", "z", "wavy3d"), - ("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 window_data(request, window_states_generator, window_metrics_generator): - states_str = request.param[0] - dir_str = request.param[1] - metric_str = request.param[2] - - return WindowData(window_states_generator(states_str), - dir_str, - window_metrics_generator(metric_str, dir_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] - 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) - - class WindowResults: def __init__(self, queue, data): pair = PairData(data) @@ -191,9 +114,110 @@ class WindowResults: self.consistent, self.dissipation_pos, self.dissipation_neg) -@pytest.fixture(scope="session") -def window_results(queue, window_data): - return WindowResults(queue, window_data) +halo = 3 + + +def parse_grid_sizes(grid_sizes_str): + return map(int, grid_sizes_str.split("x")) + + +def halo_sizes(*grid_sizes): + return map(lambda n: n + 2*halo, grid_sizes) + + +class GridData: + ndim = ref.gas.ndim + nvars = ref.gas.nvars + dirs = ref.gas.dirs + + def __init__(self, grid_states, grid_sizes_str, dir_str, grid_metrics): + self.dir_str = dir_str + self.direction = self.dirs[self.dir_str] + self.dir_internal = self.direction-1 + self.set_grid_sizes(grid_sizes_str) + + self.states = grid_states + self.metrics, self.jacobians = grid_metrics + + self.generalized_fluxes = np.empty( + (self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), + dtype=np.float64, order="F") + for i in range(self.nxhalo): + for j in range(self.nyhalo): + for k in range(self.nzhalo): + flux = ref.gas.flux(self.states[:, i, j, k]) + self.generalized_fluxes[:, i, j, k] = np.dot( + flux, + self.metrics[:, :, i, j, k]/self.jacobians[i, j, k] + )[:, self.dir_internal] + + def set_grid_sizes(self, grid_sizes_str): + self.nx, self.ny, self.nz = parse_grid_sizes(grid_sizes_str) + + self.nxhalo, self.nyhalo, self.nzhalo = halo_sizes(self.nx, self.ny, self.nz) + + self.output_dims = (self.nvars, self.nx, self.ny, self.nz) + + +class GridResults: + def __init__(self, queue, data): + self.queue = queue + self.nx = data.nx + self.ny = data.nx + self.nz = data.nx + self.nvars = data.nvars + self.ndim = data.ndim + self.output_dims = data.output_dims + self.dir_str = data.dir_str + self.direction = data.direction + self.dir_internal = data.dir_internal + + self.states = data.states + self.generalized_fluxes = data.generalized_fluxes + self.metrics = data.metrics + self.jacobians = data.jacobians + + self.generalized_flux_derivatives = np.full(self.output_dims, 0.0, + dtype=np.float64, order="F") + for i in range(self.nx): + for j in range(self.ny): + for k in range(self.nz): + flux_pos = self.local_flux( + *self.window_variables(i, j, k)) + neg_loc = ( + i - (1 if 0 == self.dir_internal else 0), + j - (1 if 1 == self.dir_internal else 0), + k - (1 if 2 == self.dir_internal else 0)) + flux_neg = self.local_flux( + *self.window_variables(*neg_loc)) + self.generalized_flux_derivatives[:, i, j, k] = ( + flux_pos - flux_neg) + + def window_variables(self, i, j, k): + def lo(n, d): + return n + halo - (2 if d == self.dir_internal else 0) + + def hi(n, d): + return n + halo + (4 if d == self.dir_internal else 1) + + xlo = lo(i, 0) + xhi = hi(i, 0) + ylo = lo(j, 1) + yhi = hi(j, 1) + zlo = lo(k, 2) + zhi = hi(k, 2) + local_states = self.states[:, xlo:xhi, ylo:yhi, zlo:zhi] + local_metrics = self.metrics[:, :, xlo:xhi, ylo:yhi, zlo:zhi] + local_jacobians = self.jacobians[xlo:xhi, ylo:yhi, zlo:zhi] + return (np.squeeze(local_states), + np.moveaxis(np.squeeze(local_metrics), -1, 0), + np.squeeze(local_jacobians)) + + def local_flux(self, state_window, metric_window, jacobian_window): + local_data = WindowData( + state_window, self.dir_str, + (metric_window, jacobian_window)) + return WindowResults(self.queue, local_data).weno_flux # vim: foldmethod=marker diff --git a/test/test_eigensystem.py b/test/test_eigensystem.py index 504979b..e96c05a 100644 --- a/test/test_eigensystem.py +++ b/test/test_eigensystem.py @@ -34,12 +34,9 @@ import logging import pytest import utilities as u -import weno_reference_implementation as ref # noqa -from data_for_test import ( # noqa - pair_data, window_data, window_states_generator, window_metrics_generator) -def test_roe_uniform_grid_ideal_gas(queue, pair_data): # noqa +def test_roe_uniform_grid_ideal_gas(queue, pair_data): data = pair_data def check_roe_identity(states, R, R_inv): diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index dd37737..426570f 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -34,185 +34,6 @@ import logging import pytest import utilities as u -import weno_reference_implementation as ref -from data_for_test import WindowResults, WindowData - - -halo = 3 - - -def parse_grid_sizes(grid_sizes_str): - return map(int, grid_sizes_str.split("x")) - - -def halo_sizes(*grid_sizes): - return map(lambda n: n + 2*halo, grid_sizes) - - -class GridData: - ndim = ref.gas.ndim - nvars = ref.gas.nvars - dirs = ref.gas.dirs - - def __init__(self, grid_states, grid_sizes_str, dir_str, grid_metrics): - self.dir_str = dir_str - self.direction = self.dirs[self.dir_str] - self.dir_internal = self.direction-1 - self.set_grid_sizes(grid_sizes_str) - - self.states = grid_states - self.metrics, self.jacobians = grid_metrics - - self.generalized_fluxes = np.empty( - (self.nvars, self.nxhalo, self.nyhalo, self.nzhalo), - dtype=np.float64, order="F") - for i in range(self.nxhalo): - for j in range(self.nyhalo): - for k in range(self.nzhalo): - flux = ref.gas.flux(self.states[:, i, j, k]) - self.generalized_fluxes[:, i, j, k] = np.dot( - flux, - self.metrics[:, :, i, j, k]/self.jacobians[i, j, k] - )[:, self.dir_internal] - - def set_grid_sizes(self, grid_sizes_str): - self.nx, self.ny, self.nz = parse_grid_sizes(grid_sizes_str) - - self.nxhalo, self.nyhalo, self.nzhalo = halo_sizes(self.nx, self.ny, self.nz) - - self.output_dims = (self.nvars, self.nx, self.ny, self.nz) - - -@pytest.fixture(scope="session") -def grid_states_generator(): - def grid_states(states_str, grid_sizes_str): - nxhalo, nyhalo, nzhalo = halo_sizes(*parse_grid_sizes(grid_sizes_str)) - if states_str == "flat": - states = np.empty( - (ref.gas.nvars, nxhalo, nyhalo, nzhalo), - dtype=np.float64, order="F") - rho = 1.0 - vel = np.array([1.0, 1.0, 1.0]) - pressure = 1.6 - state = ref.gas.conserved(rho, vel, pressure) - for i in range(nxhalo): - for j in range(nyhalo): - for k in range(nzhalo): - states[:, i, j, k] = state - else: - raise ValueError("States string {} not supported".format(states_str)) - return states - - return grid_states - - -@pytest.fixture(scope="session") -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)) - if (metric_str == "uniform"): - metrics = np.full( - (ndim, ndim, nxhalo, nyhalo, nzhalo), - 0.0, dtype=np.float64, order="F") - for i in range(nxhalo): - for j in range(nyhalo): - for k in range(nzhalo): - for m in range(ndim): - metrics[m, m, i, j, k] = 1.0 - jacobians = np.full( - (nxhalo, nyhalo, nzhalo), - 1.0, dtype=np.float64, order="F") - else: - raise ValueError("Metric string {} not supported".format(metric_str)) - return metrics, jacobians - - return grid_metrics - - -@pytest.fixture(scope="session", params=[ - ("flat", "10x10x10", "x", "uniform"), - ("flat", "10x10x10", "y", "uniform"), - ("flat", "10x10x10", "z", "uniform") - # FIXME: add remaining cases - ]) -def grid_data(request, grid_states_generator, grid_metrics_generator): - states_str = request.param[0] - grid_sizes_str = request.param[1] - dir_str = request.param[2] - metric_str = request.param[3] - - return GridData(grid_states_generator(states_str, grid_sizes_str), - grid_sizes_str, - dir_str, - grid_metrics_generator(metric_str, grid_sizes_str, dir_str)) - - -class GridResults: - def __init__(self, queue, data): - self.queue = queue - self.nx = data.nx - self.ny = data.nx - self.nz = data.nx - self.nvars = data.nvars - self.ndim = data.ndim - self.output_dims = data.output_dims - self.dir_str = data.dir_str - self.direction = data.direction - self.dir_internal = data.dir_internal - - self.states = data.states - self.generalized_fluxes = data.generalized_fluxes - self.metrics = data.metrics - self.jacobians = data.jacobians - - self.generalized_flux_derivatives = np.full(self.output_dims, 0.0, - dtype=np.float64, order="F") - for i in range(self.nx): - for j in range(self.ny): - for k in range(self.nz): - flux_pos = self.local_flux( - *self.window_variables(i, j, k)) - neg_loc = ( - i - (1 if 0 == self.dir_internal else 0), - j - (1 if 1 == self.dir_internal else 0), - k - (1 if 2 == self.dir_internal else 0)) - flux_neg = self.local_flux( - *self.window_variables(*neg_loc)) - self.generalized_flux_derivatives[:, i, j, k] = ( - flux_pos - flux_neg) - - def window_variables(self, i, j, k): - def lo(n, d): - return n + halo - (2 if d == self.dir_internal else 0) - - def hi(n, d): - return n + halo + (4 if d == self.dir_internal else 1) - - xlo = lo(i, 0) - xhi = hi(i, 0) - ylo = lo(j, 1) - yhi = hi(j, 1) - zlo = lo(k, 2) - zhi = hi(k, 2) - local_states = self.states[:, xlo:xhi, ylo:yhi, zlo:zhi] - local_metrics = self.metrics[:, :, xlo:xhi, ylo:yhi, zlo:zhi] - local_jacobians = self.jacobians[xlo:xhi, ylo:yhi, zlo:zhi] - return (np.squeeze(local_states), - np.moveaxis(np.squeeze(local_metrics), -1, 0), - np.squeeze(local_jacobians)) - - def local_flux(self, state_window, metric_window, jacobian_window): - local_data = WindowData( - state_window, self.dir_str, - (metric_window, jacobian_window)) - return WindowResults(self.queue, local_data).weno_flux - - -@pytest.fixture(scope="session") -def grid_results(queue, grid_data): - return GridResults(queue, grid_data) def test_compute_flux_derivatives(queue, grid_results): diff --git a/test/test_flux_window_ops.py b/test/test_flux_window_ops.py index a79b187..367bb77 100644 --- a/test/test_flux_window_ops.py +++ b/test/test_flux_window_ops.py @@ -34,8 +34,6 @@ import logging import pytest import utilities as u -from data_for_test import (window_results, window_data, - window_states_generator, window_metrics_generator) def test_pointwise_eigenvalues(queue, window_results): -- GitLab From 076628ec7ceaa553d31e565ddf3105a03817af0b Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Thu, 1 Oct 2020 21:26:25 -0600 Subject: [PATCH 21/25] fix even more flake complaints --- test/data_for_test.py | 1 - test/test_flux_derivatives.py | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/test/data_for_test.py b/test/data_for_test.py index 7ef126a..2dad9f1 100644 --- a/test/data_for_test.py +++ b/test/data_for_test.py @@ -21,7 +21,6 @@ THE SOFTWARE. """ import numpy as np -import utilities as u import weno_reference_implementation as ref diff --git a/test/test_flux_derivatives.py b/test/test_flux_derivatives.py index 426570f..b52c59b 100644 --- a/test/test_flux_derivatives.py +++ b/test/test_flux_derivatives.py @@ -20,7 +20,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -import numpy as np +import numpy as np # noqa: F401 import numpy.linalg as la # noqa: F401 import pyopencl as cl # noqa: F401 import pyopencl.array # noqa -- GitLab From de1cc8bdf3b28ef8755a7ccaf630d5ee1ea7b58d Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Mon, 5 Oct 2020 19:59:45 -0600 Subject: [PATCH 22/25] add tests with a single jump and fix bug in Fortran code --- WENO.F90 | 299 ++++++++++++++++++++++++----------------------- test/conftest.py | 39 ++++++- 2 files changed, 184 insertions(+), 154 deletions(-) diff --git a/WENO.F90 b/WENO.F90 index 6e02d65..5509694 100644 --- a/WENO.F90 +++ b/WENO.F90 @@ -47,160 +47,163 @@ subroutine compute_flux_derivatives(nvars, ndim, nx, ny, nz, d, & ! for the next three loops, note carefully that the inner loop indices have a different range - !$loopy begin tagged: flux_x_compute - do k=1,nz - do j=1,ny - do i=0,nx - call convert_to_generalized_frozen(nvars, ndim, & - 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, & - combined_frozen_metrics, & - generalized_states_frozen, & - generalized_fluxes_frozen) - - 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) - - call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) - - call split_characteristic_fluxes(nvars, & - generalized_states_frozen, & - generalized_fluxes_frozen, & - R_inv, & - wavespeeds, & - characteristic_fluxes_pos, & - characteristic_fluxes_neg) - - call weno_flux(nvars, & - fluxes(:, i-2:i+3, j, k), & - characteristic_fluxes_pos, characteristic_fluxes_neg, & - combined_frozen_metrics(1), R, weno_flux_tmp(:, i, j, k)) + if (d == 1) then + !$loopy begin tagged: flux_x_compute + do k=1,nz + do j=1,ny + do i=0,nx + call convert_to_generalized_frozen(nvars, ndim, & + 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, & + combined_frozen_metrics, & + generalized_states_frozen, & + generalized_fluxes_frozen) + + 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) + + call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) + + call split_characteristic_fluxes(nvars, & + generalized_states_frozen, & + generalized_fluxes_frozen, & + R_inv, & + wavespeeds, & + characteristic_fluxes_pos, & + characteristic_fluxes_neg) + + call weno_flux(nvars, & + fluxes(:, i-2:i+3, j, k), & + characteristic_fluxes_pos, characteristic_fluxes_neg, & + combined_frozen_metrics(1), R, weno_flux_tmp(:, i, j, k)) + end do + end do end do - end do - end do - !$loopy end tagged: flux_x_compute - - !$loopy begin tagged: flux_x_diff - do k=1,nz - do j=1,ny - do i=1,nx - do v=1,nvars - flux_derivatives(v, i, j, k) = & - (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i-1, j, k)) - end do + !$loopy end tagged: flux_x_compute + + !$loopy begin tagged: flux_x_diff + do k=1,nz + do j=1,ny + do i=1,nx + do v=1,nvars + flux_derivatives(v, i, j, k) = & + (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i-1, j, k)) + end do + end do + end do end do - end do - end do - !$loopy end tagged: flux_x_diff - - !$loopy begin tagged: flux_y_compute - do k=1,nz - do i=1,nx - do j=0,ny - call convert_to_generalized_frozen(nvars, ndim, & - 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, & - combined_frozen_metrics, & - generalized_states_frozen, & - generalized_fluxes_frozen) - - 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) - - call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) - - call split_characteristic_fluxes(nvars, & - generalized_states_frozen, & - generalized_fluxes_frozen, & - R_inv, & - wavespeeds, & - characteristic_fluxes_pos, & - characteristic_fluxes_neg) - - call weno_flux(nvars, & - fluxes(:, i, j-2:j+3, k), & - characteristic_fluxes_pos, characteristic_fluxes_neg, & - combined_frozen_metrics(2), R, weno_flux_tmp(:, i, j, k)) - + !$loopy end tagged: flux_x_diff + + else if (d == 2) then + !$loopy begin tagged: flux_y_compute + do k=1,nz + do i=1,nx + do j=0,ny + call convert_to_generalized_frozen(nvars, ndim, & + 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, & + combined_frozen_metrics, & + generalized_states_frozen, & + generalized_fluxes_frozen) + + 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) + + call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) + + call split_characteristic_fluxes(nvars, & + generalized_states_frozen, & + generalized_fluxes_frozen, & + R_inv, & + wavespeeds, & + characteristic_fluxes_pos, & + characteristic_fluxes_neg) + + call weno_flux(nvars, & + fluxes(:, i, j-2:j+3, k), & + characteristic_fluxes_pos, characteristic_fluxes_neg, & + combined_frozen_metrics(2), R, weno_flux_tmp(:, i, j, k)) + + end do + end do end do - end do - end do - !$loopy end tagged: flux_y_compute - - !$loopy begin tagged: flux_y_diff - do k=1,nz - do j=1,ny - do i=1,nx - do v=1,nvars - flux_derivatives(v, i, j, k) = & - (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j-1, k)) - end do + !$loopy end tagged: flux_y_compute + + !$loopy begin tagged: flux_y_diff + do k=1,nz + do j=1,ny + do i=1,nx + do v=1,nvars + flux_derivatives(v, i, j, k) = & + (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j-1, k)) + end do + end do + end do end do - end do - end do - !$loopy end tagged: flux_y_diff - - !$loopy begin tagged: flux_z_compute - do j=1,ny - do i=1,nx - do k=0,nz - call convert_to_generalized_frozen(nvars, ndim, & - 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, & - combined_frozen_metrics, & - generalized_states_frozen, & - generalized_fluxes_frozen) - - 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) - - call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) - - call split_characteristic_fluxes(nvars, & - generalized_states_frozen, & - generalized_fluxes_frozen, & - R_inv, & - wavespeeds, & - characteristic_fluxes_pos, & - characteristic_fluxes_neg) - - call weno_flux(nvars, & - fluxes(:, i, j, k-2:k+3), & - characteristic_fluxes_pos, characteristic_fluxes_neg, & - combined_frozen_metrics(3), R, weno_flux_tmp(:, i, j, k)) + !$loopy end tagged: flux_y_diff + + else + !$loopy begin tagged: flux_z_compute + do j=1,ny + do i=1,nx + do k=0,nz + call convert_to_generalized_frozen(nvars, ndim, & + 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, & + combined_frozen_metrics, & + generalized_states_frozen, & + generalized_fluxes_frozen) + + 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) + + call lax_wavespeeds(nvars, lambda_pointwise, lambda_roe, wavespeeds) + + call split_characteristic_fluxes(nvars, & + generalized_states_frozen, & + generalized_fluxes_frozen, & + R_inv, & + wavespeeds, & + characteristic_fluxes_pos, & + characteristic_fluxes_neg) + + call weno_flux(nvars, & + fluxes(:, i, j, k-2:k+3), & + characteristic_fluxes_pos, characteristic_fluxes_neg, & + combined_frozen_metrics(3), R, weno_flux_tmp(:, i, j, k)) + end do + end do end do - end do - end do - !$loopy end tagged: flux_z_compute - - !$loopy begin tagged: flux_z_diff - do k=1,nz - do j=1,ny - do i=1,nx - do v=1,nvars - flux_derivatives(v, i, j, k) = & - (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j, k-1)) - flux_derivatives(v, i, j, k) = 0.0d0 - end do + !$loopy end tagged: flux_z_compute + + !$loopy begin tagged: flux_z_diff + do k=1,nz + do j=1,ny + do i=1,nx + do v=1,nvars + flux_derivatives(v, i, j, k) = & + (weno_flux_tmp(v, i, j, k) - weno_flux_tmp(v, i, j, k-1)) + end do + end do + end do end do - end do - end do - !$loopy end tagged: flux_z_diff + !$loopy end tagged: flux_z_diff + end if end subroutine subroutine pointwise_eigenvalues(nvars, ndim, d, states, metrics, lambda_pointwise) diff --git a/test/conftest.py b/test/conftest.py index 5fb6967..355de79 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -193,12 +193,12 @@ def window_results(queue, window_data): @pytest.fixture(scope="session") def grid_states_generator(): - def grid_states(states_str, grid_sizes_str): + def grid_states(states_str, grid_sizes_str, dir_str): nxhalo, nyhalo, nzhalo = halo_sizes(*parse_grid_sizes(grid_sizes_str)) + states = np.empty( + (ref.gas.nvars, nxhalo, nyhalo, nzhalo), + dtype=np.float64, order="F") if states_str == "flat": - states = np.empty( - (ref.gas.nvars, nxhalo, nyhalo, nzhalo), - dtype=np.float64, order="F") rho = 1.0 vel = np.array([1.0, 1.0, 1.0]) pressure = 1.6 @@ -207,6 +207,30 @@ def grid_states_generator(): for j in range(nyhalo): for k in range(nzhalo): states[:, i, j, k] = state + elif states_str == "single_jump": + rho_left = 2.0 + vel_left = np.array([2.0, 2.0, 2.0]) + pressure_left = 3.2 + state_left = ref.gas.conserved(rho_left, vel_left, pressure_left) + + rho_right = 1.0 + vel_right = np.array([1.0, 1.0, 1.0]) + pressure_right = 1.6 + state_right = ref.gas.conserved(rho_right, vel_right, pressure_right) + + def is_right(i, j, k): + right_x = (dir_str == "x" and i > nxhalo/2) + right_y = (dir_str == "y" and j > nyhalo/2) + right_z = (dir_str == "z" and k > nzhalo/2) + return right_x or right_y or right_z + + for i in range(nxhalo): + for j in range(nyhalo): + for k in range(nzhalo): + if is_right(i, j, k): + states[:, i, j, k] = state_right + else: + states[:, i, j, k] = state_left else: raise ValueError("States string {} not supported".format(states_str)) return states @@ -242,7 +266,10 @@ def grid_metrics_generator(): @pytest.fixture(scope="session", params=[ ("flat", "10x10x10", "x", "uniform"), ("flat", "10x10x10", "y", "uniform"), - ("flat", "10x10x10", "z", "uniform") + ("flat", "10x10x10", "z", "uniform"), + ("single_jump", "10x10x10", "x", "uniform"), + ("single_jump", "10x10x10", "y", "uniform"), + ("single_jump", "10x10x10", "z", "uniform"), # FIXME: add remaining cases ]) def grid_data(request, grid_states_generator, grid_metrics_generator): @@ -251,7 +278,7 @@ def grid_data(request, grid_states_generator, grid_metrics_generator): dir_str = request.param[2] metric_str = request.param[3] - return GridData(grid_states_generator(states_str, grid_sizes_str), + return GridData(grid_states_generator(states_str, grid_sizes_str, dir_str), grid_sizes_str, dir_str, grid_metrics_generator(metric_str, grid_sizes_str, dir_str)) -- GitLab From 33dbab98fc3f7c49a27cfc6a3c5cebe6f0d5d448 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Mon, 19 Oct 2020 22:28:43 -0600 Subject: [PATCH 23/25] Add tests with uneven dimensions and fix test setup bug this exposed --- test/conftest.py | 3 +++ test/data_for_test.py | 4 ++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/test/conftest.py b/test/conftest.py index 355de79..2e6c728 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -270,6 +270,9 @@ def grid_metrics_generator(): ("single_jump", "10x10x10", "x", "uniform"), ("single_jump", "10x10x10", "y", "uniform"), ("single_jump", "10x10x10", "z", "uniform"), + ("single_jump", "12x14x16", "x", "uniform"), + ("single_jump", "12x14x16", "y", "uniform"), + ("single_jump", "12x14x16", "z", "uniform"), # FIXME: add remaining cases ]) def grid_data(request, grid_states_generator, grid_metrics_generator): diff --git a/test/data_for_test.py b/test/data_for_test.py index 2dad9f1..2cf9c2f 100644 --- a/test/data_for_test.py +++ b/test/data_for_test.py @@ -162,8 +162,8 @@ class GridResults: def __init__(self, queue, data): self.queue = queue self.nx = data.nx - self.ny = data.nx - self.nz = data.nx + self.ny = data.ny + self.nz = data.nz self.nvars = data.nvars self.ndim = data.ndim self.output_dims = data.output_dims -- GitLab From 039fa75bebaf032e98b8b312b5bd7b2763998147 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Mon, 5 Oct 2020 20:57:50 -0600 Subject: [PATCH 24/25] Added a jump whose strength varies in x,y,z --- test/conftest.py | 48 +++++++++++++++++++++++++++++++---------------- test/utilities.py | 2 +- 2 files changed, 33 insertions(+), 17 deletions(-) diff --git a/test/conftest.py b/test/conftest.py index 2e6c728..10c58eb 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -198,15 +198,28 @@ def grid_states_generator(): states = np.empty( (ref.gas.nvars, nxhalo, nyhalo, nzhalo), dtype=np.float64, order="F") + + def is_right(i, j, k): + right_x = (dir_str == "x" and i >= nxhalo/2) + right_y = (dir_str == "y" and j >= nyhalo/2) + right_z = (dir_str == "z" and k >= nzhalo/2) + return right_x or right_y or right_z + + def set_states_jump(state_left, state_right): + for i in range(nxhalo): + for j in range(nyhalo): + for k in range(nzhalo): + if is_right(i, j, k): + states[:, i, j, k] = state_right + else: + states[:, i, j, k] = state_left + if states_str == "flat": rho = 1.0 vel = np.array([1.0, 1.0, 1.0]) pressure = 1.6 state = ref.gas.conserved(rho, vel, pressure) - for i in range(nxhalo): - for j in range(nyhalo): - for k in range(nzhalo): - states[:, i, j, k] = state + set_states_jump(state, state) elif states_str == "single_jump": rho_left = 2.0 vel_left = np.array([2.0, 2.0, 2.0]) @@ -218,19 +231,19 @@ def grid_states_generator(): pressure_right = 1.6 state_right = ref.gas.conserved(rho_right, vel_right, pressure_right) - def is_right(i, j, k): - right_x = (dir_str == "x" and i > nxhalo/2) - right_y = (dir_str == "y" and j > nyhalo/2) - right_z = (dir_str == "z" and k > nzhalo/2) - return right_x or right_y or right_z + set_states_jump(state_left, state_right) + elif states_str == "single_jump_varying": + rho_left = 2.0 + vel_left = np.array([2.0, 4.0, 6.0]) + pressure_left = 3.2 + state_left = ref.gas.conserved(rho_left, vel_left, pressure_left) - for i in range(nxhalo): - for j in range(nyhalo): - for k in range(nzhalo): - if is_right(i, j, k): - states[:, i, j, k] = state_right - else: - states[:, i, j, k] = state_left + rho_right = 1.0 + vel_right = np.array([1.0, 2.0, 3.0]) + pressure_right = 1.6 + state_right = ref.gas.conserved(rho_right, vel_right, pressure_right) + + set_states_jump(state_left, state_right) else: raise ValueError("States string {} not supported".format(states_str)) return states @@ -273,6 +286,9 @@ def grid_metrics_generator(): ("single_jump", "12x14x16", "x", "uniform"), ("single_jump", "12x14x16", "y", "uniform"), ("single_jump", "12x14x16", "z", "uniform"), + ("single_jump_varying", "10x10x10", "x", "uniform"), + ("single_jump_varying", "10x10x10", "y", "uniform"), + ("single_jump_varying", "10x10x10", "z", "uniform"), # FIXME: add remaining cases ]) def grid_data(request, grid_states_generator, grid_metrics_generator): diff --git a/test/utilities.py b/test/utilities.py index 3f1fa07..793a407 100644 --- a/test/utilities.py +++ b/test/utilities.py @@ -35,7 +35,7 @@ import os.path # {{{ arrays def compare_arrays(a, b): - assert a == approx(b, rel=1e-12, abs=1e-14) + assert a == approx(b, rel=1e-12, abs=1e-13) def transposed_array(a): -- GitLab From b671af7678369a23b5b27e0dba8c55b9cf30e3e3 Mon Sep 17 00:00:00 2001 From: "Timothy A. Smith" <2timothy18@gmail.com> Date: Tue, 20 Oct 2020 22:02:06 -0600 Subject: [PATCH 25/25] add tests with a double hump --- test/conftest.py | 47 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 46 insertions(+), 1 deletion(-) diff --git a/test/conftest.py b/test/conftest.py index 10c58eb..994bf91 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -194,7 +194,9 @@ def window_results(queue, window_data): @pytest.fixture(scope="session") def grid_states_generator(): def grid_states(states_str, grid_sizes_str, dir_str): - 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 = (nxhalo - nx)/2 states = np.empty( (ref.gas.nvars, nxhalo, nyhalo, nzhalo), dtype=np.float64, order="F") @@ -214,6 +216,29 @@ def grid_states_generator(): else: states[:, i, j, k] = state_left + def is_left_hump(i, j, k): + hump_x = (dir_str == "x" and i >= nx/5 + halo and i <= 2*nx/5 + halo) + hump_y = (dir_str == "y" and j >= ny/5 + halo and j <= 2*ny/5 + halo) + hump_z = (dir_str == "z" and k >= nz/5 + halo and k <= 2*nz/5 + halo) + return hump_x or hump_y or hump_z + + def is_right_hump(i, j, k): + hump_x = (dir_str == "x" and i >= 3*nx/5 + halo and i <= 4*nx/5 + halo) + hump_y = (dir_str == "y" and j >= 3*ny/5 + halo and j <= 4*ny/5 + halo) + hump_z = (dir_str == "z" and k >= 3*nz/5 + halo and k <= 4*nz/5 + halo) + return hump_x or hump_y or hump_z + + def set_states_double_hump(state_base, state_left, state_right): + for i in range(nxhalo): + for j in range(nyhalo): + for k in range(nzhalo): + if is_left_hump(i, j, k): + states[:, i, j, k] = state_left + elif is_right_hump(i, j, k): + states[:, i, j, k] = state_right + else: + states[:, i, j, k] = state_base + if states_str == "flat": rho = 1.0 vel = np.array([1.0, 1.0, 1.0]) @@ -244,6 +269,23 @@ def grid_states_generator(): state_right = ref.gas.conserved(rho_right, vel_right, pressure_right) set_states_jump(state_left, state_right) + elif states_str == "double_hump": + rho_base = 1.0 + vel_base = np.array([1.0, 1.0, 1.0]) + pressure_base = 1.6 + state_base = ref.gas.conserved(rho_base, vel_base, pressure_base) + + rho_left = 2.0 + vel_left = np.array([2.0, 2.0, 2.0]) + pressure_left = 3.2 + state_left = ref.gas.conserved(rho_left, vel_left, pressure_left) + + rho_right = 3.0 + vel_right = np.array([3.0, 3.0, 3.0]) + pressure_right = 4.8 + state_right = ref.gas.conserved(rho_right, vel_right, pressure_right) + + set_states_double_hump(state_base, state_left, state_right) else: raise ValueError("States string {} not supported".format(states_str)) return states @@ -289,6 +331,9 @@ def grid_metrics_generator(): ("single_jump_varying", "10x10x10", "x", "uniform"), ("single_jump_varying", "10x10x10", "y", "uniform"), ("single_jump_varying", "10x10x10", "z", "uniform"), + ("double_hump", "25x25x25", "x", "uniform"), + ("double_hump", "25x25x25", "y", "uniform"), + ("double_hump", "25x25x25", "z", "uniform"), # FIXME: add remaining cases ]) def grid_data(request, grid_states_generator, grid_metrics_generator): -- GitLab