zhyang-dev/CFDEMcoupling-PFM
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

More rearranging.

Browse Source
This commit is contained in:
Thomas Lichtenegger
2015-10-06 16:55:13 +02:00
parent 6bb64cd937
commit 4282dcf080
10 changed files with 82 additions and 628 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
applications/solvers/cfdemSolverRecurrence/cfdemSolverRecurrence.C
View File

@ -104,14 +104,14 @@ int main(int argc, char *argv[])
<< " s\n\n" << endl;
}
// write stuff at output time
if (runTime.outputTime())
/* if (runTime.outputTime())
{
// magicParticles.writeConcField();
alpha2.write();
U2.write();
}
*/
runTime.write(); // does not work, yet

54
applications/solvers/cfdemSolverRecurrence/createFields.H
View File

@ -24,46 +24,6 @@
),
dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0), 1.0)
);
// recurrence fields
volVectorField URec
(
IOobject
(
"URec",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
dimensionedVector("URec",dimVelocity,vector(0,0,0))
);
volScalarField voidfractionRec
(
IOobject
(
"voidfractionRec",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
dimensionedScalar("voidfractionRec", dimensionSet(0, -3, 0, 0, 0), 1.0)
);
volVectorField UsRec
(
IOobject
(
"UsRec",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
dimensionedVector("UsRec",dimVelocity,vector(0,0,0))
);
// calculated fields
volVectorField U
@ -107,20 +67,8 @@
//===============================
//# include "createPhi.H"
#ifndef createPhi_H
#define createPhi_H
Info<< "Reading/calculating face flux field phi\n" << endl;
surfaceScalarField phi
(
IOobject
(
"phi",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(URec*voidfractionRec) & mesh.Sf()
);
// create phi in recModel and put reference here
#endif

374
applications/solvers/cfdemSolverRecurrence/createRecurrentFields.H
View File

@ -1,374 +0,0 @@
// to be moved into cfdemCloudRec
Foam::Info<< "Create recurrent-time\n" << Foam::endl;
// the documentation is in the code ;-)
/* Time
(
const fileName& rootPath,
const fileName& caseName,
const word& systemName = "system",
const word& constantName = "constant",
const bool enableFunctionObjects = true
);
*/
Foam::Time recTime("dataBase", "", "../system", "../constant", false);
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;
}
// all time directories including constant
instantList timeDirs = recTime.times();
if (verbose)
{
Info << "timeDirs " << timeDirs << endl;
}
// iterator for timeDirs
instantList::iterator it;
// create a data structure for the time indices
// constant will not be contained
// runTimeIndex -> continuousIndex
HashTable<label,word> timeIndexList(label(timeDirs.size())-1);
label contTimeIndex(0);
// create a data structure for the time values
// constant will not be contained
// continuousIndex -> time.value()
HashTable<label,scalar> timeValueList(label(timeDirs.size())-1);
// fill the data structure for the time indices
for (it=timeDirs.begin(); it != timeDirs.end(); ++it)
{
// set run-time
recTime.setTime(*it, it->value());
// skip constant
if (recTime.timeName() == "constant")
{
continue;
}
// insert the time name into the hash-table with a continuous second index
timeIndexList.insert(recTime.timeName(), contTimeIndex);
if (verbose)
{
Info << "current time " << recTime.timeName() << endl;
Info << "insert " << recTime.timeName() << " , " << contTimeIndex << endl;
}
// insert the time value
timeValueList.insert(contTimeIndex, recTime.timeOutputValue());
// increment continuousIndex
contTimeIndex++;
if (verbose)
{
Info << "contTimeIndex " << contTimeIndex << endl;
}
}
if (verbose)
{
Info << endl;
Info << "Found " << label(timeDirs.size()) << " time folders" << endl;
Info << "Used " << (label(timeDirs.size())-1) << " time folders" << endl;
Info << "Found " << label(timeIndexList.size()) << " time steps" << endl;
}
// check time step of provided data
scalar dtCur(0.0);
scalar dtOld(0.0);
if (verbose)
{
Info << "timeValueList : " << timeValueList << endl;
}
forAll(timeValueList, i)
{
// compute time step
if (timeDirs[i].value() == timeDirs.last().value())
{
if (verbose)
{
Info << ".. leaving loop at " << timeDirs[i] << endl;
}
// leave loop
break;
}
if (verbose)
{
Info << "timeDirs.fcIndex(i)].value(), timeDirs[i].value() : "
<< timeDirs[timeDirs.fcIndex(i)].value() << " " << timeDirs[i].value()
<< endl;
}
// the documentation is in the code ;-)
// fcIndex() - return forward circular index, i.e. the next index
dtCur = timeDirs[timeDirs.fcIndex(i)].value() - timeDirs[i].value();
if (dtOld == 0.0)
{
dtOld = dtCur;
}
// do not use dtOld != dtCur as test, you will run into numerical trouble
// e.g. dtOld = 0.1 dtCur = 0.1 delta = 1.41553e-15
if (abs(dtOld - dtCur) > SMALL)
{
Info << "dtCur, dtOld = " << dtCur << " " << dtOld << endl;
FatalError << " in setting up data" << nl
<< " non-constant time-step of provided simulation data"
<< abort(FatalError);
}
}
// set deltaT
recTime.setDeltaT(dtCur, false);
if (verbose)
{
Info << "Setting deltaT to " << dtCur << endl;
Info << "Actual recTime.deltaT = " << recTime.deltaTValue() << endl;
Info << "Actual runTime.deltaT = " << runTime.deltaTValue() << endl;
}
// allocate space for flow fields
// this is not the final stage
// see TEqn.H or other corresponding files for the fields that actually needed
PtrList<volScalarField> alpha1pl(label(timeIndexList.size()));
PtrList<volScalarField> alpha2pl(label(timeIndexList.size()));
PtrList<volVectorField> U1pl(label(timeIndexList.size()));
PtrList<volVectorField> U2pl(label(timeIndexList.size()));
Info << "Reading fields" << endl;
for (it=timeDirs.begin(); it != timeDirs.end(); ++it)
{
// set time
recTime.setTime(*it, it->value());
// skip constant
if (recTime.timeName() == "constant")
{
continue;
}
if (verbose)
{
Info << "Reading at t = " << recTime.timeName() << endl;
}
alpha1pl.set
(
timeIndexList(recTime.timeName()),
new volScalarField
(
IOobject
(
"voidfraction",
recTime.timePath(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
)
);
IOobject headerAlpha2
(
"particlefraction",
recTime.timePath(),
mesh,
IOobject::NO_READ
);
if (headerAlpha2.headerOk())
{
alpha2pl.set
(
timeIndexList(recTime.timeName()),
new volScalarField
(
IOobject
(
"particlefraction",
recTime.timePath(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
)
);
}
else
{
alpha2pl.set
(
timeIndexList(recTime.timeName()),
new volScalarField
(
IOobject
(
"particlefraction",
recTime.timePath(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
scalar(1.0) - alpha1pl[timeIndexList(recTime.timeName())]
)
);
}
U1pl.set
(
timeIndexList(recTime.timeName()),
new volVectorField
(
IOobject
(
"U",
recTime.timePath(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
)
);
U2pl.set
(
timeIndexList(recTime.timeName()),
new volVectorField
(
IOobject
(
"Us",
recTime.timePath(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
)
);
}
// go to the earliest time step
it=timeDirs.begin();
recTime.setTime(*it, it->value());
// skip constant
if (recTime.timeName() == "constant")
{
it++;
recTime.setTime(*it, it->value());
}
Info<< "\nComputing recurrence matrix\n" << endl;
// create output file
std::ostringstream str_pid;
str_pid << pid();
const std::string my_pid(str_pid.str());
OFstream matrixFile("recurrenceMatrix."+my_pid);
// create recurrence matrix
SymmetricSquareMatrix<scalar> recurrenceMatrix(timeIndexList.size(), timeIndexList.size(), scalar(-666));
scalar maxElemVal(0.0);
// compute recurrence matrix elements
forAll(timeIndexList, ti)
{
forAll(timeIndexList, tj)
{
// main diagonal
if (ti == tj)
{
recurrenceMatrix[ti][tj] = 1;
continue;
}
// skip one half of the matrix
if (ti > tj)
{
continue;
}
// compute elements
recurrenceMatrix[ti][tj]
= sumSqr(alpha1pl[ti].internalField() - alpha1pl[tj].internalField())
+ sum(magSqr(U1pl[ti].internalField() - U1pl[tj].internalField()));
recurrenceMatrix[tj][ti] = recurrenceMatrix[ti][tj];
if (maxElemVal < recurrenceMatrix[ti][tj])
{
maxElemVal = recurrenceMatrix[ti][tj];
}
}
}
// normalize matrix elements
forAll(timeIndexList, ti)
{
forAll(timeIndexList, tj)
{
recurrenceMatrix[ti][tj] /= maxElemVal;
}
}
// write recurrence matrix
matrixFile << recurrenceMatrix;
// create time-sequence output file
OFstream timeSequenceFile("timeSequence."+my_pid);
// create working fields, these will change according to the recurrent rules
volScalarField alpha1 = alpha1pl[0];
volScalarField alpha2 = alpha2pl[0];
volVectorField U1 = U1pl[0];
volVectorField U2 = U2pl[0];
// write stuff ?!
//runTime.write();
// runTime.write does not work somehow. As of now, I have no idea why.
alpha2.write();
U2.write();

112
applications/solvers/cfdemSolverRecurrence/recurrenceTreatment.H
View File

@ -1,112 +0,0 @@
// if still in current sequence
// advance time = increment virtualTimeIndex
// do nothing else
// else
// generate new sequence
// random sequence length
// reset virtualTimeIndex
// jump to beginning of sequence
// if (in sequence)
if (virtualTimeIndex < sequenceStart + sequenceLength)
{
// do nothing
}
else
{
timeSequenceFile << "new sequence : ";
int startLoop = 0;
int endLoop = 0;
// search the other half of the recurrence matrix for
// the new starting point of the next sequence
if (virtualTimeIndex < recurrenceMatrix.n()/2)
{
startLoop = recurrenceMatrix.n()/2;
endLoop = recurrenceMatrix.n();
}
else
{
startLoop = 0;
endLoop = recurrenceMatrix.n()/2-1;
}
timeSequenceFile << "vti: " << virtualTimeIndex << " st " << startLoop << " end " << endLoop;
for (label j = startLoop; j < endLoop; j++)
{
if (recurrenceMatrix[j][virtualTimeIndex] < nextBestMinimum)
{
secondBestMinimum = nextBestMinimum;
nextBestMinimum = recurrenceMatrix[j][virtualTimeIndex];
sequenceStart2 = sequenceStart;
sequenceStart = j;
continue;
}
if (recurrenceMatrix[j][virtualTimeIndex] < secondBestMinimum)
{
secondBestMinimum = recurrenceMatrix[j][virtualTimeIndex];
sequenceStart2 = j;
}
}
// use next best or second best minimum
if (ranGen.scalar01() <= 0.5)
{
if (verbose)
{
Info << " ... using second best minimum" << endl;
}
sequenceStart = sequenceStartOld;
timeSequenceFile << " second best minimum : ";
}
if (sequenceStartOld == sequenceStart)
{
if (verbose)
{
Info << " ... using second best minimum to prevent starting again at the same time" << endl;
}
sequenceStart = sequenceStart2;
}
sequenceStartOld = sequenceStart;
// write to timeSequenceFile
timeSequenceFile << sequenceStart << endl;
// set new starting point
nextBestMinimum = GREAT;
secondBestMinimum = GREAT;
// generate new sequence length
sequenceLength = ranGen.integer(lowerSeqLim, upperSeqLim);
// be wordy?
if (verbose)
{
Info << "Found new starting point at i = " << sequenceStart << endl;
Info << "new sequence length = " << sequenceLength << endl;
}
// trim new sequence length
if (sequenceStart + sequenceLength > timeIndexList.size()-1)
{
//sequenceLength = sequenceLength - (sequenceStart + sequenceLength - timeIndexList.size());
sequenceLength = timeIndexList.size()-1 - sequenceStart; // simplified version of the above line
}
// reset virtual time index
virtualTimeIndex = sequenceStart;
} // end else if(in sequence)

26
src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.C
View File

@ -67,7 +67,7 @@ void cfdemCloudRec::giveDEMdata()
if(verbose_) Info << "giveDEMdata done." << endl;
}
void cfdemCloud::setForces()
void cfdemCloudRec::setForces()
{
// resetArray(fluidVel_,numberOfParticles(),3);
// resetArray(impForces_,numberOfParticles(),3);
@ -80,25 +80,22 @@ void cfdemCloud::setForces()
// * * * write top level fields * * * //
// * * * * * * * * * * * * * * * protected Member Functions * * * * * * * * * * * * * //
void cfdemCloud::setForces()
{
// resetArray(fluidVel_,numberOfParticles(),3);
// resetArray(impForces_,numberOfParticles(),3);
// resetArray(expForces_,numberOfParticles(),3);
// resetArray(DEMForces_,numberOfParticles(),3);
// resetArray(Cds_,numberOfParticles(),1);
// for (int i=0;i<cfdemCloud::nrForceModels();i++) cfdemCloud::forceM(i).setForce();
}
// * * * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * * //
void cfdemCloudRec::initRecFields()
{
recM().initRecFields();
recModel_->initRecFields();
}
void cfdemCloudRec::updateRecFields()
{
recM().updateRecFields();
recModel_->updateRecFields();
}
void cfdemCloudRec::writeRecFields()
{
recModel_->writeRecFields();
}
bool cfdemCloudRec::evolve
@ -215,11 +212,6 @@ bool cfdemCloudRec::evolve
dataExchangeM().couple(1);
}
if(verbose_){
#include "debugInfo.H"
}
clockM().start(25,"dumpDEMdata");
// do particle IO
IOM().dumpDEMdata();

5
src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.H
View File

@ -90,8 +90,9 @@ public:
return recModel_;
}
void initRecFields(volScalarField&,volScalarField&,volVectorField&,volVectorField&);
void updateRecFields(volScalarField&,volScalarField&,volVectorField&,volVectorField&);
void initRecFields();
void updateRecFields();
void writeRecFields();
};

44
src/lagrangian/cfdemParticle/subModels/recModel/recModel/recModel.C
View File

@ -57,31 +57,33 @@ Foam::recModel::recModel
IOobject
(
"controlDict",
system(),
*this,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
sm.mesh().time().system(),
sm.mesh(),
IOobject::MUST_READ,
IOobject::NO_WRITE
)
),
verbose_(sm.verbose()),
recTime("dataBase", "", "../system", "../constant", false),
timeDirs(recTime.times()),
timeIndexList(label(timeDirs.size())-1),
timeValueList(label(timeDirs.size())-1),
numRecFields(label(timeDirs.size())),
timeIndexList(numRecFields-1),
timeValueList(numRecFields-1),
recurrenceMatrix(numRecFields,numRecFields,scalar(0)),
contTimeIndex(0),
sequenceStart(0),
sequenceEnd(0),
virtualTimeIndex(0),
virtualStartIndex(0),
startTime_(controlDict_.lookup("startTime")),
endTime_(controlDict_.lookup("endTime")),
timeStep_(controlDict_.lookup("deltaT")),
startTime_(readScalar(controlDict_.lookup("startTime"))),
endTime_(readScalar(controlDict_.lookup("endTime"))),
timeStep_(readScalar(controlDict_.lookup("deltaT"))),
virtualTimeIndexList(0),
virtualTimeIndexListPos(0),
lowerSeqLim(max(1, label(timeIndexList.size()/20))),
upperSeqLim(label(timeIndexList.size()/5))
lowerSeqLim(max(1, label(numRecFields/20))),
upperSeqLim(label(numRecFields/5))
{
if (verbose)
if (verbose_)
{
// be informative on properties of the "recTime" Time-object
Info << "recTime.rootPath() " << recTime.rootPath() << endl;
@ -103,8 +105,6 @@ Foam::recModel::~recModel()
// * * * * * * * * * * * * * public Member Functions * * * * * * * * * * * * //
virtual void Foam::recModel::updateRecFields()
{}
// * * * * * * * * * * * * * private Member Functions * * * * * * * * * * * * //
@ -128,7 +128,7 @@ void Foam::recModel::readTimeSeries()
timeIndexList.insert(recTime.timeName(), contTimeIndex);
if (verbose)
if (verbose_)
{
Info << "current time " << recTime.timeName() << endl;
Info << "insert " << recTime.timeName() << " , " << contTimeIndex << endl;
@ -141,13 +141,13 @@ void Foam::recModel::readTimeSeries()
// increment continuousIndex
contTimeIndex++;
if (verbose)
if (verbose_)
{
Info << "contTimeIndex " << contTimeIndex << endl;
}
}
if (verbose)
if (verbose_)
{
Info << endl;
Info << "Found " << label(timeDirs.size()) << " time folders" << endl;
@ -162,7 +162,7 @@ scalar Foam::recModel::checkTimeStep()
scalar dtCur(0.0);
scalar dtOld(0.0);
if (verbose)
if (verbose_)
{
Info << "timeValueList : " << timeValueList << endl;
}
@ -172,7 +172,7 @@ scalar Foam::recModel::checkTimeStep()
// compute time step
if (timeDirs[i].value() == timeDirs.last().value())
{
if (verbose)
if (verbose_)
{
Info << ".. leaving loop at " << timeDirs[i] << endl;
}
@ -180,7 +180,7 @@ scalar Foam::recModel::checkTimeStep()
break;
}
if (verbose)
if (verbose_)
{
Info << "timeDirs.fcIndex(i)].value(), timeDirs[i].value() : "
<< timeDirs[timeDirs.fcIndex(i)].value() << " " << timeDirs[i].value()
@ -208,7 +208,7 @@ scalar Foam::recModel::checkTimeStep()
// set deltaT
recTime.setDeltaT(dtCur, false);
if (verbose)
if (verbose_)
{
Info << "Setting deltaT to " << dtCur << endl;
Info << "Actual recTime.deltaT = " << recTime.deltaTValue() << endl;

9
src/lagrangian/cfdemParticle/subModels/recModel/recModel/recModel.H
View File

@ -27,7 +27,7 @@ License
#include "cfdemCloud.H"
//#include "dataExchangeModel.H"
#include "fvCFD.H"
#include "HashTable.H"
#include "labelPair.H"
#include <mpi.h>
@ -45,8 +45,9 @@ protected:
// Protected data
const dictionary& dict_;
const dictionary& controlDict_;
cfdemCloud& particleCloud_;
const IOdictionary& controlDict_;
bool verbose_;
SymmetricSquareMatrix<scalar> recurrenceMatrix;
Foam::Time recTime;
@ -75,6 +76,7 @@ protected:
scalar timeStep_;
scalar recTimeStep_;
label numRecFields;
label totRecSteps;
label virtualTimeIndex;
@ -89,7 +91,7 @@ protected:
private:
void readTimeSeries();
void checkTimeStep();
scalar checkTimeStep();
public:
@ -139,6 +141,7 @@ public:
virtual void initRecFields() = 0;
virtual void updateRecFields() = 0;
virtual void writeRecFields() = 0;
scalar recTimeStep()
{
return recTimeStep_;

65
src/lagrangian/cfdemParticle/subModels/recModel/standardRecModel/standardRecModel.C
View File

@ -23,7 +23,7 @@ License
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "IOModel.H"
#include "Random.H"
#include "standardRecModel.H"
#include "addToRunTimeSelectionTable.H"
@ -55,12 +55,12 @@ standardRecModel::standardRecModel
:
recModel(dict,sm),
propsDict_(dict.subDict(typeName + "Props")),
velFieldName_(propsDict_.lookup("velRecFieldName")),
URec_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
UFieldName_(propsDict_.lookup("velRecFieldName")),
voidfractionFieldName_(propsDict_.lookup("voidfractionRecFieldName")),
voidfractionRec_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
UsFieldName_(propsDict_.lookup("granVelRecFieldName")),
UsRec_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
voidfractionRecpl(numRecFields),
URecpl(numRecFields),
UsRecpl(numRecFields)
{
readFieldSeries();
computeRecMatrix();
@ -73,7 +73,7 @@ standardRecModel::standardRecModel
computeRecPath();
int listSizeBytes = 2*sizeof(sequenceStart)*virtualTimeIndexList.size();
MPI_Bcast(&virtualTimeIndeList[0], listSizeBytes, MPI_BYTE, root, MPI_COMM_WORLD);
MPI_Bcast(&virtualTimeIndexList[0], listSizeBytes, MPI_BYTE, root, MPI_COMM_WORLD);
sequenceStart=virtualTimeIndexList[0].first();
sequenceEnd=virtualTimeIndexList[0].second();
@ -91,11 +91,6 @@ standardRecModel::~standardRecModel()
void standardRecModel::initRecFields()
{
voidfractionRec_ = voidfractionRecpl[0];
URec = URecpl[0];
UsRec = UsRecpl[0];
correctBC?
}
void standardRecModel::updateRecFields()
@ -108,11 +103,15 @@ void standardRecModel::updateRecFields()
sequenceEnd=virtualTimeIndexList[virtualTimeIndexListPos].second();
virtualTimeIndex=sequenceStart;
}
voidfractionRec_ = voidfractionRecpl[virtualTimeIndex];
URec_ = URecpl[virtualTimeIndex];
UsRec_ = UsRecpl[virtualTimeIndex];
}
void standardRecModel::writeRecFields()
{
// need to check if this writes to the correct destination
correctBC?
// voidfractionRecpl[virtualTimeIndex].write();
// URecpl[virtualTimeIndex].write();
// UsRecpl[virtualTimeIndex].write();
}
@ -131,7 +130,7 @@ void standardRecModel::readFieldSeries()
continue;
}
if (verbose)
if (verbose_)
{
Info << "Reading at t = " << recTime.timeName() << endl;
}
@ -143,13 +142,13 @@ void standardRecModel::readFieldSeries()
(
IOobject
(
"voidfraction",
voidfractionFieldName_,
recTime.timePath(),
mesh,
particleCloud_.mesh(),
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
particleCloud_.mesh()
)
);
@ -160,13 +159,13 @@ void standardRecModel::readFieldSeries()
(
IOobject
(
"U",
UFieldName_,
recTime.timePath(),
mesh,
particleCloud_.mesh(),
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
particleCloud_.mesh()
)
);
UsRecpl.set
@ -176,23 +175,23 @@ void standardRecModel::readFieldSeries()
(
IOobject
(
"Us",
UsFieldName_,
recTime.timePath(),
mesh,
particleCloud_.mesh(),
IOobject::MUST_READ,
IOobject::NO_WRITE
),
mesh
particleCloud_.mesh()
)
);
}
}
}
void standardRecModel::computeRecMatrix()
{
Info<< "\nComputing recurrence matrix\n" << endl;
recurrenceMatrix(timeIndexList.size(), timeIndexList.size(), scalar(0));
recurrenceMatrix=0.0;
scalar maxElemVal(0.0);
// compute recurrence matrix elements
@ -215,8 +214,7 @@ void standardRecModel::computeRecMatrix()
// compute elements
recurrenceMatrix[ti][tj]
= sumSqr(alpha1pl[ti].internalField() - alpha1pl[tj].internalField())
+ sum(magSqr(U1pl[ti].internalField() - U1pl[tj].internalField()));
= sumSqr(voidfractionRecpl[ti].internalField() - voidfractionRecpl[tj].internalField());
recurrenceMatrix[tj][ti] = recurrenceMatrix[ti][tj];
@ -239,15 +237,16 @@ void standardRecModel::computeRecMatrix()
void standardRecModel::computeRecPath()
{
// random recurrence stuff
#include "Random.H" // random numbers
Random ranGen(osRandomInteger());
label virtualTimeIndex=0;
label recSteps=0;
label seqStart=0;
label seqLength=0;
labelPair seqStartEnd(0,0);
label seqLength=ranGen.integer(lowerSeqLim, upperSeqLim);
labelPair seqStartEnd(seqStart,seqStart+seqLength-1);
virtualTimeIndexList.append(seqStartEnd);
virtualTimeIndex = seqStart+seqLength-1;
recSteps+=seqLength;
while(recSteps<=totRecSteps)
{

17
src/lagrangian/cfdemParticle/subModels/recModel/standardRecModel/standardRecModel.H
View File

@ -26,6 +26,7 @@ License
#define standardRecModel_H
#include "recModel.H"
//#include "PtrList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -68,23 +69,19 @@ public:
void initRecFields();
void updateRecFields();
void writeRecFields();
private:
dictionary propsDict_;
word velFieldName_;
const volVectorField& URec_;
word UFieldName_;
word voidfractionFieldName_;
const volScalarField& voidfractionRec_;
word UsFieldName_;
const volVectorField& UsRec_;
word UsFieldName_;
PtrList<volScalarField> voidfractionRecpl(label(timeIndexList.size()));
PtrList<volVectorField> URecpl(label(timeIndexList.size()));
PtrList<volVectorField> UsRecpl(label(timeIndexList.size()));
PtrList<volScalarField> voidfractionRecpl;
PtrList<volVectorField> URecpl;
PtrList<volVectorField> UsRecpl;
};