View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0003100OpenFOAMBugpublic2018-11-06 05:312018-11-12 22:51
Reporteralbertop Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSOtherOS Version(please specify)
Product Versiondev 
Fixed in Versiondev 
Summary0003100: Calculation of r in NVDTVD.H
DescriptionThe value of r used in limiter functions is now calculated as ( https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/NVDTVD.H ):

r = 2*(gradcf/gradf) - 1

However, in some cases, this may not guarantee positivity of r, so it would be best to take the maximum of this quantity and 0:

r = max(0.0, 2*(gradcf/gradf) - 1)

I am also not sure about the rationale of using:

r = 2*1000*sign(gradcf)*sign(gradf) - 1

when mag(gradcf) >= 1000*mag(gradf).
Wouldn't it be sufficient to ensure that the denominator gradf is finite (>= small number)?

TagsNo tags attached.

Activities

albertop

2018-11-06 05:46

reporter   ~0010134

Similar considerations apply to the "r" function in NVDVTVDV.H

henry

2018-11-06 07:17

manager   ~0010135

This form of r was developed and tuned over a very long time, in particular the r = 2*1000*sign(gradcf)*sign(gradf) - 1 to avoid division overflow problems. Can you provide an alternative comparison approach which avoids the division which you are happy with?

Could you provide any examples for which the current formulation for r is inappropriate? Note that r does not assume that the scheme will be TVD formal so it is not clear that r must be positive; it is the responsibility of the particular scheme to ensure this.

albertop

2018-11-06 08:29

reporter   ~0010137

OK. I did not understand the logic behind r = 2*1000*sign(gradcf)*sign(gradf) - 1 because it returns either 2001 or 1999, depending on the sign of sign(gradcf)*sign(gradf). I was thinking about clipping (or summing/subtracting a small value depending on the sign) the denominator when calculating gradcf/gradf, but I have not compared to the current implementation.

About the calculations of r, I noticed this because I amddoing a code review for some work I am finishing. I agree most schemes take care of forcing the limiter function to be positive, but not all. For example, to mention one available in OpenFOAM, the symmetric van Albada and the OSPRE schemes have a (not strictly) positive limiter function for r <= -1 and r >= 0, so if r becomes slightly negative, the limiter function is not positive or zero.

henry

2018-11-06 08:58

manager   ~0010138

I do not see any problem with the current logic and it replaced a simpler version which added or subtracted a small amount but on some systems with some compilers this caused problems. If you can demonstrate an issue with the current implementation and an alternative with sufficient testing we can consider changing it.

It is not clear from your request if you require r to be >= 0 or >= -1 and I have not seen any issues with the operation of van Albada or OSPRE; do you believe either of these schemes needs to be changed?

albertop

2018-11-06 22:40

reporter   ~0010145

In the two schemes mentioned, which are TVD, I think the value of r passed to the limiter function should be positive so that also the limiter function is guaranteed to be positive. In most other limiters this is not necessary because the value of the limiter function is clipped to zero in any case.

henry

2018-11-06 23:01

manager   ~0010146

Unless we write separate "r" functions for each scheme this is not possible and if the scheme requires r to be limited in some manner it should do it. This avoids a lot of code duplication by keeping "r" as general as possible and different schemes either don't need to limit "r" or can do what is needed to it.

albertop

2018-11-06 23:19

reporter   ~0010147

I am not asking to write separate r functions. Clipping the r value to zero in the limiter function would be fine.
I am still not clear about what scheme would admit r < 0 (none of the TVD), but if it's needed anywhere, that's OK.

henry

2018-11-06 23:37

manager   ~0010148

You are correct that TVD schemes would require r to be positive but the current structure supports non-TVD schemes.

OK, do you know which limiter functions should be changed?

albertop

2018-11-07 15:47

reporter   ~0010154

I will finish checking the schemes implemented in OpenFOAM and share the proposed changes.

albertop

2018-11-11 04:53

reporter   ~0010173

I have checked all the schemes in the limitedSchemes directory. I think only the van Albada and OSPRE limited schemes may be affected by this, so I propose to use max(0, r) in these two.
I attach the two modified files with the proposed change.
OSPRE.H (2,795 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 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/>.

Class
    Foam::OSPRELimiter

Description
    Class with limiter function which returns the limiter for the
    OSPRE differencing scheme based on r obtained from the LimiterFunc
    class.

    Used in conjunction with the template class LimitedScheme.

SourceFiles
    OSPRE.C

\*---------------------------------------------------------------------------*/

#ifndef OSPRE_H
#define OSPRE_H

#include "vector.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
                           Class OSPRELimiter Declaration
\*---------------------------------------------------------------------------*/

template<class LimiterFunc>
class OSPRELimiter
:
    public LimiterFunc
{

public:

    OSPRELimiter(Istream&)
    {}

    scalar limiter
    (
        const scalar cdWeight,
        const scalar faceFlux,
        const typename LimiterFunc::phiType& phiP,
        const typename LimiterFunc::phiType& phiN,
        const typename LimiterFunc::gradPhiType& gradcP,
        const typename LimiterFunc::gradPhiType& gradcN,
        const vector& d
    ) const
    {
        scalar r = max
        (
            0,
            LimiterFunc::r
            (
                faceFlux, phiP, phiN, gradcP, gradcN, d
            )
        );

        scalar rrp1 = r*(r + 1);
        return 1.5*rrp1/(rrp1 + 1);
    }
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //
OSPRE.H (2,795 bytes)   
vanAlbada.H (2,797 bytes)    
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 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/>.

Class
    Foam::vanAlbadaLimiter

Description
    Class with limiter function which returns the limiter for the
    vanAlbada differencing scheme based on r obtained from the LimiterFunc
    class.

    Used in conjunction with the template class LimitedScheme.

SourceFiles
    vanAlbada.C

\*---------------------------------------------------------------------------*/

#ifndef vanAlbada_H
#define vanAlbada_H

#include "vector.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

/*---------------------------------------------------------------------------*\
                           Class vanAlbadaLimiter Declaration
\*---------------------------------------------------------------------------*/

template<class LimiterFunc>
class vanAlbadaLimiter
:
    public LimiterFunc
{

public:

    vanAlbadaLimiter(Istream&)
    {}

    scalar limiter
    (
        const scalar cdWeight,
        const scalar faceFlux,
        const typename LimiterFunc::phiType& phiP,
        const typename LimiterFunc::phiType& phiN,
        const typename LimiterFunc::gradPhiType& gradcP,
        const typename LimiterFunc::gradPhiType& gradcN,
        const vector& d
    ) const
    {
        scalar r = max
        (
            0,
            LimiterFunc::r
            (
                faceFlux, phiP, phiN, gradcP, gradcN, d
            )
        );

        return r*(r + 1)/(sqr(r) + 1);
    }
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //
vanAlbada.H (2,797 bytes)   

henry

2018-11-12 22:51

manager   ~0010177

Thanks Alberto,

Resolved by commit c3018b3f0779e68aa3e517b0d18ee13ceb9b50a3

Note that I could not commit the patches provided because they were out of date with respect to current OpenFOAM-dev, please pull the latest to ensure future patches are consistent.

Issue History

Date Modified Username Field Change
2018-11-06 05:31 albertop New Issue
2018-11-06 05:46 albertop Note Added: 0010134
2018-11-06 07:17 henry Note Added: 0010135
2018-11-06 08:29 albertop Note Added: 0010137
2018-11-06 08:58 henry Note Added: 0010138
2018-11-06 22:40 albertop Note Added: 0010145
2018-11-06 23:01 henry Note Added: 0010146
2018-11-06 23:19 albertop Note Added: 0010147
2018-11-06 23:37 henry Note Added: 0010148
2018-11-07 15:47 albertop Note Added: 0010154
2018-11-11 04:53 albertop File Added: OSPRE.H
2018-11-11 04:53 albertop File Added: vanAlbada.H
2018-11-11 04:53 albertop Note Added: 0010173
2018-11-12 22:51 henry Assigned To => henry
2018-11-12 22:51 henry Status new => resolved
2018-11-12 22:51 henry Resolution open => fixed
2018-11-12 22:51 henry Fixed in Version => dev
2018-11-12 22:51 henry Note Added: 0010177