diff --git a/python/ibmppp/ibmppp.py b/python/ibmppp/ibmppp.py index 104a8e9..b634fc6 100644 --- a/python/ibmppp/ibmppp.py +++ b/python/ibmppp/ibmppp.py @@ -356,6 +356,35 @@ class ibmppp: self.exchangeGhostCells(key) self.imposeBoundaryConditions(key) + def loadStatistics(self,filename,key): + '''Loads statistics from a h5 file.''' + if self.__rank==0: + # Construct file name + file_h5 = filename+'.h5' + # Open the file and write data + fid = h5py.File(file_h5,'r') + gid = fid[key] + # Load mean + did = gid['mean']; + fcont = did.attrs['F_CONTIGUOUS'] + dims = did.attrs['DIM'] + order = 'F' if fcont else 'C' + self.spatialMean[key] = np.reshape(did,dims,order=order) + # Load meansquare + did = gid['meansquare']; + fcont = did.attrs['F_CONTIGUOUS'] + dims = did.attrs['DIM'] + order = 'F' if fcont else 'C' + self.spatialMsqr[key] = np.reshape(did,dims,order=order) + # Load nsample + did = gid['nsample']; + fcont = did.attrs['F_CONTIGUOUS'] + dims = did.attrs['DIM'] + order = 'F' if fcont else 'C' + self.spatialNsample[key] = np.reshape(did,dims,order=order) + # Close the file + fid.close() + def saveField(self,filename,key,append=False,keepAllGhost=False,xdmf=False): '''Writes chunks in new processor layout to h5 files''' # Determine filename of output file @@ -940,6 +969,110 @@ class ibmppp: if collocate: self.shiftField(keyout,self.__getShiftDirection(keyout,'p')) + def dissipation(self,keepDerivatives=False): + '''Computes dissipation rate on pressure grid.''' + # D = 2 [dudx^2 + dvdy^2 + dwdz^2] + (dudy+dvdx)^2 + (dudz+dwdx)^2 + (dvdz+dwdy)^2 + # Viscosity is not multiplied! + # Allocate a new output field on pressure grid + self.__allocateField('D','p',shift=(0,0,0)) + # Compute derivatives on after another and add to field + # dudx*dudx + if 'ux' not in self.field: + 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.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] + self.freeField('tmp') + if not keepDerivatives: + self.freeField('ux') + # dvdy*dvdy + if 'vy' not in self.field: + 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.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] + self.freeField('tmp') + if not keepDerivatives: + self.freeField('vy') + # dwdz*dwdz + if 'wz' not in self.field: + 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.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] + self.freeField('tmp') + if not keepDerivatives: + self.freeField('wz') + # (dudy+dvdx)^2 + if 'uy' not in self.field: + self.differentiateField('u','y') + if 'vx' not in self.field: + self.differentiateField('v','x') + symm = self.__deriveSymmetry(self.__boundarySymmetries['uy'],'mult',symm2=self.__boundarySymmetries['vx']) + self.__allocateField('tmp','D',shift=(-1,-1,0),symmetry=symm) + ii1,jj1,kk1,iout,jout,kout = self.__sliceCollocated('uy','tmp') + ii2,jj2,kk2,iout,jout,kout = self.__sliceCollocated('vx','tmp') + self.field['tmp'][iout,jout,kout] += np.square(self.field['uy'][ii1,jj1,kk1]+self.field['vx'][ii2,jj2,kk2]) + self.exchangeGhostCells('tmp') + self.imposeBoundaryConditions('tmp') + 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] + self.freeField('tmp') + if not keepDerivatives: + self.freeField('uy') + self.freeField('vx') + # (dudz+dwdx)^2 + if 'uz' not in self.field: + self.differentiateField('u','z') + if 'wx' not in self.field: + self.differentiateField('w','x') + symm = self.__deriveSymmetry(self.__boundarySymmetries['uz'],'mult',symm2=self.__boundarySymmetries['wx']) + self.__allocateField('tmp','D',shift=(-1,0,-1),symmetry=symm) + ii1,jj1,kk1,iout,jout,kout = self.__sliceCollocated('uz','tmp') + ii2,jj2,kk2,iout,jout,kout = self.__sliceCollocated('wx','tmp') + self.field['tmp'][iout,jout,kout] += np.square(self.field['uz'][ii1,jj1,kk1]+self.field['wx'][ii2,jj2,kk2]) + self.exchangeGhostCells('tmp') + self.imposeBoundaryConditions('tmp') + 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] + self.freeField('tmp') + if not keepDerivatives: + self.freeField('uz') + self.freeField('wx') + # (dvdz+dwdy)^2 + if 'vz' not in self.field: + self.differentiateField('v','z') + if 'wy' not in self.field: + self.differentiateField('w','y') + symm = self.__deriveSymmetry(self.__boundarySymmetries['vz'],'mult',symm2=self.__boundarySymmetries['wy']) + self.__allocateField('tmp','D',shift=(0,-1,-1),symmetry=symm) + ii1,jj1,kk1,iout,jout,kout = self.__sliceCollocated('vz','tmp') + ii2,jj2,kk2,iout,jout,kout = self.__sliceCollocated('wy','tmp') + self.field['tmp'][iout,jout,kout] += np.square(self.field['vz'][ii1,jj1,kk1]+self.field['wy'][ii2,jj2,kk2]) + self.exchangeGhostCells('tmp') + self.imposeBoundaryConditions('tmp') + 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] + self.freeField('tmp') + if not keepDerivatives: + self.freeField('vz') + self.freeField('wy') + # Exchange ghost cells and set boundary conditions + self.exchangeGhostCells('D') + self.imposeBoundaryConditions('D') + def spatialStatistics(self,key,homogeneous=None): '''Computes spatial mean and meansquare over homogeneous directions.''' # If no homogeneous directions are provided, use periodicity @@ -1064,6 +1197,22 @@ class ibmppp: self.__comm.Send(tmpsum, dest=0,tag=2) self.__comm.Send(tmpssq, dest=0,tag=3) + def subtractMean(self,key): + '''Subtract mean from field.''' + # Check if this operation can be executed + if self.__rank==0: + if not key in spatialMean or spatialMean[key] is None: + raise ValueError('Statistics not available') + # Broadcast statistics + self.spatialMean[key] = self.__comm.bcast(self.spatialMean[key],root=0) + # Get dimensions + dims = self.spatialMean[key].shape + # Slice chunk + if dims[0]==1: + slx = slice(0,1) + else: + slx = slice(self.ib) + def collocateVelocity(self): '''Moves velocity components onto the pressure grid.''' raise NotImplementedError('TBD')