From d451c503fcefa29910a342e929d41369331df83e Mon Sep 17 00:00:00 2001 From: Ayman Khalil Date: Tue, 23 Feb 2016 12:09:21 +0300 Subject: [PATCH] Add test cases for solution downsampling --- src/petclaw/solution.py | 17 ++++++++++++++++ src/pyclaw/controller.py | 1 + src/pyclaw/solution.py | 20 ++++++++++++++++++- src/pyclaw/tests/test_io.py | 39 +++++++++++++++++++++++++++++++++++-- src/pyclaw/util.py | 22 +++++++++++++++++++++ 5 files changed, 96 insertions(+), 3 deletions(-) diff --git a/src/petclaw/solution.py b/src/petclaw/solution.py index 0e3cdd1ba..e3bc716d5 100644 --- a/src/petclaw/solution.py +++ b/src/petclaw/solution.py @@ -38,6 +38,7 @@ def _init_ds_solution(self): for i in range(self.domain.num_dim)], proc_sizes=self.state.q_da.getProcSizes()) ds_state = self._init_ds_state(self.state) self._ds_solution = pyclaw.Solution(ds_state, ds_domain) + self._ds_solution.t = self.t def _init_ds_state(self, state): """ @@ -92,6 +93,22 @@ def _downsample_global_to_local(self, da_for_ds, gVec, num_eqn, ds_domain_ranges """ from skimage.transform import downscale_local_mean + """ + downscale_local_mean downsamples n-dimensional array by local averaging. + First, it views the array as blocks of the downsampling factors, then it computes the local average of each block. + + Examples + -------- + >>> a = np.arange(15).reshape(3, 5) + >>> a + array([[ 0, 1, 2, 3, 4], + [ 5, 6, 7, 8, 9], + [10, 11, 12, 13, 14]]) + >>> downscale_local_mean(a, (2, 3)) + array([[ 3.5, 4. ], + [ 5.5, 4.5]]) + + """ # Create local array with ghost cells q_local_vec = da_for_ds.createLocalVec() diff --git a/src/pyclaw/controller.py b/src/pyclaw/controller.py index 603bfe1e7..75dc1266a 100644 --- a/src/pyclaw/controller.py +++ b/src/pyclaw/controller.py @@ -194,6 +194,7 @@ def __init__(self): self.F_file_name = 'F' r"""(string) - Name of text file containing functionals""" self.downsampling_factors = None + r"""(tuple) - A tuple of factors in each grid dimension that will be used in downsampling the solution by local averaging""" # ========== Access methods =============================================== def __str__(self): diff --git a/src/pyclaw/solution.py b/src/pyclaw/solution.py index 161fcac87..f0ce84512 100755 --- a/src/pyclaw/solution.py +++ b/src/pyclaw/solution.py @@ -388,6 +388,7 @@ def _init_ds_solution(self): domain = Domain(self.domain.grid.lower,self.domain.grid.upper,self.domain.grid.num_cells/np.array(self.downsampling_factors)) state = State(domain,self.state.num_eqn,self.state.num_aux) self._ds_solution = Solution(state, domain) + self._ds_solution.t = self.t def _init_ds_state(self, state): """ @@ -402,7 +403,24 @@ def downsample(self, write_aux,write_p): """ Returns a downsampled version of the solution by local averaging over the downsampling factors """ + from skimage.transform import downscale_local_mean + """ + downscale_local_mean downsamples n-dimensional array by local averaging. + First, it views the array as blocks of the downsampling factors, then it computes the local average of each block. + + Examples + -------- + >>> a = np.arange(15).reshape(3, 5) + >>> a + array([[ 0, 1, 2, 3, 4], + [ 5, 6, 7, 8, 9], + [10, 11, 12, 13, 14]]) + >>> downscale_local_mean(a, (2, 3)) + array([[ 3.5, 4. ], + [ 5.5, 4.5]]) + + """ for i in range(len(self.states)): state = self.states[i] @@ -416,7 +434,7 @@ def downsample(self, write_aux,write_p): else: ds_state.q = downscale_local_mean(state.q, (1,) + self.downsampling_factors) - # Dowsample aux + # Downsample aux if write_aux: ds_state.aux = downscale_local_mean(state.aux, (1,) + self.downsampling_factors) diff --git a/src/pyclaw/tests/test_io.py b/src/pyclaw/tests/test_io.py index 1703569a6..e7012bad0 100644 --- a/src/pyclaw/tests/test_io.py +++ b/src/pyclaw/tests/test_io.py @@ -1,7 +1,7 @@ import os import numpy as np from clawpack.pyclaw import Solution -from clawpack.pyclaw.util import check_solutions_are_same +from clawpack.pyclaw.util import check_solutions_are_same, check_ds_sol_is_downsampled_from_a_sol class IOTest(): @property @@ -33,6 +33,14 @@ def test_io_from_hdf5_with_aux(self): regression_dir = os.path.join(self.test_data_dir,'./advection_2d_with_aux') self.read_write_and_compare(self.file_formats,regression_dir,'hdf5',0,aux=True) + def test_ds_from_hdf5(self): + regression_dir = os.path.join(self.test_data_dir,'./Sedov_regression_hdf') + self.read_downsample_and_compare(self.file_formats,regression_dir,'hdf5',1,(2,2,2)) + + def test_ds_from_hdf5_with_aux(self): + regression_dir = os.path.join(self.test_data_dir,'./advection_2d_with_aux') + self.read_downsample_and_compare(self.file_formats,regression_dir,'hdf5',0,(2,2),aux=True) + def read_write_and_compare(self, file_formats,regression_dir,regression_format,frame_num,aux=False): r"""Test IO file formats: - Reading in an HDF file @@ -59,4 +67,31 @@ def read_write_and_compare(self, file_formats,regression_dir,regression_format,f # Compare solutions # Probably better to do this by defining __eq__ for each class for fmt, sol in s.iteritems(): - check_solutions_are_same(sol,ref_sol) \ No newline at end of file + check_solutions_are_same(sol,ref_sol) + + def read_downsample_and_compare(self, file_formats,regression_dir,regression_format,frame_num,downsampling_factors,aux=False): + r"""Test downsampling: + - Read in a solution from file + - Downsample the solution + - Check that q & aux arrays are correctly downsampled + """ + a_sol = self.solution + a_sol.read(frame_num,path=regression_dir,file_format=regression_format,read_aux=aux) + if aux: + assert (a_sol.state.aux is not None) + + # Write solution file in each format + io_test_dir = os.path.join(self.this_dir,'./io_test') + for fmt in file_formats: + a_sol.downsampling_factors = downsampling_factors + a_sol.downsample(aux, False).write(frame_num,path=io_test_dir,file_format=fmt,write_aux=aux) + + # Read solutions back in + s = {} + for fmt in file_formats: + s[fmt] = self.solution + s[fmt].read(frame_num,path=io_test_dir,file_format=fmt,write_aux=aux) + + # Compare solutions + for fmt, ds_sol in s.iteritems(): + check_ds_sol_is_downsampled_from_a_sol(ds_sol, a_sol) \ No newline at end of file diff --git a/src/pyclaw/util.py b/src/pyclaw/util.py index 43d61b66c..21fba960d 100755 --- a/src/pyclaw/util.py +++ b/src/pyclaw/util.py @@ -319,6 +319,28 @@ def check_solutions_are_same(sol_a,sol_b): for attr in ['units','on_lower_boundary','on_upper_boundary']: if hasattr(ref_dim,attr): assert getattr(dim,attr) == getattr(ref_dim,attr) + +def check_ds_sol_is_downsampled_from_a_sol(sol_ds, sol_a): + from skimage.transform import downscale_local_mean + + assert len(sol_a.states) == len(sol_ds.states) + assert sol_a.t == sol_ds.t + for state in sol_a.states: + for ds_state in sol_ds.states: + if ds_state.patch.patch_index == state.patch.patch_index: + break + + # Required state attributes + assert np.linalg.norm(downscale_local_mean(state.q, (1,) + sol_a.downsampling_factors) - ds_state.q) < 1.e-6 # Not sure why this can be so large + if state.aux is not None: + assert np.linalg.norm(downscale_local_mean(state.aux, (1,) + sol_a.downsampling_factors) - ds_state.aux) < 1.e-16 + for attr in ['t', 'num_eqn', 'num_aux']: + assert getattr(state,attr) == getattr(ds_state,attr) + # Optional state attributes + for attr in ['patch_index', 'level']: + if hasattr(ds_state,attr): + assert getattr(state,attr) == getattr(ds_state,attr) + # ============================================================================ # F2PY Utility Functions # ============================================================================