Merge pull request clawpack#528 from aymkhalil/pyclaw_hdf5_parallel_write

Pyclaw hdf5 parallel read/write

Author: ketch

.travis.yml

@@ -29,11 +29,18 @@ script:
   - nosetests --first-pkg-wins --with-doctest --exclude=limiters --exclude=sharpclaw --exclude=io --exclude=example --with-coverage --cover-package=clawpack.pyclaw
   - mv .coverage temp
   - if [[ "${TEST_PACKAGE}" == "pyclaw" ]]; then
-        nosetests -v --first-pkg-wins --exclude=limiters --exclude=sharpclaw --exclude=io --with-coverage --cover-package=clawpack.pyclaw;
+        nosetests -v --first-pkg-wins --exclude=limiters --exclude=sharpclaw --with-coverage --cover-package=clawpack.pyclaw --include=IOTest;
     fi
   - if [[ "${TEST_PACKAGE}" == "petclaw" ]]; then
        mpirun -n 4 nosetests -v --first-pkg-wins --exclude=limiters --exclude=sharpclaw --exclude=io;
     fi
+  - if [[ "${TEST_PACKAGE}" == "petclaw" ]]; then
+       cd ../petclaw/tests;
+       mpirun -n 4 nosetests -v --first-pkg-wins --include=IOTest;
+    fi
+  - if [[ "${TEST_PACKAGE}" == "petclaw" ]]; then
+       cd ../../pyclaw;
+    fi
   - mv temp .coverage.doctest
   - coverage combine
 

examples/acoustics_2d_variable/test_acoustics_2d_variable_io.py

@@ -58,7 +58,7 @@ def verify_acoustics_io(controller):
     tempdir = './_io_test_results'
     classic_tests = gen_variants(acoustics_2d_interface.setup, 
                                  verify_acoustics_io, solver_type='classic', 
-                                 outdir=tempdir, num_cells=(50, 50))
+                                 outdir=tempdir, num_cells=(50, 50), disable_petsc=True)
 
     import shutil
     from itertools import chain

src/petclaw/geometry.py

@@ -90,4 +90,10 @@ def __init__(self,geom):
         elif isinstance(geom[0],pyclaw_geometry.Dimension):
             self.patches = [Patch(geom)]
 
+class Dimension(pyclaw_geometry.Dimension):
+    def __init__(self, lower, upper, num_cells, name='x',
+                 on_lower_boundary=None,on_upper_boundary=None, units=None):
+        super(Dimension, self).__init__(lower, upper, num_cells, name,
+                 on_lower_boundary,on_upper_boundary, units)
+
 

src/petclaw/io/__init__.py

@@ -15,7 +15,7 @@
 # Check for HDF 5 support
 try:
     import h5py
-    from clawpack.pyclaw.io import hdf5
+    import hdf5
     __all__ += ['hdf5.read','hdf5.write']
 except:
     logging.debug("No hdf5 support found.")

src/petclaw/io/hdf5.py

@@ -0,0 +1,221 @@
+#!/usr/bin/env python
+# encoding: utf-8
+r"""
+Routines for reading and writing a HDF5 output file
+
+This module reads and writes hdf5 files via the following module:
+    h5py - http://code.google.com/p/h5py/
+
+To install h5py, you must also install the hdf5 library from the website:
+    http://www.hdfgroup.org/HDF5/release/obtain5.html
+"""
+
+from mpi4py import MPI
+import os
+import logging
+
+from clawpack.petclaw import geometry
+from clawpack import petclaw
+
+logger = logging.getLogger('pyclaw.io')
+
+try:
+    import h5py
+except:
+    logging.critical("Could not import h5py!")
+    error_msg = ("Could not import h5py, please install " +
+                 "either h5py.  See the doc_string for more " +
+                 "information.")
+    raise Exception(error_msg)
+
+def write(solution,frame,path,file_prefix='claw',write_aux=False,
+                options={},write_p=False):
+    r"""
+    Write out a Solution to a HDF5 file.
+    
+    :Input:
+     - *solution* - (:class:`~pyclaw.solution.Solution`) Pyclaw solution 
+       object to input into
+     - *frame* - (int) Frame number
+     - *path* - (string) Root path
+     - *file_prefix* - (string) Prefix for the file name.  ``default = 'claw'``
+     - *write_aux* - (bool) Boolean controlling whether the associated 
+       auxiliary array should be written out.  ``default = False``     
+     - *options* - (dict) Optional argument dictionary, see 
+       `HDF5 Option Table`_
+    
+    .. _`HDF5 Option Table`:
+    
+    +-----------------+------------------------------------------------------+
+    | Key             | Value                                                |
+    +=================+======================================================+
+    | compression     | (None, string ["gzip" | "lzf" | "szip"] or int 0-9)  |
+    |                 | Enable dataset compression. DEFLATE, LZF and (where  |
+    |                 | available) SZIP are supported. An integer is         |
+    |                 | interpreted as a GZIP level for backwards            |
+    |                 | compatibility.                                       |
+    +-----------------+------------------------------------------------------+
+    |compression_opts | (None, or special value) Setting for compression     |
+    |                 | filter; legal values for each filter type are:       |
+    |                 |                                                      |
+    |                 | - *gzip* - (int) 0-9                                 |
+    |                 | - *lzf* - None allowed                               |
+    |                 | - *szip* - (tuple) 2-tuple ('ec'|'nn', even integer  |
+    |                 |     0-32)                                            |
+    |                 |                                                      |
+    |                 | See the filters module for a detailed description of |
+    |                 | each of these filters.                               |
+    +-----------------+------------------------------------------------------+
+    | chunks          | (None, True or shape tuple) Store the dataset in     |
+    |                 | chunked format. Automatically selected if any of the |
+    |                 | other keyword options are given. If you don't provide|
+    |                 | a shape tuple, the library will guess one for you.   |
+    +-----------------+------------------------------------------------------+
+    | shuffle         | (True/False) Enable/disable data shuffling, which can|
+    |                 | improve compression performance. Automatically       |
+    |                 | enabled when compression is used.                    |
+    +-----------------+------------------------------------------------------+
+    | fletcher32      | (True/False) Enable Fletcher32 error detection; may  |
+    |                 | be used with or without compression.                 |
+    +-----------------+------------------------------------------------------+
+    """
+    option_defaults = {'compression':None,'compression_opts':None,
+                       'chunks':None,'shuffle':False,'fletcher32':False}
+    for (k,v) in option_defaults.iteritems():
+        options[k] = options.get(k,v)
+    
+    filename = os.path.join(path,'%s%s.hdf' % 
+                                (file_prefix,str(frame).zfill(4)))
+
+    if options['compression'] is not None:
+        err_msg = "Compression (filters) are not available for parallel h5py yet."
+        logging.critical(err_msg)
+        raise Exception(err_msg)
+
+    with h5py.File(filename,'w',driver='mpio',comm=MPI.COMM_WORLD) as f:
+        # For each patch, write out attributes
+        for state in solution.states:
+            patch = state.patch
+            # Create group for this patch
+            subgroup = f.create_group('patch%s' % patch.patch_index)
+
+            # General patch properties
+            subgroup.attrs['t'] = state.t
+            subgroup.attrs['num_eqn'] = state.num_eqn
+            subgroup.attrs['num_aux'] = state.num_aux
+            for attr in ['num_ghost','patch_index','level']:
+                if hasattr(patch,attr):
+                    if getattr(patch,attr) is not None:
+                        subgroup.attrs[attr] = getattr(patch,attr)
+
+            # Add the dimension names as a attribute
+            subgroup.attrs['dimensions'] = patch.get_dim_attribute('name')
+            # Dimension properties
+            for dim in patch.dimensions:
+                for attr in ['num_cells','lower','delta','upper',
+                             'units']:
+                    if hasattr(dim,attr):
+                        if getattr(dim,attr) is not None:
+                            attr_name = '%s.%s' % (dim.name,attr)
+                            subgroup.attrs[attr_name] = getattr(dim,attr)
+
+            if write_p:
+                q = state.p
+            else:
+                q = state.q
+            r = patch._da.getRanges()
+            globalSize = []
+            globalSize.append(q.shape[0])
+            globalSize.extend(patch.num_cells_global)
+            dset = subgroup.create_dataset('q',globalSize,dtype='float',**options)
+            to_hdf5_dataset(q, dset, len(patch.dimensions), r)
+
+            if write_aux and state.num_aux > 0:
+                r = patch._da.getRanges()
+                globalSize = []
+                globalSize.append(state.num_aux)
+                globalSize.extend(patch.num_cells_global)
+                dset = subgroup.create_dataset('aux',globalSize,dtype='float',**options)
+                to_hdf5_dataset(state.aux, dset, len(patch.dimensions), r)
+
+
+def read(solution,frame,path='./',file_prefix='claw',read_aux=True,
+                options={}):
+    r"""
+    Read in a HDF5 file into a Solution
+
+    :Input:
+     - *solution* - (:class:`~pyclaw.solution.Solution`) Pyclaw object to be
+       output
+     - *frame* - (int) Frame number
+     - *path* - (string) Root path
+     - *file_prefix* - (string) Prefix for the file name.  ``default = 'claw'``
+     - *write_aux* - (bool) Boolean controlling whether the associated
+       auxiliary array should be written out.  ``default = False``
+     - *options* - (dict) Optional argument dictionary, not used for reading.
+    """
+    filename = os.path.join(path,'%s%s.hdf' %
+                                (file_prefix,str(frame).zfill(4)))
+    patches = []
+
+    with h5py.File(filename,'r',driver='mpio',comm=MPI.COMM_WORLD) as f:
+        for patch in f.itervalues():
+            # Construct each dimension
+            dimensions = []
+            dim_names = patch.attrs['dimensions']
+            for dim_name in dim_names:
+                dim = geometry.Dimension(
+                    patch.attrs["%s.lower" % dim_name],
+                    patch.attrs["%s.upper" % dim_name],
+                    patch.attrs["%s.num_cells" % dim_name],
+                    name = dim_name)
+                # Optional attributes
+                for attr in ['units']:
+                    attr_name = "%s.%s" % (dim_name,attr)
+                    if patch.attrs.get(attr_name, None):
+                        setattr(dim,attr,patch.attrs["%s.%s" % (dim_name,attr)])
+                dimensions.append(dim)
+
+            pyclaw_patch = petclaw.Patch(dimensions)
+
+            # Fetch general patch properties
+            for attr in ['t','num_eqn','patch_index','level']:
+                setattr(pyclaw_patch,attr,patch.attrs[attr])
+
+            state = petclaw.state.State(pyclaw_patch, \
+                                        patch.attrs['num_eqn'],patch.attrs['num_aux'])
+            state.t = patch.attrs['t']
+
+            globalSize = []
+            globalSize.append(state.q.shape[0])
+            globalSize.extend(pyclaw_patch.num_cells_global)
+            r = pyclaw_patch._da.getRanges()
+
+            dset = patch['q'][:].reshape(globalSize)
+            state.q = from_hdf5_dataset(dset, len(pyclaw_patch.dimensions), r, state.q.shape)
+
+            # Read in aux if applicable
+            if read_aux and patch.get('aux',None) is not None:
+                dset = patch['aux'][:]
+                state.aux = from_hdf5_dataset(dset, len(pyclaw_patch.dimensions), r, state.aux.shape)
+
+            solution.states.append(state)
+            patches.append(pyclaw_patch)
+
+        solution.domain = geometry.Domain(patches)
+
+def to_hdf5_dataset(arr, dset, ndim, ranges):
+    if ndim == 1:
+        dset[:,ranges[0][0]:ranges[0][1]] = arr
+    elif ndim == 2:
+        dset[:,ranges[0][0]:ranges[0][1],ranges[1][0]:ranges[1][1]] = arr
+    elif ndim == 3:
+        dset[:,ranges[0][0]:ranges[0][1],ranges[1][0]:ranges[1][1],ranges[2][0]:ranges[2][1]] = arr
+
+def from_hdf5_dataset(dset, ndim, ranges, shape):
+    if ndim == 1:
+        return dset[:,ranges[0][0]:ranges[0][1]].reshape(shape,order='F')
+    elif ndim == 2:
+        return dset[:,ranges[0][0]:ranges[0][1],ranges[1][0]:ranges[1][1]].reshape(shape,order='F')
+    elif ndim == 3:
+        return dset[:,ranges[0][0]:ranges[0][1],ranges[1][0]:ranges[1][1],ranges[2][0]:ranges[2][1]].reshape(shape,order='F')

src/petclaw/solution.py

@@ -5,3 +5,23 @@ class Solution(Solution):
     """ Parallel Solution class.
     """
     __doc__ += pyclaw.util.add_parent_doc(pyclaw.Solution)
+
+    def get_read_func(self, file_format):
+        from clawpack.petclaw import io
+        if file_format == 'petsc':
+            return io.petsc.read
+        elif file_format == 'hdf5':
+            return io.hdf5.read
+        else:
+            raise ValueError("File format %s not supported." % file_format)
+
+    def get_write_func(self, file_format):
+        from clawpack.petclaw import io
+        if 'petsc' in file_format:
+            return io.petsc.write
+        elif 'hdf5' in file_format:
+            return io.hdf5.write
+        else:
+            raise ValueError("File format %s not supported." % file_format)
+
+

src/petclaw/tests/test_io.py

@@ -0,0 +1,23 @@
+from clawpack import pyclaw
+from clawpack import petclaw
+import os
+
+class PetClawIOTest(pyclaw.IOTest):
+    @property
+    def solution(self):
+        return petclaw.Solution()
+
+    @property
+    def file_formats(self):
+        return ['hdf5']
+
+    @property
+    def this_dir(self):
+        return os.path.dirname(os.path.abspath(__file__))
+
+    @property
+    def test_data_dir(self):
+        return os.path.join(self.this_dir, '../../pyclaw/tests/test_data')
+
+    def test_io_from_binary(self):
+        return

src/pyclaw/__init__.py

@@ -34,3 +34,5 @@
 
 __all__.append('BC')
 from .solver import BC
+__all__.extend('IOTest')
+from .tests.test_io import IOTest

src/pyclaw/setup.py

@@ -10,6 +10,7 @@ def configuration(parent_package='',top_path=None):
     config.add_subpackage('io')
     config.add_subpackage('limiters')
     config.add_subpackage('examples')
+    config.add_subpackage('tests')
     return config
 
 if __name__ == '__main__':