From 6f618f3d5db2931e60b79cd26c2834ed6fd6ad06 Mon Sep 17 00:00:00 2001 From: "Michael Stumpf (ifhcluster)" Date: Tue, 26 Mar 2019 12:26:12 +0000 Subject: [PATCH] interpolation tool --- matlab/interp_uvwp_ucf.m | 292 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 292 insertions(+) create mode 100644 matlab/interp_uvwp_ucf.m diff --git a/matlab/interp_uvwp_ucf.m b/matlab/interp_uvwp_ucf.m new file mode 100644 index 0000000..b4ae96b --- /dev/null +++ b/matlab/interp_uvwp_ucf.m @@ -0,0 +1,292 @@ +function [] = interp_uvwp_ucf(hucf,dout,nxpn,nypn,nzpn,nxprocsn,nyprocsn,nzprocsn,varargin) + % [] = interp_uvwp_ucf(hucf,dout,nxpn,nypn,nzpn,nxprocsn,nyprocsn,nzprocsn,varargin) + % Interpolates a given flow field onto a new grid. + % This function reads every chunk of the original field only once, while + % trying to minimize the memory footprint by writing the new chunks as + % soon as possible to disk. Memory for the new chunks is only allocated if + % they contain partial data. + % Input + % hucf UCF handle (ustar,ucfmulti) + % dout output directory + % nxpn,... number of grid points of the new pressure grid + % nxprocs,... number of processors of the new processor grid + % ? verbosity verbose output? (includes estimation of the current memory footprint) [default: 0] + % ? warnings do some sanity checks and throw warnings? [default: 0] + % ? filebase basename of the output files [default: uvwp] + % ? iseq sequence number of the output files [default: 0] + % ? time simulation time of the output files [default: 0.0] + + % Parse input + par = inputParser; + addParamValue(par,'verbosity',0,@isnumeric); + addParamValue(par,'warnings',0,@isnumeric); + addParamValue(par,'filebase','uvwp',@ischar); + addParamValue(par,'iseq',0,@isnumeric); + addParamValue(par,'time',0.0,@isnumeric); + parse(par,varargin{:}); + flag_verb = par.Results.verbosity; + flag_warn = par.Results.warning; + filebase = par.Results.filebase; + iseq = par.Results.iseq; + simtime = par.Results.time; + + % Read info from tar-archive + [xuo,yuo,zuo,xvo,yvo,zvo,xwo,ywo,zwo,xpo,ypo,zpo] = read_grid_ucf(hucf); + [params] = read_parameters_ucf(hucf); + nprocso = params.parallel.nprocs; + a = params.domain.a; + b = params.domain.b; + c = params.domain.c; + d = params.domain.d; + e = params.domain.e; + f = params.domain.f; + xperiodic = params.domain.xperiodic; + yperiodic = params.domain.yperiodic; + zperiodic = params.domain.zperiodic; + + [nxun,nyun,nzun,nxvn,nyvn,nzvn,nxwn,nywn,nzwn,... + xun,yun,zun,xvn,yvn,zvn,xwn,ywn,zwn,xpn,ypn,zpn,dxn,dyn,dzn]=... + generate_grid(a,b,c,d,e,f,nxpn,nypn,nzpn,xperiodic,yperiodic,zperiodic); + [ibegun,iendun,jbegun,jendun,kbegun,kendun,... + ibegvn,iendvn,jbegvn,jendvn,kbegvn,kendvn,... + ibegwn,iendwn,jbegwn,jendwn,kbegwn,kendwn,... + ibegpn,iendpn,jbegpn,jendpn,kbegpn,kendpn]=generate_procgrid(... + nxun,nyun,nzun,nxvn,nyvn,nzvn,nxwn,nywn,nzwn,... + nxpn,nypn,nzpn,nxprocsn,nyprocsn,nzprocsn); + + if flag_warn + if (~xperiodic && mod(nxpn,2)==0) || (~yperiodic && mod(nypn,2)==0) || (~zperiodic && mod(nzpn,2)==0) + warning('Pressure grid should contain extra point in non-periodic directions.') + end + if abs((dxn-dyn)/dxn)>1e-8 || abs((dyn-dzn)/dyn)>1e-8 + warning('New mesh is not equidistant: %f, %f, %f.',dxn,dyn,dzn); + end + if mod(nxpn-~xperiodic,64)~=0 + warning('nxpn is not a multiple of 64: MGD performance might be bad.'); + end + if mod(nypn-~yperiodic,64)~=0 + warning('nypn is not a multiple of 64: MGD performance might be bad.'); + end + if mod(nzpn-~xperiodic,64)~=0 + warning('nzpn is not a multiple of 64: MGD performance might be bad.'); + end + end + + % Cell arrays which hold interpolation results + un = cell(nxprocsn,nyprocsn,nzprocsn); + vn = cell(nxprocsn,nyprocsn,nzprocsn); + wn = cell(nxprocsn,nyprocsn,nzprocsn); + pn = cell(nxprocsn,nyprocsn,nzprocsn); + + % Cell arrays which indicate if points have been interpolated yet + un_ind = cell(nxprocsn,nyprocsn,nzprocsn); + vn_ind = cell(nxprocsn,nyprocsn,nzprocsn); + wn_ind = cell(nxprocsn,nyprocsn,nzprocsn); + pn_ind = cell(nxprocsn,nyprocsn,nzprocsn); + + % Logical arrays which indicate whether memory has been allocated + un_init = false(nxprocsn,nyprocsn,nzprocsn); + vn_init = false(nxprocsn,nyprocsn,nzprocsn); + wn_init = false(nxprocsn,nyprocsn,nzprocsn); + pn_init = false(nxprocsn,nyprocsn,nzprocsn); + + % Logical array which indicate whether a chunk has been finalized + chunk_final = false(nxprocsn,nyprocsn,nzprocsn); + + imem = 0; + maxmem = 0; + for iproco=0:nprocso-1 + if flag_verb + fprintf('Processing original chunk %5d ',iproco); + end + [uo,ibuo,jbuo,kbuo,nxuol,nyuol,nzuol,... + vo,ibvo,jbvo,kbvo,nxvol,nyvol,nzvol,... + wo,ibwo,jbwo,kbwo,nxwol,nywol,nzwol,... + po,ibpo,jbpo,kbpo,nxpol,nypol,nzpol,ighost] = ... + read_uvwp_chunk_ucf(hucf,'rank',iproco); + imem = imem+(numel(uo)+numel(vo)+numel(wo)+numel(po))*8; + if flag_verb + fprintf('(current memory footprint: %7.0f MiB)\n',imem/(1024*1024)); + end + maxmem = max(imem,maxmem); + + [xuol,yuol,zuol] = grid_chunk(xuo,yuo,zuo,ibuo,jbuo,kbuo,nxuol,nyuol,nzuol,ighost); + [xvol,yvol,zvol] = grid_chunk(xvo,yvo,zvo,ibvo,jbvo,kbvo,nxvol,nyvol,nzvol,ighost); + [xwol,ywol,zwol] = grid_chunk(xwo,ywo,zwo,ibwo,jbwo,kbwo,nxwol,nywol,nzwol,ighost); + [xpol,ypol,zpol] = grid_chunk(xpo,ypo,zpo,ibpo,jbpo,kbpo,nxpol,nypol,nzpol,ighost); + + Fu = griddedInterpolant({xuol,yuol,zuol},uo); + Fv = griddedInterpolant({xvol,yvol,zvol},vo); + Fw = griddedInterpolant({xwol,ywol,zwol},wo); + Fp = griddedInterpolant({xpol,ypol,zpol},po); + + for izprocn=0:nzprocsn-1 + for iyprocn=0:nyprocsn-1 + for ixprocn=0:nxprocsn-1 + if chunk_final(ixprocn+1,iyprocn+1,izprocn+1) + continue; + end + ibun = ibegun(ixprocn+1); nxunl = iendun(ixprocn+1)-ibegun(ixprocn+1)+1; + jbun = jbegun(iyprocn+1); nyunl = jendun(iyprocn+1)-jbegun(iyprocn+1)+1; + kbun = kbegun(izprocn+1); nzunl = kendun(izprocn+1)-kbegun(izprocn+1)+1; + ibvn = ibegvn(ixprocn+1); nxvnl = iendvn(ixprocn+1)-ibegvn(ixprocn+1)+1; + jbvn = jbegvn(iyprocn+1); nyvnl = jendvn(iyprocn+1)-jbegvn(iyprocn+1)+1; + kbvn = kbegvn(izprocn+1); nzvnl = kendvn(izprocn+1)-kbegvn(izprocn+1)+1; + ibwn = ibegwn(ixprocn+1); nxwnl = iendwn(ixprocn+1)-ibegwn(ixprocn+1)+1; + jbwn = jbegwn(iyprocn+1); nywnl = jendwn(iyprocn+1)-jbegwn(iyprocn+1)+1; + kbwn = kbegwn(izprocn+1); nzwnl = kendwn(izprocn+1)-kbegwn(izprocn+1)+1; + ibpn = ibegpn(ixprocn+1); nxpnl = iendpn(ixprocn+1)-ibegpn(ixprocn+1)+1; + jbpn = jbegpn(iyprocn+1); nypnl = jendpn(iyprocn+1)-jbegpn(iyprocn+1)+1; + kbpn = kbegpn(izprocn+1); nzpnl = kendpn(izprocn+1)-kbegpn(izprocn+1)+1; + [xunl,yunl,zunl] = grid_chunk(xun,yun,zun,ibun,jbun,kbun,nxunl,nyunl,nzunl,0); + [xvnl,yvnl,zvnl] = grid_chunk(xvn,yvn,zvn,ibvn,jbvn,kbvn,nxvnl,nyvnl,nzvnl,0); + [xwnl,ywnl,zwnl] = grid_chunk(xwn,ywn,zwn,ibwn,jbwn,kbwn,nxwnl,nywnl,nzwnl,0); + [xpnl,ypnl,zpnl] = grid_chunk(xpn,ypn,zpn,ibpn,jbpn,kbpn,nxpnl,nypnl,nzpnl,0); + + % u + ixu = find(xunl>=xuol(1) & xunl<=xuol(end)); + iyu = find(yunl>=yuol(1) & yunl<=yuol(end)); + izu = find(zunl>=zuol(1) & zunl<=zuol(end)); + if ~isempty(ixu) && ~isempty(iyu) && ~isempty(izu) + if ~un_init(ixprocn+1,iyprocn+1,izprocn+1) + un{ixprocn+1,iyprocn+1,izprocn+1} = zeros(nxunl,nyunl,nzunl); + un_ind{ixprocn+1,iyprocn+1,izprocn+1} = false(nxunl,nyunl,nzunl); + un_init(ixprocn+1,iyprocn+1,izprocn+1) = true; + imem = imem+(nxunl*nyunl*nzunl)*(8+1); + maxmem = max(imem,maxmem); + if ~xperiodic + % staggered grid goes beyond boundary in non-periodic direction: + % if field is coarsened, no interpolation point will be found! + % Thus, just leave it and mark as processed + un_ind{ixprocn+1,iyprocn+1,izprocn+1}(xunlxuo(end),:,:) = true; + end + end + un{ixprocn+1,iyprocn+1,izprocn+1}(ixu,iyu,izu) = Fu({xunl(ixu),yunl(iyu),zunl(izu)}); + un_ind{ixprocn+1,iyprocn+1,izprocn+1}(ixu,iyu,izu) = true; + end + + % v + ixv = find(xvnl>=xvol(1) & xvnl<=xvol(end)); + iyv = find(yvnl>=yvol(1) & yvnl<=yvol(end)); + izv = find(zvnl>=zvol(1) & zvnl<=zvol(end)); + if ~isempty(ixv) && ~isempty(iyv) && ~isempty(izv) + if ~vn_init(ixprocn+1,iyprocn+1,izprocn+1) + vn{ixprocn+1,iyprocn+1,izprocn+1} = zeros(nxvnl,nyvnl,nzvnl); + vn_ind{ixprocn+1,iyprocn+1,izprocn+1} = false(nxvnl,nyvnl,nzvnl); + vn_init(ixprocn+1,iyprocn+1,izprocn+1) = true; + imem = imem+(nxunl*nyunl*nzunl)*(8+1); + maxmem = max(imem,maxmem); + if ~yperiodic + % staggered grid goes beyond boundary in non-periodic direction: + % if field is coarsened, no interpolation point will be found! + % Thus, just leave it and mark as processed + vn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:,yvnlyvo(end),:) = true; + end + end + vn{ixprocn+1,iyprocn+1,izprocn+1}(ixv,iyv,izv) = Fv({xvnl(ixv),yvnl(iyv),zvnl(izv)}); + vn_ind{ixprocn+1,iyprocn+1,izprocn+1}(ixv,iyv,izv) = true; + end + + % w + ixw = find(xwnl>=xwol(1) & xwnl<=xwol(end)); + iyw = find(ywnl>=ywol(1) & ywnl<=ywol(end)); + izw = find(zwnl>=zwol(1) & zwnl<=zwol(end)); + if ~isempty(ixw) && ~isempty(iyw) && ~isempty(izw) + if ~wn_init(ixprocn+1,iyprocn+1,izprocn+1) + wn{ixprocn+1,iyprocn+1,izprocn+1} = zeros(nxwnl,nywnl,nzwnl); + wn_ind{ixprocn+1,iyprocn+1,izprocn+1} = false(nxwnl,nywnl,nzwnl); + wn_init(ixprocn+1,iyprocn+1,izprocn+1) = true; + imem = imem+(nxunl*nyunl*nzunl)*(8+1); + maxmem = max(imem,maxmem); + if ~zperiodic + % staggered grid goes beyond boundary in non-periodic direction: + % if field is coarsened, no interpolation point will be found! + % Thus, just leave it and mark as processed + wn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:,:,zwnlzwo(end)) = true; + end + end + wn{ixprocn+1,iyprocn+1,izprocn+1}(ixw,iyw,izw) = Fw({xwnl(ixw),ywnl(iyw),zwnl(izw)}); + wn_ind{ixprocn+1,iyprocn+1,izprocn+1}(ixw,iyw,izw) = true; + end + + % p + ixp = find(xpnl>=xpol(1) & xpnl<=xpol(end)); + iyp = find(ypnl>=ypol(1) & ypnl<=ypol(end)); + izp = find(zpnl>=zpol(1) & zpnl<=zpol(end)); + if ~isempty(ixp) && ~isempty(iyp) && ~isempty(izp) + if ~pn_init(ixprocn+1,iyprocn+1,izprocn+1) + pn{ixprocn+1,iyprocn+1,izprocn+1} = zeros(nxpnl,nypnl,nzpnl); + pn_ind{ixprocn+1,iyprocn+1,izprocn+1} = false(nxpnl,nypnl,nzpnl); + pn_init(ixprocn+1,iyprocn+1,izprocn+1) = true; + imem = imem+(nxunl*nyunl*nzunl)*(8+1); + maxmem = max(imem,maxmem); + end + pn{ixprocn+1,iyprocn+1,izprocn+1}(ixp,iyp,izp) = Fp({xpnl(ixp),ypnl(iyp),zpnl(izp)}); + pn_ind{ixprocn+1,iyprocn+1,izprocn+1}(ixp,iyp,izp) = true; + end + + if un_init(ixprocn+1,iyprocn+1,izprocn+1) && ... + vn_init(ixprocn+1,iyprocn+1,izprocn+1) && ... + wn_init(ixprocn+1,iyprocn+1,izprocn+1) && ... + pn_init(ixprocn+1,iyprocn+1,izprocn+1) && ... + all(un_ind{ixprocn+1,iyprocn+1,izprocn+1}(:)) && ... + all(vn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:)) && ... + all(wn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:)) && ... + all(pn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:)) + ichunkn = ixprocn*nyprocsn*nzprocsn+iyprocn*nzprocsn+izprocn; + fout = sprintf('%s/%s_%04d.%05d',dout,filebase,iseq,ichunkn); + write_uvwp_chunk_ucf(fout,... + un{ixprocn+1,iyprocn+1,izprocn+1},ibun,jbun,kbun,nxunl,nyunl,nzunl,... + vn{ixprocn+1,iyprocn+1,izprocn+1},ibvn,jbvn,kbvn,nxvnl,nyvnl,nzvnl,... + wn{ixprocn+1,iyprocn+1,izprocn+1},ibwn,jbwn,kbwn,nxwnl,nywnl,nzwnl,... + pn{ixprocn+1,iyprocn+1,izprocn+1},ibpn,jbpn,kbpn,nxpnl,nypnl,nzpnl,... + simtime,0); + un{ixprocn+1,iyprocn+1,izprocn+1} = []; + vn{ixprocn+1,iyprocn+1,izprocn+1} = []; + wn{ixprocn+1,iyprocn+1,izprocn+1} = []; + pn{ixprocn+1,iyprocn+1,izprocn+1} = []; + un_ind{ixprocn+1,iyprocn+1,izprocn+1} = []; + vn_ind{ixprocn+1,iyprocn+1,izprocn+1} = []; + wn_ind{ixprocn+1,iyprocn+1,izprocn+1} = []; + pn_ind{ixprocn+1,iyprocn+1,izprocn+1} = []; + chunk_final(ixprocn+1,iyprocn+1,izprocn+1) = true; + imem = imem-(... + (nxunl*nyunl*nzunl)+... + (nxvnl*nyvnl*nzvnl)+... + (nxwnl*nywnl*nzwnl)+... + (nxpnl*nypnl*nzpnl) ... + )*(8+1); + end + end + end + end + imem = imem-(numel(uo)+numel(vo)+numel(wo)+numel(po))*8; + end + + % Check if everything has been finalized + for izprocn=0:nzprocsn-1 + for iyprocn=0:nyprocsn-1 + for ixprocn=0:nxprocsn-1 + if ~chunk_final(ixprocn+1,iyprocn+1,izprocn+1) + ichunkn = ixprocn*nyprocsn*nzprocsn+iyprocn*nzprocsn+izprocn; + fout = sprintf('%s/%s_%04d.%05d',dout,filebase,iseq,ichunkn); + if flag_warn + warning('Chunk %2d,%2d,%2d (rank: %5d) has not been finalized! Writing anyway...',ixprocn,iyprocn,izprocn,ichunkn); + end + write_uvwp_chunk_ucf(fout,... + un{ixprocn+1,iyprocn+1,izprocn+1},ibun,jbun,kbun,nxunl,nyunl,nzunl,... + vn{ixprocn+1,iyprocn+1,izprocn+1},ibvn,jbvn,kbvn,nxvnl,nyvnl,nzvnl,... + wn{ixprocn+1,iyprocn+1,izprocn+1},ibwn,jbwn,kbwn,nxwnl,nywnl,nzwnl,... + pn{ixprocn+1,iyprocn+1,izprocn+1},ibpn,jbpn,kbpn,nxpnl,nypnl,nzpnl,... + simetime,0); + end + end + end + end + if flag_verb + fprintf('Done. (maximum memory footprint: %7.0f MiB)\n',maxmem/(1024*1024)); + end +end