293 lines
14 KiB
Matlab
293 lines
14 KiB
Matlab
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.warnings;
|
|
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.geometry.a;
|
|
b = params.geometry.b;
|
|
c = params.geometry.c;
|
|
d = params.geometry.d;
|
|
e = params.geometry.e;
|
|
f = params.geometry.f;
|
|
xperiodic = params.geometry.xperiodic;
|
|
yperiodic = params.geometry.yperiodic;
|
|
zperiodic = params.geometry.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}(xunl<xuo(1),:,:) = true;
|
|
un_ind{ixprocn+1,iyprocn+1,izprocn+1}(xunl>xuo(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}(:,yvnl<yvo(1),:) = true;
|
|
vn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:,yvnl>yvo(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}(:,:,zwnl<zwo(1)) = true;
|
|
wn_ind{ixprocn+1,iyprocn+1,izprocn+1}(:,:,zwnl>zwo(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
|