From c38449cacdc219b060ccfec9fdbb9c01e96bdffa Mon Sep 17 00:00:00 2001 From: Thomas Lichtenegger Date: Tue, 13 Oct 2015 12:55:43 +0200 Subject: [PATCH] Runs. Sometimes. --- .../derived/cfdemCloudRec/cfdemCloudRec.C | 11 ++- .../derived/cfdemCloudRec/cfdemCloudRec.H | 2 + .../forceModelRec/forceModelRec.C | 2 + .../subModels/recModel/recModel/recModel.C | 27 ++++++- .../subModels/recModel/recModel/recModel.H | 2 + .../standardRecModel/standardRecModel.C | 75 ++++++++++++------- .../standardRecModel/standardRecModel.H | 2 + 7 files changed, 88 insertions(+), 33 deletions(-) diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.C b/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.C index 3e4230fc..801b3375 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.C +++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.C @@ -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; } diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.H b/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.H index d58c31b4..02d826cb 100644 --- a/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.H +++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.H @@ -65,6 +65,8 @@ protected: void giveDEMdata(); void setForces(); + + bool coupleRecForce_; public: diff --git a/src/lagrangian/cfdemParticle/subModels/forceModelRec/forceModelRec/forceModelRec.C b/src/lagrangian/cfdemParticle/subModels/forceModelRec/forceModelRec/forceModelRec.C index 7382fccf..23fa9eba 100644 --- a/src/lagrangian/cfdemParticle/subModels/forceModelRec/forceModelRec/forceModelRec.C +++ b/src/lagrangian/cfdemParticle/subModels/forceModelRec/forceModelRec/forceModelRec.C @@ -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]; diff --git a/src/lagrangian/cfdemParticle/subModels/recModel/recModel/recModel.C b/src/lagrangian/cfdemParticle/subModels/recModel/recModel/recModel.C index bdad5279..96fb8546 100644 --- a/src/lagrangian/cfdemParticle/subModels/recModel/recModel/recModel.C +++ b/src/lagrangian/cfdemParticle/subModels/recModel/recModel/recModel.C @@ -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; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/cfdemParticle/subModels/recModel/recModel/recModel.H b/src/lagrangian/cfdemParticle/subModels/recModel/recModel/recModel.H index 8b30b109..f073153f 100644 --- a/src/lagrangian/cfdemParticle/subModels/recModel/recModel/recModel.H +++ b/src/lagrangian/cfdemParticle/subModels/recModel/recModel/recModel.H @@ -51,6 +51,7 @@ protected: Foam::Time recTime; instantList timeDirs; label numRecFields; + SymmetricSquareMatrix recurrenceMatrixLocal; SymmetricSquareMatrix recurrenceMatrix; // create a data structure for the time indices @@ -152,6 +153,7 @@ public: return recTimeStep_; } void writeRecMatrix() const; + void writeRecPath() const; }; diff --git a/src/lagrangian/cfdemParticle/subModels/recModel/standardRecModel/standardRecModel.C b/src/lagrangian/cfdemParticle/subModels/recModel/standardRecModel/standardRecModel.C index a35cb7f8..e731468a 100644 --- a/src/lagrangian/cfdemParticle/subModels/recModel/standardRecModel/standardRecModel.C +++ b/src/lagrangian/cfdemParticle/subModels/recModel/standardRecModel/standardRecModel.C @@ -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 diff --git a/src/lagrangian/cfdemParticle/subModels/recModel/standardRecModel/standardRecModel.H b/src/lagrangian/cfdemParticle/subModels/recModel/standardRecModel/standardRecModel.H index e9b8455d..9edc1659 100644 --- a/src/lagrangian/cfdemParticle/subModels/recModel/standardRecModel/standardRecModel.H +++ b/src/lagrangian/cfdemParticle/subModels/recModel/standardRecModel/standardRecModel.H @@ -89,6 +89,8 @@ private: PtrList voidfractionRecpl; PtrList URecpl; PtrList UsRecpl; + + label seqEnd(label, label&); };