zhyang-dev/CFDEMcoupling-PFM
1
0
Fork 0
mirror of https://github.com/ParticulateFlow/CFDEMcoupling-PFM.git synced 2025-12-08 06:37:44 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

recurrence model: large, incomplete data bases

Browse Source 
The standardRecModel class reads all fields of the data base at once. This
might create a problem in cases with a large number of snapshots and/or
a large mesh, as the machine's memory (RAM) would become a bottleneck.

Thus, the class gerhardsRecModel implements an incomplete data base, with
a user-defined number of slots M. In cases with a larger number of
snapshots on disk N, with N > M, the class will manage its data base
in a fashion similar to an operating system managing memory pages.

The class gerhardsRecModel implements a least-recently used (LRU)
algorithm, which vacates the least-used of the dataBase's M slots.
An integer list is used to track the slots' usage, each access to a
slot is logged, thus the least-used slot can be easily determined.

In order to fully utilize the LRU algorithm, the computation of the
recurrence matrix in sqrDiffNorm.C had to modified such, that three
nested for-loops are used, instead of two nested loops. Thus, a certain
number of additional, essentially no-op, loop iterations are expended
in order to accomodate the LRU algorithm. Keeping the two nested loops
would have reaped only part of LRU's potential gains.

In order to accomodate data base management in the classes derived
from the class recModel, some const specifiers had to be removed.
For informational purposes, a method has been added to the class.

The class gerhardsRecModel first checks the existence of all N
to-be-read fields, and then fills the data base with the first M
fields.

Further features included in this commit:

All elements of the recurrence model are initialized with the value of -1.0
Thus, some form of error checking is introduced, as negative values should not
remain within the matrix after computation has finished.

Skipping the 0 directory of the data base. This, might be useful when
using rStatAnalysis as post-processing tool.

The class multiIntervalPath was adapted to use OpenFOAM's
methods for parallel communication.

Reference:

  Modern Operating Systems
  Andrew S. Tannenbaum,
  Prentice Hall, 1992
This commit is contained in:
MarkoRamius
2018-03-08 16:31:01 +11:00
parent 708e4c465c
commit 4403f4e191
21 changed files with 1697 additions and 225 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
src/.gitignore vendored
View File

@ -5,3 +5,5 @@
log_*
log.*
*~
lnInclude

1
src/recurrence/Make/files
View File

@ -2,6 +2,7 @@ recBase/recBase.C
recModel/recModel/recModel.C
recModel/recModel/newRecModel.C
recModel/standardRecModel/standardRecModel.C
recModel/gerhardsRecModel/gerhardsRecModel.C
recNorm/recNorm/recNorm.C
recNorm/recNorm/newRecNorm.C
recNorm/sqrDiffNorm/sqrDiffNorm.C

1026
src/recurrence/recModel/gerhardsRecModel/gerhardsRecModel.C Normal file
View File

File diff suppressed because it is too large Load Diff

209
src/recurrence/recModel/gerhardsRecModel/gerhardsRecModel.H Normal file
View File

@ -0,0 +1,209 @@
/*---------------------------------------------------------------------------*\
CFDEMcoupling academic - Open Source CFD-DEM coupling
Contributing authors:
Thomas Lichtenegger, Gerhard Holzinger
Copyright (C) 2015- Johannes Kepler University, Linz
-------------------------------------------------------------------------------
License
This file is part of CFDEMcoupling academic.
CFDEMcoupling academic is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
CFDEMcoupling academic is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with CFDEMcoupling academic. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::gerhardsRecModel
Description
A recurrence model for a dataBase with a smaller, or equal number of slots
than there are snapshots on disk.
In the case of less slots than snapshots, we follow the Least Recently Used
page replacement algorithm described in computer science literature.
A label list is used to keep track of which slot was most recently used,
the least recently used slot is then vacated if a new snapshot needs to be
read from disk.
Reference:
\verbatim
"Modern Operating Systems"
Andrew S. Tannenbaum,
Prentice Hall, 1992
\endverbatim
SourceFiles
gerhardsRecModelI.H
gerhardsRecModel.C
\*---------------------------------------------------------------------------*/
#ifndef gerhardsRecModel_H
#define gerhardsRecModel_H
#include "recModel.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class gerhardsRecModel Declaration
\*---------------------------------------------------------------------------*/
class gerhardsRecModel
:
public recModel
{
protected:
dictionary propsDict_;
word dataBaseName_;
Foam::Time recTime;
instantList timeDirs;
Switch skipZero_;
Switch checkSkipZero_;
label numRecFields_;
const Switch verboseVerbose_;
const Switch verboseVacate_;
// matrix that contains the recurrence ERROR
SymmetricSquareMatrix<scalar> recurrenceMatrix_;
// create a data structure for the time indices
// constant will not be contained
// runTimeIndex -> continuousIndex
HashTable<label,word> timeIndexList_;
HashTable<word, label> revTimeIndexList_;
// create a data structure for the time values
// constant will not be contained
// continuousIndex -> time.value()
HashTable<label,scalar> timeValueList_;
label contTimeIndex;
const label lowerSeqLim_;
const label upperSeqLim_;
scalar checkTimeStep();
inline label getVolScalarFieldIndex(word, label) const;
inline label getVolVectorFieldIndex(word, label) const;
inline label getSurfaceScalarFieldIndex(word, label) const;
void readFieldSeries();
void readTimeSeries();
void fetchStateForDataBase(label);
void updateStorageUsageList(label);
void readNewSnapshot(label, label);
Switch checkSkipZero();
public:
//- Runtime type information
TypeName("gerhardsRecModel");
// Constructors
//- Construct from components
gerhardsRecModel
(
const dictionary& dict,
recBase& base
);
// Destructor
~gerhardsRecModel();
void exportVolScalarField(word, volScalarField&);
void exportVolVectorField(word, volVectorField&);
void exportSurfaceScalarField(word, surfaceScalarField&);
const volScalarField& exportVolScalarField(word, label);
const volVectorField& exportVolVectorField(word, label);
const surfaceScalarField& exportSurfaceScalarField(word, label);
// tmp<surfaceScalarField> exportAveragedSurfaceScalarField(word, scalar, label index = -1);
void exportAveragedVolVectorField(volVectorField&, word, scalar, label index = -1) const;
SymmetricSquareMatrix<scalar>& recurrenceMatrix();
const HashTable<label,word>& timeIndexList() const;
label lowerSeqLim() const;
label upperSeqLim() const;
label numRecFields() const;
label numDataBaseFields() const;
void updateRecFields();
void writeRecMatrix() const;
private:
/*
Number of fields kept in the dataBase, this number is smaller, or equal, as the
number of snapshots.
*/
const label numDataBaseFields_;
/*
As the data base holds less states than there are snapshots, we need to translate
the index of the snapshot to an appropriate index of the dataBase.
*/
List<label> storageIndex_;
/*
As we repeatedly need to replace fields in the dataBase with new snapshots from
disk, we need a way to judge, which slot to vacate for the new snapshot.
Here, we follow the lead of paging algorithms used by OS'es in memory management.
The storageUsageList_ represents the recentness of a slot's use; with 1 for the
most recently used element, and higher values for less recently used ones.
*/
List<label> storageUsageList_;
// this class does not read all the fields
List<PtrList<volScalarField> > volScalarFieldList_;
List<PtrList<volVectorField> > volVectorFieldList_;
List<PtrList<surfaceScalarField> > surfaceScalarFieldList_;
// book-keeping
label nrOfReadsFromDisk_;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "gerhardsRecModelI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

106
src/recurrence/recModel/gerhardsRecModel/gerhardsRecModelI.H Normal file
View File

@ -0,0 +1,106 @@
/*---------------------------------------------------------------------------*\
CFDEMcoupling academic - Open Source CFD-DEM coupling
Contributing authors:
Thomas Lichtenegger, Gerhard Holzinger
Copyright (C) 2015- Johannes Kepler University, Linz
-------------------------------------------------------------------------------
License
This file is part of CFDEMcoupling academic.
CFDEMcoupling academic is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
CFDEMcoupling academic is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with CFDEMcoupling academic. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline Foam::label Foam::gerhardsRecModel::getVolScalarFieldIndex(word fieldname, label index) const
{
if(index < 0 || index > numRecFields_ - 1)
{
FatalError<<"gerhardsRecModel: Index out of bounds for volScalarField with name "
<< fieldname << abort(FatalError);
}
label fieldI(-1);
for(int i=0; i<volScalarFieldNames_.size(); i++)
{
if(volScalarFieldNames_[i].match(fieldname))
{
fieldI = i;
}
}
if (fieldI == -1)
{
FatalError<<"gerhardsRecModel: Could not find volScalarField with name "
<< fieldname << abort(FatalError);
}
return fieldI;
}
inline Foam::label Foam::gerhardsRecModel::getVolVectorFieldIndex(word fieldname, label index) const
{
if(index < 0 || index > numRecFields_ - 1)
{
FatalError<<"gerhardsRecModel: Index out of bounds for volVectorField with name "
<< fieldname << abort(FatalError);
}
label fieldI(-1);
for(int i=0; i<volVectorFieldNames_.size(); i++)
{
if(volVectorFieldNames_[i].match(fieldname))
{
fieldI = i;
}
}
if (fieldI == -1)
{
FatalError<<"gerhardsRecModel: Could not find volVectorField with name "
<< fieldname << abort(FatalError);
}
return fieldI;
}
inline Foam::label Foam::gerhardsRecModel::getSurfaceScalarFieldIndex(word fieldname, label index) const
{
if(index < 0 || index > numRecFields_ - 1)
{
FatalError<<"gerhardsRecModel: Index out of bounds for surfaceScalarField with name "
<< fieldname << abort(FatalError);
}
label fieldI(-1);
for(int i=0; i<surfaceScalarFieldNames_.size(); i++)
{
if(surfaceScalarFieldNames_[i].match(fieldname))
{
fieldI = i;
}
}
if (fieldI == -1)
{
FatalError<<"gerhardsRecModel: Could not find surfaceScalarField with name "
<< fieldname << abort(FatalError);
}
return fieldI;
}

3
src/recurrence/recModel/recModel/recModel.C
View File

@ -63,7 +63,7 @@ recModel::recModel
IOobject::NO_WRITE
)
),
verbose_(false),
verbose_(dict.lookupOrDefault<Switch>("verbose", false)),
volScalarFieldNames_(recProperties_.lookup("volScalarFields")),
volVectorFieldNames_(recProperties_.lookup("volVectorFields")),
surfaceScalarFieldNames_(recProperties_.lookup("surfaceScalarFields")),
@ -80,7 +80,6 @@ recModel::recModel
virtualTimeIndexList_(0),
virtualTimeIndexListPos(0)
{
if (recProperties_.found("verbose")) verbose_=true;
recTimeStep_ = -1.0;
totRecSteps_ = -1;
}

15
src/recurrence/recModel/recModel/recModel.H
View File

@ -125,16 +125,24 @@ public:
// Member Functions
virtual void exportVolScalarField(word, volScalarField&) const = 0;
/*virtual void exportVolScalarField(word, volScalarField&) const = 0;
virtual void exportVolVectorField(word, volVectorField&) const = 0;
virtual void exportSurfaceScalarField(word, surfaceScalarField&) const = 0;
virtual const volScalarField& exportVolScalarField(word, label) const = 0;
virtual const volVectorField& exportVolVectorField(word, label) const = 0;
virtual const surfaceScalarField& exportSurfaceScalarField(word, label) const = 0;
virtual const surfaceScalarField& exportSurfaceScalarField(word, label) const = 0;*/
virtual void exportVolScalarField(word, volScalarField&) = 0;
virtual void exportVolVectorField(word, volVectorField&) = 0;
virtual void exportSurfaceScalarField(word, surfaceScalarField&) = 0;
virtual const volScalarField& exportVolScalarField(word, label) = 0;
virtual const volVectorField& exportVolVectorField(word, label) = 0;
virtual const surfaceScalarField& exportSurfaceScalarField(word, label) = 0;
// virtual tmp<surfaceScalarField> exportAveragedSurfaceScalarField(word, scalar, label index = -1) = 0;
virtual void exportAveragedVolVectorField(volVectorField&, word, scalar, label index = -1) const = 0;
//virtual void exportAveragedVolVectorField(volVectorField&, word, scalar, label index = -1) const = 0;
virtual SymmetricSquareMatrix<scalar>& recurrenceMatrix() = 0;
@ -146,6 +154,7 @@ public:
virtual label upperSeqLim() const = 0;
virtual label numRecFields() const = 0;
virtual label numDataBaseFields() const = 0;
label totRecSteps() const;

155
src/recurrence/recModel/standardRecModel/standardRecModel.C
View File

@ -26,7 +26,6 @@ License
#include "Random.H"
#include "standardRecModel.H"
#include "addToRunTimeSelectionTable.H"
#include <mpi.h>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -56,12 +55,14 @@ standardRecModel::standardRecModel
:
recModel(dict,base),
propsDict_(dict.subDict(typeName + "Props")),
//dataBaseName_("dataBase"),
dataBaseName_(propsDict_.lookupOrDefault<word>("dataBase", word("dataBase"))),
//recTime("dataBase", "", "../system", "../constant", false),
recTime(fileName(dataBaseName_), "", "../system", "../constant", false),
timeDirs(recTime.times()),
skipZero_(propsDict_.lookupOrDefault<Switch>("skipZero", Switch(false))),
numRecFields_(skipZero_ ? label(timeDirs.size())-1 : label(timeDirs.size())),
recurrenceMatrix_(numRecFields_,scalar(0.0)),
recurrenceMatrix_(numRecFields_,scalar(-1.0)),
timeIndexList_(numRecFields_-1),
timeValueList_(numRecFields_-1),
contTimeIndex(0),
@ -71,38 +72,39 @@ standardRecModel::standardRecModel
volScalarFieldList_(volScalarFieldNames_.size()),
volVectorFieldList_(volVectorFieldNames_.size()),
surfaceScalarFieldList_(surfaceScalarFieldNames_.size())
{
if (verbose_)
{
if (verbose_)
{
// be informative on properties of the "recTime" Time-object
Info << "recTime.rootPath() " << recTime.rootPath() << endl;
Info << "recTime.caseName() " << recTime.caseName() << endl;
Info << "recTime.path() " << recTime.path() << endl;
Info << "recTime.timePath() " << recTime.timePath() << endl;
Info << "recTime.timeName() " << recTime.timeName() << endl;
Info << "timeDirs " << timeDirs << endl;
Info << "ignoring 0 directory: " << skipZero_ << endl;
// be informative on properties of the "recTime" Time-object
Info << "recTime.rootPath() " << recTime.rootPath() << endl;
Info << "recTime.caseName() " << recTime.caseName() << endl;
Info << "recTime.path() " << recTime.path() << endl;
Info << "recTime.timePath() " << recTime.timePath() << endl;
Info << "recTime.timeName() " << recTime.timeName() << endl;
Info << "timeDirs " << timeDirs << endl;
Info << "ignoring 0 directory: " << skipZero_ << endl;
}
readTimeSeries();
recTimeStep_ = checkTimeStep();
totRecSteps_ = 1+static_cast<label> ((endTime_-startTime_) / recTimeStep_);
totRecSteps_ = 1 + static_cast<label>( (endTime_-startTime_) / recTimeStep_ );
for(int i=0; i<volScalarFieldNames_.size(); i++)
{
volScalarFieldList_[i].setSize(numRecFields_);
}
for(int i=0; i<volVectorFieldNames_.size(); i++)
{
volVectorFieldList_[i].setSize(numRecFields_);
}
for(int i=0; i<surfaceScalarFieldNames_.size(); i++)
{
surfaceScalarFieldList_[i].setSize(numRecFields_);
// setRecFields();
}
// setRecFields();
}
@ -126,9 +128,9 @@ scalar standardRecModel::checkTimeStep()
if (timeIndexList_.size() == 1)
{
return 1e10;
return 1.e10;
}
forAll(timeValueList_, i)
{
// skip zero
@ -200,11 +202,11 @@ void standardRecModel::readFieldSeries()
for (instantList::iterator it=timeDirs.begin(); it != timeDirs.end(); ++it)
{
if(counter >= 0.1 * percentage * size)
{
Info << "\t" << 10 * percentage << " \% done" << endl;
percentage++;
}
counter++;
{
Info << "\t" << 10 * percentage << " \% done" << endl;
percentage++;
}
counter++;
// set time
recTime.setTime(*it, it->value());
@ -232,9 +234,10 @@ void standardRecModel::readFieldSeries()
}
for(int i=0; i<volScalarFieldNames_.size(); i++)
{
volScalarFieldList_[i].set
(
timeIndexList_(recTime.timeName()),
timeIndexList_(recTime.timeName()),
new volScalarField
(
IOobject
@ -247,12 +250,14 @@ void standardRecModel::readFieldSeries()
),
base_.mesh()
)
);
);
}
for(int i=0; i<volVectorFieldNames_.size(); i++)
for(int i=0; i<volVectorFieldNames_.size(); i++)
{
volVectorFieldList_[i].set
(
timeIndexList_(recTime.timeName()),
timeIndexList_(recTime.timeName()),
new volVectorField
(
IOobject
@ -265,12 +270,14 @@ void standardRecModel::readFieldSeries()
),
base_.mesh()
)
);
);
}
for(int i=0; i<surfaceScalarFieldNames_.size(); i++)
for(int i=0; i<surfaceScalarFieldNames_.size(); i++)
{
surfaceScalarFieldList_[i].set
(
timeIndexList_(recTime.timeName()),
timeIndexList_(recTime.timeName()),
new surfaceScalarField
(
IOobject
@ -283,7 +290,8 @@ void standardRecModel::readFieldSeries()
),
base_.mesh()
)
);
);
}
}
Info << "Reading fields done" << endl;
}
@ -292,7 +300,7 @@ void standardRecModel::readFieldSeries()
void standardRecModel::readTimeSeries()
{
bool firsttime = true;
// fill the data structure for the time indices
// fill the data structure for the time indices
for (instantList::iterator it=timeDirs.begin(); it != timeDirs.end(); ++it)
{
// set run-time
@ -317,11 +325,11 @@ void standardRecModel::readTimeSeries()
}
if (firsttime)
{
firsttime = false;
recStartTime_ = recTime.value();
}
recEndTime_ = recTime.value();
{
firsttime = false;
recStartTime_ = recTime.value();
}
recEndTime_ = recTime.value();
// insert the time name into the hash-table with a continuous second index
timeIndexList_.insert(recTime.timeName(), contTimeIndex);
@ -348,48 +356,48 @@ void standardRecModel::readTimeSeries()
if (verbose_)
{
Info << endl;
Info << "Found " << label(timeDirs.size()) << " time folders" << endl;
Info << endl;
Info << "Found " << label(timeDirs.size()) << " time folders" << endl;
Info << "Found " << label(timeIndexList_.size()) << " time steps" << endl;
Info << "database start time = " << recStartTime_ << endl;
Info << "database end time = " << recEndTime_ << endl;
Info << "database start time = " << recStartTime_ << endl;
Info << "database end time = " << recEndTime_ << endl;
}
}
void standardRecModel::exportVolScalarField(word fieldname, volScalarField& field) const
void standardRecModel::exportVolScalarField(word fieldname, volScalarField& field)
{
field = exportVolScalarField(fieldname, virtualTimeIndex);
}
void standardRecModel::exportVolVectorField(word fieldname, volVectorField& field) const
void standardRecModel::exportVolVectorField(word fieldname, volVectorField& field)
{
field = exportVolVectorField(fieldname, virtualTimeIndex);
}
void standardRecModel::exportSurfaceScalarField(word fieldname, surfaceScalarField& field) const
void standardRecModel::exportSurfaceScalarField(word fieldname, surfaceScalarField& field)
{
field = exportSurfaceScalarField(fieldname, virtualTimeIndex);
}
const volScalarField& standardRecModel::exportVolScalarField(word fieldname, label index) const
const volScalarField& standardRecModel::exportVolScalarField(word fieldname, label index)
{
const label fieldI = getVolScalarFieldIndex(fieldname, index);
return volScalarFieldList_[fieldI][index];
}
const volVectorField& standardRecModel::exportVolVectorField(word fieldname, label index) const
const volVectorField& standardRecModel::exportVolVectorField(word fieldname, label index)
{
const label fieldI = getVolVectorFieldIndex(fieldname, index);
return volVectorFieldList_[fieldI][index];
}
const surfaceScalarField& standardRecModel::exportSurfaceScalarField(word fieldname, label index) const
const surfaceScalarField& standardRecModel::exportSurfaceScalarField(word fieldname, label index)
{
const label fieldI = getSurfaceScalarFieldIndex(fieldname, index);
@ -425,21 +433,36 @@ label standardRecModel::numRecFields() const
return numRecFields_;
}
label standardRecModel::numDataBaseFields() const
{
return numRecFields_;
}
void standardRecModel::updateRecFields()
{
virtualTimeIndex=virtualTimeIndexNext;
virtualTimeIndexNext++;
if (virtualTimeIndexNext>sequenceEnd)
{
virtualTimeIndexListPos++;
sequenceStart=virtualTimeIndexList_[virtualTimeIndexListPos].first();
sequenceEnd=virtualTimeIndexList_[virtualTimeIndexListPos].second();
virtualTimeIndexNext=sequenceStart;
}
if (verbose_)
Info << "\nUpdating virtual time index to " << virtualTimeIndex << ".\n" << endl;
virtualTimeIndex = virtualTimeIndexNext;
virtualTimeIndexNext++;
if (virtualTimeIndexNext > sequenceEnd)
{
virtualTimeIndexListPos++; // this is problematic with noPath
sequenceStart = virtualTimeIndexList_[virtualTimeIndexListPos].first();
sequenceEnd = virtualTimeIndexList_[virtualTimeIndexListPos].second();
virtualTimeIndexNext = sequenceStart;
if (verbose_)
{
Info << " new sequence (start/end) : " << sequenceStart << " / " << sequenceEnd << endl;
}
}
if (verbose_)
{
Info << "\nUpdating virtual time index to " << virtualTimeIndex << ".\n" << endl;
}
}

22
src/recurrence/recModel/standardRecModel/standardRecModel.H
View File

@ -20,6 +20,14 @@ License
You should have received a copy of the GNU General Public License
along with CFDEMcoupling academic. If not, see <http://www.gnu.org/licenses/>.
Description
Standard Recurrence Model
This class sets limits for the length of a recurrence-sequence. As of now,
such a sequence is between one 20-th and 5-th of the total number of data
base entries.
\*---------------------------------------------------------------------------*/
#ifndef standardRecModel_H
@ -96,13 +104,13 @@ public:
~standardRecModel();
void exportVolScalarField(word, volScalarField&) const;
void exportVolVectorField(word, volVectorField&) const;
void exportSurfaceScalarField(word, surfaceScalarField&) const;
void exportVolScalarField(word, volScalarField&);
void exportVolVectorField(word, volVectorField&);
void exportSurfaceScalarField(word, surfaceScalarField&);
const volScalarField& exportVolScalarField(word, label) const;
const volVectorField& exportVolVectorField(word, label) const;
const surfaceScalarField& exportSurfaceScalarField(word, label) const;
const volScalarField& exportVolScalarField(word, label);
const volVectorField& exportVolVectorField(word, label);
const surfaceScalarField& exportSurfaceScalarField(word, label);
// tmp<surfaceScalarField> exportAveragedSurfaceScalarField(word, scalar, label index = -1);
void exportAveragedVolVectorField(volVectorField&, word, scalar, label index = -1) const;
@ -115,6 +123,7 @@ public:
label upperSeqLim() const;
label numRecFields() const;
label numDataBaseFields() const;
void updateRecFields();
@ -125,6 +134,7 @@ private:
List<PtrList<volScalarField> > volScalarFieldList_;
List<PtrList<volVectorField> > volVectorFieldList_;
List<PtrList<surfaceScalarField> > surfaceScalarFieldList_;
};

28
src/recurrence/recModel/standardRecModel/standardRecModelI.H
View File

@ -85,22 +85,22 @@ inline Foam::label Foam::standardRecModel::getSurfaceScalarFieldIndex(word field
FatalError<<"standardRecModel: Index out of bounds for surfaceScalarField with name "
<< fieldname << abort(FatalError);
}
label fieldI(-1);
for(int i=0; i<surfaceScalarFieldNames_.size(); i++)
{
if(surfaceScalarFieldNames_[i].match(fieldname))
{
fieldI = i;
}
}
if (fieldI == -1)
{
FatalError<<"standardRecModel: Could not find surfaceScalarField with name "
<< fieldname << abort(FatalError);
}
return fieldI;
{
fieldI = i;
}
}
if (fieldI == -1)
{
FatalError<<"standardRecModel: Could not find surfaceScalarField with name "
<< fieldname << abort(FatalError);
}
return fieldI;
}

5
src/recurrence/recNorm/recNorm/recNorm.C
View File

@ -24,7 +24,7 @@ License
#include "error.H"
#include "recNorm.H"
#include <unistd.h>
//#include <unistd.h>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -52,9 +52,8 @@ recNorm::recNorm
:
base_(base),
recProperties_(dict),
verbose_(false)
verbose_(dict.lookupOrDefault<Switch>("verbose", false))
{
if (recProperties_.found("verbose")) verbose_=true;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //

2
src/recurrence/recNorm/recNorm/recNorm.H
View File

@ -26,7 +26,7 @@ License
#define recNorm_H
#include "recBase.H"
#include "fvCFD.H"
//#include "fvCFD.H"
#include "HashTable.H"
// #include "labelPair.H"
#include "OFstream.H"

175
src/recurrence/recNorm/sqrDiffNorm/sqrDiffNorm.C
View File

@ -69,71 +69,141 @@ sqrDiffNorm::~sqrDiffNorm()
void sqrDiffNorm::computeRecMatrix()
{
Info<< "\nsqrDiffNorm: computing recurrence matrix\n" << endl;
Info << nl << type() << ": computing recurrence matrix\n" << endl;
const HashTable<label,word>& timeIndexList( base_.recM().timeIndexList() );
SymmetricSquareMatrix<scalar>& recurrenceMatrix( base_.recM().recurrenceMatrix() );
scalar normIJ(0.0);
scalar maxNormIJ(0.0);
label size = timeIndexList.size();
/*
total number of computed elements: total number of matrix entries
minus the main diagonal entries, divided by two,
since the matrix is symmetric
*/
label size = (timeIndexList.size()*(timeIndexList.size()-1))/2;
label counter = 0;
label percentage = 0;
// perform un-normalized calculation for lower half of the recurrence matrix
forAll(timeIndexList, ti)
label N(this->base_.recM().numRecFields());
label M(this->base_.recM().numDataBaseFields());
if (verbose_)
{
if(counter >= 0.1 * percentage * size)
{
Info << "\t" << 10 * percentage << " \% done" << endl;
percentage++;
}
counter++;
forAll(timeIndexList, tj)
{
if(verbose_)
Info<<"\n Doing calculation for element " << ti << " " << tj << "\n" << endl;
// skip main diagonal and upper half
if (ti >= tj)
{
recurrenceMatrix[ti][tj] = 0;
continue;
}
// compute elements
if (fieldType_ == "volScalarField")
normIJ = normVSF(ti,tj);
else if (fieldType_ == "volVectorField")
normIJ = normVVF(ti,tj);
else if (fieldType_ == "surfaceScalarField")
normIJ = normSSF(ti,tj);
else
FatalError<<"sqrDiffNorm: Unknown field type " << fieldType_ << abort(FatalError);
recurrenceMatrix[ti][tj] = normIJ;
if (normIJ > maxNormIJ)
maxNormIJ = normIJ;
}
Info << " N = " << N << ", M = " << M << endl;
}
label ti(0);
label tj(0);
label tmp(-1);
for (int j=0; j<=N/(M-1); j++)
{
for (int i=0; i<timeIndexList.size(); i++)
{
if (verbose_)
{
Info << " i = " << i << ", j = " << j << endl;
}
if(counter >= 0.1 * percentage * size)
{
Info << "\t" << 10 * percentage << " \% done" << endl;
percentage++;
}
for (int k=i; k<i+(M-1); k++)
{
ti = i;
tj = j*(M-1) + k;
if (ti > tj)
{
tmp = ti;
ti = tj;
tj = tmp;
}
// skip coordinates outside the recurrence space
if (ti >= N or tj >= N)
{
continue;
}
// start
// skip main diagonal and upper half
if (ti >= tj)
{
recurrenceMatrix[ti][tj] = 0;
continue;
}
if (verbose_)
{
Info << " Doing calculation for element "
<< ti << " " << tj << endl;
}
counter++;
// compute elements
if (fieldType_ == "volScalarField")
{
normIJ = normVSF(ti,tj);
}
else if (fieldType_ == "volVectorField")
{
normIJ = normVVF(ti,tj);
}
else if (fieldType_ == "surfaceScalarField")
{
normIJ = normSSF(ti,tj);
}
else
{
FatalError
<< "sqrDiffNorm: Unknown field type " << fieldType_
<< abort(FatalError);
}
recurrenceMatrix[ti][tj] = normIJ;
if (normIJ > maxNormIJ)
{
maxNormIJ = normIJ;
}
// end
}
}
}
// normalize matrix and copy lower to upper half
if(normConstant_ > 0.0) maxNormIJ = normConstant_;
forAll(timeIndexList, ti)
{
forAll(timeIndexList, tj)
{
if (ti >= tj)
continue;
recurrenceMatrix[ti][tj] /= maxNormIJ;
recurrenceMatrix[tj][ti] = recurrenceMatrix[ti][tj];
}
forAll(timeIndexList, tj)
{
if (ti >= tj) continue;
if (recurrenceMatrix[ti][tj] < 0)
{
FatalErrorInFunction << "Error in computation of recurrence matrix!"
<< nl << "Negative elements encountered. This should not happen!"
<< abort(FatalError);
}
recurrenceMatrix[ti][tj] /= maxNormIJ;
recurrenceMatrix[tj][ti] = recurrenceMatrix[ti][tj];
}
}
Info<< "\nComputing recurrence matrix done\n" << endl;
Info << "\nComputing recurrence matrix done\n" << endl;
}
scalar sqrDiffNorm::normVSF(label ti, label tj)
@ -141,6 +211,7 @@ scalar sqrDiffNorm::normVSF(label ti, label tj)
const volScalarField& t1( base_.recM().exportVolScalarField(fieldName_,ti) );
const volScalarField& t2( base_.recM().exportVolScalarField(fieldName_,tj) );
dimensionedScalar tNorm( fvc::domainIntegrate( sqr( t1 - t2 ) ) );
return tNorm.value();
}
@ -149,6 +220,7 @@ scalar sqrDiffNorm::normVVF(label ti, label tj)
const volVectorField& t1( base_.recM().exportVolVectorField(fieldName_,ti) );
const volVectorField& t2( base_.recM().exportVolVectorField(fieldName_,tj) );
dimensionedScalar tNorm( fvc::domainIntegrate( magSqr( t1 - t2 ) ) );
return tNorm.value();
}
@ -158,6 +230,7 @@ scalar sqrDiffNorm::normSSF(label ti, label tj)
const surfaceScalarField& t2( base_.recM().exportSurfaceScalarField(fieldName_,tj) );
volVectorField t12 (fvc::reconstruct( t1-t2 ) );
dimensionedScalar tNorm( fvc::domainIntegrate( magSqr( t12 ) ) );
return tNorm.value();
}

25
src/recurrence/recPath/multiIntervalPath/multiIntervalPath.C
View File

@ -125,27 +125,26 @@ multiIntervalPath::~multiIntervalPath()
void multiIntervalPath::getRecPath()
{
label numRecIntervals = 0;
int root = 0;
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
if(rank==root)
if(Pstream::master())
{
computeRecPath();
numRecIntervals=virtualTimeIndexList_.size();
computeRecPath();
numRecIntervals = virtualTimeIndexList_.size();
}
MPI_Bcast(&numRecIntervals,sizeof(numRecIntervals),MPI_BYTE,root,MPI_COMM_WORLD);
Pstream::scatter(numRecIntervals);
if(rank!=root)
virtualTimeIndexList_.setSize(numRecIntervals);
if(not Pstream::master())
{
virtualTimeIndexList_.setSize(numRecIntervals);
}
int listSizeBytes = 2*sizeof(numRecIntervals)*numRecIntervals;
MPI_Bcast(&virtualTimeIndexList_[0], listSizeBytes, MPI_BYTE, root, MPI_COMM_WORLD);
Pstream::scatter(virtualTimeIndexList_);
if(verbose_)
Info << "\nRecurrence path communicated to all processors.\n" << endl;
{
Info << "\nRecurrence path communicated to all processors.\n" << endl;
}
}

32
src/recurrence/recPath/noPath/noPath.C
View File

@ -68,27 +68,26 @@ noPath::~noPath()
void noPath::getRecPath()
{
label numRecIntervals = 0;
int root = 0;
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
if(rank==root)
if(Pstream::master())
{
computeRecPath();
numRecIntervals=virtualTimeIndexList_.size();
computeRecPath();
numRecIntervals=virtualTimeIndexList_.size();
}
MPI_Bcast(&numRecIntervals,sizeof(numRecIntervals),MPI_BYTE,root,MPI_COMM_WORLD);
Pstream::scatter(numRecIntervals);
if(rank!=root)
virtualTimeIndexList_.setSize(numRecIntervals);
if(not Pstream::master())
{
virtualTimeIndexList_.setSize(numRecIntervals);
}
int listSizeBytes = 2*sizeof(numRecIntervals)*numRecIntervals;
MPI_Bcast(&virtualTimeIndexList_[0], listSizeBytes, MPI_BYTE, root, MPI_COMM_WORLD);
Pstream::scatter(virtualTimeIndexList_);
if(verbose_)
Info << "\nRecurrence path communicated to all processors.\n" << endl;
{
Info << "\nRecurrence path communicated to all processors.\n" << endl;
}
}
@ -96,11 +95,16 @@ void noPath::computeRecPath()
{
labelPair seqStartEnd(0,1);
virtualTimeIndexList_.append(seqStartEnd);
if (verbose_)
{
Info << " virtualTimeIndexList_ : " << virtualTimeIndexList_ << endl;
}
}
label noPath::seqEnd(label seqStart, label & seqLength)
{
return 0;
return 0;
}
void noPath::computeJumpVector()

3
src/recurrence/recPath/recPath/recPath.C
View File

@ -53,10 +53,9 @@ recPath::recPath
:
base_(base),
recProperties_(dict),
verbose_(false),
verbose_(dict.lookupOrDefault<Switch>("verbose", false)),
virtualTimeIndexList_( base_.recM().virtualTimeIndexList() )
{
if (recProperties_.found("verbose")) verbose_=true;
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //

2
src/recurrence/recPath/recPath/recPath.H
View File

@ -30,7 +30,7 @@ License
#include "HashTable.H"
#include "labelPair.H"
#include "OFstream.H"
#include <mpi.h>
//#include <mpi.h>
namespace Foam

103
src/recurrence/recPath/simpleRandomPath/simpleRandomPath.C
View File

@ -68,27 +68,26 @@ simpleRandomPath::~simpleRandomPath()
void simpleRandomPath::getRecPath()
{
label numRecIntervals = 0;
int root = 0;
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
if(rank==root)
if(Pstream::master())
{
computeRecPath();
numRecIntervals=virtualTimeIndexList_.size();
computeRecPath();
numRecIntervals = virtualTimeIndexList_.size();
}
MPI_Bcast(&numRecIntervals,sizeof(numRecIntervals),MPI_BYTE,root,MPI_COMM_WORLD);
Pstream::scatter(numRecIntervals);
if(rank!=root)
virtualTimeIndexList_.setSize(numRecIntervals);
if(not Pstream::master())
{
virtualTimeIndexList_.setSize(numRecIntervals);
}
int listSizeBytes = 2*sizeof(numRecIntervals)*numRecIntervals;
MPI_Bcast(&virtualTimeIndexList_[0], listSizeBytes, MPI_BYTE, root, MPI_COMM_WORLD);
Pstream::scatter(virtualTimeIndexList_);
if(verbose_)
Info << "\nRecurrence path communicated to all processors.\n" << endl;
{
Info << "\nRecurrence path communicated to all processors.\n" << endl;
}
}
@ -99,21 +98,21 @@ void simpleRandomPath::computeRecPath()
Random ranGen(osRandomInteger());
label virtualTimeIndex=0;
label recSteps=0;
label seqStart=0;
label virtualTimeIndex = 0;
label recSteps = 0;
label seqStart = 0;
label lowerSeqLim( base_.recM().lowerSeqLim() );
label upperSeqLim( base_.recM().upperSeqLim() );
label seqLength=ranGen.integer(lowerSeqLim, upperSeqLim);
label seqLength = ranGen.integer(lowerSeqLim, upperSeqLim);
virtualTimeIndex=seqEnd(seqStart,seqLength);
virtualTimeIndex = seqEnd(seqStart,seqLength);
labelPair seqStartEnd(seqStart,virtualTimeIndex);
virtualTimeIndexList_.append(seqStartEnd);
recSteps+=seqLength;
recSteps += seqLength;
SymmetricSquareMatrix<scalar>& recurrenceMatrix( base_.recM().recurrenceMatrix() );
label numRecFields( base_.recM().numRecFields() );
if(base_.recM().totRecSteps() == 1)
{
Info<< "\nPrimitive recurrence path with one element.\n" << endl;
@ -129,9 +128,9 @@ void simpleRandomPath::computeRecPath()
// the new starting point of the next sequence
if (virtualTimeIndex < numRecFields/2)
{
startLoop = numRecFields/2;
endLoop = numRecFields-1 - upperSeqLim/2; // avoid running into end of recurrence database too often
endLoop--; // start of next sequence one snapshot AFTER minimum position
startLoop = numRecFields/2;
endLoop = numRecFields-1 - upperSeqLim/2; // avoid running into end of recurrence database too often
endLoop--; // start of next sequence one snapshot AFTER minimum position
}
else
{
@ -151,61 +150,73 @@ void simpleRandomPath::computeRecPath()
}
seqLength = ranGen.integer(lowerSeqLim, upperSeqLim);
virtualTimeIndex=seqEnd(seqStart,seqLength);
virtualTimeIndex = seqEnd(seqStart,seqLength);
labelPair seqStartEnd(seqStart,virtualTimeIndex);
virtualTimeIndexList_.append(seqStartEnd);
recSteps+=seqLength;
recSteps += seqLength;
}
Info<< "\nComputing recurrence path done\n" << endl;
if (verbose_)
{
Info << " virtualTimeIndexList_ : " << virtualTimeIndexList_ << endl;
}
computeJumpVector();
}
label simpleRandomPath::seqEnd(label seqStart, label & seqLength)
{
label numRecFields( base_.recM().numRecFields() );
if(seqStart+seqLength > numRecFields-1)
seqLength=numRecFields-1-seqStart;
return seqStart+seqLength;
label numRecFields( base_.recM().numRecFields() );
if(seqStart+seqLength > numRecFields-1)
{
seqLength=numRecFields-1-seqStart;
}
return seqStart+seqLength;
}
void simpleRandomPath::computeJumpVector()
{
Info << "\nComputing recurrence jump vector\n" << endl;
OFstream jumpvec("rec_jump.dat");
label numRecFields( base_.recM().numRecFields() );
label startLoop = 0;
label endLoop = 0;
SymmetricSquareMatrix<scalar>& recurrenceMatrix( base_.recM().recurrenceMatrix() );
for (label i = 0; i < numRecFields; i++)
{
if (i < numRecFields/2)
if (i < numRecFields/2)
{
startLoop = numRecFields/2;
endLoop = numRecFields-1;
startLoop = numRecFields/2;
endLoop = numRecFields-1;
}
else
{
startLoop = 0;
endLoop = numRecFields/2-1;
startLoop = 0;
endLoop = numRecFields/2-1;
}
scalar nextMinimum(GREAT);
label jumpdest = 0;
label jumpdest = 0;
for (label j = startLoop; j < endLoop; j++)
{
if (recurrenceMatrix[j][i] < nextMinimum)
{
nextMinimum = recurrenceMatrix[j][i];
jumpdest = j+1;
continue;
}
if (recurrenceMatrix[j][i] < nextMinimum)
{
nextMinimum = recurrenceMatrix[j][i];
jumpdest = j+1;
continue;
}
}
jumpvec << jumpdest << endl;
}
Info<< "\nComputing recurrence jump vector done\n" << endl;
}

5
src/recurrence/recPath/simpleRandomPath/simpleRandomPath.H
View File

@ -20,6 +20,11 @@ License
You should have received a copy of the GNU General Public License
along with CFDEMcoupling academic. If not, see <http://www.gnu.org/licenses/>.
Description
Create a path based on the recurrence statistics and add some randomness.
The minimum data base size is 5 entries for this model to work.
\*---------------------------------------------------------------------------*/
#ifndef simpleRandomPath_H

2
src/recurrence/recStatAnalysis/noRecStatAnalysis/noRecStatAnalysis.C
View File

@ -27,8 +27,6 @@ License
#include "noRecStatAnalysis.H"
#include "recModel.H"
#include "addToRunTimeSelectionTable.H"
#include <mpi.h>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam

1
src/recurrence/recStatAnalysis/standardRecStatAnalysis/standardRecStatAnalysis.C
View File

@ -27,7 +27,6 @@ License
#include "standardRecStatAnalysis.H"
#include "recModel.H"
#include "addToRunTimeSelectionTable.H"
#include <mpi.h>
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //