Added statistics routines
This commit is contained in:
parent
fcca37ad93
commit
b2edd57f8b
|
|
@ -356,6 +356,35 @@ class ibmppp:
|
||||||
self.exchangeGhostCells(key)
|
self.exchangeGhostCells(key)
|
||||||
self.imposeBoundaryConditions(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):
|
def saveField(self,filename,key,append=False,keepAllGhost=False,xdmf=False):
|
||||||
'''Writes chunks in new processor layout to h5 files'''
|
'''Writes chunks in new processor layout to h5 files'''
|
||||||
# Determine filename of output file
|
# Determine filename of output file
|
||||||
|
|
@ -940,6 +969,110 @@ class ibmppp:
|
||||||
if collocate:
|
if collocate:
|
||||||
self.shiftField(keyout,self.__getShiftDirection(keyout,'p'))
|
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):
|
def spatialStatistics(self,key,homogeneous=None):
|
||||||
'''Computes spatial mean and meansquare over homogeneous directions.'''
|
'''Computes spatial mean and meansquare over homogeneous directions.'''
|
||||||
# If no homogeneous directions are provided, use periodicity
|
# 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(tmpsum, dest=0,tag=2)
|
||||||
self.__comm.Send(tmpssq, dest=0,tag=3)
|
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):
|
def collocateVelocity(self):
|
||||||
'''Moves velocity components onto the pressure grid.'''
|
'''Moves velocity components onto the pressure grid.'''
|
||||||
raise NotImplementedError('TBD')
|
raise NotImplementedError('TBD')
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue