From b2edd57f8b2ece6521d08270e3d2e3250419fa55 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Wed, 16 Sep 2020 11:55:11 +0200 Subject: [PATCH 01/12] Added statistics routines --- python/ibmppp/ibmppp.py | 149 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 149 insertions(+) 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') From 1aba520ad8e60f85f488eb61dea1bc08674b5dda Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Fri, 11 Dec 2020 14:29:10 +0100 Subject: [PATCH 02/12] update --- matlab/@ucf/ucf.m | 20 +- matlab/@ustar/ustar.m | 8 +- matlab/find_chunk.m | 18 ++ matlab/find_region_ucfmf.m | 219 ++++++++++++++++++++++ matlab/read_field_block_ucfmf.m | 108 +++++++++++ matlab/read_field_complete_ucf.m | 56 ++++++ matlab/read_grid_ucf.m | 2 +- matlab/read_procgrid_ucf.m | 2 +- matlab/read_statistics_channel_scal_ucf.m | 18 +- 9 files changed, 437 insertions(+), 14 deletions(-) create mode 100644 matlab/find_chunk.m create mode 100644 matlab/find_region_ucfmf.m create mode 100644 matlab/read_field_block_ucfmf.m create mode 100644 matlab/read_field_complete_ucf.m diff --git a/matlab/@ucf/ucf.m b/matlab/@ucf/ucf.m index 14a8451..db4ec49 100644 --- a/matlab/@ucf/ucf.m +++ b/matlab/@ucf/ucf.m @@ -372,7 +372,7 @@ classdef ucf < handle [obj.currentSetDatatype,'=>',obj.currentSetDatatype]); end function [params] = readParameters(obj,tstep,dset) - % [data,params] = obj.readParameters(tstep,dset) + % [params] = obj.readParameters(tstep,dset) % Reads a raw set of parameters without parsing. % Input % tstep index of timestep to be read @@ -551,6 +551,24 @@ classdef ucf < handle % istep step index ndset = obj.numSetPerStep(istep); end + function [dpos,dsize,dtype,dnum] = getDataSetPosition(obj,tstep,dset) + % [params] = obj.readParameters(tstep,dset) + % Reads a raw set of parameters without parsing. + % Input + % tstep index of timestep to be read + % dset index of dataset to be read + % Output + % fpos offset to beginning of dataset + fseek(obj.fileID,obj.findSet(tstep,dset),'bof'); + if obj.Debug + fprintf('Skipped to position %d\n',ftell(obj.fileID)); + end + obj.readHeaderSet(); + dpos = ftell(obj.fileID); + dsize = obj.currentSetSize; + dtype = obj.currentSetDatatype; + dnum = obj.currentSetNumElements; + end end %% ------------------------------------------------------------------------%% %% PRIVATE METHODS %% diff --git a/matlab/@ustar/ustar.m b/matlab/@ustar/ustar.m index 13d0153..7ea732c 100644 --- a/matlab/@ustar/ustar.m +++ b/matlab/@ustar/ustar.m @@ -72,7 +72,7 @@ classdef ustar < handle % indexing data. % Input % tarfile path to TAR file - % indexfile path to index file (in json format) + % indexfile path to index file (in json/msgpack/taridx format) obj.File = tarfile; obj.IndexFile = indexfile; obj.IOMode = 'read'; @@ -248,7 +248,7 @@ classdef ustar < handle error('Unable to open file: %s',obj.IndexFile); end fseek(indexfileID,0,'bof'); - jsonstr = fread(indexfileID,'char=>char')'; + jsonstr = fread(indexfileID,'schar=>char')'; fclose(indexfileID); % Parse JSON and reconstruct filenames if ~isempty(which('jsonlab.loadjson')) @@ -316,7 +316,7 @@ classdef ustar < handle tarFileSize = zeros(1,nsubfile); tarFileOffset = zeros(1,nsubfile); for isub=1:nsubfile - tarFileName{isub} = deblank(fread(indexfileID,[1,256],'char=>char')); + tarFileName{isub} = deblank(fread(indexfileID,[1,256],'schar=>char')); tarFileOffset(isub) = fread(indexfileID,1,'int64=>double'); tarFileSize(isub) = fread(indexfileID,1,'int64=>double'); end @@ -339,7 +339,7 @@ classdef ustar < handle % in 'current*' class-variables. % Input % scanMode when set to true, omit parts which are not needed during scan - header = fread(obj.fileID,[1,obj.blockSize],'char=>char'); + header = fread(obj.fileID,[1,obj.blockSize],'schar=>char'); % Extract header information name = header(1:100); mode = header(101:108); diff --git a/matlab/find_chunk.m b/matlab/find_chunk.m new file mode 100644 index 0000000..ff1a633 --- /dev/null +++ b/matlab/find_chunk.m @@ -0,0 +1,18 @@ +function [col,row,pln] = find_chunk(xq,x,ibeg,iend,a,b,periodic) + + if periodic + xq = mod(xq-a,b-a)+a; + elseif any(xqb) + error('Query points must be located inside domain.'); + end + + x = reshape(x,[],1); + xq = reshape(xq,1,[]); + [~,ii] = min(abs(bsxfun(@minus,x,xq))); + + nxq = numel(xq); + col = zeros(1,nxq); + for iq=1:nxq + col(iq) = find(ii(iq)>=ibeg & ii(iq)<=iend,1,'first'); + end +end \ No newline at end of file diff --git a/matlab/find_region_ucfmf.m b/matlab/find_region_ucfmf.m new file mode 100644 index 0000000..f3de3fa --- /dev/null +++ b/matlab/find_region_ucfmf.m @@ -0,0 +1,219 @@ +function [ii,jj,kk,xq,yq,zq] = find_region_ucfmf(fdir,iseq,field,xb,yb,zb,varargin) + % [ii,jj,kk,xq,yq,zq] = find_region_ucfmf(fdir,iseq,field,xb,yb,zb,varargin) + % Determines the global index of grid points which lie within the bounds + % xb,yb,zb specified in the parameters, i.e. the grid points for which + % xb(1) <= x <= xb(2), + % yb(1) <= y <= yb(2), + % zb(1) <= z <= zb(2). + % If the domain is periodic, the specified region can lie "outside" of the domain, + % i.e. it is automatically translated to coordinates in the inner domain. + % If xb(1)==xb(2), a slice located at the nearest neighbor point will be created. + % Input + % fdir directory with files + % iseq sequence index + % field field to be indexed + % {'u','v','w','p','s1','s2',...} + % xb,yb,zb arrays which describe interval, e.g. [0,1] for all points between 0 and 1 + % ? fuzzy Include the first point outside of the region? (default: 0) + % ? verbosity verbose output? (default: 0) + % ? debug debug output? (default: 0) + % Output + % ii,jj,kk array with field indices + % xq,yq,zq corresponding grid, monotonically increasing, even if periodic boundary is crossed + + % Parse input + par = inputParser; + addParamValue(par,'verbosity',0,@isnumeric); + addParamValue(par,'debug',0,@isnumeric); + addParamValue(par,'fuzzy',0,@isnumeric); + parse(par,varargin{:}); + iverb = par.Results.verbosity; + idebug= par.Results.debug; + ifuzzy= par.Results.fuzzy; + + % Parse bounds + xmin = xb(1); + xmax = xb(2); + ymin = yb(1); + ymax = yb(2); + zmin = zb(1); + zmax = zb(2); + + % Parse field + switch field(1) + case 'u'; ifield=1; fbase='uvwp'; cmesh='u'; + case 'v'; ifield=2; fbase='uvwp'; cmesh='v'; + case 'w'; ifield=3; fbase='uvwp'; cmesh='w'; + case 'p'; ifield=4; fbase='uvwp'; cmesh='p'; + case 's'; ifield=str2double(field(2:end)); fbase='scal'; cmesh='p'; + otherwise; error('Invalid field: %s',field); + end + + % Read auxillary files + fparam = sprintf('%s/parameters_%04d.asc',fdir,iseq); + fgrid = sprintf('%s/grid_%04d.bin',fdir,iseq); + fproc = sprintf('%s/proc_%04d.bin',fdir,iseq); + params = read_parameters_ucf(fparam,'tarmode',0,'verbosity',iverb); + [x.u,y.u,z.u,x.v,y.v,z.v,x.w,y.w,z.w,x.p,y.p,z.p] = read_grid_ucf(fgrid,'tarmode',0,'verbosity',iverb,'debug',idebug); + x = x.(cmesh); + y = y.(cmesh); + z = z.(cmesh); + [ibeg.u,iend.u,jbeg.u,jend.u,kbeg.u,kend.u,... + ibeg.v,iend.v,jbeg.v,jend.v,kbeg.v,kend.v,... + ibeg.w,iend.w,jbeg.w,jend.w,kbeg.w,kend.w,... + ibeg.p,iend.p,jbeg.p,jend.p,kbeg.p,kend.p] = read_procgrid_ucf(fproc,'tarmode',0,'verbosity',iverb,'debug',idebug); + ibeg = ibeg.(cmesh); + jbeg = jbeg.(cmesh); + kbeg = kbeg.(cmesh); + iend = iend.(cmesh); + jend = jend.(cmesh); + kend = kend.(cmesh); + nx = numel(x); + ny = numel(y); + nz = numel(z); + + % Parse domain information + a = params.geometry.a; + b = params.geometry.b; + c = params.geometry.c; + d = params.geometry.d; + e = params.geometry.e; + f = params.geometry.f; + x_periodic = params.geometry.xperiodic; + y_periodic = params.geometry.yperiodic; + z_periodic = params.geometry.zperiodic; + + % Check if bounds are valid + if xmaxb || xmaxb + error('Interval must be located inside domain.'); + else + nxloop = 0; + nxoffset = 0; + end + if y_periodic + nyloop = floor((ymax-ymin)/(d-c)); + nyoffset = floor(ymin/(d-c)); + ymin = mod(ymin-c,d-c)+c; + ymax = mod(ymax-c,d-c)+c; + elseif ymind || ymaxd + error('Interval must be located inside domain.'); + else + nyloop = 0; + nyoffset = 0; + end + if z_periodic + nzloop = floor((zmax-zmin)/(f-e)); + nzoffset = floor(zmin/(f-e)); + zmin = mod(zmin-e,f-e)+e; + zmax = mod(zmax-e,f-e)+e; + elseif zminf || zmaxf + error('Interval must be located inside domain.'); + else + nzloop = 0; + nzoffset = 0; + end + + % Find gridpoints closest to the bounds + if ifuzzy + [~,imin] = find((x-xmin)<0,1,'last'); + [~,imax] = find((x-xmax)>0,1,'first'); + [~,jmin] = find((y-ymin)<0,1,'last'); + [~,jmax] = find((y-ymax)>0,1,'first'); + [~,kmin] = find((z-zmin)<0,1,'last'); + [~,kmax] = find((z-zmax)>0,1,'first'); + else + [~,imin] = find((x-xmin)>=0,1,'first'); + [~,imax] = find((x-xmax)<=0,1,'last'); + [~,jmin] = find((y-ymin)>=0,1,'first'); + [~,jmax] = find((y-ymax)<=0,1,'last'); + [~,kmin] = find((z-zmin)>=0,1,'first'); + [~,kmax] = find((z-zmax)<=0,1,'last'); + end + if isempty(imax); imax=1; end + if isempty(jmax); jmax=1; end + if isempty(kmax); kmax=1; end + + % If a slice is requested, correct previous results + if islice + [~,imin] = min(abs(x-xmin)); + imax = imin; + end + if jslice + [~,jmin] = min(abs(y-ymin)); + jmax = jmin; + end + if kslice + [~,kmin] = min(abs(z-zmin)); + kmax = kmin; + end + + % Construct the points which are inside of the interval and reconstruct the grid + if imin<=imax + ii = [repmat(mod(imin+(0:nx-1)-1,nx)+1,[1,nxloop]),imin:imax]; + shift = [reshape((repmat([zeros(1,nx-imin+1),ones(1,imin-1)],nxloop,1)+(1:nxloop)')',1,[]),... + ones(1,imax-imin+1)+nxloop]-1; + xq = x(ii)+(shift+nxoffset)*(b-a); + else + ii = [imin:nx,repmat(1:nx,[1,nxloop]),1:imax]; + shift = [ones(1,nx-imin+1),... + reshape((repmat(ones(1,nx),nxloop,1)+(1:nxloop)')',1,[]),... + ones(1,imax)+nxloop+1]-1; + xq = x(ii)+(shift+nxoffset)*(b-a); + end + if jmin<=jmax + jj = [repmat(mod(jmin+(0:ny-1)-1,ny)+1,[1,nyloop]),jmin:jmax]; + shift = [reshape((repmat([zeros(1,ny-jmin+1),ones(1,jmin-1)],nyloop,1)+(1:nyloop)')',1,[]),... + ones(1,jmax-jmin+1)+nyloop]-1; + yq = y(jj)+(shift+nyoffset)*(d-c); + else + jj = [jmin:ny,repmat(1:ny,[1,nyloop]),1:jmax]; + shift = [ones(1,ny-jmin+1),... + reshape((repmat(ones(1,ny),nyloop,1)+(1:nyloop)')',1,[]),... + ones(1,jmax)+nyloop+1]-1; + yq = y(jj)+(shift+nyoffset)*(d-c); + end + if kmin<=kmax + kk = [repmat(mod(kmin+(0:nz-1)-1,nz)+1,[1,nzloop]),kmin:kmax]; + shift = [reshape((repmat([zeros(1,nz-kmin+1),ones(1,kmin-1)],nzloop,1)+(1:nzloop)')',1,[]),... + ones(1,kmax-kmin+1)+nzloop]-1; + zq = z(kk)+(shift+nzoffset)*(f-e); + else + kk = [kmin:nz,repmat(1:nz,[1,nzloop]),1:kmax]; + shift = [ones(1,nz-kmin+1),... + reshape((repmat(ones(1,nz),nzloop,1)+(1:nzloop)')',1,[]),... + ones(1,kmax)+nzloop+1]-1; + zq = z(kk)+(shift+nzoffset)*(f-e); + end +end \ No newline at end of file diff --git a/matlab/read_field_block_ucfmf.m b/matlab/read_field_block_ucfmf.m new file mode 100644 index 0000000..6010ba6 --- /dev/null +++ b/matlab/read_field_block_ucfmf.m @@ -0,0 +1,108 @@ +function [q] = read_field_block_ucfmf(fdir,iseq,field,ii,jj,kk,varargin) + % [q] = read_field_block_ucfmf(fdir,iseq,field,ii,jj,kk,varargin) + % Reads a specified block of a field from processor chunks which hold the data (multi-file version) + % Input + % fdir directory with files + % iseq sequence index + % field field to be read + % {'u','v','w','p','s1','s2',...} + % ii,jj,kk list of global field indices to be read + % ? step index of step to be read (default: 1) + % ? verbosity verbose output? (default: 0) + % ? debug debug output? (default: 0) + % Output + % q partial field + + % Parse optional input arguments + par = inputParser; + addParamValue(par,'step',1,@isnumeric); + addParamValue(par,'verbosity',0,@isnumeric); + addParamValue(par,'debug',0,@isnumeric); + parse(par,varargin{:}); + istep = par.Results.step; + iverb = par.Results.verbosity; + idebug= par.Results.debug; + + % Parse field + switch field(1) + case 'u'; ifield=1; fbase='uvwp'; cmesh='u'; + case 'v'; ifield=2; fbase='uvwp'; cmesh='v'; + case 'w'; ifield=3; fbase='uvwp'; cmesh='w'; + case 'p'; ifield=4; fbase='uvwp'; cmesh='p'; + case 's'; ifield=str2double(field(2:end)); fbase='scal'; cmesh='s'; + otherwise; error('Invalid field: %s',field); + end + + % Read processor grid + fproc = sprintf('%s/proc_%04d.bin',fdir,iseq); + [ibeg.u,iend.u,jbeg.u,jend.u,kbeg.u,kend.u,... + ibeg.v,iend.v,jbeg.v,jend.v,kbeg.v,kend.v,... + ibeg.w,iend.w,jbeg.w,jend.w,kbeg.w,kend.w,... + ibeg.p,iend.p,jbeg.p,jend.p,kbeg.p,kend.p,... + ibeg.s,iend.s,jbeg.s,jend.s,kbeg.s,kend.s] = ... + read_procgrid_ucf(fproc,'tarmode',0,'verbosity',iverb,'debug',idebug); + ibeg = ibeg.(cmesh); + jbeg = jbeg.(cmesh); + kbeg = kbeg.(cmesh); + iend = iend.(cmesh); + jend = jend.(cmesh); + kend = kend.(cmesh); + nxprocs = length(ibeg); + nyprocs = length(jbeg); + nzprocs = length(kbeg); + + % Determine the corresponding chunks + nxb = numel(ii); + nyb = numel(jj); + nzb = numel(kk); + col = zeros(1,nxb); + row = zeros(1,nyb); + pln = zeros(1,nzb); + for iproc=1:nxprocs + col(ii>=ibeg(iproc) & ii<=iend(iproc)) = iproc; + end + for iproc=1:nyprocs + row(jj>=jbeg(iproc) & jj<=jend(iproc)) = iproc; + end + for iproc=1:nzprocs + pln(kk>=kbeg(iproc) & kk<=kend(iproc)) = iproc; + end + + % Setup fast indexing for blocks + iib = [1:nxb]; + jjb = [1:nyb]; + kkb = [1:nzb]; + + % Print some info (if verbose) + if iverb + nchunk = numel(unique(col))*numel(unique(row))*numel(unique(pln)); + fprintf('Reading [%d x %d x %d] points from %d chunks.\n',nxb,nyb,nzb,nchunk); + end + + % Load data + q = zeros(nxb,nyb,nzb); + for thiscol=unique(col) + for thisrow=unique(row) + for thispln=unique(pln) + % Load the data + iproc = (thiscol-1)*nyprocs*nzprocs+(thisrow-1)*nzprocs+(thispln-1); + fdata = sprintf('%s/%s_%04d.%05d',fdir,fbase,iseq,iproc); + [data,ib,jb,kb,nxl,nyl,nzl,ighost] = read_field_chunk_ucf(fdata,ifield,'tarmode',0,'ghost',0,'step',istep,'verbosity',iverb,'debug',idebug); + % Determine the points which are covered here + iiq = ii(col==thiscol); + jjq = jj(row==thisrow); + kkq = kk(pln==thispln); + % Translate to local indexing of loaded data + iiql = iiq-ib+1; + jjql = jjq-jb+1; + kkql = kkq-kb+1; + % Translate to local indexing of block + iiqb = iib(col==thiscol); + jjqb = jjb(row==thisrow); + kkqb = kkb(pln==thispln); + % Add loaded data to block + q(iiqb,jjqb,kkqb) = data(iiql,jjql,kkql); + end + end + end +end diff --git a/matlab/read_field_complete_ucf.m b/matlab/read_field_complete_ucf.m new file mode 100644 index 0000000..163326d --- /dev/null +++ b/matlab/read_field_complete_ucf.m @@ -0,0 +1,56 @@ +function [q,x,y,z] = read_field_complete_ucf(hucf,field,varargin) + % [q,x,y,z] = read_field_complete_ucf(hucf,field,varargin) + % Reads a specific field from all processor chunks.(UCF tar version) + % Input + % hucf handle of UCF object (ustar,ucfmulti) + % field field to be read + % {'u','v','w','p','s1','s2',...} + % ? step index of step to be read (default: 1) + % ? verbosity verbose output? (default: 0) + % ? debug debug output? (default: 0) + % Output + % q complete field + % x,y,z corresponding grid + + % Parse optional input arguments + par = inputParser; + addParamValue(par,'step',1,@isnumeric); + addParamValue(par,'verbosity',0,@isnumeric); + addParamValue(par,'debug',0,@isnumeric); + parse(par,varargin{:}); + istep = par.Results.step; + + % Parse field + switch field(1) + case 'u'; cmesh = 'u'; + case 'v'; cmesh = 'v'; + case 'w'; cmesh = 'w'; + case 'p'; cmesh = 'p'; + case 's'; cmesh = 'p'; + otherwise; error('Invalid field: %s',field); + end + + % Read parameters + params = read_parameters_ucf(hucf); + + % Construct an array with final mesh size + nx = params.mesh.(sprintf('nx%s',cmesh)); + ny = params.mesh.(sprintf('ny%s',cmesh)); + nz = params.mesh.(sprintf('nz%s',cmesh)); + q = zeros(nx,ny,nz); + + % Read grid + [x.u,y.u,z.u,x.v,y.v,z.v,x.w,y.w,z.w,x.p,y.p,z.p] = read_grid_ucf(hucf); + x = x.(cmesh); + y = y.(cmesh); + z = z.(cmesh); + + % Get number of processors + nprocs = params.parallel.nprocs; + + % Now read chunk by chunk + for iproc=0:nprocs-1 + [data,ib,jb,kb,nxl,nyl,nzl] = read_field_chunk_ucf(hucf,field,'rank',iproc,'ghost',0,'step',istep); + q(ib:ib+nxl-1,jb:jb+nyl-1,kb:kb+nzl-1) = data; + end +end diff --git a/matlab/read_grid_ucf.m b/matlab/read_grid_ucf.m index 197a197..0e36547 100644 --- a/matlab/read_grid_ucf.m +++ b/matlab/read_grid_ucf.m @@ -33,7 +33,7 @@ function [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,xs,ys,zs] = read_grid_ucf(file,var end % Define sets to be read - isDualMesh=(obj.UCFVersion>=2); + isDualMesh=(obj.UCFVersion>=2) && (obj.NumDataset>4); if isDualMesh sets = {'u','v','w','p','s'}; else diff --git a/matlab/read_procgrid_ucf.m b/matlab/read_procgrid_ucf.m index 10ef043..8547d7e 100644 --- a/matlab/read_procgrid_ucf.m +++ b/matlab/read_procgrid_ucf.m @@ -40,7 +40,7 @@ function [ibegu,iendu,jbegu,jendu,kbegu,kendu,... end % Define sets to be read - isDualMesh=(obj.UCFVersion>=2); + isDualMesh=(obj.UCFVersion>=2) && (obj.NumDataset>4); if isDualMesh sets = {'u','v','w','p','s'}; else diff --git a/matlab/read_statistics_channel_scal_ucf.m b/matlab/read_statistics_channel_scal_ucf.m index 4e7860a..354f372 100644 --- a/matlab/read_statistics_channel_scal_ucf.m +++ b/matlab/read_statistics_channel_scal_ucf.m @@ -6,17 +6,21 @@ function [tbeg,tend,nstat,sm,ss,su,sv] = read_statistics_channel_scal_ucf(file,v addParamValue(par,'debug',0,@isnumeric); parse(par,varargin{:}); - % Define sets to be read - sets = {'um','uu','vv','ww','uv','uub','uvb'}; - nset = numel(sets); - % Open UCF file and read obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug); - obj.open(file); - tend = obj.getSimulationTime(1); + switch class(file) + case 'char' + obj.open(file); + case {'ustar','ucfmulti'} + ptr = file.pointer('statistics_scal.bin'); + obj.opentar(ptr); + otherwise + error('Input file type not supported: %s',class(file)); + end + tend = obj.getSimulationTime(); % Read parameters of first set - params = obj.readParameters(1,1); + params = obj.readParameters(1,1) % Parse parameters tbeg = typecast(params(1),'double'); From 0284e8f95422cac32d1589f256420a503b3fc9c6 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Fri, 22 Jan 2021 14:10:00 +0100 Subject: [PATCH 03/12] hawk update --- python/ucf/tools/ucf_file_diff | 64 ++++++++++++++++++ python/ucf/tools/ucf_file_diff_single_multi | 73 +++++++++++++++++++++ python/ucf/tools/ucf_file_integrity | 50 ++++++++++++++ python/ucf/ucf.py | 6 +- 4 files changed, 190 insertions(+), 3 deletions(-) create mode 100755 python/ucf/tools/ucf_file_diff create mode 100755 python/ucf/tools/ucf_file_diff_single_multi create mode 100755 python/ucf/tools/ucf_file_integrity diff --git a/python/ucf/tools/ucf_file_diff b/python/ucf/tools/ucf_file_diff new file mode 100755 index 0000000..341daa6 --- /dev/null +++ b/python/ucf/tools/ucf_file_diff @@ -0,0 +1,64 @@ +#!/usr/bin/env python3 +import ucf +import argparse +import os +import sys +import traceback +import numpy as np + +parser = argparse.ArgumentParser(description='Tests basic integrity of a UCF file') +parser.add_argument('file1', metavar='file1',help='input file #1') +parser.add_argument('file2', metavar='file2',help='input file #2') +parser.add_argument("--debug", help="enable debug output? [default: False]", action="store_true") +parser.add_argument("--stop", help="stop as soon as an error has been detected? [default: False]", action="store_true") +parser.add_argument("--tolerance", metavar='tol', help="floating point tolerance for comparison [default: 1e-12]", type=float) +parser.add_argument("-v","--verbose", help="enable verbose output? [default: False]", action="store_true") +args = parser.parse_args() + +file_in1 = args.file1 +file_in2 = args.file2 +flag_debug = args.debug +flag_verbose = args.verbose +fptol = args.tolerance +flag_stop = args.stop + +if fptol is None: + fptol = 1e-12 + +if os.path.isfile(file_in1)==False: + print('File not found:',file_in1) + sys.exit(1) +if os.path.isfile(file_in2)==False: + print('File not found:',file_in2) + sys.exit(1) + +obj1 = ucf.UCF(file=file_in1,verbosity=flag_verbose,debug=flag_debug) +obj2 = ucf.UCF(file=file_in2,verbosity=flag_verbose,debug=flag_debug) +if obj1.NumTimestep!=obj2.NumTimestep: + print("Number of time steps differ: {:d}, {:d}".format(obj1.NumTimestep,obj2.NumTimestep)) + sys.exit(2) + +errfound=False +for istep in range(0,obj1.NumTimestep): + for iset in range(0,obj1._UCF__numSetPerStep[istep]): + (data1,param1) = obj1.readSet(step=istep+1,dset=iset+1) + (data2,param2) = obj2.readSet(step=istep+1,dset=iset+1) + if not param1==param2: + print('Step #{:d}, set #{:d}: parameter mismatch. {}, {}'.format(istep,iset,param1,param2)) + if flag_stop: + sys.exit(3) + else: + errfound=True + if not np.allclose(data1,data2,rtol=fptol): + data_diff = np.abs(data1-data2) + idx = np.argmax(data_diff) + vald = data_diff.flatten()[idx] + val1 = data1.flatten()[idx] + val2 = data2.flatten()[idx] + print('Step #{:d}, set #{:d}: difference in FP data above tolerance. idx={}, diff={}, val1={}, val2={}'.format(istep,iset,idx,vald,val1,val2)) + if flag_stop: + sys.exit(4) + else: + errfound=True +if errfound: + sys.exit(5) diff --git a/python/ucf/tools/ucf_file_diff_single_multi b/python/ucf/tools/ucf_file_diff_single_multi new file mode 100755 index 0000000..4c39d3c --- /dev/null +++ b/python/ucf/tools/ucf_file_diff_single_multi @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +import ucf +import argparse +import os +import sys +import traceback +import numpy as np + +parser = argparse.ArgumentParser(description='Verifies that content of single UCF file is equivalent to multiple UCF files') +parser.add_argument('file1', metavar='file_single',help='input file #1 (single file)') +parser.add_argument('file2', metavar='filebase_multi',help='input file #2 (multiple files)') +parser.add_argument("--debug", help="enable debug output? [default: False]", action="store_true") +parser.add_argument("--stop", help="stop as soon as an error has been detected? [default: False]", action="store_true") +parser.add_argument("--tolerance", metavar='tol', help="floating point tolerance for comparison [default: 1e-12]", type=float) +parser.add_argument("-v","--verbose", help="enable verbose output? [default: False]", action="store_true") +args = parser.parse_args() + +file_single = args.file1 +fiba_multi = args.file2 +flag_debug = args.debug +flag_verbose = args.verbose +fptol = args.tolerance +flag_stop = args.stop + +if fptol is None: + fptol = 1e-12 + +if os.path.isfile(file_single)==False: + print('File not found:',os.path.isfile(file_single)) + sys.exit(1) + +objs = ucf.UCF(file=file_single,verbosity=flag_verbose,debug=flag_debug) +#if obj1.NumTimestep!=obj2.NumTimestep: + +errfound=False +for istep in range(0,objs.NumTimestep): + print('Comparing rank {:5d}'.format(istep)) + file_multi = '{}{:05d}'.format(fiba_multi,istep) + if os.path.isfile(file_multi)==False: + print('File not found:',os.path.isfile(file_multi)) + if flag_stop: + sys.exit(1) + else: + errfound=True + objm = ucf.UCF(file=file_multi,verbosity=flag_verbose,debug=flag_debug) + if objs._UCF__numSetPerStep[istep]!=objm._UCF__numSetPerStep[0]: + print("Number of datasets differs: {:d}, {:d}".format(objs._UCF__numSetPerStep[istep],objm._UCF__numSetPerStep[0])) + if flag_stop: + sys.exit(2) + else: + errfound=True + for iset in range(0,objs._UCF__numSetPerStep[istep]): + (data1,param1) = objs.readSet(step=istep+1,dset=iset+1) + (data2,param2) = objm.readSet(step=1,dset=iset+1) + if not param1==param2: + print('Step #{:d}, set #{:d}: parameter mismatch. {}, {}'.format(istep,iset,param1,param2)) + if flag_stop: + sys.exit(3) + else: + errfound=True + if not np.allclose(data1,data2,rtol=fptol): + data_diff = np.abs(data1-data2) + idx = np.argmax(data_diff) + vald = data_diff.flatten()[idx] + val1 = data1.flatten()[idx] + val2 = data2.flatten()[idx] + print('Step #{:d}, set #{:d}: difference in FP data above tolerance. idx={}, diff={}, val1={}, val2={}'.format(istep,iset,idx,vald,val1,val2)) + if flag_stop: + sys.exit(4) + else: + errfound=True +if errfound: + sys.exit(5) diff --git a/python/ucf/tools/ucf_file_integrity b/python/ucf/tools/ucf_file_integrity new file mode 100755 index 0000000..d97c929 --- /dev/null +++ b/python/ucf/tools/ucf_file_integrity @@ -0,0 +1,50 @@ +#!/usr/bin/env python3 +import ucf +import argparse +import os +import sys +import traceback + +parser = argparse.ArgumentParser(description='Tests basic integrity of a UCF file') +parser.add_argument('infile', metavar='file', nargs='+',help='input file') +parser.add_argument("--debug", help="enable debug output? [default: False]", action="store_true") +parser.add_argument("--print-error", help="print possible error messages? [default: False]", action="store_true") +parser.add_argument("--stop", help="stop as soon as a corrupted file has been found? [default: False]", action="store_true") +parser.add_argument("-v","--verbose", help="enable verbose output? [default: False]", action="store_true") +args = parser.parse_args() + +file_in = args.infile +flag_debug = args.debug +flag_verbose = args.verbose +flag_error = args.print_error +flag_stop = args.stop + +numerror=0 +for file in file_in: + print(file+": ",end="") + if os.path.isfile(file)==False: + print('not found.') + numerror+=1 + continue + try: + obj = ucf.UCF(file=file,verbosity=flag_verbose,debug=flag_debug) + for it in range(0,obj.NumTimestep): + for iset in range(0,obj._UCF__numSetPerStep[it]): + obj._UCF__findSet(it+1,iset+1) + obj._UCF__readHeaderSet() + obj.close() + print('successful ({:d} time steps, {:d} data sets)'.format(obj.NumTimestep,obj.NumDataset)) + except Exception: + print('corrupt') + if flag_error: + print(traceback.print_exc()) + if flag_stop: + sys.exit(1) + numerror+=1 + +if numerror==0: + print('Test completed with no errors.') + sys.exit(0) +else: + print("Test completed with {:d} error(s).".format(numerror)) + sys.exit(1) diff --git a/python/ucf/ucf.py b/python/ucf/ucf.py index 70eeeac..d94c13c 100644 --- a/python/ucf/ucf.py +++ b/python/ucf/ucf.py @@ -63,7 +63,7 @@ class UCF: # Scan through file to get the basic structure (steps/sets) self.__timeStep = np.zeros(self.__scanBuffSize,dtype=np.float64) - self.__posStep = np.zeros(self.__scanBuffSize,dtype=np.int32) + self.__posStep = np.zeros(self.__scanBuffSize,dtype=np.int64) self.__numSetPerStep = np.zeros(self.__scanBuffSize,dtype=np.int32) istep = 0; while self.__fileID.tell() Date: Wed, 3 Feb 2021 12:56:13 +0100 Subject: [PATCH 04/12] split mpiio ucf files while taring --- python/ucf/tools/ucf_mpiio_tar | 128 +++++++++++++++++++++++++++++++++ 1 file changed, 128 insertions(+) create mode 100755 python/ucf/tools/ucf_mpiio_tar diff --git a/python/ucf/tools/ucf_mpiio_tar b/python/ucf/tools/ucf_mpiio_tar new file mode 100755 index 0000000..55b9ba8 --- /dev/null +++ b/python/ucf/tools/ucf_mpiio_tar @@ -0,0 +1,128 @@ +#!/usr/bin/env python3 +import sys +import io +import os +import tarfile +import argparse +import numpy as np +import ucf +import configparser + +parser = argparse.ArgumentParser(description='Reads an ucf.tar archive, downsamples it and saves it to a new ucf.tar archive. Can be used as a pipe.') +parser.add_argument('indir', metavar='dirin', help='input directory') +parser.add_argument('iseq', metavar='iseq', help='sequence number') +parser.add_argument('base', metavar='base', help='filebase to be archived: "uvwp" or "scal"') +parser.add_argument("-o", "--outfile", metavar='file',nargs='?', default=None, help="name of the output file [default: snapshot_XXXX.ucf.tar]", action="store") +parser.add_argument("-v", "--verbose", help="activate verbose output", action="store_true") +args = parser.parse_args() + +dir_in = args.indir +iseq = int(args.iseq) +base = args.base +verbose = args.verbose + +if not base in ('uvwp','scal'): + raise ValueError('Invalid base key: {}',base) + +if args.outfile is None: + if base=='uvwp': + file_out = 'snapshot_{:04d}.ucf.tar'.format(iseq) + elif base=='scal': + file_out = 'snapshot_scalar_{:04d}.ucf.tar'.format(iseq) +else: + file_out = args.outfile + +# Check if all required files are available +files_req = ('parameters_{:04d}.asc','proc_{:04d}.bin','grid_{:04d}.bin','particles_{:04d}.bin',base+'_{:04d}.bin') +for files in files_req: + path_check = '{}/{}'.format(dir_in,files.format(iseq)) + if not os.path.isfile(path_check): + raise IOError('File does not exist: {}'.format(path_check)) +if verbose: + print("[x] parameters") + print("[x] grid") + print("[x] processor grid") + print("[x] particles") + if base=='uvwp': + print("[x] fluid field") + elif base=='scal': + print("[x] scalar field") + +# Open tar file for output +ftar = tarfile.open(name=file_out,mode='w',pax_headers=tarfile.USTAR_FORMAT) + +def transform_filename(filename,iseq): + return os.path.basename(file_in).replace('_{:04d}'.format(iseq),'') + +# Parse parameters to construct file headers, then add it to tar +file_in = '{}/parameters_{:04d}.asc'.format(dir_in,iseq) +config = configparser.ConfigParser() +config.read(file_in) +nxprocs = int(config['parallel']['nxprocs']) +nyprocs = int(config['parallel']['nyprocs']) +nzprocs = int(config['parallel']['nzprocs']) +nprocs = nxprocs*nyprocs*nzprocs + +print(transform_filename(file_in,iseq)) +fid = open(file_in,'rb') +info = tarfile.TarInfo(name=transform_filename(file_in,iseq)) +info.size = os.path.getsize(file_in) +ftar.addfile(info,fileobj=fid) +fid.close() + +# Add proc.bin +file_in = '{}/proc_{:04d}.bin'.format(dir_in,iseq) +print(transform_filename(file_in,iseq)) +fid = open(file_in,'rb') +info = tarfile.TarInfo(name=transform_filename(file_in,iseq)) +info.size = os.path.getsize(file_in) +ftar.addfile(info,fileobj=fid) +fid.close() + +# Add grid.bin +file_in = '{}/grid_{:04d}.bin'.format(dir_in,iseq) +print(transform_filename(file_in,iseq)) +fid = open(file_in,'rb') +info = tarfile.TarInfo(name=transform_filename(file_in,iseq)) +info.size = os.path.getsize(file_in) +ftar.addfile(info,fileobj=fid) +fid.close() + +# Add particles.bin +file_in = '{}/particles_{:04d}.bin'.format(dir_in,iseq) +print(transform_filename(file_in,iseq)) +fid = open(file_in,'rb') +info = tarfile.TarInfo(name=transform_filename(file_in,iseq)) +info.size = os.path.getsize(file_in) +ftar.addfile(info,fileobj=fid) +fid.close() + +# Split uvwp/scal files and add them +def positionFromRank(rank,nxp,nyp,nzp): + ip = rank//(nyp*nzp) + jp = (rank//nzp)%nyp + kp = rank%nzp + return (ip,jp,kp) + +file_in = '{}/{}_{:04d}.bin'.format(dir_in,base,iseq) +ucf_data = ucf.UCF(file=file_in) +if nprocs!=ucf_data.NumTimestep: + raise ValueError('Number of steps does not equal number of procs: {}, {}',ucf_data.NumTimestep,nprocs) + +for iproc in range(0,nprocs): + print(base+'.{:05d}'.format(iproc)) + # Construct a new UCF file + ucf_data.addFileHeaderToBuffer( + rank=iproc, + rankijk=positionFromRank(iproc,nxprocs,nyprocs,nzprocs), + ftype=ucf_data._UCF__typeID + ) + ucf_data.copyStepToBuffer(iproc+1,recursive=True) + ucf_bytes = ucf_data.flushBuffer() + # Write it to tar + info = tarfile.TarInfo(name=base+'.{:05d}'.format(iproc)) + info.size = len(ucf_bytes) + ftar.addfile(info,fileobj=io.BytesIO(ucf_bytes)) + +# Close tar file +ftar.close() From dda1d88a9a00abf2d634419eef7604fe54d04b07 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Wed, 3 Feb 2021 13:05:04 +0100 Subject: [PATCH 05/12] renamed file --- python/ucf/tools/{ucf_mpiio_tar => ucftar_mpiio} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename python/ucf/tools/{ucf_mpiio_tar => ucftar_mpiio} (100%) diff --git a/python/ucf/tools/ucf_mpiio_tar b/python/ucf/tools/ucftar_mpiio similarity index 100% rename from python/ucf/tools/ucf_mpiio_tar rename to python/ucf/tools/ucftar_mpiio From 2d7358c0626a8743deb7b057641e1d8d3497f966 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Wed, 3 Feb 2021 13:07:38 +0100 Subject: [PATCH 06/12] renamed file --- python/ucf/tools/{ucftar_mpiio => ucftar_mpiio_pack} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename python/ucf/tools/{ucftar_mpiio => ucftar_mpiio_pack} (100%) diff --git a/python/ucf/tools/ucftar_mpiio b/python/ucf/tools/ucftar_mpiio_pack similarity index 100% rename from python/ucf/tools/ucftar_mpiio rename to python/ucf/tools/ucftar_mpiio_pack From 7687d4f022f878099efb98328ae4f67de20ab9b5 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Wed, 3 Feb 2021 16:22:15 +0100 Subject: [PATCH 07/12] unpack files to MPIIO --- python/ucf/tools/ucftar_mpiio_unpack | 66 ++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) create mode 100755 python/ucf/tools/ucftar_mpiio_unpack diff --git a/python/ucf/tools/ucftar_mpiio_unpack b/python/ucf/tools/ucftar_mpiio_unpack new file mode 100755 index 0000000..eb9fda0 --- /dev/null +++ b/python/ucf/tools/ucftar_mpiio_unpack @@ -0,0 +1,66 @@ +#!/usr/bin/env python3 +import sys +import io +import os +import tarfile +import argparse +import numpy as np +import ucf +import configparser + +parser = argparse.ArgumentParser(description='Unpacks an UCF tar file as MPIIO input file.') +parser.add_argument('infile', metavar='filein', help='input file') +parser.add_argument('iseq', metavar='iseq', help='sequence number') +#parser.add_argument('base', metavar='base', help='filebase to be archived: "uvwp" or "scal"') +parser.add_argument("-o", "--outdir", metavar='directory',nargs='?', default=None, help="name of the output directory", action="store") +parser.add_argument("-v", "--verbose", help="activate verbose output", action="store_true") +args = parser.parse_args() + +file_in = args.infile +iseq = int(args.iseq) +#base = args.base +verbose = args.verbose + +if args.outdir is None: + dir_out = './' +else: + dir_out = args.outdir + +# Open tar file for input +ftar = tarfile.open(name=file_in,mode='r') + +# Extract files +def transform_filename(filename,iseq): + return filename.replace('.','_{:04d}.'.format(iseq)) + +files_all = ftar.getnames() +#files_aux = ('parameters.asc','proc.bin','grid.bin','particles.bin') +files_aux = [s for s in files_all if not s.startswith('uvwp') or s.startswith('scal')] +for files in files_aux: + print(files) + ucf_bytes = ftar.extractfile(files).read() + file_out = dir_out+transform_filename(files,iseq) + fidw = open(file_out,'wb') + fidw.write(ucf_bytes) + fidw.close() + +files_data = [s for s in files_all if s.startswith('uvwp') or s.startswith('scal')] +base = files_data[0].split('.')[0] +iproc = 0 +file_out = dir_out+'/'+base+'_{:04d}.bin'.format(iseq) +fidw = open(file_out,'wb') +for files in files_data: + print(files) + sproc = files.split('.')[1] + if iproc!=int(sproc): + raise ValueError('Invalid file order in tar: {}, {}'.format(iproc,files)) + ucf_bytes = ftar.extractfile(files).read() + if iproc==0: + fidw.write(ucf_bytes) + else: + fidw.write(ucf_bytes[64:]) + iproc += 1 +fidw.close() + +# Close tar file +ftar.close() From f92bd36a76ce27341e3656d64302b78ff68793b7 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Wed, 3 Feb 2021 16:22:33 +0100 Subject: [PATCH 08/12] bugfix: filename --- python/ucf/tools/ucftar_mpiio_pack | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/ucf/tools/ucftar_mpiio_pack b/python/ucf/tools/ucftar_mpiio_pack index 55b9ba8..442719a 100755 --- a/python/ucf/tools/ucftar_mpiio_pack +++ b/python/ucf/tools/ucftar_mpiio_pack @@ -8,7 +8,7 @@ import numpy as np import ucf import configparser -parser = argparse.ArgumentParser(description='Reads an ucf.tar archive, downsamples it and saves it to a new ucf.tar archive. Can be used as a pipe.') +parser = argparse.ArgumentParser(description='Packs UCF MPIIO data into a tar archive. The single output file is split into multiple files due to file size restrictions.') parser.add_argument('indir', metavar='dirin', help='input directory') parser.add_argument('iseq', metavar='iseq', help='sequence number') parser.add_argument('base', metavar='base', help='filebase to be archived: "uvwp" or "scal"') @@ -52,7 +52,7 @@ if verbose: ftar = tarfile.open(name=file_out,mode='w',pax_headers=tarfile.USTAR_FORMAT) def transform_filename(filename,iseq): - return os.path.basename(file_in).replace('_{:04d}'.format(iseq),'') + return os.path.basename(filename).replace('_{:04d}'.format(iseq),'') # Parse parameters to construct file headers, then add it to tar file_in = '{}/parameters_{:04d}.asc'.format(dir_in,iseq) From 43a3ab56f6b17c9c342d7215c1d66facdf0a04a9 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Fri, 12 Feb 2021 02:02:19 +0100 Subject: [PATCH 09/12] implemented single file saving with downsampling, field manipulation using statistics --- python/ibmppp/ibmppp.py | 149 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 143 insertions(+), 6 deletions(-) diff --git a/python/ibmppp/ibmppp.py b/python/ibmppp/ibmppp.py index b634fc6..6fd5e65 100644 --- a/python/ibmppp/ibmppp.py +++ b/python/ibmppp/ibmppp.py @@ -621,6 +621,108 @@ class ibmppp: fid.write(' \n') fid.close() + def saveFieldSingleFile(self,filename,key,append=False,downsample=False,xdmf=False,verbose=False): + '''Writes field data into a single h5 file''' + # Determine filename of output file + file_out = filename+'.h5' + # Compute global nx,ny,nz + nxg = self.__procGrid[key][1][-1]-self.__procGrid[key][0][0]+1 + nyg = self.__procGrid[key][3][-1]-self.__procGrid[key][2][0]+1 + nzg = self.__procGrid[key][5][-1]-self.__procGrid[key][4][0]+1 + # Adjust these values if downsampling is requested + if downsample: + nxg = np.ceil(nxg/2).astype('int') + nyg = np.ceil(nyg/2).astype('int') + nzg = np.ceil(nzg/2).astype('int') + # Setup the correct slicing + if not downsample: + ishift = 0 + jshift = 0 + kshift = 0 + ifb = self.__localChunkBounds[key][0] + jfb = self.__localChunkBounds[key][2] + kfb = self.__localChunkBounds[key][4] + nxl = self.__localChunkSize[key][0] + nyl = self.__localChunkSize[key][1] + nzl = self.__localChunkSize[key][2] + step = 1 + else: + ishift = 1 if self.__localChunkBounds[key][0]%2==0 else 0 + jshift = 1 if self.__localChunkBounds[key][2]%2==0 else 0 + kshift = 1 if self.__localChunkBounds[key][4]%2==0 else 0 + ifb = np.ceil((self.__localChunkBounds[key][0]+ishift)/2).astype('int') + jfb = np.ceil((self.__localChunkBounds[key][2]+jshift)/2).astype('int') + kfb = np.ceil((self.__localChunkBounds[key][4]+kshift)/2).astype('int') + nxl = np.ceil(self.__localChunkSize[key][0]/2).astype('int') + nyl = np.ceil(self.__localChunkSize[key][1]/2).astype('int') + nzl = np.ceil(self.__localChunkSize[key][2]/2).astype('int') + step = 2 + isl_ch = slice(self.__nghx+ishift,self.__nghx+self.__localChunkSize[key][0],step) + jsl_ch = slice(self.__nghy+jshift,self.__nghy+self.__localChunkSize[key][1],step) + ksl_ch = slice(self.__nghz+kshift,self.__nghz+self.__localChunkSize[key][2],step) + isl_h5 = slice(ifb-1,ifb+nxl-1) + jsl_h5 = slice(jfb-1,jfb+nyl-1) + ksl_h5 = slice(kfb-1,kfb+nzl-1) + # If append flag is set, we add a dataset to the existing file + if append: + ioflag = 'a' + else: + ioflag = 'w' + # I avoid parallel IO because h5py is not available everywhere with + # MPI support. So rank 0 does all the writing and we communicate manually. + if self.__rank==0: + if verbose: + print('[saveFieldSingleFile] rank {} / {} writing to file'.format(self.__rank,self.__nproc-1),end='\r') + fid = h5py.File(file_out,ioflag) + gid = fid.create_group('/'+key) + gid.create_dataset('nx',data=nxg) + gid.create_dataset('ny',data=nyg) + gid.create_dataset('nz',data=nzg) + did = gid.create_dataset('data',(nzg,nyg,nxg),dtype=self.__precision) + did[ksl_h5,jsl_h5,isl_h5] = np.transpose(self.field[key][isl_ch,jsl_ch,ksl_ch]) + fid.close() + self.__comm.send(True, dest=self.__rank+1, tag=1) + else: + msg = self.__comm.recv(source=self.__rank-1,tag=1) + if verbose: + print('[saveFieldSingleFile] rank {} / {} writing to file'.format(self.__rank,self.__nproc-1),end='\r') + fid = h5py.File(file_out,'a') + fid['/'+key+'/data'][ksl_h5,jsl_h5,isl_h5] = np.transpose(self.field[key][isl_ch,jsl_ch,ksl_ch]) + fid.close() + if self.__rank!=self.__nproc-1: + self.__comm.send(True, dest=self.__rank+1, tag=1) + # Create an XDMF file for the field + if xdmf and self.__rank==0: + if verbose: + print('[saveFieldSingleFile] generating XDMF',end='\r') + # Construct XDMF filename + file_xdmf = filename+'_'+key+'.xdmf' + filename_h5 = os.path.basename(file_out) + # Open file and write header + fid = open(file_xdmf,'w') + fid.write('\n') + fid.write('\n') + fid.write(' \n') + fid.write(' \n'.format(key)) + #fid.write(' \n') + fid.write(' \n') + fid.write(' \n') + fid.close() + def saveStatistics(self,filename): '''Writes all gathered statistics to a h5 file.''' if self.__rank==0: @@ -1201,17 +1303,52 @@ class ibmppp: '''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: + if not key in self.spatialMean or self.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) + # Slice statistics for each processor + ib,ie = (0,1) if dims[0]==1 else (self.__localChunkBounds[key][0]-1,self.__localChunkBounds[key][1]) + jb,je = (0,1) if dims[1]==1 else (self.__localChunkBounds[key][2]-1,self.__localChunkBounds[key][3]) + kb,ke = (0,1) if dims[2]==1 else (self.__localChunkBounds[key][4]-1,self.__localChunkBounds[key][5]) + # Subtract mean from field + self.field[key][ + self.__nghx:-self.__nghx, + self.__nghy:-self.__nghy, + self.__nghz:-self.__nghz] -= self.spatialMean[key][ib:ie,jb:je,kb:ke] + # Exchange ghost cells and set boundary conditions + self.exchangeGhostCells(key) + self.imposeBoundaryConditions(key) + + def normalizeByStandardDeviation(self,key): + '''Normalize field by its standard deviation.''' + # Check if this operation can be executed + if self.__rank==0: + if (not key in self.spatialMean or self.spatialMean[key] is None or + not key in self.spatialMsqr or self.spatialMsqr[key] is None): + raise ValueError('Statistics not available') + # Broadcast statistics + self.spatialMean[key] = self.__comm.bcast(self.spatialMean[key],root=0) + self.spatialMsqr[key] = self.__comm.bcast(self.spatialMsqr[key],root=0) + # Compute standard deviation + assert(self.spatialMean[key].shape==self.spatialMsqr[key].shape) + spatialStdDev = np.sqrt(self.spatialMsqr[key]-np.square(self.spatialMean[key])) + # Get dimensions + dims = spatialStdDev.shape + # Slice statistics for each processor + ib,ie = (0,1) if dims[0]==1 else (self.__localChunkBounds[key][0]-1,self.__localChunkBounds[key][1]) + jb,je = (0,1) if dims[1]==1 else (self.__localChunkBounds[key][2]-1,self.__localChunkBounds[key][3]) + kb,ke = (0,1) if dims[2]==1 else (self.__localChunkBounds[key][4]-1,self.__localChunkBounds[key][5]) + # Subtract mean from field + self.field[key][ + self.__nghx:-self.__nghx, + self.__nghy:-self.__nghy, + self.__nghz:-self.__nghz] /= spatialStdDev[ib:ie,jb:je,kb:ke] + # Exchange ghost cells and set boundary conditions + self.exchangeGhostCells(key) + self.imposeBoundaryConditions(key) def collocateVelocity(self): '''Moves velocity components onto the pressure grid.''' From 2e2f987f7f3497f38bee4217741fc905370d8fc6 Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Fri, 12 Feb 2021 02:11:15 +0100 Subject: [PATCH 10/12] bugfix: statistics communication for derived fields --- python/ibmppp/ibmppp.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/python/ibmppp/ibmppp.py b/python/ibmppp/ibmppp.py index 6fd5e65..62dcde7 100644 --- a/python/ibmppp/ibmppp.py +++ b/python/ibmppp/ibmppp.py @@ -1306,7 +1306,8 @@ class ibmppp: if not key in self.spatialMean or self.spatialMean[key] is None: raise ValueError('Statistics not available') # Broadcast statistics - self.spatialMean[key] = self.__comm.bcast(self.spatialMean[key],root=0) + sendbuf = self.spatialMean[key] if self.__rank==0 else None + self.spatialMean[key] = self.__comm.bcast(sendbuf,root=0) # Get dimensions dims = self.spatialMean[key].shape # Slice statistics for each processor @@ -1330,8 +1331,10 @@ class ibmppp: not key in self.spatialMsqr or self.spatialMsqr[key] is None): raise ValueError('Statistics not available') # Broadcast statistics - self.spatialMean[key] = self.__comm.bcast(self.spatialMean[key],root=0) - self.spatialMsqr[key] = self.__comm.bcast(self.spatialMsqr[key],root=0) + sendbuf = self.spatialMean[key] if self.__rank==0 else None + self.spatialMean[key] = self.__comm.bcast(sendbuf,root=0) + sendbuf = self.spatialMsqr[key] if self.__rank==0 else None + self.spatialMsqr[key] = self.__comm.bcast(sendbuf,root=0) # Compute standard deviation assert(self.spatialMean[key].shape==self.spatialMsqr[key].shape) spatialStdDev = np.sqrt(self.spatialMsqr[key]-np.square(self.spatialMean[key])) From f707360f17e4da4836e0fc5364c751f35841bbbd Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Wed, 12 May 2021 13:07:14 +0200 Subject: [PATCH 11/12] updated output --- python/ibmppp/ibmppp.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/python/ibmppp/ibmppp.py b/python/ibmppp/ibmppp.py index 62dcde7..5f5a47c 100644 --- a/python/ibmppp/ibmppp.py +++ b/python/ibmppp/ibmppp.py @@ -678,6 +678,14 @@ class ibmppp: gid.create_dataset('nx',data=nxg) gid.create_dataset('ny',data=nyg) gid.create_dataset('nz',data=nzg) + x0 = self.grid[key][0][0] + y0 = self.grid[key][1][0] + z0 = self.grid[key][2][0] + dx = self.__dx[key]*step + gid.create_dataset('x0',data=x0) + gid.create_dataset('y0',data=y0) + gid.create_dataset('z0',data=z0) + gid.create_dataset('dx',data=dx) did = gid.create_dataset('data',(nzg,nyg,nxg),dtype=self.__precision) did[ksl_h5,jsl_h5,isl_h5] = np.transpose(self.field[key][isl_ch,jsl_ch,ksl_ch]) fid.close() @@ -691,10 +699,12 @@ class ibmppp: fid.close() if self.__rank!=self.__nproc-1: self.__comm.send(True, dest=self.__rank+1, tag=1) + # Ensure that no processor is still working on the file + self.__comm.Barrier() # Create an XDMF file for the field if xdmf and self.__rank==0: if verbose: - print('[saveFieldSingleFile] generating XDMF',end='\r') + print('[saveFieldSingleFile] generating XDMF'+4*'\t',end='\r') # Construct XDMF filename file_xdmf = filename+'_'+key+'.xdmf' filename_h5 = os.path.basename(file_out) @@ -705,10 +715,6 @@ class ibmppp: fid.write(' \n') fid.write(' \n'.format(key)) #fid.write(' \n') fid.write(' \n') fid.close() + # Clear carriage return output if verbose + if self.__rank==0 and verbose: + print(10*'\t',end='\r') + return None def saveStatistics(self,filename): '''Writes all gathered statistics to a h5 file.''' From 87f65fecfd4ededb209e7e27e9409238f9f2b98c Mon Sep 17 00:00:00 2001 From: Michael Krayer Date: Wed, 12 May 2021 23:37:12 +0200 Subject: [PATCH 12/12] restructured package. implemented some abstraction layers --- python/{ucf/tools => programs}/ucf_file_diff | 0 .../ucf_file_diff_single_multi | 0 .../tools => programs}/ucf_file_integrity | 0 python/{ => programs}/ucftar_downsampler | 0 .../{ucf/tools => programs}/ucftar_mpiio_pack | 0 .../tools => programs}/ucftar_mpiio_unpack | 0 python/ucf/__init__.py | 3 + python/ucf/abstraction.py | 100 ++++++++++ python/ucf/{ucf.py => base.py} | 147 +------------- python/ucf/tools.py | 183 ++++++++++++++++++ 10 files changed, 287 insertions(+), 146 deletions(-) rename python/{ucf/tools => programs}/ucf_file_diff (100%) rename python/{ucf/tools => programs}/ucf_file_diff_single_multi (100%) rename python/{ucf/tools => programs}/ucf_file_integrity (100%) rename python/{ => programs}/ucftar_downsampler (100%) rename python/{ucf/tools => programs}/ucftar_mpiio_pack (100%) rename python/{ucf/tools => programs}/ucftar_mpiio_unpack (100%) create mode 100644 python/ucf/__init__.py create mode 100644 python/ucf/abstraction.py rename python/ucf/{ucf.py => base.py} (80%) create mode 100644 python/ucf/tools.py diff --git a/python/ucf/tools/ucf_file_diff b/python/programs/ucf_file_diff similarity index 100% rename from python/ucf/tools/ucf_file_diff rename to python/programs/ucf_file_diff diff --git a/python/ucf/tools/ucf_file_diff_single_multi b/python/programs/ucf_file_diff_single_multi similarity index 100% rename from python/ucf/tools/ucf_file_diff_single_multi rename to python/programs/ucf_file_diff_single_multi diff --git a/python/ucf/tools/ucf_file_integrity b/python/programs/ucf_file_integrity similarity index 100% rename from python/ucf/tools/ucf_file_integrity rename to python/programs/ucf_file_integrity diff --git a/python/ucftar_downsampler b/python/programs/ucftar_downsampler similarity index 100% rename from python/ucftar_downsampler rename to python/programs/ucftar_downsampler diff --git a/python/ucf/tools/ucftar_mpiio_pack b/python/programs/ucftar_mpiio_pack similarity index 100% rename from python/ucf/tools/ucftar_mpiio_pack rename to python/programs/ucftar_mpiio_pack diff --git a/python/ucf/tools/ucftar_mpiio_unpack b/python/programs/ucftar_mpiio_unpack similarity index 100% rename from python/ucf/tools/ucftar_mpiio_unpack rename to python/programs/ucftar_mpiio_unpack diff --git a/python/ucf/__init__.py b/python/ucf/__init__.py new file mode 100644 index 0000000..52782ba --- /dev/null +++ b/python/ucf/__init__.py @@ -0,0 +1,3 @@ +from .base import UCF +from .tools import read_grid, read_procgrid, read_chunk, read_particles, grid_chunk +from .abstraction import Ustar, UCFSnapshot, ChunkIterator \ No newline at end of file diff --git a/python/ucf/abstraction.py b/python/ucf/abstraction.py new file mode 100644 index 0000000..e05dec7 --- /dev/null +++ b/python/ucf/abstraction.py @@ -0,0 +1,100 @@ +class ChunkIterator: + def __init__(self,snapshot,dset=-1,keep_ghost=True): + self.snapshot = snapshot + self.dset = dset + self.keep_ghost = keep_ghost + self.iter_rank = 0 + def __iter__(self): + self.iter_rank = 0 + return self + def __next__(self): + if self.iter_rank0: - col['rank'] = ioffset; ioffset+=1 - if ihybrid>0: - col['id'] = ioffset; ioffset+=1 - col['x'] = ioffset; ioffset+=1 - col['y'] = ioffset; ioffset+=1 - col['z'] = ioffset; ioffset+=1 - col['r'] = ioffset; ioffset+=1 - col['rho']= ioffset; ioffset+=1 - col['ax'] = ioffset; ioffset+=1 - col['ay'] = ioffset; ioffset+=1 - col['az'] = ioffset; ioffset+=1 - col['u'] = ioffset; ioffset+=1 - col['v'] = ioffset; ioffset+=1 - col['w'] = ioffset; ioffset+=1 - col['ox'] = ioffset; ioffset+=1 - col['oy'] = ioffset; ioffset+=1 - col['oz'] = ioffset; ioffset+=1 - col['fx'] = ioffset; ioffset+=1 - col['fy'] = ioffset; ioffset+=1 - col['fz'] = ioffset; ioffset+=1 - col['tx'] = ioffset; ioffset+=1 - col['ty'] = ioffset; ioffset+=1 - col['tz'] = ioffset; ioffset+=1 - if idem>0: - col['fxc'] = ioffset; ioffset+=1 - col['fyc'] = ioffset; ioffset+=1 - col['fzc'] = ioffset; ioffset+=1 - col['txc'] = ioffset; ioffset+=1 - col['tyc'] = ioffset; ioffset+=1 - col['tzc'] = ioffset; ioffset+=1 - if iscal>0: - for ii in range(0,iscal): - col['s'+str(ii)] = ioffset; ioffset+=1 - col['q'+str(ii)] = ioffset; ioffset+=1 - return col + self.__currentSetNumElements = 0 \ No newline at end of file diff --git a/python/ucf/tools.py b/python/ucf/tools.py new file mode 100644 index 0000000..461e2d7 --- /dev/null +++ b/python/ucf/tools.py @@ -0,0 +1,183 @@ +def read_grid(file,verbosity=False,debug=False): + from .base import UCF + obj = UCF(file=file,verbosity=verbosity,debug=debug) + output = [] + for iset in range(0,obj.NumDataset): + (data,params) = obj.readSet(step=1,dset=iset+1) + nx = params[0] + ny = params[1] + nz = params[2] + output.append((data[0:nx],data[nx:nx+ny],data[nx+ny:nx+ny+nz])) + #if obj.UCFVersion<2: + if obj.NumDataset<5: + output.extend(output[-1:]) + obj.close() + return output + +def read_procgrid(file,verbosity=False,debug=False): + from .base import UCF + obj = UCF(file=file,verbosity=verbosity,debug=debug) + output = [] + for iset in range(0,obj.NumDataset): + (data,params) = obj.readSet(step=1,dset=iset+1) + nxp = params[0] + nyp = params[1] + nzp = params[2] + output.append(( + data[0:nxp], # ibeg + data[nxp:2*nxp], # iend + data[2*nxp:2*nxp+nyp], # jbeg + data[2*nxp+nyp:2*nxp+2*nyp], # jend + data[2*nxp+2*nyp:2*nxp+2*nyp+nzp], # kbeg + data[2*nxp+2*nyp+nzp:2*nxp+2*nyp*2*nzp] # kend + )) + #if obj.UCFVersion<2: + if obj.NumDataset<5: + output.extend(output[-1:]) + obj.close() + return output + +def read_chunk(file,step=1,dset=-1,keep_ghost=True,verbosity=False,debug=False): + from .base import UCF + obj = UCF(file=file,verbosity=verbosity,debug=debug) + if not isinstance(dset,list): + if dset==-1: + dset = range(1,obj.NumDataset+1) # fix that maybe later (this is maximum over all timesteps) + else: + dset = [dset] + output = [] + for ii in dset: + tmp = dict() + (data,params) = obj.readSet(step=step,dset=ii) + tmp['ighost'] = params[0] + tmp['ibeg'] = params[1] + tmp['jbeg'] = params[2] + tmp['kbeg'] = params[3] + tmp['nxl'] = params[4] + tmp['nyl'] = params[5] + tmp['nzl'] = params[6] + tmp['nx'] = params[7] + tmp['ny'] = params[8] + tmp['nz'] = params[9] + tmp['data'] = data.reshape((tmp['nxl']+2*tmp['ighost'], + tmp['nyl']+2*tmp['ighost'], + tmp['nzl']+2*tmp['ighost']), + order='F') + tmp['rank'] = obj.IORank[0] + tmp['rankijk']= obj.IORank[1:] + if not keep_ghost and ighost: + tmp['data'] = tmp['data'][1:-1,1:-1,1:-1] + tmp['ghost'] = 0 + output.append(tmp) + obj.close() + return output + +def read_particles(file,step=-1,verbosity=False,debug=False): + from .base import UCF + import numpy + obj = UCF(file=file,verbosity=verbosity,debug=debug) + if not isinstance(step,list): + if step==-1: + step = range(1,obj.NumTimestep+1) + else: + step = [step] + # The output will be the following: + # 1) numpy array with dimension (ncol,np,ntime) + # 2) dictionary which specifies the columns + # We read the data step by step in a list, which is then converted to a 3D array + pp = [] + for ii in step: + (data,params) = obj.readSet(step=ii,dset=1) + npart = params[0] + ncol = params[1] + ncol_rank = params[2] + ncol_hybrid = params[3] + ncol_dem = params[4] + ncol_scalar = params[5] + nscal = ncol_scalar//2 + pp.append(data.reshape((ncol,npart),order='F')) + # Close UCF obeject + obj.close() + # Convert list of 2D arrays to 3D array + pp = numpy.stack(pp,axis=2) + # Create the dictionary + col = colmap_from_flags(ncol_rank,ncol_hybrid,ncol_dem,nscal) + # Return result + return (pp,col) + +def colmap_from_flags(irank,ihybrid,idem,iscal): + '''Creates a dictionary which specifies the columns of a particle array.''' + col = {} + ioffset = 0 + if irank>0: + col['rank'] = ioffset; ioffset+=1 + if ihybrid>0: + col['id'] = ioffset; ioffset+=1 + col['x'] = ioffset; ioffset+=1 + col['y'] = ioffset; ioffset+=1 + col['z'] = ioffset; ioffset+=1 + col['r'] = ioffset; ioffset+=1 + col['rho']= ioffset; ioffset+=1 + col['ax'] = ioffset; ioffset+=1 + col['ay'] = ioffset; ioffset+=1 + col['az'] = ioffset; ioffset+=1 + col['u'] = ioffset; ioffset+=1 + col['v'] = ioffset; ioffset+=1 + col['w'] = ioffset; ioffset+=1 + col['ox'] = ioffset; ioffset+=1 + col['oy'] = ioffset; ioffset+=1 + col['oz'] = ioffset; ioffset+=1 + col['fx'] = ioffset; ioffset+=1 + col['fy'] = ioffset; ioffset+=1 + col['fz'] = ioffset; ioffset+=1 + col['tx'] = ioffset; ioffset+=1 + col['ty'] = ioffset; ioffset+=1 + col['tz'] = ioffset; ioffset+=1 + if idem>0: + col['fxc'] = ioffset; ioffset+=1 + col['fyc'] = ioffset; ioffset+=1 + col['fzc'] = ioffset; ioffset+=1 + col['txc'] = ioffset; ioffset+=1 + col['tyc'] = ioffset; ioffset+=1 + col['tzc'] = ioffset; ioffset+=1 + if iscal>0: + for ii in range(0,iscal): + col['s'+str(ii)] = ioffset; ioffset+=1 + col['q'+str(ii)] = ioffset; ioffset+=1 + return col + +def grid_chunk(grid_global,chunk): + import numpy + xg,yg,zg = grid_global + # Shift indices so that they start from zero + ib = chunk['ibeg']-1 + jb = chunk['jbeg']-1 + kb = chunk['kbeg']-1 + nxl = chunk['nxl'] + nyl = chunk['nyl'] + nzl = chunk['nzl'] + ighost = chunk['ighost'] + if ighost: + nxg = len(xg) + nyg = len(yg) + nzg = len(zg) + dx = xg[2]-xg[1] + dy = yg[2]-yg[1] + dz = zg[2]-zg[1] + xl = numpy.zeros(nxl+2) + yl = numpy.zeros(nyl+2) + zl = numpy.zeros(nzl+2) + xl[0] = xg[ib]-dx + xl[1:-1] = xg[ib:ib+nxl] + xl[-1] = xg[ib+nxl-1]+dx + yl[0] = yg[jb]-dy + yl[1:-1] = yg[jb:jb+nyl] + yl[-1] = yg[jb+nyl-1]+dy + zl[0] = zg[kb]-dz + zl[1:-1] = zg[kb:kb+nzl] + zl[-1] = zg[kb+nzl-1]+dz + else: + xl = xg[ibeg:ib+nxl-1] + yl = yg[jbeg:jb+nyl-1] + zl = zg[kbeg:kb+nzl-1] + return (xl,yl,zl) \ No newline at end of file