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

Update of random number generation. Removed deprecated multiIntervalPath model.

Browse Source
This commit is contained in:
Thomas Lichtenegger
2019-11-13 15:53:27 +01:00
parent 8f94dba403
commit fe4016eb52
12 changed files with 81 additions and 336 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
applications/solvers/rStatAnalysis/Make/options
View File

@ -23,4 +23,4 @@ EXE_LIBS = \
-lmeshTools \
-l$(CFDEM_LIB_NAME) \
$(CFDEM_ADD_LIB_PATHS) \
$(CFDEM_ADD_LIBS)
$(CFDEM_ADD_LIBS)

5
applications/solvers/rStatAnalysis/rStatAnalysis.C
View File

@ -29,11 +29,6 @@ Description
\*---------------------------------------------------------------------------*/
// #include "fvCFD.H"
// #include "singlePhaseTransportModel.H"
// #include "turbulentTransportModel.H"
// #include "fvOptions.H"
#include "recBase.H"
#include "recStatAnalysis.H"

2
src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.C
View File

@ -83,7 +83,7 @@ isotropicFluctuations::isotropicFluctuations
dimensionedScalar("D0", dimensionSet(0,0,0,0,0,0,0), D0_)
),
dtDEM_(particleCloud_.dataExchangeM().DEMts()),
ranGen_(0)
ranGen_(clock::getTime()+pid())
{
if(ignoreCellsName_ != "none")
{

3
src/recurrence/Make/files
View File

@ -12,11 +12,10 @@ recNorm/readNorm/readNorm.C
recNorm/sqrDiffNorm/sqrDiffNorm.C
recNorm/noRecNorm/noRecNorm.C
recPath/recPath/recPath.C
recPath/recPath/newRecPath.C
recPath/simpleRandomPath/simpleRandomPath.C
recPath/recPath/newRecPath.C
recPath/noPath/noPath.C
recPath/MarkovPath/MarkovPath.C
recPath/multiIntervalPath/multiIntervalPath.C
recPath/predefinedPath/predefinedPath.C
recStatAnalysis/recStatAnalysis/recStatAnalysis.C
recStatAnalysis/recStatAnalysis/newRecStatAnalysis.C

18
src/recurrence/recBase/recBase.C
View File

@ -83,14 +83,20 @@ recBase::recBase
),
couplingSubStep_(recProperties_.lookupOrDefault<label>("couplingSubStep",0))
{
recModel_ -> readFieldSeries();
recNorm_ -> computeRecMatrix();
recPath_ -> getRecPath();
recModel_ -> readFieldSeries();
if (!recStatAnalysis_->suppressMatrixAndPath())
{
recNorm_ -> computeRecMatrix();
recPath_ -> getRecPath();
}
recModel_ -> init();
recModel_ -> init();
recModel_ -> writeRecMatrix();
recModel_ -> writeRecPath();
if (!recStatAnalysis_->suppressMatrixAndPath())
{
recModel_ -> writeRecMatrix();
recModel_ -> writeRecPath();
}
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * //

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

@ -150,8 +150,7 @@ void MarkovPath::computeRecPath()
void MarkovPath::extendPath()
{
const label seed = 0;
Random ranGen(seed);
Random ranGen(clock::getTime()+pid());
SymmetricSquareMatrix<scalar>& recurrenceMatrix( base_.recM().recurrenceMatrix() );

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

@ -1,195 +0,0 @@
/*---------------------------------------------------------------------------*\
CFDEMcoupling academic - Open Source CFD-DEM coupling
Contributing authors:
Thomas Lichtenegger
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/>.
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "multiIntervalPath.H"
#include "Random.H"
#include "recModel.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(multiIntervalPath, 0);
addToRunTimeSelectionTable
(
recPath,
multiIntervalPath,
dictionary
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from components
multiIntervalPath::multiIntervalPath
(
const dictionary& dict,
recBase& base
)
:
recPath(dict, base),
propsDict_(dict.subDict(typeName + "Props")),
meanIntervalSteps_(propsDict_.lookupOrDefault<label>("meanIntervalSteps",-1)),
numIntervals_(base.recM().numIntervals()),
intervalSizes_(numIntervals_),
intervalSizesCumulative_(numIntervals_),
Pjump_(0.0),
intervalWeights_(propsDict_.lookupOrDefault<scalarList>("intervalWeights",scalarList(numIntervals_,1.0))),
intervalWeightsCumulative_(intervalWeights_)
{
if(meanIntervalSteps_<0)
{
// if no mean interval length for consecutive steps is specified, use 1/5 from first interval
meanIntervalSteps_ = static_cast<label>(0.2 * intervalSizes_[0]);
}
// normalize weights
scalar wsum = 0.0;
for(int i=0;i<numIntervals_;i++)
{
intervalSizes_[i] = base.recM().numRecFields(i);
wsum += intervalWeights_[i];
}
for(int i=0;i<numIntervals_;i++)
{
intervalWeights_[i] /= wsum;
}
for(int i=0;i<numIntervals_;i++)
{
scalar sum1 = 0.0;
scalar sum2 = 0.0;
for(int j=0;j<=i;j++)
{
sum1 += intervalWeights_[j];
sum2 += intervalSizes_[j];
}
intervalWeightsCumulative_[i] = sum1;
intervalSizesCumulative_[i] = sum2;
}
// given a jump probability of P, the probability of finding a chain of length N is
// P(N) = (1 - P)^N * P, and the mean length E(N) = (1 - P) / P
Pjump_ = 1.0 / (1 + meanIntervalSteps_);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
multiIntervalPath::~multiIntervalPath()
{}
// * * * * * * * * * * * * * protected Member Functions * * * * * * * * * * * * //
void multiIntervalPath::computeRecPath()
{
Info << "\nComputing recurrence path\n" << endl;
const label seed = 0;
Random ranGen(seed);
label virtualTimeIndex=0;
label recSteps=0;
label seqStart=0;
bool prevStepWasJump = true;
SymmetricSquareMatrix<scalar>& recurrenceMatrix( base_.recM().recurrenceMatrix() );
if(base_.recM().totRecSteps() == 1)
{
Info << "\nPrimitive recurrence path with one element.\n" << endl;
return;
}
while(recSteps <= base_.recM().totRecSteps() )
{
scalar randJump = ranGen.scalar01();
// check if current virtualTimeIndex is close to separation time
bool intervalBorder = false;
label sep = 0;
for(int i = 0;i < numIntervals_; i++)
{
sep += intervalSizes_[i];
if (sep - 1 == virtualTimeIndex) intervalBorder=true;
}
if ((randJump > Pjump_ && !intervalBorder) || prevStepWasJump)
{
virtualTimeIndex++;
prevStepWasJump = false;
}
else
{
// before jump, complete former consecutive interval
labelPair seqStartEnd(seqStart,virtualTimeIndex);
virtualTimeIndexList_.append(seqStartEnd);
recSteps += virtualTimeIndex - seqStart + 1;
// now jump
// identify interval to jump to
scalar randInterval = ranGen.scalar01();
label interval = numIntervals_-1;
for(int i = numIntervals_-2 ;i >= 0; i--)
{
if (randInterval < intervalWeightsCumulative_[i]) interval=i;
}
label startLoop = 0;
if (interval > 0) startLoop = intervalSizesCumulative_[interval-1];
label endLoop = intervalSizesCumulative_[interval] - meanIntervalSteps_;
scalar nextMinimum(GREAT);
for (label j = startLoop; j <= endLoop; j++)
{
if(abs(j - virtualTimeIndex) < meanIntervalSteps_) continue;
if (recurrenceMatrix[j][virtualTimeIndex] < nextMinimum)
{
nextMinimum = recurrenceMatrix[j][virtualTimeIndex];
seqStart = j+1;
}
}
virtualTimeIndex = seqStart;
prevStepWasJump = true;
}
}
Info << "\nComputing recurrence path done\n" << endl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// ************************************************************************* //

103
src/recurrence/recPath/multiIntervalPath/multiIntervalPath.H
View File

@ -1,103 +0,0 @@
/*---------------------------------------------------------------------------*\
CFDEMcoupling academic - Open Source CFD-DEM coupling
Contributing authors:
Thomas Lichtenegger
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/>.
Description
A recurrence database consisting of N separate intervals is assumed with separation
times t0 (start time), t1, ... tN (end time) and weights w0, ... wN-1.
\*---------------------------------------------------------------------------*/
#ifndef multiIntervalPath_H
#define multiIntervalPath_H
#include "recPath.H"
#include "scalarList.H"
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class multiIntervalPath Declaration
\*---------------------------------------------------------------------------*/
class multiIntervalPath
:
public recPath
{
protected:
// Protected data
dictionary propsDict_;
void computeRecPath();
label meanIntervalSteps_;
label numIntervals_;
labelList intervalSizes_;
labelList intervalSizesCumulative_;
scalar Pjump_;
scalarList intervalWeights_;
scalarList intervalWeightsCumulative_;
public:
//- Runtime type information
TypeName("multiIntervalPath");
// Constructors
//- Construct from components
multiIntervalPath
(
const dictionary& dict,
recBase& base
);
// Destructor
virtual ~multiIntervalPath();
// Member Functions
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

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

@ -69,8 +69,7 @@ void simpleRandomPath::computeRecPath()
{
Info << "\nComputing recurrence path\n" << endl;
const label seed = 0;
Random ranGen(seed);
Random ranGen(clock::getTime()+pid());
label virtualTimeIndex = 0;
label recSteps = 0;

67
src/recurrence/recStatAnalysis/autocorrelation/autocorrelation.C
View File

@ -73,7 +73,8 @@ autocorrelation::autocorrelation
),
base.mesh(),
dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0.0)
)
),
suppressMatrixAndPath_(propsDict_.lookupOrDefault<bool>("suppressMatrixAndPath",false))
{
if (fieldtype_ != "scalar" && fieldtype_ != "vector")
{
@ -108,31 +109,23 @@ void autocorrelation::statistics()
void autocorrelation::autocorr()
{
scalar res = 0.0;
PtrList<volScalarField> scalarFieldList;
PtrList<volVectorField> vectorFieldList;
if (fieldtype_ == "scalar") scalarFieldList.transfer(base_.recM().exportVolScalarFieldList(fieldname_));
else vectorFieldList.transfer(base_.recM().exportVolVectorFieldList(fieldname_));
if (fieldtype_ == "scalar") scalarFieldList_.transfer(base_.recM().exportVolScalarFieldList(fieldname_));
else vectorFieldList_.transfer(base_.recM().exportVolVectorFieldList(fieldname_));
label tmax = base_.recM().totRecSteps();
for (label ti = delaySteps_; ti < tmax; ti++)
{
forAll(autoCorrField_, cellI)
{
if (fieldtype_ == "scalar")
{
res = scalarFieldList[ti-delaySteps_][refCell_] * scalarFieldList[ti][cellI];
}
else
{
res = vectorFieldList[ti-delaySteps_][refCell_] & vectorFieldList[ti][cellI];
}
autoCorrField_[cellI] += res;
autoCorrField_[cellI] += autocorrSummand(ti-delaySteps_,ti,refCell_,cellI);
}
}
autoCorrField_ /= (tmax - delaySteps_);
autoCorrField_ -= meanProd();
if (normalize_)
{
volScalarField meanProd(autoCorrField_);
@ -159,17 +152,55 @@ void autocorrelation::autocorr()
dimensionSet fieldDim(0,0,0,0,0);
if (fieldtype_ == "scalar")
{
fieldDim.reset(scalarFieldList[0].dimensions());
fieldDim.reset(scalarFieldList_[0].dimensions());
}
else
{
fieldDim.reset(vectorFieldList[0].dimensions());
fieldDim.reset(vectorFieldList_[0].dimensions());
}
fieldDim = fieldDim * fieldDim;
fieldDim.reset(fieldDim * fieldDim);
autoCorrField_.dimensions().reset(fieldDim);
}
autoCorrField_.write();
}
scalar autocorrelation::autocorrSummand(label t1, label t2, label c1, label c2)
{
scalar res;
if (fieldtype_ == "scalar")
{
res = scalarFieldList_[t1][c1] * scalarFieldList_[t2][c2];
}
else
{
res = vectorFieldList_[t1][c1] & vectorFieldList_[t2][c2];
}
return res;
}
volScalarField autocorrelation::meanProd()
{
volScalarField meanProd(autoCorrField_);
if (fieldtype_ == "scalar")
{
volScalarField aveField = base_.recM().exportVolScalarFieldAve(fieldname_);
forAll(meanProd, cellI)
{
meanProd[cellI] = aveField()[cellI] * aveField()[refCell_];
}
}
else
{
volVectorField aveField = base_.recM().exportVolVectorFieldAve(fieldname_);
forAll(meanProd, cellI)
{
meanProd[cellI] = aveField()[cellI] & aveField()[refCell_];
}
}
return meanProd;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

12
src/recurrence/recStatAnalysis/autocorrelation/autocorrelation.H
View File

@ -85,9 +85,21 @@ private:
bool normalize_;
PtrList<volScalarField> scalarFieldList_;
PtrList<volVectorField> vectorFieldList_;
volScalarField autoCorrField_;
bool suppressMatrixAndPath_;
void autocorr();
scalar autocorrSummand(label, label, label, label);
volScalarField meanProd();
bool suppressMatrixAndPath() {return suppressMatrixAndPath_;}
};

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

@ -102,6 +102,8 @@ public:
virtual void statistics() = 0;
virtual bool suppressMatrixAndPath() {return false;}
};