View Issue Details

IDProjectCategoryView StatusLast Update
0001650OpenFOAMBugpublic2015-04-21 07:39
Reportertniemi Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
Summary0001650: stochasticDispersionRAS: still some problems with UTurb direction
DescriptionThis 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(1-u^2)*cos(theta)
y=sqrt(1-u^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.
TagsNo tags attached.

Activities

tniemi

2015-04-08 11:56

reporter  

StochasticDispersionRAS.C (3,816 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 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(1-u*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;
}


// ************************************************************************* //
StochasticDispersionRAS.C (3,816 bytes)   

henry

2015-04-08 13:12

manager   ~0004572

Thanks for the patch, I have included the changes in OpenFOAM-2.3.x:
commit a869754bc425db63fc23d5f3c4aafa175ee2bad4
and OpenFOAM-dev
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.

tniemi

2015-04-08 13:38

reporter   ~0004573

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.

henry

2015-04-08 14:01

manager   ~0004574

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

tniemi

2015-04-09 10:40

reporter  

cachedRandom.zip (5,428 bytes)

tniemi

2015-04-09 10:41

reporter   ~0004585

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?

henry

2015-04-09 14:22

manager   ~0004590

- 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

henry

2015-04-09 15:11

manager   ~0004592

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?

henry

2015-04-09 15:48

manager   ~0004593

I have updated OpenFOAM-dev: commit aaa3d026e8530b4de39a0c1fc58a5fbbeecbc942
Please test the changes and report back

tniemi

2015-04-13 10:47

reporter   ~0004608

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?

henry

2015-04-17 15:14

manager   ~0004614

Could you provide a patch for gradientDispersionRAS to save me having to study the code?

I will check the dispersion support for MPPIC parcels

henry

2015-04-18 11:53

manager   ~0004622

basicKinematicMPPICParcel: Reinstated dispersion models
commit ebeccfff98504c868826a1e44829a495dae1e913

tniemi

2015-04-20 12:16

reporter  

GradientDispersionRAS.C (4,211 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 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;
}


// ************************************************************************* //
GradientDispersionRAS.C (4,211 bytes)   

tniemi

2015-04-20 12:18

reporter   ~0004630

Uploaded modified gradientDispersionRAS.C.

henry

2015-04-20 17:51

manager   ~0004635

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

tniemi

2015-04-21 04:35

reporter   ~0004639

Yes, in my opinion everything is fixed now.

Issue History

Date Modified Username Field Change
2015-04-08 11:56 tniemi New Issue
2015-04-08 11:56 tniemi File Added: StochasticDispersionRAS.C
2015-04-08 13:12 henry Note Added: 0004572
2015-04-08 13:38 tniemi Note Added: 0004573
2015-04-08 14:01 henry Note Added: 0004574
2015-04-09 10:40 tniemi File Added: cachedRandom.zip
2015-04-09 10:41 tniemi Note Added: 0004585
2015-04-09 14:22 henry Note Added: 0004590
2015-04-09 15:11 henry Note Added: 0004592
2015-04-09 15:48 henry Note Added: 0004593
2015-04-13 10:47 tniemi Note Added: 0004608
2015-04-17 15:14 henry Note Added: 0004614
2015-04-18 11:53 henry Note Added: 0004622
2015-04-20 12:16 tniemi File Added: GradientDispersionRAS.C
2015-04-20 12:18 tniemi Note Added: 0004630
2015-04-20 17:51 henry Note Added: 0004635
2015-04-21 04:35 tniemi Note Added: 0004639
2015-04-21 07:39 henry Status new => resolved
2015-04-21 07:39 henry Resolution open => fixed
2015-04-21 07:39 henry Assigned To => henry