fix calculation of timeStepFraction & voidfraction interpolation

merged from PUBLIC

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -81,6 +81,7 @@ cfdemCloud::cfdemCloud
     solveFlow_(true),
     verbose_(false),
     ignore_(false),
+    allowCFDsubTimestep_(true),
     limitDEMForces_(false),
     modelType_(couplingProperties_.lookup("modelType")),
     positions_(NULL),
@@ -251,17 +252,6 @@ cfdemCloud::cfdemCloud
     else
         Info << "ignoring ddt(voidfraction)" << endl;
 
-    forceModel_ = new autoPtr<forceModel>[nrForceModels()];
-    for (int i=0;i<nrForceModels();i++)
-    {
-        forceModel_[i] = forceModel::New
-        (
-            couplingProperties_,
-            *this,
-            forceModels_[i]
-        );
-    }
-
     momCoupleModel_ = new autoPtr<momCoupleModel>[momCoupleModels_.size()];
     for (int i=0;i<momCoupleModels_.size();i++)
     {
@@ -273,6 +263,17 @@ cfdemCloud::cfdemCloud
         );
     }
 
+    forceModel_ = new autoPtr<forceModel>[nrForceModels()];
+    for (int i=0;i<nrForceModels();i++)
+    {
+        forceModel_[i] = forceModel::New
+        (
+            couplingProperties_,
+            *this,
+            forceModels_[i]
+        );
+    }
+
     // run liggghts commands from cfdem
     liggghtsCommand_ = new autoPtr<liggghtsCommandModel>[liggghtsCommandModelList_.size()];
     for (int i=0;i<liggghtsCommandModelList_.size();i++)
@@ -324,7 +325,14 @@ cfdemCloud::cfdemCloud
     {
         checkPeriodicCells_ = true;
     }
-    else if (nPatchesCyclic > 0 && nPatchesNonCyclic > 0)
+
+    //hard set checkperiodic cells if wished
+    if(this->couplingProperties().found("checkPeriodicCells"))
+    {
+        checkPeriodicCells_ = couplingProperties().lookupOrDefault<Switch>("checkPeriodicCells", checkPeriodicCells_);
+    }
+
+    if (nPatchesCyclic > 0 && nPatchesNonCyclic > 0)
     {
         if (verbose_) Info << "nPatchesNonCyclic=" << nPatchesNonCyclic << ", nPatchesCyclic=" << nPatchesCyclic << endl;
         Warning << "Periodic handing is disabled because the domain is not fully periodic!\n" << endl;
@@ -612,6 +620,10 @@ bool cfdemCloud::evolve
         //      IMPLICIT FORCE CONTRIBUTION AND SOLVER USE EXACTLY THE SAME AVERAGED
         //      QUANTITIES AT THE GRID!
         Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl;
+        if(dataExchangeM().timeStepFraction() > 1.0000001)
+        {
+            FatalError << "cfdemCloud::dataExchangeM().timeStepFraction()>1: Do not do this, since dangerous. This might be due to the fact that you used a adjustable CFD time step. Please use a fixed CFD time step." << abort(FatalError);
+        }
         clockM().start(24,"interpolateEulerFields");
 
         // update voidFractionField

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -92,6 +92,8 @@ protected:
 
     bool ignore_;
 
+    bool allowCFDsubTimestep_;
+
     bool limitDEMForces_;
 
     scalar maxDEMForce_;
@@ -226,6 +228,10 @@ public:
 // public Member Functions
 
     // Access
+       bool allowCFDsubTimestep() { return allowCFDsubTimestep_; }
+
+       void setAllowCFDsubTimestep(bool b) { allowCFDsubTimestep_ = b; }
+
        void checkCG(bool);
 
        void setPos(double **&);

src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C

@@ -328,20 +328,12 @@ void averagingModel::undoWeightFields(double**const& mask) const
 
 tmp<volVectorField> averagingModel::UsInterp() const
 {
-    if (particleCloud_.dataExchangeM().couplingStep() > 1)
-    {
-        return tmp<volVectorField>
-        (
-            new volVectorField("Us_averagingModel", (1. - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_ + particleCloud_.dataExchangeM().timeStepFraction() * UsNext_)
-        );
-    }
-    else
-    {
-        return tmp<volVectorField>
-        (
-            new volVectorField("Us_averagingModel", UsNext_)
-        );
-    }
+    const scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
+
+    return tmp<volVectorField>
+    (
+        new volVectorField("Us_averagingModel", (1. - tsf) * UsPrev_ + tsf * UsNext_)
+    );
 }
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.C

@@ -209,12 +209,7 @@ bool dataExchangeModel::couple(int i) const
 scalar dataExchangeModel::timeStepFraction() const
 {
     //return fraction between previous coupling TS and actual TS
-    //scalar DEMtime = DEMts_ * couplingInterval_;
-    //scalar frac = ( ( particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value() ) - (couplingStep_) * DEMtime) / DEMtime; //Chr 05.03.2013
-    scalar frac = ( particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value() - couplingStep_ * couplingTime() ) / couplingTime();
-    if (frac < 1e-4) frac = 1.;
-
-    return frac;
+    return ( particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value() - (couplingStep_-1) * couplingTime() ) / couplingTime();
 }
 
 int dataExchangeModel::getNumberOfParticles() const

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/dataExchangeModel/dataExchangeModel.H

@@ -198,22 +198,29 @@ public:
 
         inline void checkTSsize() const
         {
-            if(particleCloud_.mesh().time().deltaT().value() >  couplingInterval_ * DEMts_ + SMALL)
+            if (particleCloud_.mesh().time().deltaT().value() >  couplingInterval_ * DEMts_ + SMALL)
             {
                 Info << "particleCloud_.mesh().time().deltaT().value() = " << particleCloud_.mesh().time().deltaT().value() << endl;
                 Info << "couplingInterval_ = " << couplingInterval_ << endl;
                 Info << "DEMts_ = " << DEMts_ << endl;
-                FatalError<<"\nError - TS bigger than coupling interval!\n"<< abort(FatalError);
+                FatalError << "\nError - CFD time-step bigger than coupling time (= DEM time step * coupling interval)!\n" << abort(FatalError);
             }
+            if (std::fabs((round(couplingTime()/particleCloud_.mesh().time().deltaT().value())*particleCloud_.mesh().time().deltaT().value())-couplingTime()) > SMALL)
+            {
+                Info << "particleCloud_.mesh().time().deltaT().value() = " << particleCloud_.mesh().time().deltaT().value() << endl;
+                Info << "couplingInterval_ = " << couplingInterval_ << endl;
+                Info << "DEMts_ = " << DEMts_ << endl;
+                Warning << "\nWarning - Coupling time (= DEM time step * coupling interval) is not a multiple of  CFD time-step!\n" << endl;
+            }
+            if (!particleCloud_.allowCFDsubTimestep())
+                if (particleCloud_.mesh().time().deltaT().value() < couplingInterval_ * DEMts_ + SMALL)
+                    FatalError << "\nYour models require: CFD time-step = coupling interval (= DEM time step * coupling interval)! \n" << abort(FatalError);
+
+            // warn if sub-TS
+            if (particleCloud_.mesh().time().deltaT().value() < couplingTime() - SMALL)
+                Warning << "You are using sub-time-steps (i.e. CFD TS < coupling time)! Check your settings properly." << endl;
         }
 
-        /*inline bool checkExactTiming() const
-        {
-            return false;
-        }*/
-
-        //void checkNClumpTypes() const {};
-
         inline void readDEMtsfromDict(dictionary& propsDict)
         {
             DEMts_ = readScalar(propsDict.lookup("DEMts"));
@@ -222,10 +229,8 @@ public:
 
         inline bool doCoupleNow() const
         {
-            if (particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value()
-             - ((1+couplingStep_)*(DEMts_*couplingInterval_))
-             + 1e-10 > 0)                                         // Chr 27.03.2013 : first coupling after DEMts_*couplingInterval_
-//             > particleCloud_.mesh().time().deltaT().value()/2) // Chr 27.03.2013
+            if (particleCloud_.mesh().time().value()-particleCloud_.mesh().time().startTime().value()-SMALL
+                > couplingStep_*DEMts_*couplingInterval_)
             {
                 return true;
             }

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/liggghtsCommandModel/liggghtsCommandModel.C

@@ -104,6 +104,13 @@ void liggghtsCommandModel::checkTimeMode(dictionary& propsDict)
             }
         }
     }
+
+    if(verbose_){
+        Info << "runFirst = " << runFirst_ << endl;
+        Info << "runLast = " << runLast_ << endl;
+        Info << "runEveryCouplingStep = " << runEveryCouplingStep_ << endl;
+        Info << "runEveryWriteStep = " << runEveryWriteStep_ << endl;
+    }
 }
 
 void liggghtsCommandModel::checkTimeSettings(const dictionary& propsDict)
@@ -112,11 +119,12 @@ void liggghtsCommandModel::checkTimeSettings(const dictionary& propsDict)
     {
         scalar DEMts = particleCloud_.dataExchangeM().DEMts();
         scalar couplingInterval = particleCloud_.dataExchangeM().couplingInterval();
+        scalar simStartTime = particleCloud_.mesh().time().startTime().value();
 
         if(runLast_) // last run
         {
             // read time options from subdict
-            endTime_ = particleCloud_.mesh().time().endTime().value()-particleCloud_.mesh().time().startTime().value();
+            endTime_ = particleCloud_.mesh().time().endTime().value()-simStartTime;
             startTime_ = endTime_;
             timeInterval_ = -1;
 
@@ -136,9 +144,9 @@ void liggghtsCommandModel::checkTimeSettings(const dictionary& propsDict)
 
                 // calculate coupling times
                 // if this makes troubles try floor((startTime_+SMALL)/.. as above
-                firstCouplingStep_ = floor(startTime_/DEMts/couplingInterval)+1;
-                lastCouplingStep_ = floor(endTime_/DEMts/couplingInterval)+1;
-                couplingStepInterval_ = floor(timeInterval_/DEMts/couplingInterval)+1;
+                firstCouplingStep_ = floor((startTime_+SMALL-simStartTime)/DEMts/couplingInterval);
+                lastCouplingStep_ = floor((endTime_+SMALL-simStartTime)/DEMts/couplingInterval);
+                couplingStepInterval_ = floor(timeInterval_+SMALL/DEMts/couplingInterval);
             }
             else      //runEveryCouplingStep  or writeStep
             {

src/lagrangian/cfdemParticle/subModels/momCoupleModel/explicitCouple/explicitCouple.C

@@ -108,9 +108,13 @@ explicitCouple::~explicitCouple()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 tmp<volVectorField> explicitCouple::expMomSource() const
 {
-    scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
+    const scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
 
-    if (1. - tsf < 1e-4) //tsf==1
+    // calc Ksl
+
+    // update KslNext in first subTS
+    // NOTE: without following if we could update f every subTS (based on current values) and use this value
+    if(tsf < particleCloud_.mesh().time().deltaT().value()/particleCloud_.dataExchangeM().couplingTime() + 0.000001 )
     {
         // calc fNext
         forAll(fNext_,cellI)
@@ -124,18 +128,12 @@ tmp<volVectorField> explicitCouple::expMomSource() const
                 if (magF > fLimit_[i]) fNext_[cellI][i] *= fLimit_[i]/magF;
             }
         }
-        return tmp<volVectorField>
-        (
-            new volVectorField("f_explicitCouple", fPrev_)
-        );
-    }
-    else
-    {
-        return tmp<volVectorField>
-        (
-            new volVectorField("f_explicitCouple", (1. - tsf) * fPrev_ + tsf * fNext_)
-        );
     }
+
+    return tmp<volVectorField>
+    (
+        new volVectorField("f_explicitCouple", (1. - tsf) * fPrev_ + tsf * fNext_)
+    );
 }
 
 void explicitCouple::resetMomSourceField() const

src/lagrangian/cfdemParticle/subModels/momCoupleModel/implicitCouple/implicitCouple.C

@@ -122,11 +122,13 @@ implicitCouple::~implicitCouple()
 
 tmp<volScalarField> implicitCouple::impMomSource() const
 {
-    scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
+    const scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
 
     // calc Ksl
 
-    if (1. - tsf < 1e-4) //tsf==1
+    // update KslNext in first subTS
+    // NOTE: without following if we could update Ksl every subTS (based on current values) and use this value
+    if(tsf < particleCloud_.mesh().time().deltaT().value()/particleCloud_.dataExchangeM().couplingTime() + 0.000001 )
     {
         scalar Ur;
 
@@ -145,18 +147,12 @@ tmp<volScalarField> implicitCouple::impMomSource() const
             // limiter
             if (KslNext_[cellI] > KslLimit_) KslNext_[cellI] = KslLimit_;
         }
-        return tmp<volScalarField>
-        (
-            new volScalarField("Ksl_implicitCouple", KslPrev_)
-        );
-    }
-    else
-    {
-        return tmp<volScalarField>
-        (
-            new volScalarField("Ksl_implicitCouple", (1. - tsf) * KslPrev_ + tsf * KslNext_)
-        );
     }
+
+    return tmp<volScalarField>
+    (
+        new volScalarField("Ksl_implicitCouple", (1. - tsf) * KslPrev_ + tsf * KslNext_)
+    );
 }
 
 void implicitCouple::resetMomSourceField() const

src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.C

@@ -102,22 +102,12 @@ voidFractionModel::~voidFractionModel()
 // * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
 tmp<volScalarField> voidFractionModel::voidFractionInterp() const
 {
-    scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
+    const scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
 
-    if (1. - tsf < 1e-4 && particleCloud_.dataExchangeM().couplingStep() > 1) //tsf==1
-    {
-        return tmp<volScalarField>
-        (
-            new volScalarField("alpha_voidFractionModel", voidfractionPrev_)
-        );
-    }
-    else
-    {
-        return tmp<volScalarField>
-        (
-            new volScalarField("alpha_voidFractionModel", (1. - tsf) * voidfractionPrev_ + tsf * voidfractionNext_)
-        );
-    }
+    return tmp<volScalarField>
+    (
+        new volScalarField("alpha_voidFractionModel", (1. - tsf) * voidfractionPrev_ + tsf * voidfractionNext_)
+    );
 }
 
 void voidFractionModel::resetVoidFractions() const