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

Runs. Sometimes.

Browse Source
This commit is contained in:
Thomas Lichtenegger
2015-10-13 12:55:43 +02:00
parent b55ba48801
commit c38449cacd
7 changed files with 88 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@ -59,6 +59,14 @@ cfdemCloudRec::cfdemCloudRec
forceModels_[i]
);
}
if(nrForceModels()>1)
{
coupleRecForce_=true;
}
else
{
coupleRecForce_=false;
}
}
// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * //
@ -76,7 +84,8 @@ void cfdemCloudRec::getDEMdata()
void cfdemCloudRec::giveDEMdata()
{
dataExchangeM().giveData("vrec","vector-atom",fluidVel_);
dataExchangeM().giveData("dragforce","vector-atom",DEMForces_);
if(coupleRecForce_)
dataExchangeM().giveData("dragforce","vector-atom",DEMForces_);
if(verbose_) Info << "giveDEMdata done." << endl;
}

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

@ -65,6 +65,8 @@ protected:
void giveDEMdata();
void setForces();
bool coupleRecForce_;
public:

2
src/lagrangian/cfdemParticle/subModels/forceModelRec/forceModelRec/forceModelRec.C
View File

@ -64,6 +64,8 @@ void forceModelRec::partToArray
const vector& Ufluid
) const
{
// need to make sure that force models don't overwrite velocities
for(int j=0;j<3;j++)
fluidVel()[index][j]=Ufluid[j];

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

@ -67,6 +67,7 @@ Foam::recModel::recModel
recTime("dataBase", "", "../system", "../constant", false),
timeDirs(recTime.times()),
numRecFields(label(timeDirs.size())),
recurrenceMatrixLocal(numRecFields,numRecFields,scalar(0)),
recurrenceMatrix(numRecFields,numRecFields,scalar(0)),
timeIndexList(numRecFields-1),
timeValueList(numRecFields-1),
@ -218,13 +219,31 @@ scalar Foam::recModel::checkTimeStep()
void Foam::recModel::writeRecMatrix() const
{
// create output file
std::ostringstream str_pid;
// version for all processors
/* std::ostringstream str_pid;
str_pid << pid();
const std::string my_pid(str_pid.str());
OFstream matrixFile("recurrenceMatrix."+my_pid);
// write recurrence matrix
matrixFile << recurrenceMatrix;
*/
// write recurrence matrix
OFstream matrixFile("recurrenceMatrix");
matrixFile << recurrenceMatrix;
}
void Foam::recModel::writeRecPath() const
{
// version for all processors
/*
std::ostringstream str_pid;
str_pid << pid();
const std::string my_pid(str_pid.str());
OFstream listFile("recurrencePath."+my_pid);
*/
// write recurrence matrix
OFstream listFile("recurrencePath");
listFile << virtualTimeIndexList;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@ -51,6 +51,7 @@ protected:
Foam::Time recTime;
instantList timeDirs;
label numRecFields;
SymmetricSquareMatrix<scalar> recurrenceMatrixLocal;
SymmetricSquareMatrix<scalar> recurrenceMatrix;
// create a data structure for the time indices
@ -152,6 +153,7 @@ public:
return recTimeStep_;
}
void writeRecMatrix() const;
void writeRecPath() const;
};

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

@ -64,21 +64,39 @@ standardRecModel::standardRecModel
UsRecpl(numRecFields)
{
readFieldSeries();
computeRecMatrix();
// make sure each processor has the same sequence of fields
int root=0;
int rank;
// MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
computeRecMatrix();
if(verbose_)
Info << "\nSumming recurrence matrix over all processors.\n" << endl;
for(int i = 0; i < numRecFields; ++i)
MPI_Allreduce
(
&recurrenceMatrixLocal[i][0],
&recurrenceMatrix[i][0],
numRecFields,
MPI_DOUBLE,
MPI_SUM,
MPI_COMM_WORLD
);
computeRecPath();
// int listSizeBytes = 2*sizeof(sequenceStart)*virtualTimeIndexList.size();
// MPI_Bcast(&virtualTimeIndexList[0], listSizeBytes, MPI_BYTE, root, MPI_COMM_WORLD);
int listSizeBytes = 2*sizeof(sequenceStart)*virtualTimeIndexList.size();
MPI_Bcast(&virtualTimeIndexList[0], listSizeBytes, MPI_BYTE, root, MPI_COMM_WORLD);
sequenceStart=virtualTimeIndexList[0].first();
sequenceEnd=virtualTimeIndexList[0].second();
virtualTimeIndex=sequenceStart;
writeRecMatrix();
writeRecPath();
}
@ -230,13 +248,11 @@ void standardRecModel::computeRecMatrix()
forAll(timeIndexList, tj)
{
if(verbose_)
{
Info<<"\n Doing calculation for element " << ti << " " << tj << "\n" << endl;
}
Info<<"\n Doing calculation for element " << ti << " " << tj << "\n" << endl;
// main diagonal
if (ti == tj)
{
recurrenceMatrix[ti][tj] = 1;
recurrenceMatrixLocal[ti][tj] = 0;
continue;
}
@ -247,14 +263,14 @@ void standardRecModel::computeRecMatrix()
}
// compute elements
recurrenceMatrix[ti][tj]
recurrenceMatrixLocal[ti][tj]
= sumSqr(voidfractionRecpl[ti].internalField() - voidfractionRecpl[tj].internalField());
recurrenceMatrix[tj][ti] = recurrenceMatrix[ti][tj];
recurrenceMatrixLocal[tj][ti] = recurrenceMatrixLocal[ti][tj];
// if (maxElemVal < recurrenceMatrix[ti][tj])
// if (maxElemVal < recurrenceMatrixLocal[ti][tj])
// {
// maxElemVal = recurrenceMatrix[ti][tj];
// maxElemVal = recurrenceMatrixLocal[ti][tj];
// }
}
}
@ -264,7 +280,7 @@ void standardRecModel::computeRecMatrix()
//{
// forAll(timeIndexList, tj)
// {
// recurrenceMatrix[ti][tj] /= maxElemVal;
// recurrenceMatrixLocal[ti][tj] /= maxElemVal;
// }
//}
Info<< "\nComputing recurrence matrix done\n" << endl;
@ -279,24 +295,23 @@ void standardRecModel::computeRecPath()
label recSteps=0;
label seqStart=0;
label seqLength=ranGen.integer(lowerSeqLim, upperSeqLim);
labelPair seqStartEnd(seqStart,seqStart+seqLength-1);
virtualTimeIndex=seqEnd(seqStart,seqLength);
labelPair seqStartEnd(seqStart,virtualTimeIndex);
virtualTimeIndexList.append(seqStartEnd);
virtualTimeIndex = seqStart+seqLength-1;
recSteps+=seqLength;
while(recSteps<=totRecSteps)
{
label startLoop = 0;
label endLoop = 0;
seqLength = ranGen.integer(lowerSeqLim, upperSeqLim);
scalar nextMinimum(GREAT);
// 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()-1;
endLoop--; // start of next sequence one snapshot AFTER minimum position
}
else
{
@ -304,6 +319,8 @@ void standardRecModel::computeRecPath()
endLoop = recurrenceMatrix.n()/2-1;
}
scalar nextMinimum(GREAT);
for (label j = startLoop; j < endLoop; j++)
{
if (recurrenceMatrix[j][virtualTimeIndex] < nextMinimum)
@ -314,19 +331,21 @@ void standardRecModel::computeRecPath()
}
}
// trim new sequence length
if (seqStart + seqLength > timeIndexList.size()-1)
{
seqLength = timeIndexList.size()-1 - seqStart;
}
labelPair seqStartEnd(seqStart,seqStart+seqLength-1);
seqLength = ranGen.integer(lowerSeqLim, upperSeqLim);
virtualTimeIndex=seqEnd(seqStart,seqLength);
labelPair seqStartEnd(seqStart,virtualTimeIndex);
virtualTimeIndexList.append(seqStartEnd);
virtualTimeIndex = seqStart+seqLength-1;
recSteps+=seqLength;
}
Info<< "\nComputing recurrence path done\n" << endl;
}
label standardRecModel::seqEnd(label seqStart, label & seqLength)
{
if(seqStart+seqLength>numRecFields-1)
seqLength=numRecFields-1-seqStart;
return seqStart+seqLength;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

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

@ -89,6 +89,8 @@ private:
PtrList<volScalarField> voidfractionRecpl;
PtrList<volVectorField> URecpl;
PtrList<volVectorField> UsRecpl;
label seqEnd(label, label&);
};