From a0643f951bf9427ee01fbe3131416448bd87bc50 Mon Sep 17 00:00:00 2001 From: Tiago Pestana Date: Wed, 16 Sep 2020 20:06:11 +0200 Subject: [PATCH] bugfix: fixed size of the readable portion of chunk['data'] (now it works also when no ghost cells are present) --- python/ibmppp/ibmppp.py | 64 ++++++++++++++++++++++------------------- 1 file changed, 34 insertions(+), 30 deletions(-) diff --git a/python/ibmppp/ibmppp.py b/python/ibmppp/ibmppp.py index b634fc6..2bdd1b0 100644 --- a/python/ibmppp/ibmppp.py +++ b/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]=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('\n') fid.write('\n') fid.write(' \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('\n') fid.write('\n') fid.write(' \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('\n') fid.write('\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