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');