added hdf5 and netcdf write/read to petsc + usage example on 2d acous…

…tics
clawpack · Apr 30, 2014 · ec82777 · ec82777
1 parent 6b907f4
commit ec82777
Show file tree

Hide file tree

Showing 3 changed files with 46 additions and 8 deletions.
diff --git a/examples/acoustics_2d_variable/acoustics_2d_interface.py b/examples/acoustics_2d_variable/acoustics_2d_interface.py
@@ -88,6 +88,8 @@ def setup(kernel_language='Fortran',use_petsc=False,outdir='./_output',solver_ty
     claw.num_output_times = 20
     claw.write_aux_init = True
     claw.setplot = setplot
+    if use_petsc:
+        claw.output_options = {'format':'binary'}
 
     # Solve
     claw.tfinal = 0.6

diff --git a/examples/acoustics_2d_variable/test_acoustics_2d_variable_io.py b/examples/acoustics_2d_variable/test_acoustics_2d_variable_io.py
@@ -28,13 +28,15 @@ def verify_acoustics_io(controller):
         sol_0_test = Solution()
         sol_0_test.read(0,path=controller.outdir,
                         file_format=controller.output_format,
-                        file_prefix=None,read_aux=True)
+                        file_prefix=None,read_aux=True,
+                        options=controller.output_options)
         test_aux = sol_0_test.state.get_aux_global()
 
         sol_20_test = Solution()
         sol_20_test.read(20,path=controller.outdir,
                         file_format=controller.output_format,
-                        file_prefix=None,read_aux=False)
+                        file_prefix=None,read_aux=False,
+                        options=controller.output_options)
         test_q = sol_20_test.state.get_q_global()
 
 

diff --git a/src/petclaw/io/petsc.py b/src/petclaw/io/petsc.py
@@ -75,7 +75,8 @@ def write(solution,frame,path='./',file_prefix='claw',write_aux=False,
                          'num_aux':solution.num_aux,'num_dim':solution.domain.num_dim,
                          'write_aux':write_aux,
                          'problem_data' : solution.problem_data,
-                         'mapc2p': solution.state.grid.mapc2p}
+                         'mapc2p': solution.state.grid.mapc2p,
+                         'file_format':file_format}
         if write_p:
             sol_dict['num_eqn'] = solution.mp
 
@@ -99,7 +100,15 @@ def write(solution,frame,path='./',file_prefix='claw',write_aux=False,
     elif file_format == 'vtk':
         viewer = PETSc.Viewer().createASCII(viewer_filename, PETSc.Viewer.Mode.WRITE, format=PETSc.Viewer.Format.ASCII_VTK)
         if write_aux:
-            aux_viewer = PETSc.Viewer().createASCII(aux_filename, PETSc.Viewer.Mode.WRITE) 
+            aux_viewer = PETSc.Viewer().createASCII(aux_filename, PETSc.Viewer.Mode.WRITE)
+    elif file_format=='hdf5':
+        viewer = PETSc.Viewer().createHDF5(viewer_filename, PETSc.Viewer.Mode.WRITE)
+        if write_aux:
+            aux_viewer = PETSc.Viewer().createHDF5(aux_filename, PETSc.Viewer.Mode.WRITE)
+    elif file_format=='netcdf':
+        viewer = PETSc.Viewer().createHDF5(viewer_filename, PETSc.Viewer.Mode.WRITE)
+        if write_aux:
+            aux_viewer = PETSc.Viewer().createHDF5(aux_filename, PETSc.Viewer.Mode.WRITE)
     else:
         raise IOError('format type %s not supported' % file_format)
 
@@ -179,11 +188,17 @@ def read(solution,frame,path='./',file_prefix='claw',read_aux=False,options={}):
     num_aux       = value_dict['num_aux']
     num_eqn       = value_dict['num_eqn']
 
+    aux_warning = False
+
     # now set up the PETSc viewer
     if file_format == 'ascii':
         viewer = PETSc.Viewer().createASCII(viewer_filename, PETSc.Viewer.Mode.READ)
         if read_aux:
-            aux_viewer = PETSc.Viewer().createASCII(aux_viewer_filename, PETSc.Viewer.Mode.READ)
+            if os.path.exists(aux_viewer_filename):
+                aux_viewer = PETSc.Viewer().createASCII(aux_viewer_filename, PETSc.Viewer.Mode.READ)
+            else:
+                aux_warning = True
+                read_aux = False
     elif file_format == 'binary':
         if hasattr(PETSc.Viewer,'createMPIIO'):
             viewer = PETSc.Viewer().createMPIIO(viewer_filename, PETSc.Viewer.Mode.READ)
@@ -196,13 +211,32 @@ def read(solution,frame,path='./',file_prefix='claw',read_aux=False,options={}):
                 else:
                     aux_viewer = PETSc.Viewer().createBinary(aux_viewer_filename, PETSc.Viewer.Mode.READ)
             else:
-                from warnings import warn
-                aux_file_path = os.path.join(path,aux_viewer_filename)
-                warn('read_aux=True but aux file %s does not exist' % aux_file_path)
+                aux_warning = True
+                read_aux = False
+    elif file_format == 'hdf5':
+        viewer = PETSc.Viewer().createHDF5(viewer_filename, PETSc.Viewer.Mode.READ)
+        if read_aux:
+            if os.path.exists(aux_viewer_filename):
+                aux_viewer = PETSc.Viewer().createHDF5(aux_viewer_filename, PETSc.Viewer.Mode.READ)
+            else:         
+                aux_warning = True
+                read_aux=False
+    elif file_format == 'netcdf':
+        viewer = PETSc.Viewer().createNetCDF(viewer_filename, PETSc.Viewer.Mode.READ)
+        if read_aux:
+            if os.path.exists(aux_viewer_filename):
+                aux_viewer = PETSc.Viewer().createNetCDF(aux_viewer_filename, PETSc.Viewer.Mode.READ)
+            else:
+                aux_warning = True
                 read_aux=False
     else:
         raise IOError('format type %s not supported' % file_format)
 
+    if aux_warning:
+        from warnings import warn
+        aux_file_path = os.path.join(path,aux_viewer_filename)
+        warn('read_aux=True but aux file %s does not exist' % aux_file_path)
+
     patches = []
     for m in xrange(nstates):
         patch_dict = pickle.load(pickle_file)