bugfix: fixed size of the readable portion of chunk['data'] (now it works also when no ghost cells are present)

.gitignore

@@ -0,0 +1 @@
+*__pycache__

bash/ucftar_pack

@@ -72,7 +72,7 @@ set -- "${POSITIONAL[@]}"
 
 # Parse input filename
 din=$1
-seqnum=$2 
+seqnum=$2
 
 # Construct output filename if not specified
 if [ -z "$fout" ]; then
@@ -106,7 +106,8 @@ if [ ! -s ${din}/${fproc} ]; then
 fi
 if [ ! -s ${din}/${fpart} ]; then
   (>&2 echo "[Error] File not found or empty: ${din}/${fpart}")
-  exit 1
+  # exit 1
+  fpart=
 fi
 
 # Check if all velocity fields exist and are not empty
@@ -174,7 +175,7 @@ while [ $packloop -eq 1 ]; do
     snumfilecur=0
     scursize=1024   # Initialize with size of trailing zero bloc
     unset flist
-    while true; do 
+    while true; do
       sfilesize=$(wc -c ${din}/${flist_sbase[${scounter}]} | awk '{print $1}') # raw size
       sfilesize=$(((${sfilesize}+511)/512*512+512))                            # adjust to 512 byte-blocks and add header
       spredictsize=$((${scursize}+${sfilesize}))
@@ -198,8 +199,8 @@ while [ $packloop -eq 1 ]; do
         fi
         flag_splitdone=1
         break
-      fi 
-    done 
+      fi
+    done
   fi
   # Create tar archive
   tar $flagtar --directory=${din} -cf ${fout} ${flist[@]}
@@ -211,7 +212,7 @@ while [ $packloop -eq 1 ]; do
   # Compare checksums (CNC32), if flag is set
   flistx=($(echo ${flist[@]} | sed s/"_$seqnum"/""/g))
   if [ $checksum -eq 1 ]; then
-    for ii in "${!flistx[@]}"; do 
+    for ii in "${!flistx[@]}"; do
       if [ $verbose -eq 1 ]; then
         (>&2 echo "Verifying checksum: ${flist[$ii]}")
       fi

python/ibmppp/__init__.py

@@ -0,0 +1 @@
+from .ibmppp import *

python/ibmppp/ibmppp.py

@@ -146,8 +146,8 @@ class ibmppp:
           self.__localChunkSize[key][1] = self.__localChunkBounds[key][3] - self.__localChunkBounds[key][2] + 1
           self.__localChunkSize[key][2] = self.__localChunkBounds[key][5] - self.__localChunkBounds[key][4] + 1
           # Verify that local grid size is not smaller than ghost cell size
-          if (self.__localChunkSize[key][0]<self.__nghx or 
-              self.__localChunkSize[key][1]<self.__nghy or 
+          if (self.__localChunkSize[key][0]<self.__nghx or
+              self.__localChunkSize[key][1]<self.__nghy or
               self.__localChunkSize[key][2]<self.__nghz):
             raise ValueError('Local grid size must be greater than number of ghost cells in each direction!')
         # Initialize neighbor array
@@ -191,17 +191,17 @@ class ibmppp:
             (np.arange(-self.__nghx,0)*self.__dx[key]+self.grid[key][0][self.__localChunkBounds[key][0]-1],  # lower ghost cells
              self.grid[key][0][self.__localChunkBounds[key][0]-1:self.__localChunkBounds[key][1]],           # chunk grid
              np.arange(1,self.__nghx+1)*self.__dx[key]+self.grid[key][0][self.__localChunkBounds[key][1]-1]) # upper ghost cells
-            ) 
+            )
           self.localGrid[key][1] = np.concatenate(
             (np.arange(-self.__nghy,0)*self.__dx[key]+self.grid[key][1][self.__localChunkBounds[key][2]-1],  # lower ghost cells
              self.grid[key][1][self.__localChunkBounds[key][2]-1:self.__localChunkBounds[key][3]],           # chunk grid
              np.arange(1,self.__nghy+1)*self.__dx[key]+self.grid[key][1][self.__localChunkBounds[key][3]-1]) # upper ghost cells
-            ) 
+            )
           self.localGrid[key][2] = np.concatenate(
             (np.arange(-self.__nghz,0)*self.__dx[key]+self.grid[key][2][self.__localChunkBounds[key][4]-1],  # lower ghost cells
              self.grid[key][2][self.__localChunkBounds[key][4]-1:self.__localChunkBounds[key][5]],           # chunk grid
              np.arange(1,self.__nghz+1)*self.__dx[key]+self.grid[key][2][self.__localChunkBounds[key][5]-1]) # upper ghost cells
-            ) 
+            )
 
     def loadParticles(self):
         '''Reads particles from file and distributes them to the correct processor.'''
@@ -260,9 +260,9 @@ class ibmppp:
                 # Find particles which affect that chunk
                 li_part =  ((pp[col['x'],:]>=bdmin[0]-pp[col['r'],:]) &
                             (pp[col['x'],:]<=bdmax[0]+pp[col['r'],:]) &
-                            (pp[col['y'],:]>=bdmin[1]-pp[col['r'],:]) & 
+                            (pp[col['y'],:]>=bdmin[1]-pp[col['r'],:]) &
                             (pp[col['y'],:]<=bdmax[1]+pp[col['r'],:]) &
-                            (pp[col['z'],:]>=bdmin[2]-pp[col['r'],:]) & 
+                            (pp[col['z'],:]>=bdmin[2]-pp[col['r'],:]) &
                             (pp[col['z'],:]<=bdmax[2]+pp[col['r'],:]))
                 # Send them to that processor
                 sendbuf = (np.ascontiguousarray(pp[:,li_part]),col)
@@ -291,19 +291,19 @@ class ibmppp:
         # Also determine file which contains the desired field as well as the
         # dataset ID
         if key[0]=='u':
-          filebase = self.__dir_base+'uvwp.'
+          filebase = self.__dir_base+'/uvwp.'
           dset = 1
         elif key[0]=='v':
-          filebase = self.__dir_base+'uvwp.'
+          filebase = self.__dir_base+'/uvwp.'
           dset = 2
         elif key[0]=='w':
-          filebase = self.__dir_base+'uvwp.'
+          filebase = self.__dir_base+'/uvwp.'
           dset = 3
         elif key[0]=='p':
-          filebase = self.__dir_base+'uvwp.'
+          filebase = self.__dir_base+'/uvwp.'
           dset = 4
         elif key[0]=='s':
-          filebase = self.__dir_base+'scal.'
+          filebase = self.__dir_base+'/scal.'
           dset = int(key[1])
           if dset!=1:
             self.fields[key] = None
@@ -348,10 +348,14 @@ class ibmppp:
               if kb1!=chunk['kbeg'] or ke1!=chunk['kbeg']+chunk['nzl']-1:
                   raise ValueError('Loaded chunk does not correspond to processor grid (kbeg={}, nzl={} vs kbeg={},nzl={}'.format(kb1,ke1-kb1+1,chunk['kbeg'],chunk['nzl']))
               ngh1 = chunk['ighost']
+              ibeg = [i + ngh1 for i in [0,0,0]]
+              iend = [i-ngh1 for i in chunk['data'].shape]
               # No need to care about ghost cells here, we communicate them later
               self.field[key][ib1-ib+self.__nghx:ie1-ib+self.__nghx+1,
                               jb1-jb+self.__nghy:je1-jb+self.__nghy+1,
-                              kb1-kb+self.__nghz:ke1-kb+self.__nghz+1] = chunk['data'][ngh1:-ngh1,ngh1:-ngh1,ngh1:-ngh1]
+                              kb1-kb+self.__nghz:ke1-kb+self.__nghz+1] = chunk['data'][ibeg[0]:iend[0],\
+                                                                                       ibeg[1]:iend[1],\
+                                                                                       ibeg[2]:iend[2]]
         # Exchange ghost cells and set boundary conditions
         self.exchangeGhostCells(key)
         self.imposeBoundaryConditions(key)
@@ -449,7 +453,7 @@ class ibmppp:
         # Construct XDMF filename
         file_xdmf = filename+'_'+key+'.xdmf'
         # Open file and write header
-        fid = open(file_xdmf,'w') 
+        fid = open(file_xdmf,'w')
         fid.write('<?xml version="1.0"?>\n')
         fid.write('<Xdmf Version="2.0">\n')
         fid.write('   <Domain>\n')
@@ -575,7 +579,7 @@ class ibmppp:
         # Construct XDMF filename
         file_xdmf = filename+'_'+key+'.xdmf'
         # Open file and write header
-        fid = open(file_xdmf,'w') 
+        fid = open(file_xdmf,'w')
         fid.write('<?xml version="1.0"?>\n')
         fid.write('<Xdmf Version="2.0">\n')
         fid.write('   <Domain>\n')
@@ -688,7 +692,7 @@ class ibmppp:
             # Construct XDMF filename
             file_xdmf = filename+'.xdmf'
             # Open file
-            fid = open(file_xdmf,'w') 
+            fid = open(file_xdmf,'w')
             # Write header
             fid.write('<?xml version="1.0"?>\n')
             fid.write('<Xdmf Version="2.0">\n')
@@ -832,7 +836,7 @@ class ibmppp:
         self.differentiateField('u','x')
       symm = self.__deriveSymmetry(self.__boundarySymmetries['ux'],'mult',symm2=self.__boundarySymmetries['ux'])
       self.__allocateField('tmp','ux',shift=(0,0,0),symmetry=symm)
-      self.field['tmp'] += 0.5*np.square(self.field['ux']) 
+      self.field['tmp'] += 0.5*np.square(self.field['ux'])
       self.shiftField('tmp',self.__getShiftDirection('tmp','Q'))
       iin,jin,kin,iout,jout,kout = self.__sliceCollocated('tmp','Q')
       self.field['Q'][iout,jout,kout] -= self.field['tmp'][iin,jin,kin]
@@ -844,7 +848,7 @@ class ibmppp:
         self.differentiateField('v','y')
       symm = self.__deriveSymmetry(self.__boundarySymmetries['vy'],'mult',symm2=self.__boundarySymmetries['vy'])
       self.__allocateField('tmp','vy',shift=(0,0,0),symmetry=symm)
-      self.field['tmp'] += 0.5*np.square(self.field['vy']) 
+      self.field['tmp'] += 0.5*np.square(self.field['vy'])
       self.shiftField('tmp',self.__getShiftDirection('tmp','Q'))
       iin,jin,kin,iout,jout,kout = self.__sliceCollocated('tmp','Q')
       self.field['Q'][iout,jout,kout] -= self.field['tmp'][iin,jin,kin]
@@ -856,7 +860,7 @@ class ibmppp:
         self.differentiateField('w','z')
       symm = self.__deriveSymmetry(self.__boundarySymmetries['wz'],'mult',symm2=self.__boundarySymmetries['wz'])
       self.__allocateField('tmp','wz',shift=(0,0,0),symmetry=symm)
-      self.field['tmp'] += 0.5*np.square(self.field['wz']) 
+      self.field['tmp'] += 0.5*np.square(self.field['wz'])
       self.shiftField('tmp',self.__getShiftDirection('tmp','Q'))
       iin,jin,kin,iout,jout,kout = self.__sliceCollocated('tmp','Q')
       self.field['Q'][iout,jout,kout] -= self.field['tmp'][iin,jin,kin]
@@ -981,7 +985,7 @@ class ibmppp:
         self.differentiateField('u','x')
       symm = self.__deriveSymmetry(self.__boundarySymmetries['ux'],'mult',symm2=self.__boundarySymmetries['ux'])
       self.__allocateField('tmp','ux',shift=(0,0,0),symmetry=symm)
-      self.field['tmp'] += 2.0*np.square(self.field['ux']) 
+      self.field['tmp'] += 2.0*np.square(self.field['ux'])
       self.shiftField('tmp',self.__getShiftDirection('tmp','D'))
       iin,jin,kin,iout,jout,kout = self.__sliceCollocated('tmp','D')
       self.field['D'][iout,jout,kout] += self.field['tmp'][iin,jin,kin]
@@ -993,7 +997,7 @@ class ibmppp:
         self.differentiateField('v','y')
       symm = self.__deriveSymmetry(self.__boundarySymmetries['vy'],'mult',symm2=self.__boundarySymmetries['vy'])
       self.__allocateField('tmp','vy',shift=(0,0,0),symmetry=symm)
-      self.field['tmp'] += 2.0*np.square(self.field['vy']) 
+      self.field['tmp'] += 2.0*np.square(self.field['vy'])
       self.shiftField('tmp',self.__getShiftDirection('tmp','D'))
       iin,jin,kin,iout,jout,kout = self.__sliceCollocated('tmp','D')
       self.field['D'][iout,jout,kout] += self.field['tmp'][iin,jin,kin]
@@ -1005,7 +1009,7 @@ class ibmppp:
         self.differentiateField('w','z')
       symm = self.__deriveSymmetry(self.__boundarySymmetries['wz'],'mult',symm2=self.__boundarySymmetries['wz'])
       self.__allocateField('tmp','wz',shift=(0,0,0),symmetry=symm)
-      self.field['tmp'] += 2.0*np.square(self.field['wz']) 
+      self.field['tmp'] += 2.0*np.square(self.field['wz'])
       self.shiftField('tmp',self.__getShiftDirection('tmp','D'))
       iin,jin,kin,iout,jout,kout = self.__sliceCollocated('tmp','D')
       self.field['D'][iout,jout,kout] += self.field['tmp'][iin,jin,kin]
@@ -1306,7 +1310,7 @@ class ibmppp:
           self.grid[key][idir]      += shift[idir]*0.5*self.__dx[key]
           self.localGrid[key][idir] += shift[idir]*0.5*self.__dx[key]
       # Perform the interpolation
-      if subsequent: 
+      if subsequent:
         # Three (at maximum) subsequent 1D interpolations.
         # NOTE: This algorithm requires ghost cell communication for each substep, which can be avoided if we shifted in all directions simultaneously!
         # Construct index range of target
@@ -1452,14 +1456,14 @@ class ibmppp:
       self.__symmetrizeWall(key,(0,+1,0))
       self.__symmetrizeWall(key,(0,0,-1))
       self.__symmetrizeWall(key,(0,0,+1))
-      
+
     def printMemoryUsage(self):
       '''Print estimation of memory usage (max over all processors)'''
       mem = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024.
       mem = self.__comm.reduce(sendobj=mem,op=MPI.MAX,root=0)
       if self.__rank==0:
         print('Current memory usage: {:8.2f} MiB'.format(mem))
-    
+
     def printMemoryUsageField(self,key):
       '''Print estimation of memory usage per allocated field.'''
       bytesNumpy = np.dtype(self.__precision).itemsize
@@ -1478,8 +1482,8 @@ class ibmppp:
       #   v  -> -v            v  -> -v
       #   w  -> -w            w  ->  w
       #   p  ->  p            p  ->  p
-      self.__boundarySymmetries = {'u':  np.zeros((3,3,3),dtype='i'), 
-                                   'v':  np.zeros((3,3,3),dtype='i'), 
+      self.__boundarySymmetries = {'u':  np.zeros((3,3,3),dtype='i'),
+                                   'v':  np.zeros((3,3,3),dtype='i'),
                                    'w':  np.zeros((3,3,3),dtype='i'),
                                    'p':  np.zeros((3,3,3),dtype='i'),
                                    's1': np.zeros((3,3,3),dtype='i'),
@@ -1679,7 +1683,7 @@ class ibmppp:
       rank_dst = self.__nghbr[ip_dst,jp_dst,kp_dst]
       rank_src = self.__nghbr[ip_src,jp_src,kp_src]
       # [send/recv] create a tag
-      tag = ip_dst*100+jp_dst*10+kp_dst 
+      tag = ip_dst*100+jp_dst*10+kp_dst
       # [send/recv] get the indices of data to be sent/received
       (nxl,nyl,nzl) = self.__localChunkSize[key]
       if positionDst[0]==-1:
@@ -1819,7 +1823,7 @@ class ibmppp:
     @staticmethod
     def __rankFromPosition(ip,jp,kp,nxp,nyp,nzp):
       return ip*nyp*nzp+jp*nzp+kp
-  
+
     @staticmethod
     def __positionFromRank(rank,nxp,nyp,nzp):
       ip = rank//(nyp*nzp)
@@ -1851,7 +1855,7 @@ class ibmppp:
 def gaussianFilterStencilWidth(sigma=(1,1,1),truncate=4.0):
   if len(sigma)!=3:
         raise ValueError('Expected one value of sigma for each direction, but only got {}.'.format(len(sigma)))
-  # Compute stencil width 
+  # Compute stencil width
   stencil = [0,0,0]
   for ii in range(0,3):
     stencil[ii] = 2*int(truncate*sigma[ii]+0.5)+1

python/ucf/__init__.py

@@ -0,0 +1 @@
+from .ucf import *