From fe4016eb52ab68202eb804a5a6bc77ad18ef9087 Mon Sep 17 00:00:00 2001
From: Thomas Lichtenegger <thomas.lichtenegger@jku.at>
Date: Wed, 13 Nov 2019 15:53:27 +0100
Subject: [PATCH] Update of random number generation. Removed deprecated
 multiIntervalPath model.

---
 .../solvers/rStatAnalysis/Make/options        |   2 +-
 .../solvers/rStatAnalysis/rStatAnalysis.C     |   5 -
 .../isotropicFluctuations.C                   |   2 +-
 src/recurrence/Make/files                     |   3 +-
 src/recurrence/recBase/recBase.C              |  20 +-
 .../recPath/MarkovPath/MarkovPath.C           |   3 +-
 .../multiIntervalPath/multiIntervalPath.C     | 195 ------------------
 .../multiIntervalPath/multiIntervalPath.H     | 103 ---------
 .../simpleRandomPath/simpleRandomPath.C       |   3 +-
 .../autocorrelation/autocorrelation.C         |  67 ++++--
 .../autocorrelation/autocorrelation.H         |  12 ++
 .../recStatAnalysis/recStatAnalysis.H         |   2 +
 12 files changed, 81 insertions(+), 336 deletions(-)
 delete mode 100644 src/recurrence/recPath/multiIntervalPath/multiIntervalPath.C
 delete mode 100644 src/recurrence/recPath/multiIntervalPath/multiIntervalPath.H

diff --git a/applications/solvers/rStatAnalysis/Make/options b/applications/solvers/rStatAnalysis/Make/options
index f91d7947..0d54afbc 100644
--- a/applications/solvers/rStatAnalysis/Make/options
+++ b/applications/solvers/rStatAnalysis/Make/options
@@ -23,4 +23,4 @@ EXE_LIBS = \
     -lmeshTools \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
-$(CFDEM_ADD_LIBS)
+     $(CFDEM_ADD_LIBS)
diff --git a/applications/solvers/rStatAnalysis/rStatAnalysis.C b/applications/solvers/rStatAnalysis/rStatAnalysis.C
index 38356182..72523bfd 100644
--- a/applications/solvers/rStatAnalysis/rStatAnalysis.C
+++ b/applications/solvers/rStatAnalysis/rStatAnalysis.C
@@ -29,11 +29,6 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-// #include "fvCFD.H"
-// #include "singlePhaseTransportModel.H"
-// #include "turbulentTransportModel.H"
-// #include "fvOptions.H"
-
 #include "recBase.H"
 #include "recStatAnalysis.H"
 
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.C b/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.C
index d67b1391..ea8ec320 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.C
@@ -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")
     {
diff --git a/src/recurrence/Make/files b/src/recurrence/Make/files
index 31c09b23..eabb090e 100644
--- a/src/recurrence/Make/files
+++ b/src/recurrence/Make/files
@@ -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
diff --git a/src/recurrence/recBase/recBase.C b/src/recurrence/recBase/recBase.C
index 53b3fdcb..0b030bc3 100644
--- a/src/recurrence/recBase/recBase.C
+++ b/src/recurrence/recBase/recBase.C
@@ -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_        ->  writeRecMatrix();
-  recModel_        ->  writeRecPath();
+    recModel_        ->  init();
+  
+    if (!recStatAnalysis_->suppressMatrixAndPath())
+    {
+        recModel_        ->  writeRecMatrix();
+        recModel_        ->  writeRecPath();
+    }
 }
 
 // * * * * * * * * * * * * * * * * Destructors  * * * * * * * * * * * * * * //
diff --git a/src/recurrence/recPath/MarkovPath/MarkovPath.C b/src/recurrence/recPath/MarkovPath/MarkovPath.C
index 2dcd427d..88c77875 100644
--- a/src/recurrence/recPath/MarkovPath/MarkovPath.C
+++ b/src/recurrence/recPath/MarkovPath/MarkovPath.C
@@ -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() );
 
diff --git a/src/recurrence/recPath/multiIntervalPath/multiIntervalPath.C b/src/recurrence/recPath/multiIntervalPath/multiIntervalPath.C
deleted file mode 100644
index 15483e7f..00000000
--- a/src/recurrence/recPath/multiIntervalPath/multiIntervalPath.C
+++ /dev/null
@@ -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
-
-// ************************************************************************* //
diff --git a/src/recurrence/recPath/multiIntervalPath/multiIntervalPath.H b/src/recurrence/recPath/multiIntervalPath/multiIntervalPath.H
deleted file mode 100644
index 06c1e691..00000000
--- a/src/recurrence/recPath/multiIntervalPath/multiIntervalPath.H
+++ /dev/null
@@ -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
-
-// ************************************************************************* //
diff --git a/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.C b/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.C
index e31d8a66..687ef55c 100644
--- a/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.C
+++ b/src/recurrence/recPath/simpleRandomPath/simpleRandomPath.C
@@ -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;
diff --git a/src/recurrence/recStatAnalysis/autocorrelation/autocorrelation.C b/src/recurrence/recStatAnalysis/autocorrelation/autocorrelation.C
index 0da4e408..1a3c253f 100644
--- a/src/recurrence/recStatAnalysis/autocorrelation/autocorrelation.C
+++ b/src/recurrence/recStatAnalysis/autocorrelation/autocorrelation.C
@@ -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;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/recurrence/recStatAnalysis/autocorrelation/autocorrelation.H b/src/recurrence/recStatAnalysis/autocorrelation/autocorrelation.H
index ae08ad51..c25e5284 100644
--- a/src/recurrence/recStatAnalysis/autocorrelation/autocorrelation.H
+++ b/src/recurrence/recStatAnalysis/autocorrelation/autocorrelation.H
@@ -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_;}
 };
 
 
diff --git a/src/recurrence/recStatAnalysis/recStatAnalysis/recStatAnalysis.H b/src/recurrence/recStatAnalysis/recStatAnalysis/recStatAnalysis.H
index accf1b38..ee86c9d5 100644
--- a/src/recurrence/recStatAnalysis/recStatAnalysis/recStatAnalysis.H
+++ b/src/recurrence/recStatAnalysis/recStatAnalysis/recStatAnalysis.H
@@ -102,6 +102,8 @@ public:
 
     virtual void statistics() = 0;
 
+    virtual bool suppressMatrixAndPath() {return false;}
+
 
 };