read arbitrary field from legacy format

This commit is contained in:
Michael Stumpf (ifhcluster) 2019-04-05 11:08:38 +02:00
parent 78b67be465
commit 356beef2d9
1 changed files with 157 additions and 0 deletions

View File

@ -0,0 +1,157 @@
function [data,ib,jb,kb,nxl,nyl,nzl,ighost] = read_field_chunk_legacy(file,fproc,ighost,field,varargin)
% [data,ib,jb,kb,nxl,nyl,nzl,ighost] = read_field_chunk_legacy(file,fproc,ighost,field,varargin)
% Reads arbitrary field from one processor chunk.
% Input
% file path to chunk
% fproc path to 'current_proc_grid.tmp'
% ighost data written with ghost cells?
% field field to be read
% {'u','v','w','p','s1','s2',...}
% ? ghost keep ghost cells? (default: 1)
% ? verbosity verbose output? (default: 0)
% ? nparam number of parameters in header (default: 10)
% ? precision precision of data (default: 'float64')
% 'float32'
% 'float64'
% ? endian endianess of file (default: 'a')
% ? reclen FORTRAN record length in bytes (default: 4)
% Output
% data partial field with dim(nxl+2*ighost,nyl+2*ighost,nzl+2*ighost)
% ib,jb,kb global index of first grid point of the partial field w/o ghost cells
% nxl,nyl,nzl local field size w/o ghost cells
% ighost ghost cell flag
% Parse optional input arguments
par = inputParser;
addParamValue(par,'ghost',1,@isnumeric);
addParamValue(par,'verbosity',0,@isnumeric);
addParamValue(par,'nparam',10,@isnumeric);
addParamValue(par,'precision','float64',@ischar);
addParamValue(par,'endian','a',@ischar);
addParamValue(par,'reclen',4,@isnumeric);
parse(par,varargin{:});
keepghost = par.Results.ghost;
verbosity = par.Results.verbosity;
nparam = par.Results.nparam;
precision = par.Results.precision;
endian = par.Results.endian;
reclen = par.Results.reclen;
% Parse field
if ischar(field)
switch field(1)
case 'u'; ifield=1; fbase='uvwp';
case 'v'; ifield=2; fbase='uvwp';
case 'w'; ifield=3; fbase='uvwp';
case 'p'; ifield=4; fbase='uvwp';
case 's'; ifield=str2double(field(2:end)); fbase='scal';
end
else
error('field must be of type char');
end
% First read processor grid from tmp file
[ibegu,iendu,jbegu,jendu,kbegu,kendu,...
ibegv,iendv,jbegv,jendv,kbegv,kendv,...
ibegw,iendw,jbegw,jendw,kbegw,kendw,...
ibegp,iendp,jbegp,jendp,kbegp,kendp,...
nxprocs,nyprocs,nzprocs] = read_procgrid_legacy(fproc);
% Determine processor rank from filename
[idxbeg,idxend] = regexp(file,'[\d]+$');
if isempty(idxbeg)
error('Invalid file: does not contain rank. %s',fuvwp);
end
iproc = str2double(file(idxbeg:idxend));
% Determine local array sizes
ixproc = floor(iproc/(nyprocs*nzprocs));
iyproc = mod(floor(iproc/nzprocs),nyprocs);
izproc = mod(iproc,nzprocs);
ibu = ibegu(ixproc+1);
jbu = jbegu(iyproc+1);
kbu = kbegu(izproc+1);
nxul = iendu(ixproc+1)-ibu+1;
nyul = jendu(iyproc+1)-jbu+1;
nzul = kendu(izproc+1)-kbu+1;
ibv = ibegv(ixproc+1);
jbv = jbegv(iyproc+1);
kbv = kbegv(izproc+1);
nxvl = iendv(ixproc+1)-ibv+1;
nyvl = jendv(iyproc+1)-jbv+1;
nzvl = kendv(izproc+1)-kbv+1;
ibw = ibegw(ixproc+1);
jbw = jbegw(iyproc+1);
kbw = kbegw(izproc+1);
nxwl = iendw(ixproc+1)-ibw+1;
nywl = jendw(iyproc+1)-jbw+1;
nzwl = kendw(izproc+1)-kbw+1;
ibp = ibegp(ixproc+1);
jbp = jbegp(iyproc+1);
kbp = kbegp(izproc+1);
nxpl = iendp(ixproc+1)-ibp+1;
nypl = jendp(iyproc+1)-jbp+1;
nzpl = kendp(izproc+1)-kbp+1;
% Convert 'precision' into bytes
switch precision
case 'float32'; nprecision = 4;
case 'float64'; nprecision = 8;
otherwise; error('Invalid precision: %s',precision)
end
% Open uvwp file
fid = fopen(fuvwp,'r',endian);
if fid<0
error('File not found: %s',fuvwp);
end
% Determine header size
nheader = 4*reclen+(4+nparam)*nprecision;
% Skip datasets until correct one is reached
for iset=1:ifield
switch fbase
case 'uvwp'
if iset==1; nxl=nxul; nyl=nyul; nzl=nzul; end
if iset==2; nxl=nxvl; nyl=nyvl; nzl=nzvl; end
if iset==3; nxl=nxwl; nyl=nywl; nzl=nzwl; end
if iset==4; nxl=nxpl; nyl=nypl; nzl=nzpl; end
case 'scal'
nxl=nxpl; nyl=nypl; nzl=nzpl;
end
ndata = (nxl+2*ighost)*(nyl+2*ighost)*(nzl+2*ighost);
if iset~=ifield
nskip = nheader+2*reclen+ndata*nprecision;
fseek(fid,nskip,'cof');
else
fseek(fid,nheader+reclen,'cof');
data = fread(fid,ndata,precision);
data = reshape(data,nxl,nyl,nzl);
end
end
% Get processor bounds
switch fbase
case 'uvwp'
if ifield==1; ib=ibu; jb=jbu; kb=kbu; nxl=nxul; nyl=nyul; nzl=nzul; end
if ifield==2; ib=ibv; jb=jbv; kb=kbv; nxl=nxvl; nyl=nyvl; nzl=nzvl; end
if ifield==3; ib=ibw; jb=jbw; kb=kbw; nxl=nxwl; nyl=nywl; nzl=nzwl; end
if ifield==4; ib=ibp; jb=jbp; kb=kbp; nxl=nxpl; nyl=nypl; nzl=nzpl; end
case 'scal'
ib=ibp; jb=jbp; kb=kbp; nxl=nxpl; nyl=nypl; nzl=nzpl;
end
% Close file
fclose(fid);
% Remove ghosts if necessary
if ighost && ~keepghost
data = data(2:end-1,2:end-1,2:end-1);
ighost = 0;
end
end