Fixing generate_xdmf because paraview now requires 2D data to be 3D

Author: mikaem

mpi4py_fft/__init__.py

@@ -22,4 +22,5 @@
 from .distarray import DistArray, newDistArray, Function
 from .mpifft import PFFT
 from . import fftw
+from .fftw import fftlib
 from .io import HDF5File, NCFile, generate_xdmf

mpi4py_fft/distarray.py

@@ -56,9 +56,9 @@ class DistArray(np.ndarray):
 
     """
     def __new__(cls, global_shape, subcomm=None, val=None, dtype=float,
-                buffer=None, alignment=None, rank=0):
+                buffer=None, strides=None, alignment=None, rank=0):
         if len(global_shape[rank:]) < 2: # 1D case
-            obj = np.ndarray.__new__(cls, global_shape, dtype=dtype, buffer=buffer)
+            obj = np.ndarray.__new__(cls, global_shape, dtype=dtype, buffer=buffer, strides=strides)
             if buffer is None and isinstance(val, Number):
                 obj.fill(val)
             obj._rank = rank

mpi4py_fft/fftw/__init__.py

@@ -1,3 +1,3 @@
 from .xfftn import *
 from .factory import get_planned_FFT, export_wisdom, import_wisdom, \
-    forget_wisdom, cleanup, set_timelimit, get_fftw_lib
+    forget_wisdom, cleanup, set_timelimit, get_fftw_lib, fftlib

mpi4py_fft/io/generate_xdmf.py

@@ -38,13 +38,16 @@ def get_geometry(kind=0, dim=2):
           </DataItem>
         </Geometry>"""
 
-        return """<Geometry Type="VXVY">
+        return """<Geometry Type="VXVYVZ">
           <DataItem Format="HDF" NumberType="Float" Precision="{0}" Dimensions="{1}">
             {3}:{6}/mesh/{4}
           </DataItem>
           <DataItem Format="HDF" NumberType="Float" Precision="{0}" Dimensions="{2}">
             {3}:{6}/mesh/{5}
           </DataItem>
+          <DataItem Format="XML" NumberType="Float" Precision="8" Dimensions="1">
+            0
+          </DataItem>
         </Geometry>"""
 
     if dim == 3:
@@ -74,7 +77,7 @@ def get_topology(dims, kind=0):
     assert len(dims) in (2, 3)
     co = 'Co' if kind == 0 else ''
     if len(dims) == 2:
-        return """<Topology Dimensions="{0} {1}" Type="2D{2}RectMesh"/>""".format(dims[0], dims[1], co)
+        return """<Topology Dimensions="1 {0} {1}" Type="3D{2}RectMesh"/>""".format(dims[0], dims[1], co)
     if len(dims) == 3:
         return """<Topology Dimensions="{0} {1} {2}" Type="3D{3}RectMesh"/>""".format(dims[0], dims[1], dims[2], co)
 
@@ -83,7 +86,7 @@ def get_attribute(attr, h5filename, dims, prec):
     assert len(dims) in (2, 3)
     if len(dims) == 2:
         return """<Attribute Name="{0}" Center="Node">
-          <DataItem Format="HDF" NumberType="Float" Precision="{5}" Dimensions="{1} {2}">
+          <DataItem Format="HDF" NumberType="Float" Precision="{5}" Dimensions="1 {1} {2}">
             {3}:/{4}
           </DataItem>
         </Attribute>
@@ -96,7 +99,7 @@ def get_attribute(attr, h5filename, dims, prec):
         </Attribute>
         """.format(name, dims[0], dims[1], dims[2], h5filename, attr, prec)
 
-def generate_xdmf(h5filename, periodic=True, order='paraview'):
+def generate_xdmf(h5filename, periodic=True, order='visit'):
     """Generate XDMF-files
 
     Parameters

mpi4py_fft/mpifft.py

@@ -363,7 +363,7 @@ def shape(self, forward_output=True):
         """
         if forward_output is not True:
             return self.forward.input_pencil.subshape
-        return self.backward.input_pencil.subshape
+        return self.forward.output_array.shape
 
     def local_slice(self, forward_output=True):
         """The local view into the global data

tests/test_io.py

@@ -21,6 +21,12 @@ def remove_if_exists(filename):
     except OSError:
         pass
 
+def cleanup():
+    import glob
+    files = glob.glob('*.h5')+glob.glob('*.xdmf')+glob.glob('*.nc')
+    for f in files:
+        remove_if_exists(f)
+
 def test_2D(backend, forward_output):
     if backend == 'netcdf4':
         assert forward_output is False
@@ -162,6 +168,7 @@ def test_4D(backend, forward_output):
 
 if __name__ == '__main__':
     #pylint: disable=unused-import
+    cleanup()
     skip = {'hdf5': False, 'netcdf4': False}
     try:
         import h5py
@@ -181,3 +188,4 @@ def test_4D(backend, forward_output):
                 test_2D(bnd, kind)
                 if bnd == 'hdf5':
                     test_4D(bnd, kind)
+    cleanup()

tests/test_speed.py

@@ -3,11 +3,14 @@
 import pyfftw
 import scipy.fftpack as sp
 from mpi4py_fft import fftw
+import pickle
 
 try:
-    fftw.import_wisdom('wisdom.dat')
-except AssertionError:
-    pass
+    #fftw.import_wisdom('wisdom.dat')
+    pyfftw.import_wisdom(pickle.load(open('pyfftw.wisdom', 'rb')))
+    print('Wisdom imported')
+except:
+    print('Wisdom not imported')
 
 N = (64, 64, 64)
 loops = 50