diff --git a/bash/ucf_archive_pack b/bash/ucf_archive_pack new file mode 100755 index 0000000..8f19b94 --- /dev/null +++ b/bash/ucf_archive_pack @@ -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 diff --git a/bash/ucftar_pack b/bash/ucftar_pack new file mode 100755 index 0000000..6545bac --- /dev/null +++ b/bash/ucftar_pack @@ -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[@]} diff --git a/matlab/@ini/ini.m b/matlab/@ini/ini.m new file mode 100644 index 0000000..edc3c9e --- /dev/null +++ b/matlab/@ini/ini.m @@ -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.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.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.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 diff --git a/matlab/@ustar/ustar.m b/matlab/@ustar/ustar.m new file mode 100644 index 0000000..7d19fa7 --- /dev/null +++ b/matlab/@ustar/ustar.m @@ -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 diff --git a/matlab/colmap_from_flags.m b/matlab/colmap_from_flags.m new file mode 100644 index 0000000..57c1629 --- /dev/null +++ b/matlab/colmap_from_flags.m @@ -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 \ No newline at end of file diff --git a/matlab/convert_colmap_matlab2paraview.m b/matlab/convert_colmap_matlab2paraview.m new file mode 100644 index 0000000..a113353 --- /dev/null +++ b/matlab/convert_colmap_matlab2paraview.m @@ -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 \ No newline at end of file diff --git a/matlab/convert_fluid_hdf2ucf.m b/matlab/convert_fluid_hdf2ucf.m new file mode 100644 index 0000000..1699a7e --- /dev/null +++ b/matlab/convert_fluid_hdf2ucf.m @@ -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 diff --git a/matlab/convert_particles_array2struct.m b/matlab/convert_particles_array2struct.m new file mode 100644 index 0000000..2cf62fb --- /dev/null +++ b/matlab/convert_particles_array2struct.m @@ -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 \ No newline at end of file diff --git a/matlab/convert_particles_cell2struct.m b/matlab/convert_particles_cell2struct.m new file mode 100644 index 0000000..af6f1fc --- /dev/null +++ b/matlab/convert_particles_cell2struct.m @@ -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 \ No newline at end of file diff --git a/matlab/convert_particles_hdf2ucf.m b/matlab/convert_particles_hdf2ucf.m new file mode 100644 index 0000000..11fda7e --- /dev/null +++ b/matlab/convert_particles_hdf2ucf.m @@ -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 diff --git a/matlab/flags_from_colmap.m b/matlab/flags_from_colmap.m new file mode 100644 index 0000000..ad2c6b2 --- /dev/null +++ b/matlab/flags_from_colmap.m @@ -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 \ No newline at end of file diff --git a/matlab/ncol_from_colmap.m b/matlab/ncol_from_colmap.m new file mode 100644 index 0000000..8b7e16e --- /dev/null +++ b/matlab/ncol_from_colmap.m @@ -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 \ No newline at end of file diff --git a/matlab/ncol_from_flags.m b/matlab/ncol_from_flags.m new file mode 100644 index 0000000..82729c5 --- /dev/null +++ b/matlab/ncol_from_flags.m @@ -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 \ No newline at end of file diff --git a/matlab/read_field_chunk_ucf.m b/matlab/read_field_chunk_ucf.m new file mode 100644 index 0000000..c33ddbc --- /dev/null +++ b/matlab/read_field_chunk_ucf.m @@ -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 diff --git a/matlab/read_field_complete_ucf.m b/matlab/read_field_complete_ucf.m new file mode 100644 index 0000000..c20be9c --- /dev/null +++ b/matlab/read_field_complete_ucf.m @@ -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 diff --git a/matlab/read_grid_ucf.m b/matlab/read_grid_ucf.m new file mode 100644 index 0000000..b54b015 --- /dev/null +++ b/matlab/read_grid_ucf.m @@ -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 diff --git a/matlab/read_parameters_ucf.m b/matlab/read_parameters_ucf.m new file mode 100644 index 0000000..df2f4be --- /dev/null +++ b/matlab/read_parameters_ucf.m @@ -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 diff --git a/matlab/read_particle_balancing_ucf.m b/matlab/read_particle_balancing_ucf.m new file mode 100644 index 0000000..c628be3 --- /dev/null +++ b/matlab/read_particle_balancing_ucf.m @@ -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 diff --git a/matlab/read_particles_ucf.m b/matlab/read_particles_ucf.m new file mode 100644 index 0000000..32550c7 --- /dev/null +++ b/matlab/read_particles_ucf.m @@ -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 diff --git a/matlab/read_procgrid_ucf.m b/matlab/read_procgrid_ucf.m new file mode 100644 index 0000000..cba6031 --- /dev/null +++ b/matlab/read_procgrid_ucf.m @@ -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 diff --git a/matlab/read_scal_chunk_ucf.m b/matlab/read_scal_chunk_ucf.m new file mode 100644 index 0000000..8e11d26 --- /dev/null +++ b/matlab/read_scal_chunk_ucf.m @@ -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 diff --git a/matlab/read_scal_complete_ucf.m b/matlab/read_scal_complete_ucf.m new file mode 100644 index 0000000..9f0accb --- /dev/null +++ b/matlab/read_scal_complete_ucf.m @@ -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 diff --git a/matlab/read_statistics_channel_scal_ucf.m b/matlab/read_statistics_channel_scal_ucf.m new file mode 100644 index 0000000..af7ab02 --- /dev/null +++ b/matlab/read_statistics_channel_scal_ucf.m @@ -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 \ No newline at end of file diff --git a/matlab/read_statistics_channel_ucf.m b/matlab/read_statistics_channel_ucf.m new file mode 100644 index 0000000..e1e2cef --- /dev/null +++ b/matlab/read_statistics_channel_ucf.m @@ -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 \ No newline at end of file diff --git a/matlab/read_uvwp_chunk_ucf.m b/matlab/read_uvwp_chunk_ucf.m new file mode 100644 index 0000000..af5e27d --- /dev/null +++ b/matlab/read_uvwp_chunk_ucf.m @@ -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 diff --git a/matlab/read_uvwp_complete_ucf.m b/matlab/read_uvwp_complete_ucf.m new file mode 100644 index 0000000..be0b6f1 --- /dev/null +++ b/matlab/read_uvwp_complete_ucf.m @@ -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 diff --git a/matlab/write_particles_ucf.m b/matlab/write_particles_ucf.m new file mode 100644 index 0000000..21c876b --- /dev/null +++ b/matlab/write_particles_ucf.m @@ -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 diff --git a/matlab/write_scal_chunk_ucf.m b/matlab/write_scal_chunk_ucf.m new file mode 100644 index 0000000..d15ba7e --- /dev/null +++ b/matlab/write_scal_chunk_ucf.m @@ -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 diff --git a/matlab/write_scal_complete_ucf.m b/matlab/write_scal_complete_ucf.m new file mode 100644 index 0000000..f65703c --- /dev/null +++ b/matlab/write_scal_complete_ucf.m @@ -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 diff --git a/matlab/write_uvwp_chunk_ucf.m b/matlab/write_uvwp_chunk_ucf.m new file mode 100644 index 0000000..37e67f7 --- /dev/null +++ b/matlab/write_uvwp_chunk_ucf.m @@ -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 diff --git a/matlab/write_uvwp_complete_ucf.m b/matlab/write_uvwp_complete_ucf.m new file mode 100644 index 0000000..e2f257c --- /dev/null +++ b/matlab/write_uvwp_complete_ucf.m @@ -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 \ No newline at end of file diff --git a/matlab/write_uvwp_newmesh_ucf.m b/matlab/write_uvwp_newmesh_ucf.m new file mode 100644 index 0000000..bd9ff9e --- /dev/null +++ b/matlab/write_uvwp_newmesh_ucf.m @@ -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