View Issue Details

IDProjectCategoryView StatusLast Update
0000676OpenFOAMBugpublic2013-02-12 09:08
Reporteruser534Assigned Touser2 
PriorityimmediateSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
PlatformLinuxOSUbuntuOS Version10.04
Summary0000676: All extrapolated convection schemes on all coupled boundaries use wrong neighbour gradient
Descriptionfvc and all related gradient operators do not call correctBoundaryConditions on result. On coupled boundaries, this means that patch neighbour field takes the current values or the values from patch calculus, instead of correct patchNeighbourField - which would require communication.

Therefore, when a scheme uses coupled patchNeighbourField gradient, the values are wrong and inconsistent across the processor patch: conservation error.

Further to this, it is unclear that correctBoundaryConditions has been called on a main field, but since such fields are usually solved for, the issue is conceptual only. For gradients (where correctBoundaryConditions is never called), this is a major bug.
Additional InformationTo fix, add evaluateCoupled into all fvc operators, OR explicitly call correctBoundaryConditions on all schemes which use gradients on coupled boundaries.
TagsNo tags attached.

Activities

user4

2012-11-05 14:57

  ~0001763

Hi Hrvoje,

I'm trying to replicate this by running linearUpwind in parallel but I don't see any difference on the neighbouring processors. Attached a hacked version of Test-fvc.C.

Mattijs

user534

2012-11-05 15:04

  ~0001764

Well, I did it with a coupled compressible density-based solver, so I really cannot give you the test code - sorry. This is what you do:

- make a 2-processor case for anything: scalar transport will do
- write patchInternalField() and patchNeighbourField() from proc1 and proc2 at the processor boundary and compare. patchInternalField from proc1 must be equal to patchNeighbourField on proc2
- try doing gradPhi.correctBoundaryConditions() and see how it changes.

I don't see the hacked file -am I missing something?

Please shout if you need more - I can call if this helps.

Hrv

user4

2012-11-05 15:08

 

Test-fvc.C (4,494 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/>.

Application
    test

Description
    Finite volume method test code.

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

#include "fvCFD.H"
#include "syncTools.H"
#include "Random.H"

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

template<class Type>
void checkSwap
(
    const GeometricField<Type, fvPatchField, volMesh>& fld
)
{
    const fvMesh& mesh = fld.mesh();

    List<Type> bVals(mesh.nFaces()-mesh.nInternalFaces(), pTraits<Type>::zero);
    forAll(fld.boundaryField(), patchI)
    {
        const polyPatch& pp = mesh.boundaryMesh()[patchI];
        const fvPatchField<Type>& pfld = fld.boundaryField()[patchI];

        const tmp<Field<Type> > pintFld = pfld.patchInternalField();

        forAll(pfld, i)
        {
            if (pp.coupled())
            {
                label meshFaceI = pp.start()+i;
                bVals[meshFaceI-mesh.nInternalFaces()] = pintFld()[i];
            }
        }
    }

    syncTools::swapBoundaryFaceList(mesh, bVals);

    forAll(fld.boundaryField(), patchI)
    {
        const polyPatch& pp = mesh.boundaryMesh()[patchI];
        const fvPatchField<Type>& pfld = fld.boundaryField()[patchI];

        Pout<< "Patch:" << pp.name() << endl;

        forAll(pfld, i)
        {
            if (pp.coupled())
            {
                label meshFaceI = pp.start()+i;
                const Type& bVal = bVals[meshFaceI-mesh.nInternalFaces()];

                if (bVal != pfld[i])
                {
                    Pout<< "face:" << meshFaceI << nl
                        << "    val:" << pfld[i] << nl
                        << "    swp:" << bVal << nl
                        << endl;
                }
            }
        }
    }
}


int main(int argc, char *argv[])
{

#   include "setRootCase.H"

#   include "createTime.H"
#   include "createMesh.H"

    volScalarField fx("fx", pow(mesh.C().component(vector::X), 2));
    fx.write();


    // Get randomised velocity field

    volVectorField U
    (
        IOobject
        (
            "U2",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedVector("one", dimVelocity, vector::zero)
    );

    Random rndGen(0);
    forAll(U.internalField(), cellI)
    {
        rndGen.randomise(U[cellI]);
        U[cellI] -= vector(0.5, 0.5, 0.5);
    }
    U.correctBoundaryConditions();
    U.write();


    // Calculate flux

    surfaceScalarField phi
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        linearInterpolate(U) & mesh.Sf()
    );


    // Calculate divergence

    volScalarField fld(fvc::div(phi, fx));
    checkSwap(fld);



//    volScalarField gradx4(fvc::grad(fx)().component(vector::X));
//    gradx4.write();

    //volVectorField curlC(fvc::curl(1.0*mesh.C()));
    //curlC.write();

    /*
    surfaceScalarField xf(mesh.Cf().component(vector::X));
    surfaceScalarField xf4(pow(xf, 4));

    for (int i=1; i<xf4.size()-1; i++)
    {
        scalar gradx4a = (xf4[i] - xf4[i-1])/(xf[i] - xf[i-1]);
        Info<< (gradx4a - gradx4[i])/gradx4a << endl;
    }
    */

    Info<< "end" << endl;
}


// ************************************************************************* //
Test-fvc.C (4,494 bytes)   

user4

2012-11-05 15:49

 

scalarTransportFoam.C (3,623 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/>.

Application
    scalarTransportFoam

Description
    Solves a transport equation for a passive scalar

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

#include "fvCFD.H"
#include "simpleControl.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createFields.H"

    simpleControl simple(mesh);

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

    Info<< "\nCalculating scalar transport\n" << endl;

    #include "CourantNo.H"

    while (simple.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        while (simple.correctNonOrthogonal())
        {
            {
                tmp<volVectorField> tgradT(fvc::grad(T));
                volVectorField& gradT = tgradT();
                const volVectorField::GeometricBoundaryField& bfld =
                    gradT.boundaryField();

                Pout<< "--- Coupled patchFields before ---" << endl;
                forAll(bfld, patchI)
                {
                    if (bfld[patchI].coupled())
                    {
                        Pout<< "Patch:" << bfld[patchI].patch().name() << nl
                            << "internal:" << bfld[patchI].patchInternalField()
                            << nl
                            << "neighb:" << bfld[patchI].patchNeighbourField()
                            << endl;
                    }
                }




                gradT.correctBoundaryConditions();

                Pout<< "--- Coupled patchFields after ---" << endl;
                forAll(bfld, patchI)
                {
                    if (bfld[patchI].coupled())
                    {
                        Pout<< "Patch:" << bfld[patchI].patch().name() << nl
                            << "internal:" << bfld[patchI].patchInternalField()
                            << nl
                            << "neighb:" << bfld[patchI].patchNeighbourField()
                            << endl;
                    }
                }
                Pout<< nl << endl;
            }


            solve
            (
                fvm::ddt(T)
              + fvm::div(phi, T)
              - fvm::laplacian(DT, T)
            );
        }

        runTime.write();
    }

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //
scalarTransportFoam.C (3,623 bytes)   

user4

2012-11-05 15:51

  ~0001765

I've attached scalarTransportFoam.C with a gradient calc and a simple testcase. I cannot see any difference in grad(T) on either side, i.e. patchNeighbourField == other-side's patchInternalField.

user534

2012-12-04 14:22

  ~0001806

I have asked Dominik to prepare a case that demonstrates the problem - it will take some time (apologies).

Hrv

dhora

2013-02-11 12:30

reporter   ~0001894

Some time ago a problem has been reported in the CFD-Online forum. I'm not sure but maybe it goes into the same direction.

http://www.cfd-online.com/Forums/openfoam-bugs/94085-strange-processor-boundary-behavior-linearupwindv.html

It can be reproduced with the pitzDaily test case.

Regards
David

dhora

2013-02-11 12:30

reporter  

user534

2013-02-11 13:01

  ~0001895

Dominik has produced a trivial code that demonstrates the error. I have fixed this a while back in 1.6-ext and a number of problems we have had just went away. Linear upwinding is one like that.

If you want a hot-fix, go to linear upwind and do

    // Note: in order for the patchNeighbourField to be correct on coupled
    // boundaries, correctBoundaryConditions needs to be called.
    // The call shall be moved into the library fvc operators
    gradVf.correctBoundaryConditions();

before the grad is used.

dhora

2013-02-11 16:01

reporter   ~0001897

Unfortunately the problem still persists, tested with 1.6-ext.

David

dhora

2013-02-11 19:28

reporter   ~0001898

Now fixed in 2.1.x (commit 3e32ea0c3432da6f1235f3879101782e30bd4291)

Thanks

user2

2013-02-12 09:08

  ~0001902

linearUpwindV scheme corrected for parallel running by commit 3e32ea

Issue History

Date Modified Username Field Change
2012-11-01 11:34 user534 New Issue
2012-11-05 14:57 user4 Note Added: 0001763
2012-11-05 15:04 user534 Note Added: 0001764
2012-11-05 15:08 user4 File Added: Test-fvc.C
2012-11-05 15:49 user4 File Added: scalarTransportFoam.C
2012-11-05 15:51 user4 Note Added: 0001765
2012-12-04 14:22 user534 Note Added: 0001806
2013-02-11 12:30 dhora Note Added: 0001894
2013-02-11 12:30 dhora File Added: pitzDaily_parallel-bug.tar.gz
2013-02-11 13:01 user534 Note Added: 0001895
2013-02-11 16:01 dhora Note Added: 0001897
2013-02-11 19:28 dhora Note Added: 0001898
2013-02-12 09:08 user2 Note Added: 0001902
2013-02-12 09:08 user2 Status new => resolved
2013-02-12 09:08 user2 Resolution open => fixed
2013-02-12 09:08 user2 Assigned To => user2