precompute ghost cells: ndSparse as prerequisite

This commit is contained in:
Michael Stumpf (ifhcluster) 2019-03-29 09:10:25 +01:00
parent 6f618f3d5d
commit 57785b9cb9
1 changed files with 370 additions and 0 deletions

View File

@ -0,0 +1,370 @@
function [] = precompute_ghosts_uvwp_ucf(hucf,dout,nghost,varargin)
% [] = precompute_ghosts_uvwp_ucf(hucf,dout,nghost,varargin)
% Constructs a given number of ghost cells per processor and saves them to
% mat-files.
% The ghost cell generation is performed with a single read operation (input sweep),
% which extracts the points needed by the neighbors, and a subsequent processing
% sweep, which exchanges the ghost cell data. The data is hold in memory using
% a sparse data structure (c.f. https://de.mathworks.com/matlabcentral/fileexchange/29832-n-dimensional-sparse-arrays)
% and is also written to file in that way, so the 'ndSparse' class is a prerequisite.
% This way, the ghost cells can be added to a chunk of data with a simple add
% operation while avoiding large consumption of disk space.
% Input
% hucf UCF handle (ustar,ucfmulti)
% dout output directory
% nghost number of ghost cells to compute
% periodic: take data from neighboring processor
% non-periodic: duplicate last datapoint
% ? verbosity verbose output? (includes estimation of the memory footprint) [default: 0]
% ? filebase basename of the output files [default: ghost<nghost>]
% ? iseq sequence number of the output files [default: 0]
% Parse input
par = inputParser;
addParamValue(par,'verbosity',0,@isnumeric);
addParamValue(par,'filebase',sprintf('ghost%d',nghost),@ischar);
addParamValue(par,'iseq',0,@isnumeric);
parse(par,varargin{:});
flag_verb = par.Results.verbosity;
filebase = par.Results.filebase;
iseq = par.Results.iseq;
% Read info from tar-archive
[xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp] = read_grid_ucf(hucf);
[ibegu,iendu,jbegu,jendu,kbegu,kendu,...
ibegv,iendv,jbegv,jendv,kbegv,kendv,...
ibegw,iendw,jbegw,jendw,kbegw,kendw,...
ibegp,iendp,jbegp,jendp,kbegp,kendp] = read_procgrid_ucf(hucf);
[params] = read_parameters_ucf(hucf);
nxprocs = params.parallel.nxprocs;
nyprocs = params.parallel.nyprocs;
nzprocs = params.parallel.nzprocs;
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;
% Buffer for boundary data of each chunk
ubd = cell(nxprocs,nyprocs,nzprocs);
vbd = cell(nxprocs,nyprocs,nzprocs);
wbd = cell(nxprocs,nyprocs,nzprocs);
pbd = cell(nxprocs,nyprocs,nzprocs);
% Estimate memory requirement (conservative simplification)
if flag_verb
nxl_estm = floor(numel(xp)/nxprocs);
nyl_estm = floor(numel(yp)/nyprocs);
nzl_estm = floor(numel(zp)/nzprocs);
imem = 0;
for ii=0:nghost-1
imem = imem + 4*8*(...
2*(nxl_estm-ii)*(nyl_estm-ii) + ...
2*(nyl_estm-ii)*(nzl_estm-ii) + ...
2*(nxl_estm-ii)*(nzl_estm-ii) * ...
(nxprocs*nyprocs*nzprocs));
end
fprintf('Estimated memory requirement: %7.0f MiB\n',imem/(1024*1024));
end
% Read boundary data of each chunk into memory
for ixproc=0:nxprocs-1
for iyproc=0:nyprocs-1
for izproc=0:nzprocs-1
ichunk = ixproc*nyprocs*nzprocs+iyproc*nzprocs+izproc;
if flag_verb
fprintf('Processing original chunk %5d\n',ichunk);
end
[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_ucf(hucf,'rank',ichunk,'ghost',0);
% Set interior to zero
u(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxul-2*nghost,nyul-2*nghost,nzul-2*nghost);
v(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxvl-2*nghost,nyvl-2*nghost,nzvl-2*nghost);
w(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxwl-2*nghost,nywl-2*nghost,nzwl-2*nghost);
p(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxpl-2*nghost,nypl-2*nghost,nzpl-2*nghost);
% Save ghost points as sparse matrix
ubd{ixproc+1,iyproc+1,izproc+1} = ndSparse(u);
vbd{ixproc+1,iyproc+1,izproc+1} = ndSparse(v);
wbd{ixproc+1,iyproc+1,izproc+1} = ndSparse(w);
pbd{ixproc+1,iyproc+1,izproc+1} = ndSparse(p);
end
end
end
% Construct ghost cells for each chunk
for ixproc=0:nxprocs-1
nxul = iendu(ixproc+1)-ibegu(ixproc+1)+1;
nxvl = iendv(ixproc+1)-ibegv(ixproc+1)+1;
nxwl = iendw(ixproc+1)-ibegw(ixproc+1)+1;
nxpl = iendp(ixproc+1)-ibegp(ixproc+1)+1;
for iyproc=0:nyprocs-1
nyul = jendu(iyproc+1)-jbegu(iyproc+1)+1;
nyvl = jendv(iyproc+1)-jbegv(iyproc+1)+1;
nywl = jendw(iyproc+1)-jbegw(iyproc+1)+1;
nypl = jendp(iyproc+1)-jbegp(iyproc+1)+1;
for izproc=0:nzprocs-1
nzul = kendu(izproc+1)-kbegu(izproc+1)+1;
nzvl = kendv(izproc+1)-kbegv(izproc+1)+1;
nzwl = kendw(izproc+1)-kbegw(izproc+1)+1;
nzpl = kendp(izproc+1)-kbegp(izproc+1)+1;
ichunk = ixproc*nyprocs*nzprocs+iyproc*nzprocs+izproc;
if flag_verb
fprintf('Constructing ghost cells for chunk %5d\n',ichunk);
end
% Create sparse arrays which hold the results
ugh = ones(nxul+2*nghost,nyul+2*nghost,nzul+2*nghost);
ugh(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxul,nyul,nzul);
ugh = ndSparse(ugh);
vgh = ones(nxvl+2*nghost,nyvl+2*nghost,nzvl+2*nghost);
vgh(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxvl,nyvl,nzvl);
vgh = ndSparse(vgh);
wgh = ones(nxwl+2*nghost,nywl+2*nghost,nzwl+2*nghost);
wgh(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxwl,nywl,nzwl);
wgh = ndSparse(wgh);
pgh = ones(nxpl+2*nghost,nypl+2*nghost,nzpl+2*nghost);
pgh(1+nghost:end-nghost,1+nghost:end-nghost,1+nghost:end-nghost) = zeros(nxpl,nypl,nzpl);
pgh = ndSparse(pgh);
% Fill them with ghost cell data
% u
% 1st faces (6)
ugh(1:nghost, 1+nghost:end-nghost,1+nghost:end-nghost) = ubd{mod(ixproc-1,nxprocs)+1,iyproc+1,izproc+1}(end-nghost+1:end,:, : );
ugh(end-nghost+1:end, 1+nghost:end-nghost,1+nghost:end-nghost) = ubd{mod(ixproc+1,nxprocs)+1,iyproc+1,izproc+1}(1:nghost, :, : );
ugh(1+nghost:end-nghost,1:nghost, 1+nghost:end-nghost) = ubd{ixproc+1,mod(iyproc-1,nyprocs)+1,izproc+1}(:, end-nghost+1:end,: );
ugh(1+nghost:end-nghost,end-nghost+1:end, 1+nghost:end-nghost) = ubd{ixproc+1,mod(iyproc+1,nyprocs)+1,izproc+1}(:, 1:nghost, : );
ugh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost ) = ubd{ixproc+1,iyproc+1,mod(izproc-1,nzprocs)+1}(:, :, end-nghost+1:end);
ugh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end ) = ubd{ixproc+1,iyproc+1,mod(izproc+1,nzprocs)+1}(:, :, 1:nghost );
% 2nd edges (12)
ugh(1:nghost, 1:nghost, 1+nghost:end-nghost) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1 }(end-nghost+1:end,end-nghost+1:end,:);
ugh(1:nghost, 1+nghost:end-nghost,1:nghost ) = ubd{mod(ixproc-1,nxprocs)+1,iyproc+1, mod(izproc-1,nzprocs)+1}(end-nghost+1:end,:, end-nghost+1:end);
ugh(1+nghost:end-nghost,1:nghost, 1:nghost ) = ubd{ixproc+1, mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:, end-nghost+1:end,end-nghost+1:end);
ugh(end-nghost+1:end, end-nghost+1:end, 1+nghost:end-nghost) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1 }(1:nghost, 1:nghost, : );
ugh(end-nghost+1:end, 1+nghost:end-nghost,end-nghost+1:end ) = ubd{mod(ixproc+1,nxprocs)+1,iyproc+1, mod(izproc+1,nzprocs)+1}(1:nghost, :, 1:nghost );
ugh(1+nghost:end-nghost,end-nghost+1:end, end-nghost+1:end ) = ubd{ixproc+1, mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:, 1:nghost, 1:nghost );
ugh(1:nghost, end-nghost+1:end, 1+nghost:end-nghost) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1 }(end-nghost+1:end,1:nghost, : );
ugh(1:nghost, 1+nghost:end-nghost,end-nghost+1:end ) = ubd{mod(ixproc-1,nxprocs)+1,iyproc+1, mod(izproc+1,nzprocs)+1}(end-nghost+1:end,:, 1:nghost );
ugh(1+nghost:end-nghost,1:nghost, end-nghost+1:end ) = ubd{ixproc+1, mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:, end-nghost+1:end,1:nghost );
ugh(end-nghost+1:end, 1:nghost, 1+nghost:end-nghost) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1 }(1:nghost, end-nghost+1:end,: );
ugh(end-nghost+1:end, 1+nghost:end-nghost,1:nghost ) = ubd{mod(ixproc+1,nxprocs)+1,iyproc+1, mod(izproc-1,nzprocs)+1}(1:nghost, :, end-nghost+1:end);
ugh(1+nghost:end-nghost,end-nghost+1:end, 1:nghost ) = ubd{ixproc+1, mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:, 1:nghost, end-nghost+1:end);
% 3rd corners (8)
ugh(1:nghost, 1:nghost, 1:nghost ) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end);
ugh(end-nghost+1:end,1:nghost, 1:nghost ) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost, end-nghost+1:end,end-nghost+1:end);
ugh(1:nghost, end-nghost+1:end,1:nghost ) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,1:nghost, end-nghost+1:end);
ugh(1:nghost, 1:nghost, end-nghost+1:end) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,1:nghost );
ugh(end-nghost+1:end,end-nghost+1:end,1:nghost ) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost, 1:nghost, end-nghost+1:end);
ugh(end-nghost+1:end,1:nghost, end-nghost+1:end) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost, end-nghost+1:end,1:nghost );
ugh(1:nghost, end-nghost+1:end,end-nghost+1:end) = ubd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,1:nghost, 1:nghost );
ugh(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end) = ubd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost, 1:nghost, 1:nghost );
% Correct if not periodic
if ~xperiodic && ixproc==0
ugh(1:nghost, :, : ) = repmat(ugh(nghost+1,:,:),[nghost,1,1]);
ugh(1:nghost, 1+nghost:end-nghost,1+nghost:end-nghost) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(1,:,:),[nghost,1,1]);
end
if ~xperiodic && ixproc==nxprocs-1
ugh(end-nghost+1:end,:, : ) = repmat(ugh(end-nghost,:,:),[nghost,1,1]);
ugh(end-nghost+1:end,1+nghost:end-nghost,1+nghost:end-nghost) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(end,:,:),[nghost,1,1]);
end
if ~yperiodic && iyproc==0
ugh(:, 1:nghost, : ) = repmat(ugh(:,nghost+1,:),[1,nghost,1]);
ugh(1+nghost:end-nghost,1:nghost, 1+nghost:end-nghost) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(:,1,:),[1,nghost,1]);
end
if ~yperiodic && iyproc==nyprocs-1
ugh(:, end-nghost+1:end,: ) = repmat(ugh(:,end-nghost,:),[1,nghost,1]);
ugh(1+nghost:end-nghost,end-nghost+1:end,1+nghost:end-nghost) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(:,end,:),[1,nghost,1]);
end
if ~zperiodic && izproc==0
ugh(:, :, 1:nghost ) = repmat(ugh(:,:,nghost+1),[1,1,nghost]);
ugh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost ) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(:,:,1),[1,1,nghost]);
end
if ~zperiodic && izproc==nzprocs-1
ugh(:, :, end-nghost+1:end) = repmat(ugh(:,:,end-nghost),[1,1,nghost]);
ugh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end) = repmat(ubd{ixproc+1,iyproc+1,izproc+1}(:,:,end),[1,1,nghost]);
end
% Fill them with ghost cell data
% v
% 1st faces (6)
vgh(1:nghost, 1+nghost:end-nghost,1+nghost:end-nghost) = vbd{mod(ixproc-1,nxprocs)+1,iyproc+1,izproc+1}(end-nghost+1:end,:, : );
vgh(end-nghost+1:end, 1+nghost:end-nghost,1+nghost:end-nghost) = vbd{mod(ixproc+1,nxprocs)+1,iyproc+1,izproc+1}(1:nghost, :, : );
vgh(1+nghost:end-nghost,1:nghost, 1+nghost:end-nghost) = vbd{ixproc+1,mod(iyproc-1,nyprocs)+1,izproc+1}(:, end-nghost+1:end,: );
vgh(1+nghost:end-nghost,end-nghost+1:end, 1+nghost:end-nghost) = vbd{ixproc+1,mod(iyproc+1,nyprocs)+1,izproc+1}(:, 1:nghost, : );
vgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost ) = vbd{ixproc+1,iyproc+1,mod(izproc-1,nzprocs)+1}(:, :, end-nghost+1:end);
vgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end ) = vbd{ixproc+1,iyproc+1,mod(izproc+1,nzprocs)+1}(:, :, 1:nghost );
% 2nd edges (12)
vgh(1:nghost, 1:nghost, 1+nghost:end-nghost) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1 }(end-nghost+1:end,end-nghost+1:end,:);
vgh(1:nghost, 1+nghost:end-nghost,1:nghost ) = vbd{mod(ixproc-1,nxprocs)+1,iyproc+1, mod(izproc-1,nzprocs)+1}(end-nghost+1:end,:, end-nghost+1:end);
vgh(1+nghost:end-nghost,1:nghost, 1:nghost ) = vbd{ixproc+1, mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:, end-nghost+1:end,end-nghost+1:end);
vgh(end-nghost+1:end, end-nghost+1:end, 1+nghost:end-nghost) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1 }(1:nghost, 1:nghost, : );
vgh(end-nghost+1:end, 1+nghost:end-nghost,end-nghost+1:end ) = vbd{mod(ixproc+1,nxprocs)+1,iyproc+1, mod(izproc+1,nzprocs)+1}(1:nghost, :, 1:nghost );
vgh(1+nghost:end-nghost,end-nghost+1:end, end-nghost+1:end ) = vbd{ixproc+1, mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:, 1:nghost, 1:nghost );
vgh(1:nghost, end-nghost+1:end, 1+nghost:end-nghost) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1 }(end-nghost+1:end,1:nghost, : );
vgh(1:nghost, 1+nghost:end-nghost,end-nghost+1:end ) = vbd{mod(ixproc-1,nxprocs)+1,iyproc+1, mod(izproc+1,nzprocs)+1}(end-nghost+1:end,:, 1:nghost );
vgh(1+nghost:end-nghost,1:nghost, end-nghost+1:end ) = vbd{ixproc+1, mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:, end-nghost+1:end,1:nghost );
vgh(end-nghost+1:end, 1:nghost, 1+nghost:end-nghost) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1 }(1:nghost, end-nghost+1:end,: );
vgh(end-nghost+1:end, 1+nghost:end-nghost,1:nghost ) = vbd{mod(ixproc+1,nxprocs)+1,iyproc+1, mod(izproc-1,nzprocs)+1}(1:nghost, :, end-nghost+1:end);
vgh(1+nghost:end-nghost,end-nghost+1:end, 1:nghost ) = vbd{ixproc+1, mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:, 1:nghost, end-nghost+1:end);
% 3rd corners (8)
vgh(1:nghost, 1:nghost, 1:nghost ) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end);
vgh(end-nghost+1:end,1:nghost, 1:nghost ) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost, end-nghost+1:end,end-nghost+1:end);
vgh(1:nghost, end-nghost+1:end,1:nghost ) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,1:nghost, end-nghost+1:end);
vgh(1:nghost, 1:nghost, end-nghost+1:end) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,1:nghost );
vgh(end-nghost+1:end,end-nghost+1:end,1:nghost ) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost, 1:nghost, end-nghost+1:end);
vgh(end-nghost+1:end,1:nghost, end-nghost+1:end) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost, end-nghost+1:end,1:nghost );
vgh(1:nghost, end-nghost+1:end,end-nghost+1:end) = vbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,1:nghost, 1:nghost );
vgh(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end) = vbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost, 1:nghost, 1:nghost );
% Correct if not periodic
if ~xperiodic && ixproc==0
vgh(1:nghost, :, : ) = repmat(vgh(nghost+1,:,:),[nghost,1,1]);
vgh(1:nghost, 1+nghost:end-nghost,1+nghost:end-nghost) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(1,:,:),[nghost,1,1]);
end
if ~xperiodic && ixproc==nxprocs-1
vgh(end-nghost+1:end,:, : ) = repmat(vgh(end-nghost,:,:),[nghost,1,1]);
vgh(end-nghost+1:end,1+nghost:end-nghost,1+nghost:end-nghost) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(end,:,:),[nghost,1,1]);
end
if ~yperiodic && iyproc==0
vgh(:, 1:nghost, : ) = repmat(vgh(:,nghost+1,:),[1,nghost,1]);
vgh(1+nghost:end-nghost,1:nghost, 1+nghost:end-nghost) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(:,1,:),[1,nghost,1]);
end
if ~yperiodic && iyproc==nyprocs-1
vgh(:, end-nghost+1:end,: ) = repmat(vgh(:,end-nghost,:),[1,nghost,1]);
vgh(1+nghost:end-nghost,end-nghost+1:end,1+nghost:end-nghost) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(:,end,:),[1,nghost,1]);
end
if ~zperiodic && izproc==0
vgh(:, :, 1:nghost ) = repmat(vgh(:,:,nghost+1),[1,1,nghost]);
vgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost ) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(:,:,1),[1,1,nghost]);
end
if ~zperiodic && izproc==nzprocs-1
vgh(:, :, end-nghost+1:end) = repmat(vgh(:,:,end-nghost),[1,1,nghost]);
vgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end) = repmat(vbd{ixproc+1,iyproc+1,izproc+1}(:,:,end),[1,1,nghost]);
end
% Fill them with ghost cell data
% w
% 1st faces (6)
wgh(1:nghost, 1+nghost:end-nghost,1+nghost:end-nghost) = wbd{mod(ixproc-1,nxprocs)+1,iyproc+1,izproc+1}(end-nghost+1:end,:, : );
wgh(end-nghost+1:end, 1+nghost:end-nghost,1+nghost:end-nghost) = wbd{mod(ixproc+1,nxprocs)+1,iyproc+1,izproc+1}(1:nghost, :, : );
wgh(1+nghost:end-nghost,1:nghost, 1+nghost:end-nghost) = wbd{ixproc+1,mod(iyproc-1,nyprocs)+1,izproc+1}(:, end-nghost+1:end,: );
wgh(1+nghost:end-nghost,end-nghost+1:end, 1+nghost:end-nghost) = wbd{ixproc+1,mod(iyproc+1,nyprocs)+1,izproc+1}(:, 1:nghost, : );
wgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost ) = wbd{ixproc+1,iyproc+1,mod(izproc-1,nzprocs)+1}(:, :, end-nghost+1:end);
wgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end ) = wbd{ixproc+1,iyproc+1,mod(izproc+1,nzprocs)+1}(:, :, 1:nghost );
% 2nd edges (12)
wgh(1:nghost, 1:nghost, 1+nghost:end-nghost) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1 }(end-nghost+1:end,end-nghost+1:end,:);
wgh(1:nghost, 1+nghost:end-nghost,1:nghost ) = wbd{mod(ixproc-1,nxprocs)+1,iyproc+1, mod(izproc-1,nzprocs)+1}(end-nghost+1:end,:, end-nghost+1:end);
wgh(1+nghost:end-nghost,1:nghost, 1:nghost ) = wbd{ixproc+1, mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:, end-nghost+1:end,end-nghost+1:end);
wgh(end-nghost+1:end, end-nghost+1:end, 1+nghost:end-nghost) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1 }(1:nghost, 1:nghost, : );
wgh(end-nghost+1:end, 1+nghost:end-nghost,end-nghost+1:end ) = wbd{mod(ixproc+1,nxprocs)+1,iyproc+1, mod(izproc+1,nzprocs)+1}(1:nghost, :, 1:nghost );
wgh(1+nghost:end-nghost,end-nghost+1:end, end-nghost+1:end ) = wbd{ixproc+1, mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:, 1:nghost, 1:nghost );
wgh(1:nghost, end-nghost+1:end, 1+nghost:end-nghost) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1 }(end-nghost+1:end,1:nghost, : );
wgh(1:nghost, 1+nghost:end-nghost,end-nghost+1:end ) = wbd{mod(ixproc-1,nxprocs)+1,iyproc+1, mod(izproc+1,nzprocs)+1}(end-nghost+1:end,:, 1:nghost );
wgh(1+nghost:end-nghost,1:nghost, end-nghost+1:end ) = wbd{ixproc+1, mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:, end-nghost+1:end,1:nghost );
wgh(end-nghost+1:end, 1:nghost, 1+nghost:end-nghost) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1 }(1:nghost, end-nghost+1:end,: );
wgh(end-nghost+1:end, 1+nghost:end-nghost,1:nghost ) = wbd{mod(ixproc+1,nxprocs)+1,iyproc+1, mod(izproc-1,nzprocs)+1}(1:nghost, :, end-nghost+1:end);
wgh(1+nghost:end-nghost,end-nghost+1:end, 1:nghost ) = wbd{ixproc+1, mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:, 1:nghost, end-nghost+1:end);
% 3rd corners (8)
wgh(1:nghost, 1:nghost, 1:nghost ) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end);
wgh(end-nghost+1:end,1:nghost, 1:nghost ) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost, end-nghost+1:end,end-nghost+1:end);
wgh(1:nghost, end-nghost+1:end,1:nghost ) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,1:nghost, end-nghost+1:end);
wgh(1:nghost, 1:nghost, end-nghost+1:end) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,1:nghost );
wgh(end-nghost+1:end,end-nghost+1:end,1:nghost ) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost, 1:nghost, end-nghost+1:end);
wgh(end-nghost+1:end,1:nghost, end-nghost+1:end) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost, end-nghost+1:end,1:nghost );
wgh(1:nghost, end-nghost+1:end,end-nghost+1:end) = wbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,1:nghost, 1:nghost );
wgh(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end) = wbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost, 1:nghost, 1:nghost );
% Correct if not periodic
if ~xperiodic && ixproc==0
wgh(1:nghost, :, : ) = repmat(wgh(nghost+1,:,:),[nghost,1,1]);
wgh(1:nghost, 1+nghost:end-nghost,1+nghost:end-nghost) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(1,:,:),[nghost,1,1]);
end
if ~xperiodic && ixproc==nxprocs-1
wgh(end-nghost+1:end,:, : ) = repmat(wgh(end-nghost,:,:),[nghost,1,1]);
wgh(end-nghost+1:end,1+nghost:end-nghost,1+nghost:end-nghost) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(end,:,:),[nghost,1,1]);
end
if ~yperiodic && iyproc==0
wgh(:, 1:nghost, : ) = repmat(wgh(:,nghost+1,:),[1,nghost,1]);
wgh(1+nghost:end-nghost,1:nghost, 1+nghost:end-nghost) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(:,1,:),[1,nghost,1]);
end
if ~yperiodic && iyproc==nyprocs-1
wgh(:, end-nghost+1:end,: ) = repmat(wgh(:,end-nghost,:),[1,nghost,1]);
wgh(1+nghost:end-nghost,end-nghost+1:end,1+nghost:end-nghost) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(:,end,:),[1,nghost,1]);
end
if ~zperiodic && izproc==0
wgh(:, :, 1:nghost ) = repmat(wgh(:,:,nghost+1),[1,1,nghost]);
wgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost ) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(:,:,1),[1,1,nghost]);
end
if ~zperiodic && izproc==nzprocs-1
wgh(:, :, end-nghost+1:end) = repmat(wgh(:,:,end-nghost),[1,1,nghost]);
wgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end) = repmat(wbd{ixproc+1,iyproc+1,izproc+1}(:,:,end),[1,1,nghost]);
end
% Fill them with ghost cell data
% p
% 1st faces (6)
pgh(1:nghost, 1+nghost:end-nghost,1+nghost:end-nghost) = pbd{mod(ixproc-1,nxprocs)+1,iyproc+1,izproc+1}(end-nghost+1:end,:, : );
pgh(end-nghost+1:end, 1+nghost:end-nghost,1+nghost:end-nghost) = pbd{mod(ixproc+1,nxprocs)+1,iyproc+1,izproc+1}(1:nghost, :, : );
pgh(1+nghost:end-nghost,1:nghost, 1+nghost:end-nghost) = pbd{ixproc+1,mod(iyproc-1,nyprocs)+1,izproc+1}(:, end-nghost+1:end,: );
pgh(1+nghost:end-nghost,end-nghost+1:end, 1+nghost:end-nghost) = pbd{ixproc+1,mod(iyproc+1,nyprocs)+1,izproc+1}(:, 1:nghost, : );
pgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost ) = pbd{ixproc+1,iyproc+1,mod(izproc-1,nzprocs)+1}(:, :, end-nghost+1:end);
pgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end ) = pbd{ixproc+1,iyproc+1,mod(izproc+1,nzprocs)+1}(:, :, 1:nghost );
% 2nd edges (12)
pgh(1:nghost, 1:nghost, 1+nghost:end-nghost) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1 }(end-nghost+1:end,end-nghost+1:end,:);
pgh(1:nghost, 1+nghost:end-nghost,1:nghost ) = pbd{mod(ixproc-1,nxprocs)+1,iyproc+1, mod(izproc-1,nzprocs)+1}(end-nghost+1:end,:, end-nghost+1:end);
pgh(1+nghost:end-nghost,1:nghost, 1:nghost ) = pbd{ixproc+1, mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:, end-nghost+1:end,end-nghost+1:end);
pgh(end-nghost+1:end, end-nghost+1:end, 1+nghost:end-nghost) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1 }(1:nghost, 1:nghost, : );
pgh(end-nghost+1:end, 1+nghost:end-nghost,end-nghost+1:end ) = pbd{mod(ixproc+1,nxprocs)+1,iyproc+1, mod(izproc+1,nzprocs)+1}(1:nghost, :, 1:nghost );
pgh(1+nghost:end-nghost,end-nghost+1:end, end-nghost+1:end ) = pbd{ixproc+1, mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:, 1:nghost, 1:nghost );
pgh(1:nghost, end-nghost+1:end, 1+nghost:end-nghost) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,izproc+1 }(end-nghost+1:end,1:nghost, : );
pgh(1:nghost, 1+nghost:end-nghost,end-nghost+1:end ) = pbd{mod(ixproc-1,nxprocs)+1,iyproc+1, mod(izproc+1,nzprocs)+1}(end-nghost+1:end,:, 1:nghost );
pgh(1+nghost:end-nghost,1:nghost, end-nghost+1:end ) = pbd{ixproc+1, mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(:, end-nghost+1:end,1:nghost );
pgh(end-nghost+1:end, 1:nghost, 1+nghost:end-nghost) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,izproc+1 }(1:nghost, end-nghost+1:end,: );
pgh(end-nghost+1:end, 1+nghost:end-nghost,1:nghost ) = pbd{mod(ixproc+1,nxprocs)+1,iyproc+1, mod(izproc-1,nzprocs)+1}(1:nghost, :, end-nghost+1:end);
pgh(1+nghost:end-nghost,end-nghost+1:end, 1:nghost ) = pbd{ixproc+1, mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(:, 1:nghost, end-nghost+1:end);
% 3rd corners (8)
pgh(1:nghost, 1:nghost, 1:nghost ) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end);
pgh(end-nghost+1:end,1:nghost, 1:nghost ) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost, end-nghost+1:end,end-nghost+1:end);
pgh(1:nghost, end-nghost+1:end,1:nghost ) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(end-nghost+1:end,1:nghost, end-nghost+1:end);
pgh(1:nghost, 1:nghost, end-nghost+1:end) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,end-nghost+1:end,1:nghost );
pgh(end-nghost+1:end,end-nghost+1:end,1:nghost ) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc-1,nzprocs)+1}(1:nghost, 1:nghost, end-nghost+1:end);
pgh(end-nghost+1:end,1:nghost, end-nghost+1:end) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc-1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost, end-nghost+1:end,1:nghost );
pgh(1:nghost, end-nghost+1:end,end-nghost+1:end) = pbd{mod(ixproc-1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(end-nghost+1:end,1:nghost, 1:nghost );
pgh(end-nghost+1:end,end-nghost+1:end,end-nghost+1:end) = pbd{mod(ixproc+1,nxprocs)+1,mod(iyproc+1,nyprocs)+1,mod(izproc+1,nzprocs)+1}(1:nghost, 1:nghost, 1:nghost );
% Correct if not periodic
if ~xperiodic && ixproc==0
pgh(1:nghost, :, : ) = repmat(pgh(nghost+1,:,:),[nghost,1,1]);
pgh(1:nghost, 1+nghost:end-nghost,1+nghost:end-nghost) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(1,:,:),[nghost,1,1]);
end
if ~xperiodic && ixproc==nxprocs-1
pgh(end-nghost+1:end,:, : ) = repmat(pgh(end-nghost,:,:),[nghost,1,1]);
pgh(end-nghost+1:end,1+nghost:end-nghost,1+nghost:end-nghost) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(end,:,:),[nghost,1,1]);
end
if ~yperiodic && iyproc==0
pgh(:, 1:nghost, : ) = repmat(pgh(:,nghost+1,:),[1,nghost,1]);
pgh(1+nghost:end-nghost,1:nghost, 1+nghost:end-nghost) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(:,1,:),[1,nghost,1]);
end
if ~yperiodic && iyproc==nyprocs-1
pgh(:, end-nghost+1:end,: ) = repmat(pgh(:,end-nghost,:),[1,nghost,1]);
pgh(1+nghost:end-nghost,end-nghost+1:end,1+nghost:end-nghost) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(:,end,:),[1,nghost,1]);
end
if ~zperiodic && izproc==0
pgh(:, :, 1:nghost ) = repmat(pgh(:,:,nghost+1),[1,1,nghost]);
pgh(1+nghost:end-nghost,1+nghost:end-nghost,1:nghost ) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(:,:,1),[1,1,nghost]);
end
if ~zperiodic && izproc==nzprocs-1
pgh(:, :, end-nghost+1:end) = repmat(pgh(:,:,end-nghost),[1,1,nghost]);
pgh(1+nghost:end-nghost,1+nghost:end-nghost,end-nghost+1:end) = repmat(pbd{ixproc+1,iyproc+1,izproc+1}(:,:,end),[1,1,nghost]);
end
fout = sprintf('%s/%s_%04d.%05d.mat',dout,filebase,iseq,ichunk);
save(fout,'ugh','vgh','wgh','pgh');
end
end
end
end