diff --git a/matlab/read_domain_legacy.m b/matlab/read_domain_legacy.m new file mode 100644 index 0000000..61c2b73 --- /dev/null +++ b/matlab/read_domain_legacy.m @@ -0,0 +1,16 @@ +function [a,b,c,d,e,f] = read_domain_legacy(file) + % [a,b,c,d,e,f] = read_domain_legacy(file) + % Reads domain boundaries from 'current_domain.tmp' + % Input + % file path to 'current_domain.tmp' + % Output + % a,b,c,d,e,f domain boundaries + + border=load(file); + a=border(1); + b=border(2); + c=border(3); + d=border(4); + e=border(5); + f=border(6); +end diff --git a/matlab/read_grid_legacy.m b/matlab/read_grid_legacy.m new file mode 100644 index 0000000..17300c1 --- /dev/null +++ b/matlab/read_grid_legacy.m @@ -0,0 +1,21 @@ +function [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,... + nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,nxp,nyp,nzp] = read_grid_legacy(file,a,b,c,d,e,f,x_periodic,y_periodic,z_periodic) + % [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,... + % nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,nxp,nyp,nzp] = read_grid_legacy(file,a,b,c,d,e,f,x_periodic,y_periodic,z_periodic) + % Reconstructs grid from 'current_points.tmp'. + % Input + % file path to 'current_points.tmp' + % a,b,c,d,e,f domain boundaries + % x/y/z_periodic periodicity + % Output + % xu,yu,zu,... grid vectors + % nxu,nyu,nzu,... number of points + + nn=load(file); + nxp=nn(1); + nyp=nn(2); + nzp=nn(3); + [nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,... + xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp]=... + generate_grid(a,b,c,d,e,f,nxp,nyp,nzp,x_periodic,y_periodic,z_periodic); +end diff --git a/matlab/read_procgrid_legacy.m b/matlab/read_procgrid_legacy.m new file mode 100644 index 0000000..266b997 --- /dev/null +++ b/matlab/read_procgrid_legacy.m @@ -0,0 +1,83 @@ +function [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(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(file) + % Reads processor grids from 'current_proc_grid.tmp'. + % Input + % file path to 'current_proc_grid.tmp' + % Output + % ibegu,iendu,jbegu,jendu,kbegu,kendu processor grid u + % ibegv,iendv,jbegv,jendv,kbegv,kendv processor grid v + % ibegw,iendw,jbegw,jendw,kbegw,kendw processor grid w + % ibegp,iendp,jbegp,jendp,kbegp,kendp processor grid p + % nxprocs,nyprocs,nzprocs number of processors + + proc_grid=load(file); + nxprocs = proc_grid(1,1); + nyprocs = proc_grid(1,2); + nzprocs = proc_grid(1,3); + + ibegu = zeros(nxprocs,1); + iendu = zeros(nxprocs,1); + ibegv = zeros(nxprocs,1); + iendv = zeros(nxprocs,1); + ibegw = zeros(nxprocs,1); + iendw = zeros(nxprocs,1); + ibegp = zeros(nxprocs,1); + iendp = zeros(nxprocs,1); + + jbegu = zeros(nyprocs,1); + jendu = zeros(nyprocs,1); + jbegv = zeros(nyprocs,1); + jendv = zeros(nyprocs,1); + jbegw = zeros(nyprocs,1); + jendw = zeros(nyprocs,1); + jbegp = zeros(nyprocs,1); + jendp = zeros(nyprocs,1); + + kbegu = zeros(nzprocs,1); + kendu = zeros(nzprocs,1); + kbegv = zeros(nzprocs,1); + kendv = zeros(nzprocs,1); + kbegw = zeros(nzprocs,1); + kendw = zeros(nzprocs,1); + kbegp = zeros(nzprocs,1); + kendp = zeros(nzprocs,1); + + for ii=1:nxprocs + ibegu(ii)=proc_grid(1+ii,1); + iendu(ii)=proc_grid(1+ii,2); + ibegv(ii)=proc_grid(1+ii,3); + iendv(ii)=proc_grid(1+ii,4); + ibegw(ii)=proc_grid(1+ii,5); + iendw(ii)=proc_grid(1+ii,6); + ibegp(ii)=proc_grid(1+ii,7); + iendp(ii)=proc_grid(1+ii,8); + end + for ii=1:nyprocs + jbegu(ii)=proc_grid(1+nxprocs+ii,1); + jendu(ii)=proc_grid(1+nxprocs+ii,2); + jbegv(ii)=proc_grid(1+nxprocs+ii,3); + jendv(ii)=proc_grid(1+nxprocs+ii,4); + jbegw(ii)=proc_grid(1+nxprocs+ii,5); + jendw(ii)=proc_grid(1+nxprocs+ii,6); + jbegp(ii)=proc_grid(1+nxprocs+ii,7); + jendp(ii)=proc_grid(1+nxprocs+ii,8); + end + for ii=1:nzprocs + kbegu(ii)=proc_grid(1+nxprocs+nyprocs+ii,1); + kendu(ii)=proc_grid(1+nxprocs+nyprocs+ii,2); + kbegv(ii)=proc_grid(1+nxprocs+nyprocs+ii,3); + kendv(ii)=proc_grid(1+nxprocs+nyprocs+ii,4); + kbegw(ii)=proc_grid(1+nxprocs+nyprocs+ii,5); + kendw(ii)=proc_grid(1+nxprocs+nyprocs+ii,6); + kbegp(ii)=proc_grid(1+nxprocs+nyprocs+ii,7); + kendp(ii)=proc_grid(1+nxprocs+nyprocs+ii,8); + end +end diff --git a/matlab/read_uvwp_chunk_legacy.m b/matlab/read_uvwp_chunk_legacy.m new file mode 100644 index 0000000..a90fbfa --- /dev/null +++ b/matlab/read_uvwp_chunk_legacy.m @@ -0,0 +1,149 @@ +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