View Issue Details
ID  Project  Category  View Status  Date Submitted  Last Update 

0001650  OpenFOAM  Bug  public  20150408 11:56  20150421 07:39 
Reporter  tniemi  Assigned To  henry  
Priority  normal  Severity  minor  Reproducibility  always 
Status  resolved  Resolution  fixed  
Summary  0001650: stochasticDispersionRAS: still some problems with UTurb direction  
Description  This bug report is related to http://www.openfoam.org/mantisbt/view.php?id=1568. Unfortunately, the applied fix (using polar coords) is not correct, as it will produce more directions close to the poles of the sphere, for discussion see eg http://mathworld.wolfram.com/SpherePointPicking.html. The "correct" way of using polar coordinates would be x=sqrt(1u^2)*cos(theta) y=sqrt(1u^2)*sin(theta) z=u where u is uniform(1,1) and theta is uniform(0,2pi). I have attached a patch with corrected algo. In addition, the patch should be added to both small t and big T turbulent models in 2.3.x, as the previous commit only modified small t model. In addition, for some extra performance gain, the cachedRandom could be modified to provide access to GaussNormal(). Then the fac calculation could be changed to use it. The GaussNormal() should be 50 % faster, as it caches the second random number, which is now discarded when computing fac.  
Tags  No tags attached.  

StochasticDispersionRAS.C (3,816 bytes)
/**\ =========  \\ / F ield  OpenFOAM: The Open Source CFD Toolbox \\ / O peration  \\ / A nd  Copyright (C) 20112015 OpenFOAM Foundation \\/ M anipulation   License This file is part of OpenFOAM. OpenFOAM 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. OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>. \**/ #include "StochasticDispersionRAS.H" #include "constants.H" using namespace Foam::constant::mathematical; // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template<class CloudType> Foam::StochasticDispersionRAS<CloudType>::StochasticDispersionRAS ( const dictionary& dict, CloudType& owner ) : DispersionRASModel<CloudType>(dict, owner) {} template<class CloudType> Foam::StochasticDispersionRAS<CloudType>::StochasticDispersionRAS ( StochasticDispersionRAS<CloudType>& dm ) : DispersionRASModel<CloudType>(dm) {} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template<class CloudType> Foam::StochasticDispersionRAS<CloudType>::~StochasticDispersionRAS() {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class CloudType> Foam::vector Foam::StochasticDispersionRAS<CloudType>::update ( const scalar dt, const label cellI, const vector& U, const vector& Uc, vector& UTurb, scalar& tTurb ) { cachedRandom& rnd = this>owner().rndGen(); const scalar cps = 0.16432; const scalar k = this>kPtr_>internalField()[cellI]; const scalar epsilon = this>epsilonPtr_>internalField()[cellI] + ROOTVSMALL; const scalar UrelMag = mag(U  Uc  UTurb); const scalar tTurbLoc = min(k/epsilon, cps*pow(k, 1.5)/epsilon/(UrelMag + SMALL)); // Parcel is perturbed by the turbulence if (dt < tTurbLoc) { tTurb += dt; if (tTurb > tTurbLoc) { tTurb = 0; const scalar sigma = sqrt(2*k/3.0); // Calculate a random direction dir distributed uniformly // in spherical coordinates const scalar theta = rnd.sample01<scalar>()*twoPi; const scalar u = rnd.sample01<scalar>()*2  1; const scalar a = sqrt(1u*u); const vector dir(a*cos(theta), a*sin(theta), u); // Numerical Recipes... Ch. 7. Random Numbers... scalar x1 = 0; scalar x2 = 0; scalar rsq = 10; while ((rsq > 1)  (rsq == 0)) { x1 = 2*rnd.sample01<scalar>()  1; x2 = 2*rnd.sample01<scalar>()  1; rsq = x1*x1 + x2*x2; } scalar fac = sqrt(2*log(rsq)/rsq); fac *= mag(x1); UTurb = sigma*fac*dir; } } else { tTurb = GREAT; UTurb = vector::zero; } return Uc + UTurb; } // ************************************************************************* // 

Thanks for the patch, I have included the changes in OpenFOAM2.3.x: commit a869754bc425db63fc23d5f3c4aafa175ee2bad4 and OpenFOAMdev commit 782763997ad8394d04260c90b2ac40295166e09b Currently cachedRandom does not provide the GaussNormal function which is in the Random class but I think it would be MUCH better if the cachedRandom class reused and provided all the functionality of that class. However it is not clear how this would provide a speed improvement. 

I was thinking that the speed improvement would come from the fact, that now we are discarding the x2 value when computing fac, while using GaussNormal we would get every second normal random variable for free. (the algo gives two random variables at one pass and gaussnormal uses static variables to store the other) Of course it would be possible to do a similar static variable caching also in the dispersion model, but it would be cleaner if we could just call a GaussNormal function. 

I agree, would you like to have a go at adding the GaussNormal function to the cachedRandom class? 



Okay, I have now added a modified cachedRandom with GaussNormal. It should be a straightforward modification, but I haven't yet tested it much other than it compiles. Hopefully it is useful. List of changes:  Added bool hasGaussSample_ and scalar gaussSample_ member variables. I used members instead of static variables, since static would be shared with every instance Added normal01, which is the actual algorithm copied from random Added GaussNormal local and global functions with template specializations for label and scalar used scatter instead of reduce in GaussNormal global functions. This reminded me of the discussion regarding global randoms in http://www.openfoam.org/mantisbt/view.php?id=1556. cachedRandom uses reduce everywhere, but scatter would be better I guess? Also I realized, that since there is a globalPosition function, it could be directly used in the injection model? 

 Added bool hasGaussSample_ and scalar gaussSample_ member variables. I used members instead of static variables, since static would be shared with every instance Agreed, static variable should be avoided and in this case can be stored as private data Added normal01, which is the actual algorithm copied from random Fine Added GaussNormal local and global functions with template specializations for label and scalar OK used scatter instead of reduce in GaussNormal global functions. This reminded me of the discussion regarding global randoms in http://www.openfoam.org/mantisbt/view.php?id=1556. [^] cachedRandom uses reduce everywhere, but scatter would be better I guess? In principle yes but currently Pstream does not map scatter through to the corresponding MPI call but does so for reductions. Nevertheless we should use scatter where appropriate and update the Pstream scatter functions. > Also I realized, that since there is a globalPosition function, it could be directly used in the injection model? Agreed 

I was about to replace the code in StochasticDispersionRAS with rnd.GaussNormal<scalar> but there is a difference in the implementations: scalar fac = sqrt(2*log(rsq)/rsq); fac *= mag(x1); vs scalar fac = sqrt(2*log(rsq)/rsq); gaussSample_ = v1*fac; i.e. one has a mag. Which is correct? 

I have updated OpenFOAMdev: commit aaa3d026e8530b4de39a0c1fc58a5fbbeecbc942 Please test the changes and report back 

Okay, I have now run some (artificial) tests and everything appears to work as intended. I checked that the GaussNormal produces correct distribution and the directions are spherically uniform. The mag is used to ensure, that fac does not flip direction, but I guess you already figured it out. Also the gradientDispersionRAS could be changed to use the new formulation. BTW! I noticed that dev does not have dispersion support for MPPIC parcels. The MPPIC support was only in the big T turbulence models in the 2.3.x, so it should probably be added to dev? 

Could you provide a patch for gradientDispersionRAS to save me having to study the code? I will check the dispersion support for MPPIC parcels 

basicKinematicMPPICParcel: Reinstated dispersion models commit ebeccfff98504c868826a1e44829a495dae1e913 

GradientDispersionRAS.C (4,211 bytes)
/**\ =========  \\ / F ield  OpenFOAM: The Open Source CFD Toolbox \\ / O peration  \\ / A nd  Copyright (C) 20112015 OpenFOAM Foundation \\/ M anipulation   License This file is part of OpenFOAM. OpenFOAM 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. OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>. \**/ #include "GradientDispersionRAS.H" #include "demandDrivenData.H" #include "fvcGrad.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template<class CloudType> Foam::GradientDispersionRAS<CloudType>::GradientDispersionRAS ( const dictionary& dict, CloudType& owner ) : DispersionRASModel<CloudType>(dict, owner), gradkPtr_(NULL), ownGradK_(false) {} template<class CloudType> Foam::GradientDispersionRAS<CloudType>::GradientDispersionRAS ( const GradientDispersionRAS<CloudType>& dm ) : DispersionRASModel<CloudType>(dm), gradkPtr_(dm.gradkPtr_), ownGradK_(dm.ownGradK_) { dm.ownGradK_ = false; } // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // template<class CloudType> Foam::GradientDispersionRAS<CloudType>::~GradientDispersionRAS() { cacheFields(false); } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // template<class CloudType> void Foam::GradientDispersionRAS<CloudType>::cacheFields(const bool store) { DispersionRASModel<CloudType>::cacheFields(store); if (store) { gradkPtr_ = fvc::grad(*this>kPtr_).ptr(); ownGradK_ = true; } else { if (ownGradK_) { deleteDemandDrivenData(gradkPtr_); gradkPtr_ = NULL; ownGradK_ = false; } } } template<class CloudType> Foam::vector Foam::GradientDispersionRAS<CloudType>::update ( const scalar dt, const label cellI, const vector& U, const vector& Uc, vector& UTurb, scalar& tTurb ) { cachedRandom& rnd = this>owner().rndGen(); const scalar cps = 0.16432; const scalar k = this>kPtr_>internalField()[cellI]; const scalar epsilon = this>epsilonPtr_>internalField()[cellI] + ROOTVSMALL; const vector& gradk = this>gradkPtr_>internalField()[cellI]; const scalar UrelMag = mag(U  Uc  UTurb); const scalar tTurbLoc = min(k/epsilon, cps*pow(k, 1.5)/epsilon/(UrelMag + SMALL)); // Parcel is perturbed by the turbulence if (dt < tTurbLoc) { tTurb += dt; if (tTurb > tTurbLoc) { tTurb = 0.0; scalar sigma = sqrt(2.0*k/3.0); vector dir = gradk/(mag(gradk) + SMALL); scalar fac = 0.0; // In 2D calculations the grad(k) is always // away from the axis of symmetry // This creates a 'hole' in the spray and to // prevent this we let fac be both negative/positive if (this>owner().mesh().nSolutionD() == 2) { fac = rnd.GaussNormal<scalar>(); } else { fac = mag(rnd.GaussNormal<scalar>()); } UTurb = sigma*fac*dir; } } else { tTurb = GREAT; UTurb = vector::zero; } return Uc + UTurb; } // ************************************************************************* // 

Uploaded modified gradientDispersionRAS.C. 

Thanks, I have just pushed this change. Shall I close this report now? 

Yes, in my opinion everything is fixed now. 
Date Modified  Username  Field  Change 

20150408 11:56  tniemi  New Issue  
20150408 11:56  tniemi  File Added: StochasticDispersionRAS.C  
20150408 13:12  henry  Note Added: 0004572  
20150408 13:38  tniemi  Note Added: 0004573  
20150408 14:01  henry  Note Added: 0004574  
20150409 10:40  tniemi  File Added: cachedRandom.zip  
20150409 10:41  tniemi  Note Added: 0004585  
20150409 14:22  henry  Note Added: 0004590  
20150409 15:11  henry  Note Added: 0004592  
20150409 15:48  henry  Note Added: 0004593  
20150413 10:47  tniemi  Note Added: 0004608  
20150417 15:14  henry  Note Added: 0004614  
20150418 11:53  henry  Note Added: 0004622  
20150420 12:16  tniemi  File Added: GradientDispersionRAS.C  
20150420 12:18  tniemi  Note Added: 0004630  
20150420 17:51  henry  Note Added: 0004635  
20150421 04:35  tniemi  Note Added: 0004639  
20150421 07:39  henry  Status  new => resolved 
20150421 07:39  henry  Resolution  open => fixed 
20150421 07:39  henry  Assigned To  => henry 