update from taifun

matlab/read_grid_ucf.m

@@ -1,6 +1,5 @@
-function [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,...
+function [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,xs,ys,zs] = read_grid_ucf(file,varargin)
-  nxu,nyu,nzu,nxv,nyv,nzv,nxw,nyw,nzw,nxp,nyp,nzp] = read_grid_ucf(file,varargin)
+    % [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,xs,ys,zs] = read_grid_ucf(file,varargin)
-    % [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp] = read_grid_ucf(file,varargin)
     % Reads staggered grid from ucf file.
     % Input
     %   file              file to be read
@@ -11,7 +10,8 @@ function [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,...
     %   xu,yu,zu          velocity grid u
     %   xv,yv,zv          velocity grid v
     %   xw,yw,zw          velocity grid w
-    %   xp,yp,zp          pressure grid / scalar grid
+    %   xp,yp,zp          pressure grid
+    %   xs,ys,zs          scalar grid
 
     % Parse optional input arguments
     par = inputParser;
@@ -20,10 +20,6 @@ function [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,...
     addParamValue(par,'tarmode',0,@isnumeric);  % deprecated
     parse(par,varargin{:});
     
-    % Define sets to be read
-    sets = {'u','v','w','p'};
-    nset = numel(sets);
-    
     % Open UCF file and read
     obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
     switch class(file)
@@ -36,6 +32,15 @@ function [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,...
       error('Input file type not supported: %s',class(file));
     end
     
+    % Define sets to be read
+    isDualMesh=(obj.UCFVersion>=2);
+    if isDualMesh
+      sets = {'u','v','w','p','s'};
+    else
+      sets = {'u','v','w','p'};
+    end
+    nset = numel(sets);
+  
     % Read raw data
     x = cell(1,nset);
     y = cell(1,nset);
@@ -59,31 +64,30 @@ function [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp,...
     xu = x{iset};
     yu = y{iset};
     zu = z{iset};
-    nxu = nx{iset};
-    nyu = ny{iset};
-    nzu = nz{iset};
     
     iset = find(strcmp(sets,'v'));
     xv = x{iset};
     yv = y{iset};
     zv = z{iset};
-    nxv = nx{iset};
-    nyv = ny{iset};
-    nzv = nz{iset};
     
     iset = find(strcmp(sets,'w'));
     xw = x{iset};
     yw = y{iset};
     zw = z{iset};
-    nxw = nx{iset};
-    nyw = ny{iset};
-    nzw = nz{iset};
     
     iset = find(strcmp(sets,'p'));
     xp = x{iset};
     yp = y{iset};
     zp = z{iset};
-    nxp = nx{iset};
+    
-    nyp = ny{iset};
+    if isDualMesh
-    nzp = nz{iset};
+      iset = find(strcmp(sets,'s'));
+      xs = x{iset};
+      ys = y{iset};
+      zs = z{iset};
+    else
+      xs = xp;
+      ys = yp;
+      zs = zp;
+    end
 end

matlab/read_particles_ucf.m

@@ -12,6 +12,7 @@ function [pp,col,stime] = read_particles_ucf(file,varargin)
     % ? verbosity         verbose output? (default: no)
     % ? debug             debug output? (default: no)
     % ? tarmode           read from tar-file? (default: 0)
+    % ? append            read an append file? Only matters for tar archives (default: 0)
     % Output
     %   pp                particle data in specified format
     %   col               container map 'char'->'double' which maps quantity names to column index
@@ -24,9 +25,11 @@ function [pp,col,stime] = read_particles_ucf(file,varargin)
     addParamValue(par,'verbosity',0,@isnumeric);
     addParamValue(par,'debug',0,@isnumeric);
     addParamValue(par,'tarmode',0,@isnumeric);  % deprecated (automatically checked now)
+    addParamValue(par,'append',0,@isnumeric);
     parse(par,varargin{:});
     istep     = par.Results.step;
     mlformat  = par.Results.format;
+    isappend  = par.Results.append;
     
     % Open file
     obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
@@ -34,7 +37,11 @@ function [pp,col,stime] = read_particles_ucf(file,varargin)
     case 'char'
       obj.open(file);
     case {'ustar','ucfmulti'}
+      if isappend
+        ptr = file.pointer('append_particles.bin');
+      else
         ptr = file.pointer('particles.bin');
+      end
       obj.opentar(ptr);
     otherwise
       error('Input file type not supported: %s',class(file));

matlab/read_procgrid_ucf.m

@@ -2,12 +2,12 @@ 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_ucf(file,varargin)
+          ibegs,iends,jbegs,jends,kbegs,kends] = read_procgrid_ucf(file,varargin)
     % [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_ucf(file,varargin)
+    %  ibegs,iends,jbegs,jends,kbegs,kends] = read_procgrid_ucf(file,varargin)
     % Reads processor grids.
     % Input
     %   file              file name (if tar-mode: ustar handle)
@@ -18,7 +18,7 @@ function [ibegu,iendu,jbegu,jendu,kbegu,kendu,...
     %  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
+    %  ibegs,iends,jbegs,jends,kbegs,kends    processor grid s  
     
     % Parse optional input arguments
     par = inputParser;
@@ -27,10 +27,6 @@ function [ibegu,iendu,jbegu,jendu,kbegu,kendu,...
     addParamValue(par,'tarmode',0,@isnumeric);
     parse(par,varargin{:});
     
-    % Define sets to be read
-    sets = {'u','v','w','p'};
-    nset = numel(sets);
-    
     % Open file
     obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
     switch class(file)
@@ -43,6 +39,15 @@ function [ibegu,iendu,jbegu,jendu,kbegu,kendu,...
       error('Input file type not supported: %s',class(file));
     end
     
+    % Define sets to be read
+    isDualMesh=(obj.UCFVersion>=2);
+    if isDualMesh
+      sets = {'u','v','w','p','s'};
+    else
+      sets = {'u','v','w','p'};
+    end
+    nset = numel(sets);
+    
     % Read raw data
     ibeg = cell(1,nset);
     iend = cell(1,nset);
@@ -99,4 +104,21 @@ function [ibegu,iendu,jbegu,jendu,kbegu,kendu,...
     jendp = cast(jend{iset},'double');
     kbegp = cast(kbeg{iset},'double');
     kendp = cast(kend{iset},'double');
+    
+    if isDualMesh
+      iset = find(strcmp(sets,'s'));
+      ibegs = cast(ibeg{iset},'double');
+      iends = cast(iend{iset},'double');
+      jbegs = cast(jbeg{iset},'double');
+      jends = cast(jend{iset},'double');
+      kbegs = cast(kbeg{iset},'double');
+      kends = cast(kend{iset},'double');
+    else
+      ibegs = ibegp;
+      iends = iendp;
+      jbegs = jbegp;
+      jends = jendp;
+      kbegs = kbegp;
+      kends = kendp;
+    end
 end

matlab/read_scal_complete_ucf.m

@@ -62,7 +62,7 @@ function [s] = read_scal_complete_ucf(file,varargin)
     % Reassign content from loop and create global arrays
     s = zeros(nx,ny,nz,nscal,class(q{1}));
     for iscal=1:nscal
-      s(ib:ib+nxl-1,jb:jb+nyl-1,kb:kb+nzl-1,iscal) = q{iset};
+      s(ib:ib+nxl-1,jb:jb+nyl-1,kb:kb+nzl-1,iscal) = q{iscal};
     end
 
     % Close the first file
@@ -113,7 +113,7 @@ function [s] = read_scal_complete_ucf(file,varargin)
 
       % Reassign content from loop and create global arrays
       for iscal=1:nscal
-        s(ib:ib+nxl-1,jb:jb+nyl-1,kb:kb+nzl-1,iscal) = q{iset};
+        s(ib:ib+nxl-1,jb:jb+nyl-1,kb:kb+nzl-1,iscal) = q{iscal};
       end
 
       % Close file
@@ -123,9 +123,9 @@ function [s] = read_scal_complete_ucf(file,varargin)
       ifile = ifile+1;
       switch class(file)
       case 'char'
-        fname = sprintf('scal.%05d',ifile);
-      case {'ustar','ucfmulti'}
         fname = sprintf('%s/%s.%05d',fdir,fbase,ifile);
+      case {'ustar','ucfmulti'}
+        fname = sprintf('scal.%05d',ifile);
       end
     end
 end

matlab/read_statistics_channel_ucf.m

@@ -1,10 +1,29 @@
 function [tbeg,tend,nstat,um,uu,vv,ww,uv,uub,uvb,varargout] = read_statistics_channel_ucf(file,varargin)
+    % [tbeg,tend,nstat,um,uu,vv,ww,uv,uub,uvb,varargout] = read_statistics_channel_ucf(file,varargin)
+    % Reads statistics file for channel data in UCF format.
+    % Input
+    %   file              file name (if tar-mode: ustar handle)
+    % ? verbosity         verbose output? (default: no)
+    % ? debug             debug output? (default: no)
+    % ? pure              read an statistics_pure file? Only matters for tar archives (default: 0)
+    % Output
+    %   tbeg              time at beginning of sampling
+    %   tend              time at end of sampling
+    %   nstat             number of samples
+    %   um                mean streamwise velocity
+    %   uu,vv,ww          <u'u'>, <v'v'>, <w'w'>
+    %   uv                <u'v'>
+    %   uub               <uu>
+    %   uvb               <uv>
+    % ? fxp,fyp,fzp       particle force in Eulerian frame
 
     % Parse optional input arguments
     par = inputParser;
     addParamValue(par,'verbosity',0,@isnumeric);
     addParamValue(par,'debug',0,@isnumeric);
+    addParamValue(par,'pure',0,@isnumeric);
     parse(par,varargin{:});
+    ispure  = par.Results.pure;
 
     % Define sets to be read
     sets = {'um','uu','vv','ww','uv','uub','uvb'};
@@ -12,7 +31,19 @@ function [tbeg,tend,nstat,um,uu,vv,ww,uv,uub,uvb,varargout] = read_statistics_ch
 
     % Open UCF file and read
     obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
+    switch class(file)
+    case 'char'
       obj.open(file);
+    case {'ustar','ucfmulti'}
+      if ispure
+        ptr = file.pointer('statistics_pure.bin');
+      else
+        ptr = file.pointer('statistics.bin');
+      end
+      obj.opentar(ptr);
+    otherwise
+      error('Input file type not supported: %s',class(file));
+    end
     tend = obj.getSimulationTime();
 
     % Read first dataset

matlab/write_uvwp_complete_ucf.m

@@ -62,35 +62,35 @@ function [] = write_uvwp_complete_ucf(file,...
           ighost = 0;
 
           % Write u
-          nxul = size(u,1)-2*ighost;
-          nyul = size(u,2)-2*ighost;
-          nzul = size(u,3)-2*ighost;
-          params = int64([ighost,ibu,jbu,kbu,nxul,nyul,nzul,nxu,nyu,nzu]);
           data = u(ibu:ieu,jbu:jeu,kbu:keu);
+          nxul = size(data,1)-2*ighost;
+          nyul = size(data,2)-2*ighost;
+          nzul = size(data,3)-2*ighost;
+          params = int64([ighost,ibu,jbu,kbu,nxul,nyul,nzul,nxu,nyu,nzu]);
           obj.appendSet(data,params);
 
           % Write v
-          nxvl = size(v,1)-2*ighost;
-          nyvl = size(v,2)-2*ighost;
-          nzvl = size(v,3)-2*ighost;
-          params = int64([ighost,ibv,jbv,kbv,nxvl,nyvl,nzvl,nxv,nyv,nzv]);
           data = v(ibv:iev,jbv:jev,kbv:kev);
+          nxvl = size(data,1)-2*ighost;
+          nyvl = size(data,2)-2*ighost;
+          nzvl = size(data,3)-2*ighost;
+          params = int64([ighost,ibv,jbv,kbv,nxvl,nyvl,nzvl,nxv,nyv,nzv]);
           obj.appendSet(data,params);
 
           % Write w
-          nxwl = size(w,1)-2*ighost;
-          nywl = size(w,2)-2*ighost;
-          nzwl = size(w,3)-2*ighost;
-          params = int64([ighost,ibw,jbw,kbw,nxwl,nywl,nzwl,nxw,nyw,nzw]);
           data = w(ibw:iew,jbw:jew,kbw:kew);
+          nxwl = size(data,1)-2*ighost;
+          nywl = size(data,2)-2*ighost;
+          nzwl = size(data,3)-2*ighost;
+          params = int64([ighost,ibw,jbw,kbw,nxwl,nywl,nzwl,nxw,nyw,nzw]);
           obj.appendSet(data,params);
 
           % Write p
-          nxpl = size(p,1)-2*ighost;
-          nypl = size(p,2)-2*ighost;
-          nzpl = size(p,3)-2*ighost;
-          params = int64([ighost,ibp,jbp,kbp,nxpl,nypl,nzpl,nxp,nyp,nzp]);
           data = p(ibp:iep,jbp:jep,kbp:kep);
+          nxpl = size(data,1)-2*ighost;
+          nypl = size(data,2)-2*ighost;
+          nzpl = size(data,3)-2*ighost;
+          params = int64([ighost,ibp,jbp,kbp,nxpl,nypl,nzpl,nxp,nyp,nzp]);
           obj.appendSet(data,params);
 
           % Close chunk file

matlab/write_uvwp_newmesh_ucf.m

@@ -75,28 +75,51 @@ function [] = write_uvwp_newmesh_ucf(file,...
           obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
           obj.create(fname,'type','field','rank',[iproc,ixproc,iyproc,izproc],'endian',par.Results.endian);
           obj.appendStep(time);
+          ighost = 0;
 
           % Interpolate to refined chunk and write
-          tmp = giu({xun(ibu:ieu),yun(jbu:jeu),zun(kbu:keu)});
+          data = giu({xun(ibu:ieu),yun(jbu:jeu),zun(kbu:keu)});
-          if any(isnan(tmp))
+          if any(isnan(data))
-            warning('NaN detected: u, proc #%d, %d occurences',iproc,sum(isnan(tmp(:))));
+            warning('NaN detected: u, proc #%d, %d occurences',iproc,sum(isnan(data(:))));
           end
-          obj.appendField(tmp,0,ibu,jbu,kbu,nxu,nyu,nzu);
+          nxul = size(data,1)-2*ighost;
-          tmp = giv({xvn(ibv:iev),yvn(jbv:jev),zvn(kbv:kev)});
+          nyul = size(data,2)-2*ighost;
-          if any(isnan(tmp))
+          nzul = size(data,3)-2*ighost;
-            warning('NaN detected: v, proc #%d, %d occurences',iproc,sum(isnan(tmp(:))));
+          params = int64([ighost,ibu,jbu,kbu,nxul,nyul,nzul,nxu,nyu,nzu]);
+          obj.appendSet(data,params);
+                    
+          % v
+          data = giv({xvn(ibv:iev),yvn(jbv:jev),zvn(kbv:kev)});
+          if any(isnan(data))
+            warning('NaN detected: v, proc #%d, %d occurences',iproc,sum(isnan(data(:))));
           end
-          obj.appendField(tmp,0,ibv,jbv,kbv,nxv,nyv,nzv);
+          nxvl = size(data,1)-2*ighost;
-          tmp = giw({xwn(ibw:iew),ywn(jbw:jew),zwn(kbw:kew)});
+          nyvl = size(data,2)-2*ighost;
-          if any(isnan(tmp))
+          nzvl = size(data,3)-2*ighost;
-            warning('NaN detected: w, proc #%d, %d occurences',iproc,sum(isnan(tmp(:))));
+          params = int64([ighost,ibv,jbv,kbv,nxvl,nyvl,nzvl,nxv,nyv,nzv]);
+          obj.appendSet(data,params);
+
+          % w
+          data = giw({xwn(ibw:iew),ywn(jbw:jew),zwn(kbw:kew)});
+          if any(isnan(data))
+            warning('NaN detected: w, proc #%d, %d occurences',iproc,sum(isnan(data(:))));
           end
-          obj.appendField(tmp,0,ibw,jbw,kbw,nxw,nyw,nzw);
+          nxwl = size(data,1)-2*ighost;
-          tmp = gip({xpn(ibp:iep),ypn(jbp:jep),zpn(kbp:kep)});
+          nywl = size(data,2)-2*ighost;
-          if any(isnan(tmp))
+          nzwl = size(data,3)-2*ighost;
-            warning('NaN detected: p, proc #%d, %d occurences',iproc,sum(isnan(tmp(:))));
+          params = int64([ighost,ibw,jbw,kbw,nxwl,nywl,nzwl,nxw,nyw,nzw]);
+          obj.appendSet(data,params);
+          
+          % p
+          data = gip({xpn(ibp:iep),ypn(jbp:jep),zpn(kbp:kep)});
+          if any(isnan(data))
+            warning('NaN detected: p, proc #%d, %d occurences',iproc,sum(isnan(data(:))));
           end
-          obj.appendField(tmp,0,ibp,jbp,kbp,nxp,nyp,nzp);
+          nxpl = size(data,1)-2*ighost;
+          nypl = size(data,2)-2*ighost;
+          nzpl = size(data,3)-2*ighost;
+          params = int64([ighost,ibp,jbp,kbp,nxpl,nypl,nzpl,nxp,nyp,nzp]);
+          obj.appendSet(data,params);
 
           obj.close();
         end