View Issue Details

IDProjectCategoryView StatusLast Update
0002634OpenFOAMBugpublic2017-07-27 16:57
ReporterShorty Assigned Towill  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version16.04
Product Versiondev 
Summary0002634: limitVelocity function
DescriptionHi all,

I guess the limitVelocity function has some bug. As far as I understand the function, it limit the magnitude of the velocity field to the one we specify with the max variable, right? However, the way it is implemented does not provide what I expect. I expect the following:

- If the magnitude of the velocity of a cell is higher the one we specified, change the velocity in a way to get the magnitude value we specified (max).

Right now it is as follow: If the magSqr of the velocity is larger than the sqr of max, we modify the velocity (reduce it) by the ratio of sqr(max)/magSqr(Ui). If we do so, the magnitude of the new vector is much smaller than the one we specified with the max keyword (example below). In addition, the actual implementation reduces the magnitude of the vector even more, the higher its magnitude is. But reducing to the one we specified will never happen.

Again:
If I got it right, we should do it as follow. If the magnitude of the velocity is larger than the one we specified, reduce the velocity in a way that the magnitude of the new vector is equal to the one we specified with the max keyword, right?

If this is correct, then we can use the simple patch I attached. Otherwise I would be happy to understand the actual imeplementation because even with a discussion with my colleague today, we did not get the point.

The example below demonstrate the above text:
Steps To ReproduceExample given:

U = (2 1 4) // of the cell we check out
max = 2;

mag(U) = 4.58258 --> bigger than max so reduce it.

Actual calculation will give us a new velocity vector as:
=========================================================
sqrMax = sqr(2) = 4
magSqrUi = magSqr(U) = 21
Unew = U * sqrMax / magSqrUi = 4 / 21 = (0.190476 0.0952381 0.380952)

However, the magnitude of the new vector is:

mag(Unew) = 0.436436 (I expect it to be 2)



Patched version calculates the new velocity vector as:
=========================================================
max = 2
magUi = mag(U) = 4.58258
Unew = U * max / magUi = (0.872872 0.436436 1.74574)

Now we get the magnitude of the new vector exact to be 2 (as we specified with the max keyword)
mag(Unew) = 2

TagslimitVelocity

Activities

Shorty

2017-07-26 23:21

reporter  

limitVelocityPatch.C (3,254 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2016 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 "limitVelocity.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
namespace fv
{
    defineTypeNameAndDebug(limitVelocity, 0);
    addToRunTimeSelectionTable
    (
        option,
        limitVelocity,
        dictionary
    );
}
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::fv::limitVelocity::limitVelocity
(
    const word& name,
    const word& modelType,
    const dictionary& dict,
    const fvMesh& mesh
)
:
    cellSetOption(name, modelType, dict, mesh),
    UName_(coeffs_.lookupOrDefault<word>("U", "U")),
    max_(readScalar(coeffs_.lookup("max")))
{
    fieldNames_.setSize(1, UName_);
    applied_.setSize(1, false);
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

bool Foam::fv::limitVelocity::read(const dictionary& dict)
{
    if (cellSetOption::read(dict))
    {
        coeffs_.lookup("max") >> max_;

        return true;
    }
    else
    {
        return false;
    }
}


void Foam::fv::limitVelocity::correct(volVectorField& U)
{
    const scalar maxU = max_;

    vectorField& Uif = U.primitiveFieldRef();

    forAll(cells_, i)
    {
        const label celli = cells_[i];

        const scalar magUi = mag(Uif[celli]);

        if (magUi > maxU)
        {
            Uif[celli] *= maxU/magUi;
        }
    }

    // handle boundaries in the case of 'all'
    if (selectionMode_ == smAll)
    {
        volVectorField::Boundary& Ubf = U.boundaryFieldRef();

        forAll(Ubf, patchi)
        {
            fvPatchVectorField& Up = Ubf[patchi];

            if (!Up.fixesValue())
            {
                forAll(Up, facei)
                {
                    const scalar magUi = mag(Up[facei]);

                    if (magUi > maxU)
                    {
                        Up[facei] *= maxU/magUi;
                    }
                }
            }
        }
    }
}


// ************************************************************************* //
limitVelocityPatch.C (3,254 bytes)   

will

2017-07-27 15:21

manager   ~0008463

Thanks for the report

Fixed in dev by commit 16b559c1093019a4a85aa8c810c9d837dfa8fa4c

Fixed in 5.x by commit 3ece3e9e81ef05c8336bb7fa5e9a8e0bd78aa6b2

Shorty

2017-07-27 16:04

reporter   ~0008465

Just one remark to the patch. Taking the sqrt() of the prefaktor is correct but we could reduce the call of functions as given in my patch :)

Well I know, this will not influence the speed of the code too much and is tweak in the micro level :)

By the way, you are welcomed Will.
Nice to see the patch already there.

will

2017-07-27 16:11

manager   ~0008466

Your patch calls mag (which contains sqrt) on every velocity value. My fix only calls sqrt for the velocity values that are larger than the maximum. Sqrt is the expensive bit, so mine will be significantly quicker.

Shorty

2017-07-27 16:54

reporter   ~0008468

Oh yes you are right :)

Thanks Will for the note. Then everything is fine.

Issue History

Date Modified Username Field Change
2017-07-26 23:21 Shorty New Issue
2017-07-26 23:21 Shorty File Added: limitVelocityPatch.C
2017-07-26 23:21 Shorty Tag Attached: limitVelocity
2017-07-27 15:21 will Assigned To => will
2017-07-27 15:21 will Status new => resolved
2017-07-27 15:21 will Resolution open => fixed
2017-07-27 15:21 will Note Added: 0008463
2017-07-27 16:04 Shorty Status resolved => feedback
2017-07-27 16:04 Shorty Resolution fixed => reopened
2017-07-27 16:04 Shorty Note Added: 0008465
2017-07-27 16:11 will Status feedback => resolved
2017-07-27 16:11 will Resolution reopened => fixed
2017-07-27 16:11 will Note Added: 0008466
2017-07-27 16:54 Shorty Status resolved => feedback
2017-07-27 16:54 Shorty Resolution fixed => reopened
2017-07-27 16:54 Shorty Note Added: 0008468
2017-07-27 16:57 will Status feedback => resolved
2017-07-27 16:57 will Resolution reopened => fixed