From b6b2e7d1f453e6537d5331f44d5ee7fc8423bab9 Mon Sep 17 00:00:00 2001
From: Thomas Lichtenegger <thomas.lichtenegger@jku.at>
Date: Thu, 15 Sep 2016 15:40:34 +0200
Subject: [PATCH] Different treatment for velocity fluctuations.

---
 .../forceModel/freeStreaming/freeStreaming.C  |  1 +
 .../isotropicFluctuations.C                   | 36 ++++++++++++++-----
 .../isotropicFluctuations.H                   |  4 +++
 3 files changed, 33 insertions(+), 8 deletions(-)

diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/freeStreaming/freeStreaming.C b/src/lagrangian/cfdemParticle/subModels/forceModel/freeStreaming/freeStreaming.C
index ecde6bca..d9910b41 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/freeStreaming/freeStreaming.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/freeStreaming/freeStreaming.C
@@ -123,6 +123,7 @@ void freeStreaming::setForce() const
 		}
 		else
 		{
+		    position = particleCloud_.position(index);
 		    radius = particleCloud_.radius(index);
                     mass = 4.188790205*radius*radius*radius * particleDensity_;
 		    grav = mass*gravAcc_;
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.C b/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.C
index d21ac700..d8f570b4 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.C
@@ -62,11 +62,12 @@ isotropicFluctuations::isotropicFluctuations
     voidfractionRec_(sm.mesh().lookupObject<volScalarField> (voidfractionRecFieldName_)),
     critVoidfraction_(propsDict_.lookupOrDefault<scalar>("critVoidfraction", 1.0)),
     D0_(readScalar(propsDict_.lookup("D0"))),
-    ranGen_(osRandomInteger())
+    ranGen_(osRandomInteger()),
+    vfluc_(NULL)
 {
     dt_=particleCloud_.mesh().time().deltaTValue();
     
-    forceSubModels_.setSize(1, "recU");
+  /*  forceSubModels_.setSize(1, "recU");
     delete[] forceSubModel_;
     forceSubModel_ = new autoPtr<forceSubModel>[nrForceSubModels()];
     Info << "nrForceSubModels()=" << nrForceSubModels() << endl;
@@ -80,19 +81,34 @@ isotropicFluctuations::isotropicFluctuations
             forceSubModels_[i]
         );
     }
+  */
+    allocateMyArrays();
 }
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 isotropicFluctuations::~isotropicFluctuations()
-{}
+{
+    delete vfluc_;
+}
 
 
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+void isotropicFluctuations::allocateMyArrays() const
+{
+    // get memory for 2d arrays
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(vfluc_,initVal,3); 
+}
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void isotropicFluctuations::setForce() const
 {
+     // realloc the arrays
+    allocateMyArrays();
+    
     vector position(0,0,0);
     scalar voidfraction(0.0);
     scalar voidfractionRec(0.0);
@@ -112,7 +128,7 @@ void isotropicFluctuations::setForce() const
 	    voidfraction=0.0;
 	    voidfractionRec=0.0;
 	    deltaVoidfrac=0.0;
-            if (cellI > -1) // particle Found
+            if (cellI > -1) // particle found
             {
 	        // particles in empty regions follow trajectories subject to gravity
 		if(voidfractionRec_[cellI] < critVoidfraction_)
@@ -135,11 +151,15 @@ void isotropicFluctuations::setForce() const
 		    if(deltaVoidfrac>0)
 		    {
                         flucU=unitRndVec()*fluctuationMag(relVolfractionExcess);
-			forceSubM(0).partToArray(index,flucU,vector::zero);
+	//		forceSubM(0).partToArray(index,flucU,vector::zero);
 		    }              
+		    for(int i = 0; i < 3; i++)
+		        vfluc_[index][i]=flucU[i];
 		}
 	    }
     }
+    
+    particleCloud_.dataExchangeM().giveData("vfluc","vector-atom", vfluc_);   
 }
 
 scalar isotropicFluctuations::fluctuationMag(const scalar relVolfractionExcess) const
@@ -174,9 +194,9 @@ vector isotropicFluctuations::unitRndVec() const
 	s=v1*v1+v2*v2;
     }
     s2=Foam::sqrt(1-s);
-    rvec[1]=2*v1*s2;
-    rvec[2]=2*v2*s2;
-    rvec[3]=1-2*s;
+    rvec[0]=2*v1*s2;
+    rvec[1]=2*v2*s2;
+    rvec[2]=1-2*s;
     return rvec;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.H b/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.H
index 06a236f4..3ef31375 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/isotropicFluctuations/isotropicFluctuations.H
@@ -66,6 +66,10 @@ private:
     scalar fluctuationMag(const scalar) const;
     
     mutable Random ranGen_;
+    
+    mutable double **vfluc_;              // Lagrangian array
+
+    void allocateMyArrays() const;
 
 public: