diff --git a/matlab/precompute_ghosts_uvwp_ucf.m b/matlab/precompute_ghosts_uvwp_ucf.m new file mode 100644 index 0000000..5115370 --- /dev/null +++ b/matlab/precompute_ghosts_uvwp_ucf.m @@ -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] + % ? 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 \ No newline at end of file