mwtkrayer
/
ucftools
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

initial commit

This commit is contained in:
Michael Stumpf 2018-11-27 09:48:44 +01:00
parent 1feeb75588
commit 00e165cdc2
33 changed files with 3579 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

153
bash/ucf_archive_pack Executable file
View File

@ -0,0 +1,153 @@
#!/bin/bash
display_help() {
(>&2 echo "Usage: $(basename "$0") [-hv] [-f outfile] filename")
(>&2 echo "UCF archive packer")
(>&2 echo)
(>&2 echo " filename path to arbitrary file with correct sequence number")
(>&2 echo " -f, --file output file (default: archive_XXXX.ucf)")
(>&2 echo " -h, --help display this help message")
(>&2 echo " -v, --verbose verbose output")
}
exit_script() {
#trap - SIGINT SIGTERM # clear the trap
#kill -- -$$ # Sends SIGTERM to child/sub processes
(>&2 echo "SIGINT/SIGTERM received: removing archive")
rm $fout
}
# Parse command line arguments
if [ $# -eq 0 ]; then
display_help
exit -1
fi
fout=""
verbose=0
POSITIONAL=()
while [[ $# -gt 0 ]]
do
key="$1"
case $key in
-h|--help)
display_help
exit 0
shift # past argument
;;
-f|--file)
fout="$2"
shift # past argument
shift # past value
;;
-v|--verbose)
verbose=1
shift # past argument
;;
*) # unknown option
POSITIONAL+=("$1") # save it in an array for later
shift # past argument
;;
esac
done
set -- "${POSITIONAL[@]}"
# Parse input filename
din=$(dirname $1)
fin=$(basename $1)
seqnum=$(echo $fin | grep -o '_[0-9]\+')
if [ -z "$seqnum" ]; then
(>&2 echo "[Error] No sequence number found in input filename.")
exit 1;
fi
# Construct output filename if not specified
if [ -z "$fout" ]; then
fout="archive${seqnum}.ucf"
fi
# Display verbose info
if [ $verbose -eq 1 ]; then
(>&2 echo "Creating archive: $fout")
fi
# Check input files
fparam="parameters${seqnum}.asc"
fgrid="grid${seqnum}.bin"
fproc="proc${seqnum}.bin"
fpart="particles${seqnum}.bin"
fbuvwp="uvwp"
fbscal="scal"
# Check if obligatory files are present
if [ ! -s ${din}/${fparam} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${fparam}")
exit 1
fi
if [ ! -s ${din}/${fgrid} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${fgrid}")
exit 1
fi
if [ ! -s ${din}/${fproc} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${fproc}")
exit 1
fi
if [ ! -s ${din}/${fpart} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${fpart}")
exit 1
fi
# Check if all velocity fields exist and are not empty
nproc=$(cat ${din}/${fparam} | grep "nprocs" | awk '{print $NF}')
if [ $verbose -eq 1 ]; then
(>&2 echo "Number of processors: $nproc")
fi
fuvwp=()
for (( ii=0; ii<$nproc; ii++ )); do
ftmp=$(printf uvwp${seqnum}.%05d $ii)
if [ ! -s ${din}/${ftmp} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${ftmp}")
exit 1
fi
fuvwp+=($ftmp)
done
# Check if all scalar fields exist (if any) and are not empty
fscal=()
if test -n "$(shopt -s nullglob; echo scal${seqnum}*)"; then
for (( ii=0; ii<$nproc; ii++ )); do
ftmp=$(printf scal${seqnum}.%05d $ii)
if [ ! -s ${din}/${ftmp} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${ftmp}")
exit 1
fi
fscal+=($ftmp)
done
fi
# Verbose info
if [ $verbose -eq 1 ]; then
(>&2 echo "[x] parameters")
(>&2 echo "[x] grid")
(>&2 echo "[x] processor grid")
(>&2 echo "[x] particles")
(>&2 echo "[x] fluid fields")
if [ ${#fscal[@]} -ne 0 ]; then
(>&2 echo "[x] scalar fields")
else
(>&2 echo "[ ] scalar fields")
fi
fi
# Now tar them (remove seqence number from file names)
if [ $verbose -eq 1 ]; then
flagtar="-v"
fi
trap exit_script SIGINT SIGTERM
tar $flagtar --format ustar --transform="flags=r;s|$seqnum||" --directory=${din} -cf ${fout} ${fparam} ${fgrid} ${fproc} ${fpart} ${fuvwp[@]} ${fscal[@]}
tarexit=$?
# Set exit status accoring to tar
if [ $tarexit -ne 0 ]; then
(>&2 echo "tar failed with exit code $tarexit")
exit 254
fi
exit 0

138
bash/ucftar_pack Executable file
View File

@ -0,0 +1,138 @@
#!/bin/bash
display_help() {
(>&2 echo "Usage: $(basename "$0") [-hv] [-f outfile] filename")
(>&2 echo "UCF archive packer v1.0")
(>&2 echo)
(>&2 echo " filename path to arbitrary file with correct sequence number")
(>&2 echo " -f, --file output file (default: archive_XXXX.ucf)")
(>&2 echo " -h, --help display this help message")
(>&2 echo " -v, --verbose verbose output")
}
# Parse command line arguments
if [ $# -eq 0 ]; then
display_help
exit -1
fi
fout=""
verbose=0
POSITIONAL=()
while [[ $# -gt 0 ]]
do
key="$1"
case $key in
-h|--help)
display_help
exit 0
shift # past argument
;;
-f|--file)
fout="$2"
shift # past argument
shift # past value
;;
-v|--verbose)
verbose=1
shift # past argument
;;
*) # unknown option
POSITIONAL+=("$1") # save it in an array for later
shift # past argument
;;
esac
done
set -- "${POSITIONAL[@]}"
# Parse input filename
din=$(dirname $1)
fin=$(basename $1)
seqnum=$(echo $fin | grep -o '_[0-9]\+')
if [ -z "$seqnum" ]; then
(>&2 echo "[Error] No sequence number found in input filename.")
exit 1;
fi
# Construct output filename if not specified
if [ -z "$fout" ]; then
fout="archive${seqnum}.ucf"
fi
# Display verbose info
if [ $verbose -eq 1 ]; then
(>&2 echo "Creating archive: $fout")
fi
# Check input files
fparam="parameters${seqnum}.asc"
fgrid="grid${seqnum}.bin"
fproc="proc${seqnum}.bin"
fpart="particles${seqnum}.bin"
fbuvwp="uvwp"
fbscal="scal"
# Check if obligatory files are present
if [ ! -s ${din}/${fparam} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${fparam}")
exit 1
fi
if [ ! -s ${din}/${fgrid} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${fgrid}")
exit 1
fi
if [ ! -s ${din}/${fproc} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${fproc}")
exit 1
fi
if [ ! -s ${din}/${fpart} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${fpart}")
exit 1
fi
# Check if all velocity fields exist and are not empty
nproc=$(cat ${din}/${fparam} | grep "nprocs" | awk '{print $NF}')
if [ $verbose -eq 1 ]; then
(>&2 echo "Number of processors: $nproc")
fi
fuvwp=()
for (( ii=0; ii<$nproc; ii++ )); do
ftmp=$(printf uvwp${seqnum}.%05d $ii)
if [ ! -s ${din}/${ftmp} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${ftmp}")
exit 1
fi
fuvwp+=($ftmp)
done
# Check if all scalar fields exist (if any) and are not empty
if grep -q "scalar" "${din}/${fparam}"; then
fscal=()
for (( ii=0; ii<$nproc; ii++ )); do
ftmp=$(printf scal${seqnum}.%05d $ii)
if [ ! -s ${din}/${ftmp} ]; then
(>&2 echo "[Error] File not found or empty: ${din}/${ftmp}")
exit 1
fi
fscal+=($ftmp)
done
fi
# Verbose info
if [ $verbose -eq 1 ]; then
(>&2 echo "[x] parameters")
(>&2 echo "[x] grid")
(>&2 echo "[x] processor grid")
(>&2 echo "[x] particles")
(>&2 echo "[x] fluid fields")
if [ ${#fscal[@]} -ne 0 ]; then
(>&2 echo "[x] scalar fields")
else
(>&2 echo "[ ] scalar fields")
fi
fi
# Now tar them (remove seqence number from file names)
if [ $verbose -eq 1 ]; then
flagtar="-v"
fi
tar $flagtar --format ustar --transform="flags=r;s|$seqnum||" --directory=${din} -cf ${fout} ${fparam} ${fgrid} ${fproc} ${fpart} ${fuvwp[@]} ${fscal[@]}

240
matlab/@ini/ini.m Normal file
View File

@ -0,0 +1,240 @@
classdef ini < handle
% Low-level utilities for INI files.
properties (Access = public)
File % file name
IOMode % file opened in read-only or read-write mode?
FileSize % file size
NumSection % maximum number of datasets in this file (over all time steps)
Verbosity % verbose output?
end
%%
properties (Access = private)
% File info
fileID
fileBeg
fileEnd
ioflag
tarflag
content
end
%% ------------------------------------------------------------------------%%
%% CONSTRUCORS/DESTRUCTORS %%
%% ------------------------------------------------------------------------%%
methods (Access = public)
function obj = ini(varargin)
% [obj] = ini(varargin)
% Default contructor
% Input
% ? verbosity verbose output? (default: no)
par = inputParser;
addOptional(par,'verbosity',0,@isnumeric);
parse(par,varargin{:});
obj.resetPublicProperties();
obj.resetPrivateProperties();
obj.setVerbosity(par.Results.verbosity);
end
function delete(obj)
% obj.delete()
% Object destructor
obj.close();
end
end
%% ------------------------------------------------------------------------%%
%% INITIALIZATION METHODS %%
%% ------------------------------------------------------------------------%%
methods(Access=public)
function open(obj,fname)
% obj.open()
% Opens an INI file.
% Input
% fname file name
obj.File = fname;
obj.IOMode = 'read';
obj.ioflag = 'r';
obj.tarflag = false;
obj.fileID = fopen(obj.File,obj.ioflag);
obj.fileBeg = 0;
fseek(obj.fileID,0,'eof');
obj.fileEnd = ftell(obj.fileID);
obj.FileSize = obj.fileEnd-obj.fileBeg;
obj.parse();
end
function opentar(obj,ptr)
% obj.opentar(ptr)
% Opens a subfile from tar in read-only mode.
% Input
% ptr pointer: [fid,first byte,number of bytes]
obj.File = 'tar-mode'; % Set generic file name
obj.IOMode = 'read'; % abstract IO mode
obj.ioflag = 'r'; % internal IO mode
obj.tarflag = true; % reading from tar archive?
% Set file ID to tar file ID
obj.fileID = ptr(1);
if obj.fileID<0
error('Unable to access file: %s',obj.File);
end
obj.fileBeg = ptr(2);
obj.fileEnd = ptr(2)+ptr(3);
obj.FileSize = ptr(3);
obj.parse();
end
end
methods(Access=public)
function close(obj)
% obj.close()
% Closes a file.
tarmode = obj.tarflag;
obj.resetPublicProperties();
obj.resetPrivateProperties();
if obj.fileID<0 || tarmode
return;
end
status = fclose(obj.fileID);
if status<0
warning('Error while closing: %d',status);
return;
end
obj.fileID = -1;
end
function setVerbosity(obj,flag)
% obj.setVerbosity(flag)
% Sets verbosity
% Input
% flag '1' -> verbose output
obj.Verbosity = flag;
end
function [data] = getContent(obj)
% [data] = obj.getContent()
% Get the entire content of the INI file.
% Output
% data entire content of file as structure
data = obj.content;
end
function [data] = getSection(obj,section)
% [data] = obj.getSection(section)
% Get the content of a single section of the INI file.
% Input
% section name of the requested section
% Output
% data content of one section of file as structure
if isfield(obj.content,section)
data = obj.content.(section);
else
error('Section not found: %s',name);
end
end
function [data] = getEntry(obj,section,entry)
% [data] = obj.getEntry(section,entry)
% Get the value of a single entry within the INI file.
% Input
% section name of the requested section
% entry name of the entry within requested section
% Output
% data value of the requested entry
if ~isfield(obj.content,section)
error('Section not found: %s',section);
elseif ~isfield(obj.content.(section),entry)
error('Entry not found: %s',entry);
else
data = obj.content.(section).(entry);
end
end
function [data] = listSections(obj)
% [data] = obj.listSections()
% Get a list of section names.
% Output
% data cell array with section names
data = fieldnames(obj.content);
end
function [data] = listEntries(obj,section)
% [data] = obj.listEntries(section)
% Get a list of entry names within a section..
% Input
% section name of the requested section
% Output
% data cell array with section names
data = fieldnames(obj.content.(section));
end
end
%% ------------------------------------------------------------------------%%
%% PRIVATE STATIC METHODS %%
%% ------------------------------------------------------------------------%%
methods(Access=private)
function parse(obj)
fseek(obj.fileID,obj.fileBeg,'bof');
% Construct a struct with content of file
obj.content = [];
curSection = '';
while ftell(obj.fileID)<obj.fileEnd
str = strtrim(fgetl(obj.fileID));
% Check for empty lines and comments (; or #)
if isempty(str)
continue;
end
if (str(1)==';')
continue;
end
if (str(1)=='#')
continue;
end
% Check for section
if (str(1)=='[') && (str(end)==']')
% This is a section: create a structure
curSection = ini.getUniqueVarname(str(2:end-1));
obj.content.(curSection) = [];
else
% This is not a section: create structure entry
[curMember,curValue] = strtok(str,'=');
curValue = ini.trimValue(curValue);
curMember = ini.getUniqueVarname(curMember);
% We now have the field name and the corresponding value as strings.
% Convert the value to a numeric value if possible.
[tmp,flag] = str2num(curValue);
if flag
curValue = tmp;
end
% Add this structure to the current section
if ~isempty(curSection)
obj.content.(curSection).(curMember) = curValue;
else
obj.content.(curMember) = curValue;
end
end
end
obj.NumSection = numel(fieldnames(obj.content));
end
function resetPublicProperties(obj)
obj.File = [];
obj.FileSize = [];
obj.IOMode = [];
obj.NumSection = [];
obj.Verbosity = [];
end
function resetPrivateProperties(obj)
obj.fileBeg = [];
obj.fileEnd = [];
obj.ioflag = [];
obj.tarflag = [];
obj.content = [];
end
end
%% ------------------------------------------------------------------------%%
%% PRIVATE STATIC METHODS %%
%% ------------------------------------------------------------------------%%
methods(Access=private,Static)
function str = trimValue(str)
str = strtrim(str);
if strcmpi(str(1),'=')
str(1)=[];
end
str = strtrim(str);
end
function varname = getUniqueVarname(varname)
%tmp = matlab.lang.makeValidName(varname);
%varname = matlab.lang.makeUniqueStrings(tmp,who,namelengthmax);
varname = strrep(varname,'-','_');
varname = genvarname(varname);
end
end
end

796
matlab/@ucf/ucf.m Normal file
View File

@ -0,0 +1,796 @@
classdef ucf < handle
% Low-level utilities for UCF files.
properties (Access = public)
File % file name
Type % file type
Class % file class
Endian % endianess
CodeVersion % version of the simulation code
UCFVersion % version of the data format ("unified container format")
NumDataset % maximum number of datasets in this file (over all time steps)
NumTimestep % number of time steps in this file
FileSize % file size
CreationTime % time of creation
IOMode % file opened in read-only or read-write mode?
IORank % rank of processor + col,row,pln
Verbosity % verbose output?
Debug % debug information?
end
properties (Access = private)
% File info
fileID
fileBeg
fileEnd
typeID
ioflag
tarflag
creationTimeUnix
versionMajor
versionMinor
versionPatch
versionFile
% Step pointers and info
posStep
numSetPerStep
timeStep
% Write-mode
isFileHeaderWritten
isStepHeaderWritten
% Current step information (should be resetable by resetCurrentStep())
currentStep
currentStepPosHeader
currentStepPosData
currentStepSize
currentStepTime
currentStepNumSet
% Current set information (should be resetable by resetCurrentSet())
currentSet
currentSetPosHeader
currentSetPosData
currentSetSize
currentSetDatatype
currentSetDatatypeNumeric
currentSetSizeof
currentSetNumParams
currentSetParams
currentSetNumElements
% Constants
magicFile = int64(81985529216486895);
magicStep = int64(11944304052957);
magicSet = int64(240217520921210);
nHeaderFile = 8;
nHeaderStep = 4;
nHeaderSet = 4;
nByteHeaderFile = 8;
nByteHeaderStep = 8;
nByteHeaderSet = 8;
nSetParamsField = 10;
nSetParamsParticle = 16;
factorMajor = 1000000000;
factorMinor = 1000000;
factorPatch = 1000;
factorTypeIDClass = 1000;
factorTypeIDKind = 10;
typeIDmatlabField = 1999;
typeIDmatlabParticle = 2999;
scanBuffSize = 4096;
end
%% ------------------------------------------------------------------------%%
%% CONSTRUCORS/DESTRUCTORS %%
%% ------------------------------------------------------------------------%%
methods(Access=public)
function obj = ucf(varargin)
% obj = ucf(varargin)
% Default contructor
% Input
% ? verbosity verbose output? (default: false)
% ? debug debug output? (default: false)
par = inputParser;
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
obj.resetPublicProperties();
obj.resetPrivateProperties();
obj.resetCurrentStep();
obj.resetCurrentSet();
obj.setVerbosity(par.Results.verbosity);
obj.setDebug(par.Results.debug);
end
function delete(obj)
% obj.delete()
% Default destructor
obj.close();
end
end
%% ------------------------------------------------------------------------%%
%% INITIALIZATION METHODS %%
%% ------------------------------------------------------------------------%%
methods(Access=public)
function open(obj,file)
% obj.open(file)
% Opens a file in read-only mode
obj.File = file; % add filename as object property
obj.IOMode = 'read'; % abstract IO mode
obj.ioflag = 'r'; % internal IO mode
obj.tarflag = false; % data is within a tar file?
% Try to open the file
obj.fileID = fopen(obj.File,obj.ioflag);
if obj.fileID<0
error('Unable to open file: %s',obj.File);
end
obj.fileBeg = ftell(fid);
fseek(obj.fileID,0,'eof');
obj.fileEnd = ftell(obj.fileID);
fseek(obj.fileID,0,'bof');
obj.FileSize = obj.fileEnd-obj.fileBeg;
% Read the file header
obj.readHeaderFile()
% Scan through file to get the basic structure (steps/sets)
obj.timeStep = zeros(1,obj.scanBuffSize,'double');
obj.posStep = zeros(1,obj.scanBuffSize,'double');
obj.numSetPerStep = zeros(1,obj.scanBuffSize,'double');
istep = 0;
while ftell(obj.fileID)<obj.FileSize
obj.readHeaderStep();
istep = istep+1;
obj.timeStep(istep) = obj.currentStepTime;
obj.posStep(istep) = ftell(obj.fileID);
obj.numSetPerStep(istep) = obj.currentStepNumSet;
if obj.currentStepSize==-1
break;
else
fseek(obj.fileID,obj.currentStepSize,'cof');
end
end
nstep = istep;
% Truncate buffered arrays
if nstep>obj.scanBuffSize
warning('Buffer overflow detected: increase scanBuffSize.');
end
obj.timeStep = obj.timeStep(1:nstep);
obj.posStep = obj.posStep(1:nstep);
obj.numSetPerStep = obj.numSetPerStep(1:nstep);
% Set some internal variables
obj.NumDataset = max(obj.numSetPerStep);
obj.NumTimestep = nstep;
obj.isFileHeaderWritten = true;
obj.isStepHeaderWritten = true;
end
function edit(obj,file)
% obj.edit(file)
% Opens a file in read-write mode
obj.File = file; % add filename as object property
obj.IOMode = 'edit'; % abstract IO mode
obj.ioflag = 'r+'; % internal IO mode
obj.tarflag = false; % data is within a tar file?
% Try to open the file
obj.fileID = fopen(obj.File,obj.ioflag);
if obj.fileID<0
error('Unable to open file: %s',obj.File);
end
% Read the file header
obj.readHeaderFile()
% Scan through file to get the basic structure (steps/sets)
obj.timeStep = zeros(1,obj.scanBuffSize,'double');
obj.posStep = zeros(1,obj.scanBuffSize,'double');
obj.numSetPerStep = zeros(1,obj.scanBuffSize,'double');
istep = 0;
while ftell(obj.fileID)<obj.FileSize
obj.readHeaderStep();
istep = istep+1;
obj.timeStep(istep) = obj.currentStepTime;
obj.posStep(istep) = ftell(obj.fileID);
obj.numSetPerStep(istep) = obj.currentStepNumSet;
if obj.currentStepSize==-1
break;
else
fseek(obj.fileID,obj.currentStepSize,'cof');
end
end
nstep = istep;
% Truncate buffered arrays
if nstep>obj.scanBuffSize
warning('Buffer overflow detected: increase scanBuffSize.');
end
obj.timeStep = obj.timeStep(1:nstep);
obj.posStep = obj.posStep(1:nstep);
obj.numSetPerStep = obj.numSetPerStep(1:nstep);
% Set some internal variables
obj.NumDataset = max(obj.numSetPerStep);
obj.NumTimestep = nstep;
obj.isFileHeaderWritten = true;
obj.isStepHeaderWritten = true;
end
function create(obj,file,varargin)
% obj.create(file,varargin)
% Creates a new file and writes file header
% Input
% file file name/path
% ? type file type {'field' (default),'particle'}
% ? rank proc rank (default: [0,0,0,0])
% ? endian endianess {'n' (dafault),'a','l','s','b'}
par = inputParser;
addOptional(par,'type','field',@ischar);
addOptional(par,'rank',[0,0,0,0],@(x)isnumeric(x)&&numel(x)==4);
addOptional(par,'endian','n',@ischar);
parse(par,varargin{:});
%
obj.File = file; % add filename as object property
obj.IOMode = 'create'; % abstract IO mode
obj.ioflag = 'w+'; % internal IO mode
obj.tarflag = false; % data is within a tar file?
obj.isFileHeaderWritten = false;
obj.isStepHeaderWritten = false;
% Try to create the file
obj.fileID = fopen(obj.File,obj.ioflag);
if obj.fileID<0
error('Unable to create file: %s',obj.File);
end
%
switch par.Results.type
case 'field'
obj.typeID = obj.typeIDmatlabField;
obj.Class = 'field';
case 'particle'
obj.typeID = obj.typeIDmatlabParticle;
obj.Class = 'particle';
otherwise
error('Unknown file type: %s',par.Results.type);
end
obj.IORank = par.Results.rank;
obj.Endian = par.Results.endian;
obj.versionMajor = 0;
obj.versionMinor = 0;
obj.versionPatch = 0;
obj.versionFile = 1;
obj.writeHeaderFile();
end
function opentar(obj,ptr)
% obj.opentar(ptr)
% Opens a subfile from tar in read-only mode
obj.File = 'tar-mode'; % Set generic file name
obj.IOMode = 'read'; % abstract IO mode
obj.ioflag = 'r'; % internal IO mode
obj.tarflag = true; % data is within a tar file?
% Set file ID to tar file ID
obj.fileID = ptr(1);
if obj.fileID<0
error('Unable to access file: %s',obj.File);
end
obj.fileBeg = ptr(2);
obj.fileEnd = ptr(2)+ptr(3);
obj.FileSize = ptr(3);
% Read the file header
obj.readHeaderFile()
% Scan through file to get the basic structure (steps/sets)
obj.timeStep = zeros(1,obj.scanBuffSize,'double');
obj.posStep = zeros(1,obj.scanBuffSize,'double');
obj.numSetPerStep = zeros(1,obj.scanBuffSize,'double');
istep = 0;
while ftell(obj.fileID)<obj.fileEnd
obj.readHeaderStep();
istep = istep+1;
obj.timeStep(istep) = obj.currentStepTime;
obj.posStep(istep) = ftell(obj.fileID);
obj.numSetPerStep(istep) = obj.currentStepNumSet;
if obj.currentStepSize==-1
break;
else
fseek(obj.fileID,obj.currentStepSize,'cof');
end
end
nstep = istep;
% Truncate buffered arrays
if nstep>obj.scanBuffSize
warning('Buffer overflow detected: increase scanBuffSize.');
end
obj.timeStep = obj.timeStep(1:nstep);
obj.posStep = obj.posStep(1:nstep);
obj.numSetPerStep = obj.numSetPerStep(1:nstep);
% Set some internal variables
obj.NumDataset = max(obj.numSetPerStep);
obj.NumTimestep = nstep;
obj.isFileHeaderWritten = true;
obj.isStepHeaderWritten = true;
end
end
%% ------------------------------------------------------------------------%%
%% PUBLIC METHODS %%
%% ------------------------------------------------------------------------%%
methods(Access=public)
function close(obj)
% obj.close()
% Closes a file
tarmode = obj.tarflag;
obj.resetPublicProperties();
obj.resetPrivateProperties();
obj.resetCurrentStep();
obj.resetCurrentSet();
if obj.fileID<0 || tarmode
return;
end
status = fclose(obj.fileID);
if status<0
warning('Unable to close file (exit code: %d)',status);
return;
end
obj.fileID = -1;
end
function [isType] = validateType(obj,type)
switch type
case 'grid'
isType = (obj.typeID==0);
case 'procgrid'
isType = (obj.typeID==10);
case 'field'
isType = (floor(obj.typeID/obj.factorTypeIDClass)==1);
case 'particle'
isType = (floor(obj.typeID/obj.factorTypeIDClass)==2);
case 'snapshot'
isType = (mod(obj.typeID,obj.factorTypeIDKind)==0);
case 'append'
isType = (mod(obj.typeID,obj.factorTypeIDKind)==1);
case {'uvwp','fluid'}
isType = (obj.typeID==1000);
case 'scalar'
isType = (obj.typeID==1010);
case 'matlab'
isType = (mod(obj.typeID,obj.factorTypeIDKind)==9);
otherwise
error('Unknown type: %s',type);
end
end
% ------------------------------------------------------------------------%
% Read %
% ------------------------------------------------------------------------%
function [data,params] = readSet(obj,tstep,dset)
% [data,params] = obj.readSet(tstep,dset)
% Reads a raw dataset from file (no reshaping or parsing of parameters)
% Input
% tstep index of timestep to be read
% dset index of dataset to be read
% Output
% data raw data with dim(1,nelements)
% params raw parameters with dim(1,nparams)
fseek(obj.fileID,obj.findSet(tstep,dset),'bof');
obj.readHeaderSet();
params = obj.currentSetParams;
data = fread(obj.fileID,[1,obj.currentSetNumElements],...
[obj.currentSetDatatype,'=>',obj.currentSetDatatype]);
end
function [params] = readParameters(obj,tstep,dset)
% [data,params] = obj.readParameters(tstep,dset)
% Reads a raw set of parameters without parsing.
% Input
% tstep index of timestep to be read
% dset index of dataset to be read
% Output
% params raw parameters with dim(1,nparams)
fseek(obj.fileID,obj.findSet(tstep,dset),'bof');
if obj.Debug
fprintf('Skipped to position %d\n',ftell(obj.fileID));
end
obj.readHeaderSet();
params = obj.currentSetParams;
end
% ------------------------------------------------------------------------%
% Append %
% ------------------------------------------------------------------------%
function appendStep(obj,time)
% [] = obj.appendStep(time)
% Creates a new step at the end of opened file.
% Input
% time simulation time
if ~obj.verifyWriteMode()
error('File must be opened in create/edit mode.');
end
%
obj.resetCurrentStep();
obj.resetCurrentSet();
%
fseek(obj.fileID,0,'eof');
obj.currentStep = obj.NumTimestep+1;
obj.currentStepPosHeader = ftell(obj.fileID);
obj.currentStepSize = 0;
obj.currentStepTime = time;
obj.currentStepNumSet = 0;
obj.writeHeaderStep();
obj.currentStepPosData = ftell(obj.fileID);
obj.currentSet = 0;
%
obj.NumTimestep = obj.currentStep;
obj.posStep(obj.currentStep) = obj.currentStepPosData;
obj.numSetPerStep(obj.currentStep) = obj.currentStepNumSet;
obj.timeStep(obj.currentStep) = obj.currentStepTime;
%
if obj.Verbosity
fprintf('Created timestep #%d with time = %f\n',obj.NumTimestep,obj.currentStepTime);
end
end
function appendSet(obj,data,params)
% [] = obj.appendSet(data,params)
% Creates a new dataset at the end of opened file and writes data.
% Input
% data data to append
% parameters dataset header parameters
if ~obj.verifyWriteMode()
error('File must be opened in create/edit mode.');
end
% Evaluate datatype
obj.currentSetDatatype = class(data);
switch obj.currentSetDatatype
case 'int32'; obj.currentSetDatatypeNumeric = 11; obj.currentSetSizeof = 4;
case 'int64'; obj.currentSetDatatypeNumeric = 12; obj.currentSetSizeof = 8;
case 'single'; obj.currentSetDatatypeNumeric = 21; obj.currentSetSizeof = 4;
case 'double'; obj.currentSetDatatypeNumeric = 22; obj.currentSetSizeof = 8;
otherwise
error('Datatype not implemented: %s',obj.currentSetDatatype);
end
% Append a new set header
fseek(obj.fileID,0,'eof');
obj.currentSet = obj.numSetPerStep(obj.currentStep)+1;
obj.currentSetPosHeader = ftell(obj.fileID);
obj.currentSetSize = numel(data)*obj.currentSetSizeof;
switch obj.Class
case 'field'; obj.currentSetNumParams = obj.nSetParamsField;
case 'particle'; obj.currentSetNumParams = obj.nSetParamsParticle;
otherwise; error('Invalid file class: %s',obj.Class);
end
obj.currentSetParams = zeros(1,obj.currentSetNumParams,'int64');
obj.currentSetParams(1:numel(params)) = params;
obj.writeHeaderSet();
% Append actual data
obj.currentSetPosData = ftell(obj.fileID);
fwrite(obj.fileID,data(:),obj.currentSetDatatype,0,obj.Endian);
% Rewrite step header
obj.currentStepNumSet = obj.currentStepNumSet+1;
obj.currentStepSize = obj.currentStepSize+...
(4+obj.currentSetNumParams)*obj.nByteHeaderStep+...
obj.currentSetSize;
obj.writeHeaderStep();
% Update internal variables
obj.numSetPerStep(obj.currentStep) = obj.currentSet;
obj.NumDataset = max(obj.currentSet,obj.NumDataset);
% Verbose output
if obj.Verbosity
fprintf('Created dataset %d at timestep #%d\n',obj.currentSet,obj.currentStep);
end
end
function exchangeSet(obj,istep,iset,data)
% [] = obj.exchangeSet(istep,iset,data)
% Exchanges the data stored in dataset (istep,iset) with the data provided.
% Input
% data data, which replaces dataset (must be same size and type)
if ~obj.verifyWriteMode()
error('File must be opened in create/edit mode.');
end
% Evaluate datatype and datasize
dtype = class(data);
switch dtype
case 'int32'; dtypeNumeric = 11; dtypeSizeof = 4;
case 'int64'; dtypeNumeric = 12; dtypeSizeof = 8;
case 'single'; dtypeNumeric = 21; dtypeSizeof = 4;
case 'double'; dtypeNumeric = 22; dtypeSizeof = 8;
otherwise
error('Datatype not implemented: %s',dtype);
end
dsize = numel(data)*dtypeSizeof;
% Navigate to correct dataset
fseek(obj.fileID,obj.findSet(istep,iset),'bof');
obj.readHeaderSet();
% Check if datatype/datasize match
if dtypeNumeric~=obj.currentSetDatatypeNumeric
error('Datatype mismatch: %s,%s',dtype,obj.currentSetDatatype)
end
if dsize~=obj.currentSetSize
error('Datasize mismatch: %d,%d',dsize,obj.currentSetSize)
end
% Write data
fwrite(obj.fileID,data(:),obj.currentSetDatatype,0,obj.Endian);
% Verbose output
if obj.Verbosity
fprintf('Edited dataset %d at timestep #%d\n',obj.currentSet,obj.currentStep);
end
end
% ------------------------------------------------------------------------%
% Setter %
% ------------------------------------------------------------------------%
function setVerbosity(obj,flag)
% obj.setVerbosity(flag)
% Sets verbosity
% Input
% flag '1' -> verbose output
obj.Verbosity = flag;
end
function setDebug(obj,flag)
% obj.setDebug(flag)
% Sets debug flag
% Input
% flag '1' -> debug output
obj.Debug = flag;
end
% ------------------------------------------------------------------------%
% Getter %
% ------------------------------------------------------------------------%
function [nstep] = getNumTimesteps(obj)
% [nstep] = obj.getNumTimesteps()
% Gets total number of time steps.
nstep = obj.NumTimestep;
end
function [stime] = getSimulationTime(obj,varargin)
% [stime] = obj.getSimulationTime(varargin)
% Gets simulation time for one/all time steps
% Input
% ? istep step index (default: -1, all steps)
par = inputParser;
addOptional(par,'istep',-1,@isnumeric);
parse(par,varargin{:});
if par.Results.istep==-1
stime = obj.timeStep;
else
stime = obj.timeStep(par.Results.istep);
end
end
function [ndset] = getNumDatasets(obj,istep)
% [ndset] = obj.getNumDatasets(istep)
% Gets number datasets for a given time step.
% Input
% istep step index
ndset = obj.numSetPerStep(istep);
end
end
%% ------------------------------------------------------------------------%%
%% PRIVATE METHODS %%
%% ------------------------------------------------------------------------%%
methods (Access = private)
function readHeaderFile(obj)
% Skip to beginning of file
fseek(obj.fileID,obj.fileBeg,'bof');
% Determine endianess
mfmt = 'aslb';
for fmt=mfmt
currentMagic = fread(obj.fileID,1,'int64=>int64',0,fmt);
fseek(obj.fileID,obj.fileBeg,'bof');
if currentMagic==obj.magicFile
break;
end
end
if currentMagic~=obj.magicFile
error('Magic mismatch: invalid file header. %d',currentMagic);
end
obj.Endian = fmt;
% Read header
header = fread(obj.fileID,[1,obj.nHeaderFile],'int64=>int64',0,obj.Endian);
if obj.Debug
fprintf('Read the following file header at %d bytes\n',0);
fprintf('%d,%d,%d,%d,%d,%d,%d,%d\n',header);
end
% Parse version
obj.versionMajor = floor(double(header(2))/obj.factorMajor);
obj.versionMinor = floor(mod(double(header(2)),obj.factorMajor)/obj.factorMinor);
obj.versionPatch = floor(mod(double(header(2)),obj.factorMinor)/obj.factorPatch);
obj.CodeVersion = sprintf('%d.%d.%d',obj.versionMajor,obj.versionMinor,obj.versionPatch);
obj.UCFVersion = mod(double(header(2)),obj.factorPatch);
% Parse time stamp (UTC)
obj.creationTimeUnix = double(header(3));
obj.CreationTime = datestr(obj.creationTimeUnix/86400 + datenum(1970,1,1));
% Parse file type
obj.typeID = double(header(4));
switch obj.typeID
case 0000; obj.Type = 'grid';
case 0010; obj.Type = 'processor grid';
case 1000; obj.Type = 'fluid snapshot';
case 1010; obj.Type = 'scalar snapshot';
case 1999; obj.Type = 'matlab field data';
case 2000; obj.Type = 'particle snapshot';
case 2001; obj.Type = 'particle append';
case 2011; obj.Type = 'particle lagrange';
case 2021; obj.Type = 'particle balancing';
case 2999; obj.Type = 'matlab particle data';
case 3000; obj.Type = 'statistics fluid';
case 3010; obj.Type = 'statistics fluid pure';
case 3020; obj.Type = 'statistics scalar';
otherwise; obj.Type = 'unknown';
end
% Parse file class
switch floor(obj.typeID/obj.factorTypeIDClass)
case 1; obj.Class = 'field';
case 2; obj.Class = 'particle';
case 3; obj.Class = 'statistics';
otherwise; obj.Class = 'unknown';
end
% Parse IO rank
obj.IORank = double(header(5:8));
end
function readHeaderStep(obj)
% Read and parse
obj.currentStepPosHeader = ftell(obj.fileID);
header = fread(obj.fileID,[1,obj.nHeaderStep],'int64=>int64',0,obj.Endian);
obj.currentStepPosData = ftell(obj.fileID);
currentMagic = header(1);
obj.currentStepSize = double(header(2));
obj.currentStepTime = typecast(header(3),'double');
obj.currentStepNumSet = double(header(4));
if obj.Debug
fprintf('Read the following step header at %d bytes\n',obj.currentStepPosHeader);
fprintf('%d,%d,%f,%d\n',header(1),header(2),typecast(header(3),'double'),header(4));
end
% Check if magic is correct
if currentMagic~=obj.magicStep
error('Magic mismatch: invalid step header. %d',currentMagic);
end
end
function readHeaderSet(obj)
% Read and parse
obj.currentSetPosHeader = ftell(obj.fileID);
header = fread(obj.fileID,[1,obj.nHeaderSet],'int64=>int64',0,obj.Endian);
obj.currentSetPosData = ftell(obj.fileID);
currentMagic = header(1);
obj.currentSetSize = double(header(2));
obj.currentSetDatatypeNumeric = double(header(3));
switch obj.currentSetDatatypeNumeric
case 11; obj.currentSetSizeof = 4; obj.currentSetDatatype = 'int32';
case 12; obj.currentSetSizeof = 8; obj.currentSetDatatype = 'int64';
case 21; obj.currentSetSizeof = 4; obj.currentSetDatatype = 'single';
case 22; obj.currentSetSizeof = 8; obj.currentSetDatatype = 'double';
otherwise
error('Unknown datatype: %d',obj.currentSetDatatypeNumeric);
end
obj.currentSetNumParams = double(header(4));
obj.currentSetNumElements = round(obj.currentSetSize/obj.currentSetSizeof);
if obj.Debug
fprintf('Read the following set header at %d bytes\n',obj.currentSetPosHeader);
fprintf('%d,%d,%d,%d\n',header(1),header(2),header(3),header(4));
end
% Check if magic is correct
if currentMagic~=obj.magicSet
error('Magic mismatch: invalid dataset header. %d',currentMagic);
end
% Read variable number of parameters
obj.currentSetParams = fread(obj.fileID,[1,obj.currentSetNumParams],'int64=>int64',0,obj.Endian);
if obj.Debug
fprintf('with parameters:')
for ii=1:obj.currentSetNumParams
fprintf(' %d',obj.currentSetParams(ii));
end
fprintf('\n');
end
end
function [flag] = verifyWriteMode(obj)
flag = false;
switch obj.IOMode
case {'edit','create'}
flag = true;
end
end
function writeHeaderFile(obj)
obj.creationTimeUnix = int64(java.lang.System.currentTimeMillis/1000); % using java for compatability before Matlab2014a
header = zeros(1,obj.nHeaderFile,'int64');
header(1) = obj.magicFile;
header(2) = obj.versionMajor*obj.factorMajor + ...
obj.versionMinor*obj.factorMinor + ...
obj.versionPatch*obj.factorPatch + ...
obj.versionFile;
header(3) = obj.creationTimeUnix;
header(4) = int64(obj.typeID);
header(5:8) = int64(obj.IORank);
fseek(obj.fileID,obj.fileBeg,'bof');
if obj.Debug
fprintf('Writing the following file header at %d bytes\n',ftell(obj.fileID));
fprintf('%d,%d,%d,%d,%d,%d,%d,%d\n',header);
end
fwrite(obj.fileID,header,'int64',0,obj.Endian);
obj.isFileHeaderWritten = true;
end
function writeHeaderStep(obj)
if ~obj.isFileHeaderWritten
error('No file header has been written yet.');
end
fseek(obj.fileID,obj.currentStepPosHeader,'bof');
if obj.Debug
fprintf('Writing the following step header at %d bytes\n',ftell(obj.fileID));
fprintf('%d,%d,%f,%d\n',obj.magicStep,obj.currentStepSize,obj.currentStepTime,obj.currentStepNumSet);
end
fwrite(obj.fileID,obj.magicStep,'int64',0,obj.Endian);
fwrite(obj.fileID,obj.currentStepSize,'int64',0,obj.Endian);
fwrite(obj.fileID,obj.currentStepTime,'double',0,obj.Endian);
fwrite(obj.fileID,obj.currentStepNumSet,'int64',0,obj.Endian);
obj.isStepHeaderWritten = true;
end
function writeHeaderSet(obj)
if ~obj.isStepHeaderWritten
error('No step header has been written yet.');
end
fseek(obj.fileID,obj.currentSetPosHeader,'bof');
header = zeros(1,obj.nHeaderSet,'int64');
header(1) = obj.magicSet;
header(2) = obj.currentSetSize;
header(3) = obj.currentSetDatatypeNumeric;
header(4) = obj.currentSetNumParams;
obj.currentSetParams = obj.currentSetParams(1:obj.currentSetNumParams);
if obj.Debug
fprintf('Writing the following set header at %d bytes\n',ftell(obj.fileID));
fprintf('%d,%d,%d,%d\n',header);
fprintf('with parameters:')
for ii=1:obj.currentSetNumParams
fprintf(' %d',obj.currentSetParams(ii));
end
fprintf('\n');
end
fwrite(obj.fileID,header,'int64',0,obj.Endian);
fwrite(obj.fileID,obj.currentSetParams,'int64',0,obj.Endian);
end
function [posHeader] = findSet(obj,tstep,dset)
% Check input
if tstep>obj.NumTimestep
error('Out of bounds: timestep. %d, %d',tstep,obj.NumTimestep);
end
if dset>obj.numSetPerStep(tstep)
error('Out of bounds: dataset. %d, %d',dset,obj.NumDataset);
end
% Navigate to correct set
fseek(obj.fileID,obj.posStep(tstep),'bof');
for iset=1:dset-1
obj.readHeaderSet();
fseek(obj.fileID,obj.currentSetSize,'cof');
end
posHeader = ftell(obj.fileID);
if obj.Debug
fprintf('Found step #%d, set #%d at position %d\n',tstep,dset,posHeader)
end
end
function resetPublicProperties(obj)
obj.File = [];
obj.Type = [];
obj.Class = [];
obj.Endian = [];
obj.CodeVersion = [];
obj.UCFVersion = [];
obj.NumDataset = [];
obj.NumTimestep = [];
obj.FileSize = [];
obj.CreationTime = [];
obj.IOMode = [];
obj.IORank = [];
obj.Verbosity = [];
obj.Debug = [];
end
function resetPrivateProperties(obj)
obj.fileBeg = [];
obj.fileEnd = [];
obj.typeID = [];
obj.ioflag = [];
obj.tarflag = [];
obj.creationTimeUnix = [];
obj.versionMajor = [];
obj.versionMinor = [];
obj.versionPatch = [];
obj.versionFile = [];
obj.posStep = [];
obj.numSetPerStep = [];
obj.timeStep = [];
obj.isFileHeaderWritten = [];
obj.isStepHeaderWritten = [];
end
function resetCurrentStep(obj)
obj.currentStep = [];
obj.currentStepPosHeader = [];
obj.currentStepPosData = [];
obj.currentStepSize = [];
obj.currentStepTime = [];
obj.currentStepNumSet = [];
end
function resetCurrentSet(obj)
obj.currentSet = [];
obj.currentSetPosHeader = [];
obj.currentSetPosData = [];
obj.currentSetSize = [];
obj.currentSetDatatype = [];
obj.currentSetDatatypeNumeric = [];
obj.currentSetSizeof = [];
obj.currentSetNumParams = [];
obj.currentSetParams = [];
obj.currentSetNumElements = [];
end
end
end

305
matlab/@ustar/ustar.m Normal file
View File

@ -0,0 +1,305 @@
classdef ustar < handle
% Low-level utilities for UNIX standard tar files.
properties (Access = public)
File % file name
IOMode % file opened in read-only or read-write mode?
NumberOfSubfiles % number of subfiles
end
properties (Access = private)
% File info
fileID
ioflag
subFile
subFileBeg
subFileSize
% Current subfile information
currentFile
currentMode
currentUID
currentGID
currentFileSize
currentModtime
currentLink
currentLinkname
currentUsername
currentGroupname
currentDevMajor
currentDevMinor
currentFileBeg
% Constants
scanBuffSize = 2^17; % buffer size of scanner (max. number of files in tar)
extrBuffSize = 4194304; % buffer size of extracter
blockSize = 512; % ustar block size (do not change)
end
%% ------------------------------------------------------------------------%%
%% CONSTRUCORS/DESTRUCTORS %%
%% ------------------------------------------------------------------------%%
methods(Access=public)
function obj = ustar()
% obj = ucf()
% Default contructor
obj.resetPublicProperties();
obj.resetPrivateProperties();
obj.resetCurrent();
end
function delete(obj)
% obj.delete()
% Default destructor
obj.close();
end
end
%% ------------------------------------------------------------------------%%
%% INITIALIZATION METHODS %%
%% ------------------------------------------------------------------------%%
methods(Access=public)
function open(obj,file)
% obj.open(file)
% Opens a file in read-only mode
obj.File = file;
obj.IOMode = 'read';
obj.ioflag = 'r';
obj.fileID = fopen(obj.File,obj.ioflag);
if obj.fileID<0
error('Unable to open file: %s',obj.File);
end
obj.scanArchive();
obj.resetCurrent();
end
end
%% ------------------------------------------------------------------------%%
%% PUBLIC METHODS %%
%% ------------------------------------------------------------------------%%
methods(Access=public)
function close(obj)
% obj.close()
% Closes a file
if obj.fileID<0
return;
end
status = fclose(obj.fileID);
if status<0
warning('Unable to close file (exit code: %d)',status);
return;
end
obj.resetPublicProperties();
obj.resetPrivateProperties();
obj.resetCurrent();
obj.fileID = -1;
end
function [ptr] = pointer(obj,fname)
% [ptr] = obj.pointer(fname)
% Returns a 'pointer' to the requested file within the tar-ball
% which can be used to read the data without extracting.
% Input
% fname file name of subfile within tar-ball
% Output
% ptr pointer: [fid,first byte,number of bytes]
idx = obj.findSubfile(fname);
ptr = [obj.fileID,obj.subFileBeg(idx),obj.subFileSize(idx)];
end
function [fname,fsize] = list(obj)
% [fname,fsize] = obj.list()
% Returns a list of name/size of all subfiles within the tar-ball
% Output
% fname cell array with filenames
% fsize array with file sizes in bytes
fname = obj.subFile;
fsize = obj.subFileSize;
end
function extract(obj,fname)
% obj.extract(fname)
% Extracts the requested subfile to a standalone file.
% Input
% fname name of subfile
idx = obj.findSubfile(fname);
fbeg = obj.subFileBeg(idx);
fsize = obj.subFileSize(idx);
fidw = fopen(fname,'w');
fseek(obj.fileID,fbeg,'bof');
% Chunk the file
nchunk = ceil(fsize/obj.extrBuffSize);
nchunkFull = floor(fsize/obj.extrBuffSize);
nchunkPart = nchunk-nchunkFull;
for ichunk=1:nchunkFull
buff = fread(obj.fileID,[1,obj.extrBuffSize],'char=>char');
fwrite(fidw,buff,'char');
end
if nchunkPart>0
sizeChunkPart = mod(fsize,obj.extrBuffSize);
buff = fread(obj.fileID,[1,sizeChunkPart],'char=>char');
fwrite(fidw,buff,'char');
end
fclose(fidw);
end
end
%% ------------------------------------------------------------------------%%
%% PRIVATE METHODS %%
%% ------------------------------------------------------------------------%%
methods(Access=private)
function scanArchive(obj)
% obj.scanArchive()
% Scans the tar-ball for subfiles and stores meta-data in class variables.
obj.subFile = cell(obj.scanBuffSize,1);
obj.subFileBeg = zeros(obj.scanBuffSize,1);
obj.subFileSize = zeros(obj.scanBuffSize,1);
% Jump to start of file
fseek(obj.fileID,0,'bof');
% Loop over (unknown) number of subfiles and evaluate header
ii = 0;
while ~obj.checkEOF()
ii = ii+1;
obj.readHeader(true);
obj.subFile{ii} = obj.currentFile;
obj.subFileSize(ii) = obj.currentFileSize;
obj.subFileBeg(ii) = obj.currentFileBeg;
nblock = ceil(obj.currentFileSize/obj.blockSize);
fseek(obj.fileID,nblock*obj.blockSize,'cof');
end
% Truncate preallocated arrays
obj.NumberOfSubfiles = ii;
obj.subFile = obj.subFile(1:ii);
obj.subFileSize = obj.subFileSize(1:ii);
obj.subFileBeg = obj.subFileBeg(1:ii);
if obj.NumberOfSubfiles>obj.scanBuffSize
warning('Number of subfiles exceeds scanBuffSize.');
end
obj.resetCurrent();
end
function readHeader(obj,scanMode)
% obj.readHeader(scanMode)
% Reads header data of a subfile in tar-ball and stores information
% in 'current*' class-variables.
% Input
% scanMode when set to true, omit parts which are not needed during scan
header = fread(obj.fileID,[1,obj.blockSize],'char=>char');
% Extract header information
name = header(1:100);
mode = header(101:108);
uid = header(109:116);
gid = header(117:124);
fsize = header(125:136);
mtime = header(137:148);
chksum = header(149:156);
link = header(157);
linkname = header(158:257);
magic = header(258:263);
version = header(264:265);
uname = header(266:297);
gname = header(298:329);
devmajor = header(330:337);
devminor = header(338:345);
prefix = header(346:500);
% Evaluate checksum
chksum1 = ustar.computeChecksum(header);
chksum2 = ustar.parseOctalStr(chksum);
if chksum1~=chksum2
error('Checksum mismatch! %d,%d',chksum1,chksum2);
end
% Evaluate magic
if ~strcmp(ustar.parseStr(magic),'ustar')
error(' Not a UNIX standard tar file.')
end
% Parse header information
obj.currentFile = ustar.parseStr([prefix,name]);
obj.currentFileBeg = ftell(obj.fileID);
obj.currentFileSize = ustar.parseOctalStr(fsize);
if ~scanMode
obj.currentMode = ustar.parseStr(mode);
obj.currentUID = ustar.parseOctalStr(uid);
obj.currentGID = ustar.parseOctalStr(gid);
obj.currentModtime = datestr(ustar.parseOctalStr(mtime)/86400+datenum(1970,1,1));
obj.currentLink = ustar.parseOctalStr(link);
obj.currentLinkname = ustar.parseStr(linkname);
obj.currentUsername = ustar.parseStr(uname);
obj.currentGroupname = ustar.parseStr(gname);
obj.currentDevMajor = ustar.parseOctalStr(devmajor);
obj.currentDevMinor = ustar.parseOctalStr(devminor);
end
end
function [isEOF] = checkEOF(obj)
% [isEOF] = obj.checkEOF()
% Checks if end-of-file is reached (two blocks of binary zeros).
% Output
% isEOF flag which indicates end-of-file
isEOF = false;
curPosition = ftell(obj.fileID);
blockref = zeros(1,obj.blockSize,'int8');
blockcur = fread(obj.fileID,[1,obj.blockSize],'int8=>int8');
if isequal(blockcur,blockref)
blockcur = fread(obj.fileID,[1,obj.blockSize],'int8=>int8');
if isequal(blockcur,blockref)
isEOF = true;
return;
end
end
fseek(obj.fileID,curPosition,'bof');
end
function [idx] = findSubfile(obj,fname)
% [idx] = obj.findSubfile(fname)
% Get index of requested subfile
% Input
% fname name of subfile
% Output
% idx index of subfile
isReqFile = ismember(obj.subFile,fname);
switch sum(isReqFile)
case 0; error('File not found: %s',fname);
case 1;
otherwise; warning('More than one matching file found.');
end
idx = find(isReqFile);
end
function resetPublicProperties(obj)
obj.File = [];
obj.IOMode = [];
obj.NumberOfSubfiles = [];
end
function resetPrivateProperties(obj)
obj.ioflag = [];
obj.subFile = [];
obj.subFileBeg = [];
obj.subFileSize = [];
end
function resetCurrent(obj)
obj.currentFile = [];
obj.currentMode = [];
obj.currentUID = [];
obj.currentGID = [];
obj.currentFileSize = [];
obj.currentModtime = [];
obj.currentLink = [];
obj.currentLinkname = [];
obj.currentUsername = [];
obj.currentGroupname = [];
obj.currentDevMajor = [];
obj.currentDevMinor = [];
obj.currentFileBeg = [];
end
end
%% ------------------------------------------------------------------------%%
%% PRIVATE STATIC METHODS %%
%% ------------------------------------------------------------------------%%
methods(Access=private,Static)
function [chksum] = computeChecksum(block)
block(149:156) = ' '; % checksum is computed with spaces in check sum field
chksum = sum(block);
end
function [str] = parseStr(str)
charZero = cast(0,'char');
str = strrep(str,charZero,'');
end
function [num] = parseOctalStr(str)
num = ustar.oct2dec_long(str2double(ustar.parseStr(str)));
end
function [dec] = oct2dec_long(oct)
dec = 0;
ii = 1;
while floor(oct/10^(ii-1))~=0
cbase = 8^(ii-1);
cfact = floor(mod(oct,10^ii)/10^(ii-1));
dec = dec + cfact*cbase;
ii = ii+1;
end
end
end
end

58
matlab/colmap_from_flags.m Normal file
View File

@ -0,0 +1,58 @@
function [col] = colmap_from_flags(irank,ihybrid,idem,iscal)
% [col] = colmap_from_flags(irank,ihybrid,idem,iscal)
% Create a containers.Map object with particle column order.
% Input
% irank rank written?
% ihybrid hybrid written?
% idem DEM written?
% iscal scalar written? (number of scalars)
% Output
% col column map which can be indexed by e.g. col('x')
col = containers.Map('KeyType','char','ValueType','double');
ioffset = 0;
if irank
col('rank') = ioffset+1;
ioffset = ioffset+1;
end
if ihybrid
col('id') = ioffset+1;
col('x') = ioffset+2;
col('y') = ioffset+3;
col('z') = ioffset+4;
col('r') = ioffset+5;
col('rho')= ioffset+6;
col('ax') = ioffset+7;
col('ay') = ioffset+8;
col('az') = ioffset+9;
col('u') = ioffset+10;
col('v') = ioffset+11;
col('w') = ioffset+12;
col('ox') = ioffset+13;
col('oy') = ioffset+14;
col('oz') = ioffset+15;
col('fx') = ioffset+16;
col('fy') = ioffset+17;
col('fz') = ioffset+18;
col('tx') = ioffset+19;
col('ty') = ioffset+20;
col('tz') = ioffset+21;
ioffset = ioffset+21;
end
if idem
col('fxc') = ioffset+1;
col('fyc') = ioffset+2;
col('fzc') = ioffset+3;
col('txc') = ioffset+4;
col('tyc') = ioffset+5;
col('tzc') = ioffset+6;
ioffset = ioffset+6;
end
if iscal
for ii=1:iscal
col(['s',sprintf('%d',ii)]) = ioffset+1;
col(['q',sprintf('%d',ii)]) = ioffset+2;
ioffset = ioffset+2;
end
end
end

41
matlab/convert_colmap_matlab2paraview.m Normal file
View File

@ -0,0 +1,41 @@
function [col] = convert_colmap_matlab2paraview(col)
% [col] = convert_colmap_matlab2paraview(col)
% Renames fieldnames of vector quantities, so that paraview
% recognizes them as vectors
colval = cell2mat(col.values);
keysold = col.keys;
keysnew = keysold;
keysnew = regexprep(keysnew,'^id$' ,'ID' );
keysnew = regexprep(keysnew,'^x$' ,'Coords_0' );
keysnew = regexprep(keysnew,'^y$' ,'Coords_1' );
keysnew = regexprep(keysnew,'^z$' ,'Coords_2' );
keysnew = regexprep(keysnew,'^r$' ,'Radius' );
keysnew = regexprep(keysnew,'^rho$' ,'Density' );
keysnew = regexprep(keysnew,'^ax$' ,'AxCoords_0' );
keysnew = regexprep(keysnew,'^ay$' ,'AxCoords_1' );
keysnew = regexprep(keysnew,'^az$' ,'AxCoords_2' );
keysnew = regexprep(keysnew,'^u$' ,'Velocity_0' );
keysnew = regexprep(keysnew,'^v$' ,'Velocity_1' );
keysnew = regexprep(keysnew,'^w$' ,'Velocity_2' );
keysnew = regexprep(keysnew,'^ox$' ,'AxVelocity_0' );
keysnew = regexprep(keysnew,'^oy$' ,'AxVelocity_1' );
keysnew = regexprep(keysnew,'^oz$' ,'AxVelocity_2' );
keysnew = regexprep(keysnew,'^fx$' ,'Force_0' );
keysnew = regexprep(keysnew,'^fy$' ,'Force_1' );
keysnew = regexprep(keysnew,'^fz$' ,'Force_2' );
keysnew = regexprep(keysnew,'^tx$' ,'Torque_0' );
keysnew = regexprep(keysnew,'^ty$' ,'Torque_1' );
keysnew = regexprep(keysnew,'^tz$' ,'Torque_2' );
keysnew = regexprep(keysnew,'^fxc$' ,'CollisionForce_0' );
keysnew = regexprep(keysnew,'^fyc$' ,'CollisionForce_1' );
keysnew = regexprep(keysnew,'^fzc$' ,'CollisionForce_2' );
keysnew = regexprep(keysnew,'^txc$' ,'CollisionTorque_0');
keysnew = regexprep(keysnew,'^tyc$' ,'CollisionTorque_1');
keysnew = regexprep(keysnew,'^tzc$' ,'CollisionTorque_2');
keysnew = regexprep(keysnew,'^s(?=\d*)','ScalarVal' );
keysnew = regexprep(keysnew,'^q(?=\d*)','ScalarFlux' );
col = containers.Map(keysnew,colval);
end

211
matlab/convert_fluid_hdf2ucf.m Normal file
View File

@ -0,0 +1,211 @@
function [] = convert_fluid_hdf2ucf(fhdf,fucf,fproc)
% [] = convert_particles_hdf2ucf(fhdf,fucf)
% Converts a H5Part file containing particle data to a file in UCF format.
% Input
% fhdf file path to HDF file (input file)
% fucf file path to UCF file (output file)
% Get basename for UCF chunks
[fdir,fbase] = fileparts(fucf);
fbasepath = fullfile(fdir,fbase);
% Loop through timesteps
fprintf('[HDF read] %s\n',fhdf);
% Read target proc grid from file
[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(fproc);
nxprocs = numel(ibegu);
nyprocs = numel(jbegu);
nzprocs = numel(kbegu);
nprocs = nxprocs*nyprocs*nzprocs;
% Open HDF
id_hdf = H5F.open(fhdf,'H5F_ACC_RDONLY','H5P_DEFAULT');
% Read simulation time from HDF and create UCF timestep
id_dset = H5D.open(id_hdf,'Time');
time = H5D.read(id_dset);
H5D.close(id_dset);
% Read periodicity for reconstruction of ghost cells
id_group = H5G.open(id_hdf,'/Domain');
id_dset = H5D.open(id_group,'PeriodicityX'); x_periodic = H5D.read(id_dset); H5D.close(id_dset);
id_dset = H5D.open(id_group,'PeriodicityY'); y_periodic = H5D.read(id_dset); H5D.close(id_dset);
id_dset = H5D.open(id_group,'PeriodicityZ'); z_periodic = H5D.read(id_dset); H5D.close(id_dset);
H5G.close(id_group);
% Read full fields from HDF
id_group = H5G.open(id_hdf,'/Fluid');
id_dset = H5D.open(id_group,'VelocityX'); u = H5D.read(id_dset); H5D.close(id_dset);
id_dset = H5D.open(id_group,'VelocityY'); v = H5D.read(id_dset); H5D.close(id_dset);
id_dset = H5D.open(id_group,'VelocityZ'); w = H5D.read(id_dset); H5D.close(id_dset);
id_dset = H5D.open(id_group,'Pressure'); p = H5D.read(id_dset); H5D.close(id_dset);
H5G.close(id_group);
% Close HDF
H5F.close(id_hdf);
clear id_hdf
% Determine global dimensions
nxu = size(u,1); nxv = size(v,1); nxw = size(w,1); nxp = size(p,1);
nyu = size(u,2); nyv = size(v,2); nyw = size(w,2); nyp = size(p,2);
nzu = size(u,3); nzv = size(v,3); nzw = size(w,3); nzp = size(p,3);
if nxp~=iendp(end) || nyp~=jendp(end) || nzp~=kendp(end)
error('Processor grid does not match data. %d,%d,%d <-> %d,%d,%d',nxp,nyp,nzp,iendp(end),jendp(end),kendp(end));
end
for ixproc=0:nxprocs-1
for iyproc=0:nyprocs-1
for izproc=0:nzprocs-1
[iproc] = locfun_proc_id(ixproc,iyproc,izproc,nxprocs,nyprocs,nzprocs);
fname = sprintf('%s.%05d',fbasepath,iproc);
fprintf('[UCF write] %s (out of %d)\n',fname,nprocs);
obj = ucf(fname,'create');
obj.setFileHeader('field','rank',[iproc,ixproc,iyproc,izproc],'endian','a');
obj.appendStep(time);
% Get processor bounadries for this file
ibu = ibegu(ixproc+1); ieu = iendu(ixproc+1);
jbu = jbegu(iyproc+1); jeu = jendu(iyproc+1);
kbu = kbegu(izproc+1); keu = kendu(izproc+1);
ibv = ibegv(ixproc+1); iev = iendv(ixproc+1);
jbv = jbegv(iyproc+1); jev = jendv(iyproc+1);
kbv = kbegv(izproc+1); kev = kendv(izproc+1);
ibw = ibegw(ixproc+1); iew = iendw(ixproc+1);
jbw = jbegw(iyproc+1); jew = jendw(iyproc+1);
kbw = kbegw(izproc+1); kew = kendw(izproc+1);
ibp = ibegp(ixproc+1); iep = iendp(ixproc+1);
jbp = jbegp(iyproc+1); jep = jendp(iyproc+1);
kbp = kbegp(izproc+1); kep = kendp(izproc+1);
% Determine local grid size
nxul = ieu-ibu+1; nxvl = iev-ibv+1; nxwl = iew-ibw+1; nxpl = iep-ibp+1;
nyul = jeu-jbu+1; nyvl = jev-jbv+1; nywl = jew-jbw+1; nypl = jep-jbp+1;
nzul = keu-kbu+1; nzvl = kev-kbv+1; nzwl = kew-kbw+1; nzpl = kep-kbp+1;
% Reconstruct ghost cells
if x_periodic
if ibu==1; ium1=nxu; else; ium1=ibu-1; end
if ibv==1; ivm1=nxv; else; ivm1=ibv-1; end
if ibw==1; iwm1=nxw; else; iwm1=ibw-1; end
if ibp==1; ipm1=nxp; else; ipm1=ibp-1; end
if ieu==nxu; iup1=1; else; iup1=ieu+1; end
if iev==nxv; ivp1=1; else; ivp1=iev+1; end
if iew==nxw; iwp1=1; else; iwp1=iew+1; end
if iep==nxp; ipp1=1; else; ipp1=iep+1; end
else
if ibu==1; ium1=1; else; ium1=ibu-1; end
if ibv==1; ivm1=1; else; ivm1=ibv-1; end
if ibw==1; iwm1=1; else; iwm1=ibw-1; end
if ibp==1; ipm1=1; else; ipm1=ibp-1; end
if ieu==nxu; iup1=nxu; else; iup1=ieu+1; end
if iev==nxv; ivp1=nxv; else; ivp1=iev+1; end
if iew==nxw; iwp1=nxw; else; iwp1=iew+1; end
if iep==nxp; ipp1=nxp; else; ipp1=iep+1; end
end
if y_periodic
if jbu==1; jum1=nyu; else; jum1=jbu-1; end
if jbv==1; jvm1=nyv; else; jvm1=jbv-1; end
if jbw==1; jwm1=nyw; else; jwm1=jbw-1; end
if jbp==1; jpm1=nyp; else; jpm1=jbp-1; end
if jeu==nyu; jup1=1; else; jup1=jeu+1; end
if jev==nyv; jvp1=1; else; jvp1=jev+1; end
if jew==nyw; jwp1=1; else; jwp1=jew+1; end
if jep==nyp; jpp1=1; else; jpp1=jep+1; end
else
if jbu==1; jum1=1; else; jum1=jbu-1; end
if jbv==1; jvm1=1; else; jvm1=jbv-1; end
if jbw==1; jwm1=1; else; jwm1=jbw-1; end
if jbp==1; jpm1=1; else; jpm1=jbp-1; end
if jeu==nyu; jup1=nyu; else; jup1=jeu+1; end
if jev==nyv; jvp1=nyv; else; jvp1=jev+1; end
if jew==nyw; jwp1=nyw; else; jwp1=jew+1; end
if jep==nyp; jpp1=nyp; else; jpp1=jep+1; end
end
if z_periodic
if kbu==1; kum1=nzu; else; kum1=kbu-1; end
if kbv==1; kvm1=nzv; else; kvm1=kbv-1; end
if kbw==1; kwm1=nzw; else; kwm1=kbw-1; end
if kbp==1; kpm1=nzp; else; kpm1=kbp-1; end
if keu==nzu; kup1=1; else; kup1=keu+1; end
if kev==nzv; kvp1=1; else; kvp1=kev+1; end
if kew==nzw; kwp1=1; else; kwp1=kew+1; end
if kep==nzp; kpp1=1; else; kpp1=kep+1; end
else
if kbu==1; kum1=1; else; kum1=kbu-1; end
if kbv==1; kvm1=1; else; kvm1=kbv-1; end
if kbw==1; kwm1=1; else; kwm1=kbw-1; end
if kbp==1; kpm1=1; else; kpm1=kbp-1; end
if keu==nzu; kup1=nzu; else; kup1=keu+1; end
if kev==nzv; kvp1=nzv; else; kvp1=kev+1; end
if kew==nzw; kwp1=nzw; else; kwp1=kew+1; end
if kep==nzp; kpp1=nzp; else; kpp1=kep+1; end
end
% Write u
%fprintf('%d,%d,%d,%d - %d,%d,%d,%d - %d,%d,%d,%d\n',ium1,ibu,ieu,iup1,jum1,jbu,jeu,jup1,kum1,kbu,keu,kup1);
%fprintf('%d,%d,%d,%d - %d,%d,%d,%d - %d,%d,%d,%d\n',ivm1,ibv,iev,ivp1,jvm1,jbv,jev,jvp1,kvm1,kbv,kev,kvp1);
%fprintf('%d,%d,%d,%d - %d,%d,%d,%d - %d,%d,%d,%d\n',iwm1,ibw,iew,iwp1,jwm1,jbw,jew,jwp1,kwm1,kbw,kew,kwp1);
%fprintf('%d,%d,%d,%d - %d,%d,%d,%d - %d,%d,%d,%d\n',ipm1,ibp,iep,ipp1,jpm1,jbp,jep,jpp1,kpm1,kbp,kep,kpp1);
tmp = zeros(nxul+2,nyul+2,nzul+2);
tmp(2:end-1,2:end-1,2:end-1) = u(ibu:ieu,jbu:jeu,kbu:keu);
tmp(1,:,:) = u(ium1,[jum1,jbu:jeu,jup1],[kum1,kbu:keu,kup1]);
tmp(end,:,:) = u(iup1,[jum1,jbu:jeu,jup1],[kum1,kbu:keu,kup1]);
tmp(:,1,:) = u([ium1,ibu:ieu,iup1],jum1,[kum1,kbu:keu,kup1]);
tmp(:,end,:) = u([ium1,ibu:ieu,iup1],jup1,[kum1,kbu:keu,kup1]);
tmp(:,:,1) = u([ium1,ibu:ieu,iup1],[jum1,jbu:jeu,jup1],kum1);
tmp(:,:,end) = u([ium1,ibu:ieu,iup1],[jum1,jbu:jeu,jup1],kup1);
obj.appendField(tmp,1,ibu,jbu,kbu,nxu,nyu,nzu);
% Write v
tmp = zeros(nxvl+2,nyvl+2,nzvl+2);
tmp(2:end-1,2:end-1,2:end-1) = v(ibv:iev,jbv:jev,kbv:kev);
tmp(1,:,:) = v(ivm1,[jvm1,jbv:jev,jvp1],[kvm1,kbv:kev,kvp1]);
tmp(end,:,:) = v(ivp1,[jvm1,jbv:jev,jvp1],[kvm1,kbv:kev,kvp1]);
tmp(:,1,:) = v([ivm1,ibv:iev,ivp1],jvm1,[kvm1,kbv:kev,kvp1]);
tmp(:,end,:) = v([ivm1,ibv:iev,ivp1],jvp1,[kvm1,kbv:kev,kvp1]);
tmp(:,:,1) = v([ivm1,ibv:iev,ivp1],[jvm1,jbv:jev,jvp1],kvm1);
tmp(:,:,end) = v([ivm1,ibv:iev,ivp1],[jvm1,jbv:jev,jvp1],kvp1);
obj.appendField(tmp,1,ibv,jbv,kbv,nxv,nyv,nzv);
% Write w
tmp = zeros(nxwl+2,nywl+2,nzwl+2);
tmp(2:end-1,2:end-1,2:end-1) = w(ibw:iew,jbw:jew,kbw:kew);
tmp(1,:,:) = w(iwm1,[jwm1,jbw:jew,jwp1],[kwm1,kbw:kew,kwp1]);
tmp(end,:,:) = w(iwp1,[jwm1,jbw:jew,jwp1],[kwm1,kbw:kew,kwp1]);
tmp(:,1,:) = w([iwm1,ibw:iew,iwp1],jwm1,[kwm1,kbw:kew,kwp1]);
tmp(:,end,:) = w([iwm1,ibw:iew,iwp1],jwp1,[kwm1,kbw:kew,kwp1]);
tmp(:,:,1) = w([iwm1,ibw:iew,iwp1],[jwm1,jbw:jew,jwp1],kwm1);
tmp(:,:,end) = w([iwm1,ibw:iew,iwp1],[jwm1,jbw:jew,jwp1],kwp1);
obj.appendField(tmp,1,ibw,jbw,kbw,nxw,nyw,nzw);
% Write p
tmp = zeros(nxpl+2,nypl+2,nzpl+2);
tmp(2:end-1,2:end-1,2:end-1) = p(ibp:iep,jbp:jep,kbp:kep);
tmp(1,:,:) = p(ipm1,[jpm1,jbp:jep,jpp1],[kpm1,kbp:kep,kpp1]);
tmp(end,:,:) = p(ipp1,[jpm1,jbp:jep,jpp1],[kpm1,kbp:kep,kpp1]);
tmp(:,1,:) = p([ipm1,ibp:iep,ipp1],jpm1,[kpm1,kbp:kep,kpp1]);
tmp(:,end,:) = p([ipm1,ibp:iep,ipp1],jpp1,[kpm1,kbp:kep,kpp1]);
tmp(:,:,1) = p([ipm1,ibp:iep,ipp1],[jpm1,jbp:jep,jpp1],kpm1);
tmp(:,:,end) = p([ipm1,ibp:iep,ipp1],[jpm1,jbp:jep,jpp1],kpp1);
obj.appendField(tmp,1,ibp,jbp,kbp,nxp,nyp,nzp);
obj.close();
clear obj
end
end
end
end
function [ind] = locfun_proc_id(ii,jj,kk,nx,ny,nz)
% local version of 'sub2ind_zero_row'
ind=ii*ny*nz+jj*nz+kk;
end

22
matlab/convert_particles_array2struct.m Normal file
View File

@ -0,0 +1,22 @@
function [ppstruct] = convert_particles_array2struct(pp,col)
% [ppstruct] = convert_particles_array2struct(pp,col)
% Converts MATLAB particle data from 'array' format to 'struct'
% Input
% pp particle data in 'array' format
% col column map
% Output
% ppstruct particle data in 'struct' format
% Convert map to (cell) arrays and sort by ascending array indices
[colvals,idx] = sort(cell2mat(col.values));
colkeys = col.keys;
colkeys = {colkeys{idx}};
nkeys = col.Count;
nstep = size(pp,3);
for istep=nstep:-1:1
for ikey=1:nkeys
ppstruct(istep).(colkeys{ikey}) = pp(colvals(ikey),:,istep);
end
end
end

22
matlab/convert_particles_cell2struct.m Normal file
View File

@ -0,0 +1,22 @@
function [ppstruct] = convert_particles_cell2struct(pp,col)
% [ppstruct] = convert_particles_cell2struct(pp,col)
% Converts MATLAB particle data from 'cell' format to 'struct'
% Input
% pp particle data in 'cell' format
% col column map
% Output
% ppstruct particle data in 'struct' format
% Convert map to (cell) arrays and sort by ascending array indices
[colvals,idx] = sort(cell2mat(col.values));
colkeys = col.keys;
colkeys = {colkeys{idx}};
nkeys = col.Count;
nstep = numel(pp);
for istep=nstep:-1:1
for ikey=1:nkeys
ppstruct(istep).(colkeys{ikey}) = pp{istep}(colvals(ikey),:);
end
end
end

87
matlab/convert_particles_hdf2ucf.m Normal file
View File

@ -0,0 +1,87 @@
function [] = convert_particles_hdf2ucf(fhdf,fucf)
% [] = convert_particles_hdf2ucf(fhdf,fucf)
% Converts a H5Part file containing particle data to a file in UCF format.
% Input
% fhdf file path to HDF file (input file)
% fucf file path to UCF file (output file)
% Check if file contains DEM data
try
h5info(fhdf,'/Parameters/DEM');
isDEM = true;
catch
isDEM = false;
end
% Create column maps
colml = ucf.partColmapFromFlags(false,true,isDEM,false);
colpv = convert_colmap_matlab2paraview(colml);
% Determine final length of set
ncol = colpv.length;
% Open HDF file ro and UCF file w+
id_file = H5F.open(fhdf,'H5F_ACC_RDONLY','H5P_DEFAULT');
obj = ucf(fucf,'create');
obj.setFileHeader('particle');
% Read number of steps from HDF5 file
id_dset = H5D.open(id_file,'/NumberOfSteps');
nstep = H5D.read(id_dset);
H5D.close(id_dset);
% Loop through timesteps
fprintf('[HDF read] %s\n',fhdf);
fprintf('[UCF write] %s\n',fucf)
for istep=1:nstep
% Open group
groupname = sprintf('Step#%d',istep-1);
id_group = H5G.open(id_file,groupname);
% Read simulation time from HDF and create UCF timestep
id_att = H5A.open(id_group,'TimeValue');
time = H5A.read(id_att);
H5A.close(id_att);
obj.appendStep(time);
% Read number of particles from attribute
id_att = H5A.open(id_group,'NumberOfParticles');
np = H5A.read(id_att);
H5A.close(id_att);
% Prepare array to store current timestep data
id_dset = H5D.open(id_group,'Radius');
id_type = H5D.get_type(id_dset);
type_size = H5T.get_size(id_type);
if type_size==4
pp = zeros(ncol,np,'single');
else
pp = zeros(ncol,np,'double');
end
H5T.close(id_type);
H5D.close(id_dset);
colval = cell2mat(colpv.values);
[colval,sortidx] = sort(colval);
colkey = colpv.keys;
colkey = colkey(sortidx);
% Read data
for icol=1:ncol
id_dset = H5D.open(id_group,colkey{icol});
tmp = H5D.read(id_dset);
H5D.close(id_dset);
pp(colval(icol),:) = tmp;
end
% Add data to UCF file
obj.appendParticle(pp,colml);
% Close current timestep
H5G.close(id_group);
end
% Close files
H5F.close(id_file);
obj.close();
end

22
matlab/flags_from_colmap.m Normal file
View File

@ -0,0 +1,22 @@
function [irank,ihybrid,idem,iscal] = flags_from_colmap(col)
% [irank,ihybrid,idem,iscal] = partFlagsFromColmap(col)
% Extracts flags from containers.Map object
% Input
% col column map which can be indexed by e.g. col('x')
% Output
% irank rank written?
% ihybrid hybrid written?
% idem DEM written?
% iscal scalar written? (number of scalars)
irank = 0;
ihybrid = 0;
idem = 0;
iscal = 0;
if col.isKey('rank'); irank = 1; end
if col.isKey('fx'); ihybrid = 1; end
if col.isKey('fxc'); idem = 1; end
while col.isKey(sprintf('s%d',iscal+1))
iscal = iscal+1;
end
end

26
matlab/ncol_from_colmap.m Normal file
View File

@ -0,0 +1,26 @@
function [ncol,ncol_rank,ncol_hybrid,ncol_dem,ncol_scal] = ncol_from_colmap(col)
% [ncol,ncol_rank,ncol_hybrid,ncol_dem,ncol_scal] = ncol_from_flags(irank,ihybrid,idem,iscal)
% Get number of columns from containers.Map object
% Input
% col column map which can be indexed by e.g. col('x')
% Output
% ncol total number of columns
% ncol_rank number of columns (rank)
% ncol_hybrid number of columns (hybrid)
% ncol_dem number of columns (DEM)
% ncol_scal number of columns (scalar)
ncol_rank = 0;
ncol_hybrid = 0;
ncol_dem = 0;
ncol_scal = 0;
if col.isKey('rank'); ncol_rank = 1; end
if col.isKey('fx'); ncol_hybrid = 21; end
if col.isKey('fxc'); ncol_dem = 6; end
ii = 1;
while col.isKey(sprintf('s%d',ii))
ncol_scal = ncol_scal+2;
ii = ii+1;
end
ncol = ncol_rank+ncol_hybrid+ncol_dem+ncol_scal;
end

25
matlab/ncol_from_flags.m Normal file
View File

@ -0,0 +1,25 @@
function [ncol,ncol_rank,ncol_hybrid,ncol_dem,ncol_scal] = ncol_from_flags(irank,ihybrid,idem,iscal)
% [ncol,ncol_rank,ncol_hybrid,ncol_dem,ncol_scal] = ncol_from_flags(irank,ihybrid,idem,iscal)
% Get number of columns from flags
% Input
% irank rank written?
% ihybrid hybrid written?
% idem DEM written?
% iscal scalar written? (number of scalars)
% Output
% ncol total number of columns
% ncol_rank number of columns (rank)
% ncol_hybrid number of columns (hybrid)
% ncol_dem number of columns (DEM)
% ncol_scal number of columns (scalar)
ncol_rank = 0;
ncol_hybrid = 0;
ncol_dem = 0;
ncol_scal = 0;
if irank; ncol_rank = 1; end
if ihybrid; ncol_hybrid = 21; end
if idem; ncol_dem = 6; end
if iscal; ncol_scal = 2*iscal; end
ncol = ncol_rank+ncol_hybrid+ncol_dem+ncol_scal;
end

57
matlab/read_field_chunk_ucf.m Normal file
View File

@ -0,0 +1,57 @@
function [data,ib,jb,kb,nxl,nyl,nzl,ighost] = read_field_chunk_ucf(file,ifield,varargin)
% [data,ib,jb,kb,nxl,nyl,nzl,ighost] = read_field_chunk_ucf(file,ifield,varargin)
% Reads single field from one processor chunk.
% Input
% file file to be read
% ifield field to be read
% 'uvwp': 1:u, 2:v, 3:w, 4:p
% 'scal': 1:s1, 2:s2, ...
% Optional input (key/value pair)
% ghost keep ghost cells? (default: yes)
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Output
% data partial field with dim(nxl+2*ighost,nyl+2*ighost,nzl+2*ighost)
% ib,jb,kb global index of first grid point of the partial field w/o ghost cells
% nxl,nyl,nzl local field size w/o ghost cells
% ighost ghost cell flag
% Parse optional input arguments
par = inputParser;
addOptional(par,'timestep',1,@isnumeric);
addOptional(par,'ghost',1,@isnumeric);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
istep = par.Results.timestep;
keepghost = par.Results.ghost;
% Open file
obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.open(file);
% Read raw data
if ~obj.validateType('field')
error('read error: no field data.');
end
[data,params] = obj.readSet(istep,ifield);
params = cast(params,'double');
ighost = params(1);
ib = params(2);
jb = params(3);
kb = params(4);
nxl = params(5);
nyl = params(6);
nzl = params(7);
nx = params(8);
ny = params(9);
nz = params(10);
data = reshape(data,nxl+2*ighost,nyl+2*ighost,nzl+2*ighost);
if ighost && ~keepghost
data = data(2:end-1,2:end-1,2:end-1);
ighost = 0;
end
% Close UCF file
obj.close();
end

96
matlab/read_field_complete_ucf.m Normal file
View File

@ -0,0 +1,96 @@
function [q] = read_field_complete_ucf(file,ifield,varargin)
% [data] = read_field_complete_ucf(file,ifield,varargin)
% Reads a specific field from all processor chunks.
% Input
% file file to be read
% ifield field to be read
% 'uvwp': 1:u, 2:v, 3:w, 4:p
% 'scal': 1:s1, 2:s2, ...
% Optional input (key/value pair)
% timestep index of timestep to be read (default: 1)
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Output
% data partial field with dim(nxl+2*ighost,nyl+2*ighost,nzl+2*ighost)
% Parse optional input arguments
par = inputParser;
addOptional(par,'timestep',1,@isnumeric);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
istep = par.Results.timestep;
% Parse filename
[fdir,fbase,fext] = fileparts(file);
fname = sprintf('%s/%s.%05d',fdir,fbase,0);
% Open first file
obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.open(fname);
if ~obj.validateType('field')
error('read error: no field data.');
end
% Read raw data
[data,params] = obj.readSet(istep,ifield);
params = cast(params,'double');
ighost = params(1);
ib = params(2);
jb = params(3);
kb = params(4);
nxl = params(5);
nyl = params(6);
nzl = params(7);
nx = params(8);
ny = params(9);
nz = params(10);
data = reshape(data,nxl+2*ighost,nyl+2*ighost,nzl+2*ighost);
if ighost
data = data(2:end-1,2:end-1,2:end-1);
end
% Create global array and assign chunk
q = zeros(nx,ny,nz,class(data));
q(ib:ib+nxl-1,jb:jb+nyl-1,kb:kb+nzl-1) = data;
% Close the first file
obj.close();
% Now loop consecutively through files
ifile = 1;
fname = sprintf('%s/%s.%05d',fdir,fbase,ifile);
while exist(fname,'file')
% Open file
obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.open(fname);
if ~obj.validateType('field')
error('read error: no field data.');
end
% Read raw data
[data,params] = obj.readSet(istep,ifield);
params = cast(params,'double');
ighost = params(1);
ib = params(2);
jb = params(3);
kb = params(4);
nxl = params(5);
nyl = params(6);
nzl = params(7);
data = reshape(data,nxl+2*ighost,nyl+2*ighost,nzl+2*ighost);
if ighost
data = data(2:end-1,2:end-1,2:end-1);
end
% Assign chunk
q(ib:ib+nxl-1,jb:jb+nyl-1,kb:kb+nzl-1) = data;
% Close the first file
obj.close();
% Move on to next file
ifile = ifile+1;
fname = sprintf('%s/%s.%05d',fdir,fbase,ifile);
end
end

67
matlab/read_grid_ucf.m Normal file
View File

@ -0,0 +1,67 @@
function [xu,yu,zu,xv,yv,zv,xw,yw,zw,xp,yp,zp] = 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
% Optional input (key/value pair)
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Output
% 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
% Parse optional input arguments
par = inputParser;
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
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);
obj.open(file);
% Read raw data
x = cell(1,nset);
y = cell(1,nset);
z = cell(1,nset);
for iset=1:nset
[data,params] = obj.readSet(1,iset);
params = cast(params,'double');
nx = params(1);
ny = params(2);
nz = params(3);
x{iset} = data(1:nx);
y{iset} = data(nx+1:nx+ny);
z{iset} = data(nx+ny+1:nx+ny+nz);
end
% Close UCF file
obj.close();
% Reassign content from loop
iset = find(strcmp(sets,'u'));
xu = x{iset};
yu = y{iset};
zu = z{iset};
iset = find(strcmp(sets,'v'));
xv = x{iset};
yv = y{iset};
zv = z{iset};
iset = find(strcmp(sets,'w'));
xw = x{iset};
yw = y{iset};
zw = z{iset};
iset = find(strcmp(sets,'p'));
xp = x{iset};
yp = y{iset};
zp = z{iset};
end

21
matlab/read_parameters_ucf.m Normal file
View File

@ -0,0 +1,21 @@
function [params] = read_parameters_ucf(file,varargin)
% [params] = read_parameters_ucf(file)
% Reads the entire parameters_XXXX.asc as written by UCF file format
% into a struct.
% Input
% file file name
% Optional input
% verbosity verbose screen output? (default: no)
% Output
% params struct with content of parameters_XXXX.asc
% Parse optional input arguments
par = inputParser;
addOptional(par,'verbosity',0,@isnumeric);
parse(par,varargin{:});
obj = ini('verbosity',par.Results.verbosity);
obj.open(file);
params = obj.getContent();
obj.close();
end

47
matlab/read_particle_balancing_ucf.m Normal file
View File

@ -0,0 +1,47 @@
function [iproc,nmaster,nslave,stime] = read_particle_balancing_ucf(file,varargin)
% [pp,col] = read_particles_ucf(file,varargin)
% Reads particle data from UCF file
% Input
% file file to be read
% Optional input (key/value pair)
% timestep timestep to be read (default: all)
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Output
% pp particle balancing data
% stime simulation time with dim(1,nt)
% Parse optional input arguments
par = inputParser;
addOptional(par,'timestep',1,@isnumeric);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
istep = par.Results.timestep;
% Open file
obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.open(file);
% Get range of steps to be read
if istep>obj.getNumTimesteps()
error('Timestep out of range.');
end
% Get simulation time
stime = obj.getSimulationTime(istep);
% Read data
[data,params] = obj.readSet(istep,1);
params = cast(params,'double');
nxprocs = params(1);
nyprocs = params(2);
nzprocs = params(3);
data = cast(reshape(data,3,nxprocs*nyprocs*nzprocs),'double');
iproc = data(1,:);
nmaster = data(2,:);
nslave = data(3,:);
% Close file
obj.close();
end

84
matlab/read_particles_ucf.m Normal file
View File

@ -0,0 +1,84 @@
function [pp,col,stime] = read_particles_ucf(file,varargin)
% [pp,col] = read_particles_ucf(file,varargin)
% Reads particle data from UCF file
% Input
% file file to be read
% Optional input (key/value pair)
% timestep timestep to be read (default: all)
% format MATLAB format (default: 'array')
% 'array' plain array with dim(ncol.np,nt)
% 'struct' structure array with short fieldnames
% 'paraview' structure array with H5Part fieldnames
% 'cell' cell array with dim(1,nt) and plain arrays with dim(ncol,np) inside
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Output
% pp particle data in specified format
% col container map 'char'->'double' which maps quantity names to column index
% stime simulation time with dim(1,nt)
% Parse optional input arguments
par = inputParser;
addOptional(par,'timestep',0,@isnumeric);
addOptional(par,'format','array',@ischar);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
istep = par.Results.timestep;
mlformat = par.Results.format;
% Open file
obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.open(file);
if ~obj.validateType('particle')
error('read error: no particle data.');
end
% Get range of steps to be read
if istep==0
steps = uint32(1:obj.getNumTimesteps());
else
steps = uint32(istep);
end
nstepout = numel(steps);
% Get simulation time
stime = obj.getSimulationTime();
stime = stime(steps);
% Read raw data
for istep=nstepout:-1:1
[data,params] = obj.readSet(steps(istep),1);
params = cast(params,'double');
np = params(1);
ncol = params(2);
ncol_rank = params(3);
ncol_hybrid = params(4);
ncol_dem = params(5);
ncol_scalar = params(6);
data = reshape(data,ncol,np);
nscal = ncol_scalar/2;
col = colmap_from_flags(ncol_rank,ncol_hybrid,ncol_dem,nscal);
% Parse data to specified format
switch mlformat
case 'struct'
pp(istep) = convert_particles_array2struct(data,col);
case 'paraview'
col = convert_colmap_matlab2paraview(col);
pp(istep) = convert_particles_array2struct(data,col);
case {'cell','array'}
pp{istep} = data;
otherwise
error('Unknown format: %s',mlformat);
end
end
% Convert cell array to 3-dim array, if mlformat=='array'
if strcmp(mlformat,'array')
pp = cat(3,pp{:});
end
% Close file
obj.close();
end

91
matlab/read_procgrid_ucf.m Normal file
View File

@ -0,0 +1,91 @@
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] = 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] = read_procgrid_ucf(file,varargin)
% Reads processor grids.
% Input
% file file to be read
% Optional input (key/value pair)
% verbosity verbose output? (default: no)
%
% Output
% ibegu,iendu,jbegu,jendu,kbegu,kendu processor grid u
% 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
% Parse optional input arguments
par = inputParser;
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',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);
obj.open(file);
% Read raw data
ibeg = cell(1,nset);
iend = cell(1,nset);
jbeg = cell(1,nset);
jend = cell(1,nset);
kbeg = cell(1,nset);
kend = cell(1,nset);
for iset=1:nset
[data,params] = obj.readSet(1,iset);
params = cast(params,'double');
nxprocs = params(1);
nyprocs = params(2);
nzprocs = params(3);
ibeg{iset} = data(1:nxprocs);
iend{iset} = data(nxprocs+1:2*nxprocs);
jbeg{iset} = data(2*nxprocs+1:2*nxprocs+nyprocs);
jend{iset} = data(2*nxprocs+nyprocs+1:2*nxprocs+2*nyprocs);
kbeg{iset} = data(2*nxprocs+2*nyprocs+1:2*nxprocs+2*nyprocs+nzprocs);
kend{iset} = data(2*nxprocs+2*nyprocs+nzprocs+1:2*nxprocs+2*nyprocs+2*nzprocs);
end
% Close UCF file
obj.close();
% Reassign content from loop
iset = find(strcmp(sets,'u'));
ibegu = ibeg{iset};
iendu = iend{iset};
jbegu = jbeg{iset};
jendu = jend{iset};
kbegu = kbeg{iset};
kendu = kend{iset};
iset = find(strcmp(sets,'v'));
ibegv = ibeg{iset};
iendv = iend{iset};
jbegv = jbeg{iset};
jendv = jend{iset};
kbegv = kbeg{iset};
kendv = kend{iset};
iset = find(strcmp(sets,'w'));
ibegw = ibeg{iset};
iendw = iend{iset};
jbegw = jbeg{iset};
jendw = jend{iset};
kbegw = kbeg{iset};
kendw = kend{iset};
iset = find(strcmp(sets,'p'));
ibegp = ibeg{iset};
iendp = iend{iset};
jbegp = jbeg{iset};
jendp = jend{iset};
kbegp = kbeg{iset};
kendp = kend{iset};
end

59
matlab/read_scal_chunk_ucf.m Normal file
View File

@ -0,0 +1,59 @@
function [s,ib,jb,kb,nxl,nyl,nzl,ighost] = read_scal_chunk_ucf(file,varargin)
% [s,ib,jb,kb,nxl,nyl,nzl,ighost] = read_scal_chunk_ucf(file,varargin)
% Reads scalar fields from one processor chunk.
% Input
% file file to be read
% Optional input (key/value pair)
% timestep index of timestep to be read (default: 1)
% ghost keep ghost cells? (default: yes)
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Output
% s partial fields with dim(nxl+2*ighost,nyl+2*ighost,nzl+2*ighost,nscal)
% ib,jb,kb,... global index of first grid point of the partial field w/o ghost cells
% nxl,nyl,nzl,... local field size w/o ghost cells
% ighost ghost cell flag
% Parse optional input arguments
par = inputParser;
addOptional(par,'timestep',1,@isnumeric);
addOptional(par,'ghost',1,@isnumeric);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
istep = par.Results.timestep;
keepghost = par.Results.ghost;
% Open file
obj = ucf('verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.open(file);
if ~obj.validateType('field')
error('read error: no field data.');
end
nscal = obj.getNumDatasets(istep);
% Read raw data
s = cell(1,nscal);
for iset=1:nscal
[data,params] = obj.readSet(istep,iset);
params = cast(params,'double');
ighost = params(1);
ib = params(2);
jb = params(3);
kb = params(4);
nxl = params(5);
nyl = params(6);
nzl = params(7);
s{iset} = reshape(data,nxl+2*ighost,nyl+2*ighost,nzl+2*ighost);
if ighost && ~keepghost
s{iset} = s{iset}(2:end-1,2:end-1,2:end-1);
ighost = 0;
end
end
% Close UCF file
obj.close();
% Convert cell array to 4-D array
s = cat(4,s{:});
end

102
matlab/read_scal_complete_ucf.m Normal file
View File

@ -0,0 +1,102 @@
function [s] = read_scal_complete_ucf(file,varargin)
% [s] = read_scal_complete_ucf(file,varargin)
% Reads u,v,w,p from all processor chunks and combines them to a single field.
% Input
% file file to be read
% Optional input (key/value pair)
% timestep index of timestep to be read (default: 1)
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Output
% s complete scalar fields with dim(nx,ny,nz,nscal)
% Parse optional input arguments
par = inputParser;
addOptional(par,'timestep',1,@isnumeric);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
istep = par.Results.timestep;
% Parse filename
[fdir,fbase,fext] = fileparts(file);
fname = sprintf('%s/%s.%05d',fdir,fbase,0);
% Open first file
obj = ucf(fname,'read','verbosity',par.Results.verbosity,'debug',par.Results.debug);
if ~obj.validateType('field')
error('read error: no field data.');
end
nscal = obj.getNumDatasets(istep);
% Read raw data
q = cell(1,nscal);
for iset=1:nscal
[data,params] = obj.readSet(istep,iset);
params = cast(params,'double');
ighost = params(1);
ib = params(2);
jb = params(3);
kb = params(4);
nxl = params(5);
nyl = params(6);
nzl = params(7);
nx = params(8);
ny = params(9);
nz = params(10);
q{iset} = reshape(data,nxl+2*ighost,nyl+2*ighost,nzl+2*ighost);
if ighost
q{iset} = q{iset}(2:end-1,2:end-1,2:end-1);
end
end
% 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};
end
% Close the first file
obj.close();
% Now loop consecutively through files
ifile = 1;
fname = sprintf('%s/%s.%05d',fdir,fbase,ifile);
while exist(fname,'file')
% Open file
obj = ucf(fname,'read','verbosity',par.Results.verbosity,'debug',par.Results.debug);
if ~obj.validateType('field')
error('read error: no field data.');
end
% Read raw data
q = cell(1,nscal);
for iset=1:nscal
[data,params] = obj.readSet(istep,iset);
params = cast(params,'double');
ighost = params(1);
ib = params(2);
jb = params(3);
kb = params(4);
nxl = params(5);
nyl = params(6);
nzl = params(7);
q{iset} = reshape(data,nxl+2*ighost,nyl+2*ighost,nzl+2*ighost);
if ighost
q{iset} = q{iset}(2:end-1,2:end-1,2:end-1);
end
end
% 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};
end
% Close file
obj.close();
% Move on to next file
ifile = ifile+1;
fname = sprintf('%s/%s.%05d',fdir,fbase,ifile);
end
end

40
matlab/read_statistics_channel_scal_ucf.m Normal file
View File

@ -0,0 +1,40 @@
function [tbeg,tend,nstat,sm,ss,su,sv] = read_statistics_channel_scal_ucf(file,varargin)
% Parse optional input arguments
par = inputParser;
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
% Define sets to be read
sets = {'um','uu','vv','ww','uv','uub','uvb'};
nset = numel(sets);
% Open UCF file and read
obj = ucf(file,'read','verbosity',par.Results.verbosity,'debug',par.Results.debug);
tend = obj.getSimulationTime(1);
% Read parameters of first set
params = obj.readParameters(1,1);
% Parse parameters
tbeg = typecast(params(1),'double');
nstat = cast(params(2),'double');
ny = cast(params(3),'double');
nscal = cast(params(4),'double'); % number of scal
nset = cast(params(5),'double'); % number of stats
sm = zeros(ny,nscal);
ss = zeros(ny,nscal);
su = zeros(ny,nscal);
sv = zeros(ny,nscal);
for iscal=1:nscal
sm(:,iscal) = obj.readSet(1,(iscal-1)*4+1)/nstat;
ss(:,iscal) = obj.readSet(1,(iscal-1)*4+2)/nstat;
su(:,iscal) = obj.readSet(1,(iscal-1)*4+3)/nstat;
sv(:,iscal) = obj.readSet(1,(iscal-1)*4+4)/nstat;
end
% Close file
obj.close();
end

52
matlab/read_statistics_channel_ucf.m Normal file
View File

@ -0,0 +1,52 @@
function [tbeg,tend,nstat,um,uu,vv,ww,uv,uub,uvb,varargout] = read_statistics_channel_ucf(file,varargin)
% Parse optional input arguments
par = inputParser;
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
% Define sets to be read
sets = {'um','uu','vv','ww','uv','uub','uvb'};
nset = numel(sets);
% Open UCF file and read
obj = ucf(file,'read','verbosity',par.Results.verbosity,'debug',par.Results.debug);
tend = obj.getSimulationTime(1);
% Read first dataset
[um,params] = obj.readSet(1,1);
% Parse parameters
tbeg = typecast(params(1),'double');
nstat = cast(params(2),'double');
nset1 = cast(params(4),'double'); % channel stats
nset2 = cast(params(5),'double'); % force eulerian
if nset1~=7
error('Invalid file.');
end
um = um/nstat;
uu = obj.readSet(1,2)/nstat;
vv = obj.readSet(1,3)/nstat;
ww = obj.readSet(1,4)/nstat;
uv = obj.readSet(1,5)/nstat;
uub = obj.readSet(1,6)/nstat;
uvb = obj.readSet(1,7)/nstat;
% If eulerian force data, read that too
if nset2==3
if par.Results.verbosity
fprintf('Eulerian force data found.\n');
end
fxp = obj.readSet(1,8);
fyp = obj.readSet(1,9);
fzp = obj.readSet(1,10);
varargout{1} = fxp;
varargout{2} = fyp;
varargout{3} = fzp;
end
% Close file
obj.close();
end

107
matlab/read_uvwp_chunk_ucf.m Normal file
View File

@ -0,0 +1,107 @@
function [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(file,varargin)
% [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(file,varargin)
% Reads u,v,w,p from one processor chunk.
% Input
% file file to be read
% Optional input (key/value pair)
% timestep index of timestep to be read (default: 1)
% ghost keep ghost cells? (default: yes)
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Output
% u,v,w,p partial fields with dim(nxl,nyl,nzl)+2*ighost
% ibu,jbu,kbu,... global index of first grid point of the partial field w/o ghost cells
% nxul,nyul,nzul,... local field size w/o ghost cells
% ighost ghost cell flag
% Parse optional input arguments
par = inputParser;
addOptional(par,'timestep',1,@isnumeric);
addOptional(par,'ghost',1,@isnumeric);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
istep = par.Results.timestep;
keepghost = par.Results.ghost;
% Define sets to be read
sets = {'u','v','w','p'};
nset = numel(sets);
% Open file
obj = ucf(file,'read','verbosity',par.Results.verbosity,'debug',par.Results.debug);
if ~obj.validateType('field')
error('read error: no field data.');
end
% Read raw data
ib = cell(1,nset);
jb = cell(1,nset);
kb = cell(1,nset);
nxl = cell(1,nset);
nyl = cell(1,nset);
nzl = cell(1,nset);
q = cell(1,nset);
for iset=1:nset
[data,params] = obj.readSet(istep,iset);
params = cast(params,'double');
ighost = params(1);
ib{iset} = params(2);
jb{iset} = params(3);
kb{iset} = params(4);
nxl{iset} = params(5);
nyl{iset} = params(6);
nzl{iset} = params(7);
q{iset} = reshape(data,nxl{iset}+2*ighost,nyl{iset}+2*ighost,nzl{iset}+2*ighost);
if ighost && ~keepghost
q{iset} = q{iset}(2:end-1,2:end-1,2:end-1);
ighost = 0;
end
end
% Close UCF file
obj.close();
% Reassign content from loop
iset = find(strcmp(sets,'u'));
u = q{iset};
ibu = ib{iset};
jbu = jb{iset};
kbu = kb{iset};
nxul = nxl{iset};
nyul = nyl{iset};
nzul = nzl{iset};
iset = find(strcmp(sets,'v'));
v = q{iset};
ibv = ib{iset};
jbv = jb{iset};
kbv = kb{iset};
nxvl = nxl{iset};
nyvl = nyl{iset};
nzvl = nzl{iset};
iset = find(strcmp(sets,'w'));
w = q{iset};
ibw = ib{iset};
jbw = jb{iset};
kbw = kb{iset};
nxwl = nxl{iset};
nywl = nyl{iset};
nzwl = nzl{iset};
iset = find(strcmp(sets,'p'));
p = q{iset};
ibp = ib{iset};
jbp = jb{iset};
kbp = kb{iset};
nxpl = nxl{iset};
nypl = nyl{iset};
nzpl = nzl{iset};
end

155
matlab/read_uvwp_complete_ucf.m Normal file
View File

@ -0,0 +1,155 @@
function [u,v,w,p] = read_uvwp_complete_ucf(file,varargin)
% [u,v,w,p] = read_uvwp_complete_ucf(file,varargin)
% Reads u,v,w,p from all processor chunks and combines them to a single field.
% Input
% file file to be read
% Optional input (key/value pair)
% timestep index of timestep to be read (default: 1)
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Output
% u,v,w,p complete fields with dim(nx,ny,nz)
% Parse optional input arguments
par = inputParser;
addOptional(par,'timestep',1,@isnumeric);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
istep = par.Results.timestep;
% Define sets to be read
sets = {'u','v','w','p'};
nset = numel(sets);
% Parse filename
[fdir,fbase,fext] = fileparts(file);
fname = sprintf('%s/%s.%05d',fdir,fbase,0);
% Open first file
obj = ucf(fname,'read','verbosity',par.Results.verbosity,'debug',par.Results.debug);
if ~obj.validateType('field')
error('read error: no field data.');
end
% Read raw data
q = cell(1,nset);
ib = cell(1,nset);
jb = cell(1,nset);
kb = cell(1,nset);
nxl = cell(1,nset);
nyl = cell(1,nset);
nzl = cell(1,nset);
nx = cell(1,nset);
ny = cell(1,nset);
nz = cell(1,nset);
for iset=1:nset
[data,params] = obj.readSet(istep,iset);
params = cast(params,'double');
ighost = params(1);
ib{iset} = params(2);
jb{iset} = params(3);
kb{iset} = params(4);
nxl{iset}= params(5);
nyl{iset}= params(6);
nzl{iset}= params(7);
nx{iset} = params(8);
ny{iset} = params(9);
nz{iset} = params(10);
q{iset} = reshape(data,nxl{iset}+2*ighost,nyl{iset}+2*ighost,nzl{iset}+2*ighost);
if ighost
q{iset} = q{iset}(2:end-1,2:end-1,2:end-1);
end
end
% Reassign content from loop and create global arrays
iset = find(strcmp(sets,'u'));
u = zeros(nx{iset},ny{iset},nz{iset},class(q{iset}));
u(ib{iset}:ib{iset}+nxl{iset}-1,...
jb{iset}:jb{iset}+nyl{iset}-1,...
kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};
iset = find(strcmp(sets,'v'));
v = zeros(nx{iset},ny{iset},nz{iset},class(q{iset}));
v(ib{iset}:ib{iset}+nxl{iset}-1,...
jb{iset}:jb{iset}+nyl{iset}-1,...
kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};
iset = find(strcmp(sets,'w'));
w = zeros(nx{iset},ny{iset},nz{iset},class(q{iset}));
w(ib{iset}:ib{iset}+nxl{iset}-1,...
jb{iset}:jb{iset}+nyl{iset}-1,...
kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};
iset = find(strcmp(sets,'p'));
p = zeros(nx{iset},ny{iset},nz{iset},class(q{iset}));
p(ib{iset}:ib{iset}+nxl{iset}-1,...
jb{iset}:jb{iset}+nyl{iset}-1,...
kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};
% Close the first file
obj.close();
% Now loop consecutively through files
ifile = 1;
fname = sprintf('%s/%s.%05d',fdir,fbase,ifile);
while exist(fname,'file')
% Open file
obj = ucf(fname,'read','verbosity',par.Results.verbosity,'debug',par.Results.debug);
if ~obj.validateType('field')
error('read error: no field data.');
end
% Read raw data
q = cell(1,nset);
ib = cell(1,nset);
jb = cell(1,nset);
kb = cell(1,nset);
nxl = cell(1,nset);
nyl = cell(1,nset);
nzl = cell(1,nset);
for iset=1:nset
[data,params] = obj.readSet(istep,iset);
params = cast(params,'double');
ighost = params(1);
ib{iset} = params(2);
jb{iset} = params(3);
kb{iset} = params(4);
nxl{iset}= params(5);
nyl{iset}= params(6);
nzl{iset}= params(7);
q{iset} = reshape(data,nxl{iset}+2*ighost,nyl{iset}+2*ighost,nzl{iset}+2*ighost);
if ighost
q{iset} = q{iset}(2:end-1,2:end-1,2:end-1);
end
end
% Reassign content from loop and create global arrays
iset = find(strcmp(sets,'u'));
u(ib{iset}:ib{iset}+nxl{iset}-1,...
jb{iset}:jb{iset}+nyl{iset}-1,...
kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};
iset = find(strcmp(sets,'v'));
v(ib{iset}:ib{iset}+nxl{iset}-1,...
jb{iset}:jb{iset}+nyl{iset}-1,...
kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};
iset = find(strcmp(sets,'w'));
w(ib{iset}:ib{iset}+nxl{iset}-1,...
jb{iset}:jb{iset}+nyl{iset}-1,...
kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};
iset = find(strcmp(sets,'p'));
p(ib{iset}:ib{iset}+nxl{iset}-1,...
jb{iset}:jb{iset}+nyl{iset}-1,...
kb{iset}:kb{iset}+nzl{iset}-1) = q{iset};
% Close file
obj.close();
% Move on to next file
ifile = ifile+1;
fname = sprintf('%s/%s.%05d',fdir,fbase,ifile);
end
end

48
matlab/write_particles_ucf.m Normal file
View File

@ -0,0 +1,48 @@
function [] = write_particles_ucf(file,pp,col,time,varargin)
% [] = write_particles_ucf(file,pp,col,time,varargin)
% Writes particles data into a UCF file.
% Input
% file file to be written (the extension will be replaced by the proc rank)
% pp particle data in 'array' format with dim(ncol,np,nt)
% col particle column map. CAUTION: data is not rearranged according to colmap yet!!!
% time simulation time with dim(nt)
% Optional input (key/value pair)
% endian endianess of the file as used by fopen (default: 'n')
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Parse variable input
par = inputParser;
addOptional(par,'endian','n',@ischar);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
ncol = size(pp,1);
np = size(pp,2);
[ncolref,ncol_rank,ncol_hybrid,ncol_dem,ncol_scal] = ncol_from_colmap(col);
nt = size(pp,3);
% Check for consistency
if ncol~=col.Count || ncol~=ncolref
error('Particle data does not match column map.');
end
if nt~=numel(time)
error('Number of time steps must match the number of time values given. %d,%d',nt,numel(time));
end
% Create file
obj = ucf(file,'create','verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.setFileHeader('particle','endian',par.Results.endian);
% Add particle data step by step
for it=1:nt
obj.appendStep(time(it));
params = int64([np,ncol,ncol_rank,ncol_hybrid,ncol_dem,ncol_scal]);
data = pp(:,:,it);
obj.appendSet(data,params);
end
% Close file
obj.close();
end

44
matlab/write_scal_chunk_ucf.m Normal file
View File

@ -0,0 +1,44 @@
function [] = write_scal_chunk_ucf(file,s,ib,jb,kb,nx,ny,nz,time,ighost,varargin)
% [] = write_scal_chunk_ucf(file,s,ibp,jbp,kbp,nxp,nyp,nzp,time,ighost,varargin)
% Writes a partial uvwp field to UCF chunk.
% Input
% file file to be written (the extension will be replaced by the proc rank)
% s partial scalar fields with dim(nx,ny,nz,nscal)
% ib,jb,kb global index of first data point of field w/o ghost cells
% nx,ny,nz global grid size
% time simulation time
% ighost does partial field contain ghost cells?
% Optional input (key/value pair)
% rank processor rank as array [myid,my_col,my_row,my_pln]
% endian endianess of the file as used by fopen (default: 'n')
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Parse variable input
par = inputParser;
addOptional(par,'rank',[0,0,0,0],@(x)isnumeric(x)&&numel(x)==4);
addOptional(par,'endian','n',@ischar);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
% Determine number of scalar fields
nscal = size(s,4);
% Create file and write to it
obj = ucf(file,'create','verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.setFileHeader('field','rank',par.Results.rank,'endian',par.Results.endian);
obj.appendStep(time);
% Write scalar field by field
for iscal=1:nscal
nxl = size(s,1)-2*ighost;
nyl = size(s,2)-2*ighost;
nzl = size(s,3)-2*ighost;
params = int64([ighost,ib,jb,kb,nxl,nyl,nzl,nx,ny,nz]);
obj.appendSet(s(:,:,:,iscal),params);
end
% Close file
obj.close();
end

75
matlab/write_scal_complete_ucf.m Normal file
View File

@ -0,0 +1,75 @@
function [] = write_scal_complete_ucf(file,s,ibegp,iendp,jbegp,jendp,kbegp,kendp,time,varargin)
% [] = write_scal_complete_ucf(file,s,ibegp,iendp,jbegp,jendp,kbegp,kendp,time,varargin)
% Writes a complete scalar field to UCF chunks.
% Input
% file file to be written (the extension will be replaced by the proc rank)
% s full scalar field with dim(nx,ny,nz,nscal)
% ibegp,jbegp,kbegp arrays with length(n.proc) with processor grid
% iendp,jendp,kendp """
% time simulation time
% Optional input (key/value pair)
% endian endianess of the file as used by fopen (default: 'n')
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Parse variable input
par = inputParser;
addOptional(par,'endian','n',@ischar);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
% Parse filename
[fdir,fbase,fext] = fileparts(file);
fbasepath = fullfile(fdir,fbase);
nxprocs = numel(ibegp);
nyprocs = numel(jbegp);
nzprocs = numel(kbegp);
nx = size(s,1);
ny = size(s,2);
nz = size(s,3);
nscal = size(s,4);
for ixproc=0:nxprocs-1
for iyproc=0:nyprocs-1
for izproc=0:nzprocs-1
[iproc] = locfun_proc_id(ixproc,iyproc,izproc,nxprocs,nyprocs,nzprocs);
fname = sprintf('%s.%05d',fbasepath,iproc);
% Extract local indices
ib = ibegp(ixproc+1);
ie = iendp(ixproc+1);
jb = jbegp(iyproc+1);
je = jendp(iyproc+1);
kb = kbegp(izproc+1);
ke = kendp(izproc+1);
% Create a new chunk file
obj = ucf(fname,'create','verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.setFileHeader('field','rank',[iproc,ixproc,iyproc,izproc],'endian',par.Results.endian);
obj.appendStep(time);
% Write data field by field
for iscal=1:nscal
ighost = 0;
nxl = size(s,1)-2*ighost;
nyl = size(s,2)-2*ighost;
nzl = size(s,3)-2*ighost;
params = int64([ighost,ib,jb,kb,nxl,nyl,nzl,nx,ny,nz]);
data = s(ib:ie,jb:je,kb:ke,iscal);
obj.appendSet(data,params);
end
% Close chunk file
obj.close();
end
end
end
end
function [ind] = locfun_proc_id(ii,jj,kk,nx,ny,nz)
% local version of 'sub2ind_zero_row'
ind=ii*ny*nz+jj*nz+kk;
end

70
matlab/write_uvwp_chunk_ucf.m Normal file
View File

@ -0,0 +1,70 @@
function [] = write_uvwp_chunk_ucf(file,...
u,ibu,jbu,kbu,nxu,nyu,nzu,...
v,ibv,jbv,kbv,nxv,nyv,nzv,...
w,ibw,jbw,kbw,nxw,nyw,nzw,...
p,ibp,jbp,kbp,nxp,nyp,nzp,...
time,ighost,varargin)
% [] = write_uvwp_chunk_ucf(file,...
% u,ibu,jbu,kbu,nxu,nyu,nzu,,...
% v,ibv,jbv,kbv,nxv,nyv,nzv,,...
% w,ibw,jbw,kbw,nxw,nyw,nzw,,...
% p,ibp,jbp,kbp,nxp,nyp,nzp,,...
% time,ighost,varargin)
% Writes a partial uvwp field to UCF chunk.
% Input
% file file to be written (the extension will be replaced by the proc rank)
% u,v,w,p partial u,v,w,p fields with dim(nx,ny,nz)
% ib,jb,kb global index of first data point of field w/o ghost cells
% nx,ny,nz global grid size
% time simulation time
% ighost does partial field contain ghost cells?
% Optional input (key/value pair)
% rank processor rank as array [myid,my_col,my_row,my_pln]
% endian endianess of the file as used by fopen (default: 'n')
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Parse variable input
par = inputParser;
addOptional(par,'rank',[0,0,0,0],@(x)isnumeric(x)&&numel(x)==4);
addOptional(par,'endian','n',@ischar);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
% Create file and write to it
obj = ucf(file,'create','verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.setFileHeader('field','rank',par.Results.rank,'endian',par.Results.endian);
obj.appendStep(time);
% 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]);
obj.appendSet(u,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]);
obj.appendSet(v,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]);
obj.appendSet(w,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]);
obj.appendSet(p,params);
% Close file
obj.close();
end

106
matlab/write_uvwp_complete_ucf.m Normal file
View File

@ -0,0 +1,106 @@
function [] = write_uvwp_complete_ucf(file,...
u,ibegu,jbegu,kbegu,iendu,jendu,kendu,...
v,ibegv,jbegv,kbegv,iendv,jendv,kendv,...
w,ibegw,jbegw,kbegw,iendw,jendw,kendw,...
p,ibegp,jbegp,kbegp,iendp,jendp,kendp,...
time,varargin)
% [] = write_uvwp_complete_ucf(file,...
% u,ibegu,jbegu,kbegu,iendu,jendu,kendu,...
% v,ibegv,jbegv,kbegv,iendv,jendv,kendv,...
% w,ibegw,jbegw,kbegw,iendw,jendw,kendw,...
% p,ibegp,jbegp,kbegp,iendp,jendp,kendp,...
% time,varargin)
% Writes a complete uvwp field to UCF chunks.
% Input
% file file to be written (the extension will be replaced by the proc rank)
% u,v,w,p full u,v,w,p fields with dim(nx,ny,nz)
% ibeg,jbeg,kbeg arrays with length(n.proc) with processor grid
% iend,jend,kend """
% time simulation time
% Optional input (key/value pair)
% endian endianess of the file as used by fopen (default: 'n')
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Parse variable input
par = inputParser;
addOptional(par,'endian','n',@ischar);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
% Parse filename
[fdir,fbase,fext] = fileparts(file);
fbasepath = fullfile(fdir,fbase);
nxprocs = numel(ibegu);
nyprocs = numel(jbegu);
nzprocs = numel(kbegu);
nxu = size(u,1); nxv = size(v,1); nxw = size(w,1); nxp = size(p,1);
nyu = size(u,2); nyv = size(v,2); nyw = size(w,2); nyp = size(p,2);
nzu = size(u,3); nzv = size(v,3); nzw = size(w,3); nzp = size(p,3);
for ixproc=0:nxprocs-1
for iyproc=0:nyprocs-1
for izproc=0:nzprocs-1
[iproc] = locfun_proc_id(ixproc,iyproc,izproc,nxprocs,nyprocs,nzprocs);
fname = sprintf('%s.%05d',fbasepath,iproc);
% Extract local indices
ibu = ibegu(ixproc+1); ibv = ibegv(ixproc+1); ibw = ibegw(ixproc+1); ibp = ibegp(ixproc+1);
ieu = iendu(ixproc+1); iev = iendv(ixproc+1); iew = iendw(ixproc+1); iep = iendp(ixproc+1);
jbu = jbegu(iyproc+1); jbv = jbegv(iyproc+1); jbw = jbegw(iyproc+1); jbp = jbegp(iyproc+1);
jeu = jendu(iyproc+1); jev = jendv(iyproc+1); jew = jendw(iyproc+1); jep = jendp(iyproc+1);
kbu = kbegu(izproc+1); kbv = kbegv(izproc+1); kbw = kbegw(izproc+1); kbp = kbegp(izproc+1);
keu = kendu(izproc+1); kev = kendv(izproc+1); kew = kendw(izproc+1); kep = kendp(izproc+1);
% Create a new chunk file
obj = ucf(fname,'create','verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.setFileHeader('field','rank',[iproc,ixproc,iyproc,izproc],'endian',par.Results.endian);
obj.appendStep(time);
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);
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);
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);
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);
obj.appendSet(data,params);
% Close chunk file
obj.close();
end
end
end
end
function [ind] = locfun_proc_id(ii,jj,kk,nx,ny,nz)
% local version of 'sub2ind_zero_row'
ind=ii*ny*nz+jj*nz+kk;
end

112
matlab/write_uvwp_newmesh_ucf.m Normal file
View File

@ -0,0 +1,112 @@
function [] = write_uvwp_newmesh_ucf(file,...
u,v,w,p,...
xuo,yuo,zuo,xvo,yvo,zvo,xwo,ywo,zwo,xpo,ypo,zpo,...
xun,yun,zun,xvn,yvn,zvn,xwn,ywn,zwn,xpn,ypn,zpn,...
ibegu,jbegu,kbegu,iendu,jendu,kendu,...
ibegv,jbegv,kbegv,iendv,jendv,kendv,...
ibegw,jbegw,kbegw,iendw,jendw,kendw,...
ibegp,jbegp,kbegp,iendp,jendp,kendp,...
time,varargin)
% [] = write_uvwp_newmesh_ucf(file,...
% u,v,w,p,...
% xuo,yuo,zuo,xvo,yvo,zvo,xwo,ywo,zwo,xpo,ypo,zpo,...
% xun,yun,zun,xvn,yvn,zvn,xwn,ywn,zwn,xpn,ypn,zpn,...
% ibegu,jbegu,kbegu,iendu,jendu,kendu,...
% ibegv,jbegv,kbegv,iendv,jendv,kendv,...
% ibegw,jbegw,kbegw,iendw,jendw,kendw,...
% ibegp,jbegp,kbegp,iendp,jendp,kendp,...
% time,varargin)
% Writes a complete uvwp field to UCF chunks with a different mesh.
% Caution: extrapolation might be performed (periodicity is not checked)
% Input
% file output file
% u,v,w,p complete uvwp fields on old grid
% xuo,yuo,zuo,... old mesh
% xun,yun,zun,... new mesh
% ibegu,jbegu,kbegu,... new processor grid
% time simulation time
% Optional input (key/value pair)
% method interpolation method (default: 'cubic')
% endian endianess of the file as used by fopen (default: 'n')
% verbosity verbose output? (default: no)
% debug debug output? (default: no)
% Parse variable input
par = inputParser;
addOptional(par,'method','cubic',@ischar);
addOptional(par,'endian','n',@ischar);
addOptional(par,'verbosity',0,@isnumeric);
addOptional(par,'debug',0,@isnumeric);
parse(par,varargin{:});
% Parse filename
[fdir,fbase,fext] = fileparts(file);
fbasepath = fullfile(fdir,fbase);
nxprocs = numel(ibegu);
nyprocs = numel(jbegu);
nzprocs = numel(kbegu);
nxu = numel(xun); nxv = numel(xvn); nxw = numel(xwn); nxp = numel(xpn);
nyu = numel(yun); nyv = numel(yvn); nyw = numel(ywn); nyp = numel(ypn);
nzu = numel(zun); nzv = numel(zvn); nzw = numel(zwn); nzp = numel(zpn);
% Create interpolants
giu = griddedInterpolant({xuo,yuo,zuo},u,par.Results.method);
giv = griddedInterpolant({xvo,yvo,zvo},v,par.Results.method);
giw = griddedInterpolant({xwo,ywo,zwo},w,par.Results.method);
gip = griddedInterpolant({xpo,ypo,zpo},p,par.Results.method);
for ixproc=0:nxprocs-1
for iyproc=0:nyprocs-1
for izproc=0:nzprocs-1
[iproc] = locfun_proc_id(ixproc,iyproc,izproc,nxprocs,nyprocs,nzprocs);
fname = sprintf('%s.%05d',fbasepath,iproc);
% Get processor boundaries
ibu = ibegu(ixproc+1); ibv = ibegv(ixproc+1); ibw = ibegw(ixproc+1); ibp = ibegp(ixproc+1);
ieu = iendu(ixproc+1); iev = iendv(ixproc+1); iew = iendw(ixproc+1); iep = iendp(ixproc+1);
jbu = jbegu(iyproc+1); jbv = jbegv(iyproc+1); jbw = jbegw(iyproc+1); jbp = jbegp(iyproc+1);
jeu = jendu(iyproc+1); jev = jendv(iyproc+1); jew = jendw(iyproc+1); jep = jendp(iyproc+1);
kbu = kbegu(izproc+1); kbv = kbegv(izproc+1); kbw = kbegw(izproc+1); kbp = kbegp(izproc+1);
keu = kendu(izproc+1); kev = kendv(izproc+1); kew = kendw(izproc+1); kep = kendp(izproc+1);
% Create output file
obj = ucf(fname,'create','verbosity',par.Results.verbosity,'debug',par.Results.debug);
obj.setFileHeader('field','rank',[iproc,ixproc,iyproc,izproc],'endian',par.Results.endian);
obj.appendStep(time);
% Interpolate to refined chunk and write
tmp = giu({xun(ibu:ieu),yun(jbu:jeu),zun(kbu:keu)});
if any(isnan(tmp))
warning('NaN detected: u, proc #%d, %d occurences',iproc,sum(isnan(tmp(:))));
end
obj.appendField(tmp,0,ibu,jbu,kbu,nxu,nyu,nzu);
tmp = giv({xvn(ibv:iev),yvn(jbv:jev),zvn(kbv:kev)});
if any(isnan(tmp))
warning('NaN detected: v, proc #%d, %d occurences',iproc,sum(isnan(tmp(:))));
end
obj.appendField(tmp,0,ibv,jbv,kbv,nxv,nyv,nzv);
tmp = giw({xwn(ibw:iew),ywn(jbw:jew),zwn(kbw:kew)});
if any(isnan(tmp))
warning('NaN detected: w, proc #%d, %d occurences',iproc,sum(isnan(tmp(:))));
end
obj.appendField(tmp,0,ibw,jbw,kbw,nxw,nyw,nzw);
tmp = gip({xpn(ibp:iep),ypn(jbp:jep),zpn(kbp:kep)});
if any(isnan(tmp))
warning('NaN detected: p, proc #%d, %d occurences',iproc,sum(isnan(tmp(:))));
end
obj.appendField(tmp,0,ibp,jbp,kbp,nxp,nyp,nzp);
obj.close();
clear obj
fclose all;
end
end
end
end
function [ind] = locfun_proc_id(ii,jj,kk,nx,ny,nz)
% local version of 'sub2ind_zero_row'
ind=ii*ny*nz+jj*nz+kk;
end