ucftools/matlab/read_uvwp_chunk_legacy.m

150 lines
5.1 KiB
Matlab

function [u,ibu,jbu,kbu,nxul,nyul,nzul,...
v,ibv,jbv,kbv,nxvl,nyvl,nzvl,...
w,ibw,jbw,kbw,nxwl,nywl,nzwl,...
p,ibp,jbp,kbp,nxpl,nypl,nzpl,ighost] = read_uvwp_chunk_legacy(fuvwp,fproc,ighost,varargin)
% [u,ibu,jbu,kbu,nxul,nyul,nzul,...
% v,ibv,jbv,kbv,nxvl,nyvl,nzvl,...
% w,ibw,jbw,kbw,nxwl,nywl,nzwl,...
% p,ibp,jbp,kbp,nxpl,nypl,nzpl,ighost] = read_uvwp_chunk_legacy(fuvwp,fproc,ighost,varargin)
% Reads u,v,w,p from one processor chunk.
% Input
% fuvwp path to uvwp chunk
% fproc path to 'current_proc_grid.tmp'
% ighost data written with ghost cells?
% ? 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
% u,v,w,p partial fields with dim(nxl,nyl,nzl)+2*ighost
% ibu,jbu,kbu,... global index of first grid point of the partial field w/o ghost cells
% nxul,nyul,nzul,... 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;
% 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(fuvwp,'[\d]+$');
if isempty(idxbeg)
error('Invalid file: does not contain rank. %s',fuvwp);
end
iproc = str2double(fuvwp(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
% Skip headers and read fields
nheader = 4*reclen+(4+nparam)*nprecision;
fseek(fid,nheader+reclen,'bof');
nx = (nxul+2*ighost);
ny = (nyul+2*ighost);
nz = (nzul+2*ighost);
u = fread(fid,nx*ny*nz,precision);
u = reshape(u,nx,ny,nz);
fseek(fid,nheader+2*reclen,'cof');
nx = (nxvl+2*ighost);
ny = (nyvl+2*ighost);
nz = (nzvl+2*ighost);
v = fread(fid,nx*ny*nz,precision);
v = reshape(v,nx,ny,nz);
fseek(fid,nheader+2*reclen,'cof');
nx = (nxwl+2*ighost);
ny = (nywl+2*ighost);
nz = (nzwl+2*ighost);
w = fread(fid,nx*ny*nz,precision);
w = reshape(w,nx,ny,nz);
fseek(fid,nheader+2*reclen,'cof');
nx = (nxpl+2*ighost);
ny = (nypl+2*ighost);
nz = (nzpl+2*ighost);
p = fread(fid,nx*ny*nz,precision);
p = reshape(p,nx,ny,nz);
fseek(fid,reclen,'cof');
% Test for eof
fread(fid,1,'char');
if ~feof(fid)
warning('End-of-file not reached: data might be corrupt.');
end
% Close file
fclose(fid);
% Remove ghosts if necessary
if ighost && ~keepghost
u = u(2:end-1,2:end-1,2:end-1);
v = v(2:end-1,2:end-1,2:end-1);
w = w(2:end-1,2:end-1,2:end-1);
p = p(2:end-1,2:end-1,2:end-1);
ighost = 0;
end
end