View Issue Details

IDProjectCategoryView StatusLast Update
0001160OpenFOAMBugpublic2014-03-07 10:08
Reporteruser528Assigned Touser4 
PrioritynormalSeveritymajorReproducibilityalways
Status resolvedResolutionfixed 
PlatformLinuxOSSUSE Linux Enterprise ServerOS Version11.1
Summary0001160: fvc::reconstruct produces different results in parallel and serial
Descriptionfvc::reconstruct results differ in parallel and serial mode on non-uniform meshes with non-uniform surfaceScalarField
Steps To Reproduce1.1 Run the attached solver on the attached case:
$ ReconstructTest
1.2 Note that testRec field becomes
internalField nonuniform List<vector> 2((1 0 0) (1.6 0 0));

2.1 Decompose the case (decomposeParDict is already in place):
$ decomposePar
2.2 Run the test in parallel
$ mpirun -np 2 ReconstructTest -parallel
2.3 Reconstruct parallel case
$ reconstructPar
2.4 Note that testRec field becomes
internalField nonuniform List<vector> 2((1 0 0) (1.5 0 0));

So, the reconstruct results are substantially different, although I expected the same results.
Additional InformationTh solver just reads "test" surfaceScalarField and calls fvc::reconstruct to reconstruct it to vector field. If the mesh is a regular rectangular mesh, everything is ok. If the "test" surfaceScalarField is uniform, everything is ok.

However, in this case, the mesh consists of two cells: a prism and a cube attached to one cell of the prism, and the surfaceScalarField is non-uniform (specifically, value on the right boundary is different from internal value), and the results of fvc::reconstruct are different.

I understand that this might not pose a problem in real applications, because one might say that if the values of a surfaceScalarField are different on opposing parallel sides of a cell, than the reconstructed value inside can be anything within a certain range. However, if I have a cube with edge length equal to one, and I have different values of the surfaceScalarField on the opposing faces, I expect the corresponding component of the reconstructed vector to be exactly the average between the two values. This is so in the parallel run, but not in serial.

Anyway, above all I expect the results of serial and parallel run to be the same.

I tried to find the reasons for this behavior, and looking into fvc::reconstruct code and adding debug output there, I found that it starts differing (serial run vs processor1) after faceVols are corrected by 1-weights, and weights differ because in serial run the face is owned by the left cell, while in parallel run the face on processor 1 is owned by the only cell on this processor, that is, by the right cell of the global mesh.

However, I am not sure whether the problem is that mesh.weights differ, or should this difference be compensated somehow by fvc::reconstruct code.

TagsNo tags attached.

Activities

user528

2014-02-17 15:35

 

ReconstructTest.zip (12,162 bytes)

user528

2014-02-18 15:02

  ~0002846

I have studied this a bit more.

Firstly, it seems to be not a problem of parallel vs. serial, but a problem of owner vs neighbor cell. For example, if in the test case I have attached, I change the order of cells in blockMeshDict, thus inverting owner/neighbor of the only internal face of the mesh, then the results do change once again.

I have looked into fvc::reconstruct code, and, although I still do not understand it completely, I think I have a correction that should be applied.

This is the code from fvc::reconstruct that defines the result:
 inv(surfaceSum(mesh.Sf()*faceVols))&surfaceSum(faceVols*ssf)

So, I understand that for each cell two values are constructed: tensor T and (assuming Type=scalar) vector A:
T=sum over faces of cell of (mesh.Sf()*faceVols)
A=sum over faces of cell of (faceVols*ssf)
and the answer is T^{-1}&A.

In this form it seems that all faces of a cell should be treated in a similar way while computing the sums. Such similarity is seen in mesh.Sf() and ssf, but on a closer inspection faceVols do not present such a similarity.

Indeed, faceVols is firstly defined as
mesh.Sf()/(mesh.magSf()*mesh.nonOrthDeltaCoeffs()),
but then it is corrected by (1.0-weights) factor.

mesh.weights is the ratio of "neighbor cell center to face center" distance to "neighbor cell center to owner cell center" distance (here for simplicity I assume mesh to be orthogonal). So (1.0-weights) is the ration of "owner cell center to face center" distance to "neighbor cell center to owner cell center".

Now assume we are calculating T and A for a particular cell. For some faces of this cell, this cell is the owner, and (1.0-weights) factor will be "distance from this cell center to face center"/"distance from this cell center to other (neighbor) cell center". But for other faces this cell will be a "neighbor" cell and (1.0-weights) factor will be "distance from other cell center to face center"/"distance from other cell center to this cell center".

So, the physical meaning of faceVols factor inside the sums over faces is definetely different depending on whether this cell is an owner of a face or a neighbor.

I would assume it to be sensible to either always take "distance from this cell center to face center" either always take "distance from other cell center to face center". In particular, this means that the (1.0-weights) factor should be replaced by different factors for different sides of a face.

Judging that for a boundary face a cell is always an owner, and that (1.0-weight) is used in that case, I can assume that (1.0-weight) should be used for the owners and therefore simply (weight) should be used for neighbors.

I have rewritten the fvc::reconstruct code to do this correction, see an uploaded file. This version produces consistent results for me, and satisfies the "unary cube" test described in the bug report, that is, for an unary cube cell the components of fvc::reconstruct result are averages between the values of field on the faces.

user528

2014-02-18 15:06

 

fvcReconstruct.C (5,389 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 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 "fvcReconstruct.H"
#include "fvMesh.H"
#include "zeroGradientFvPatchFields.H"

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

namespace Foam
{

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

namespace fvc
{

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


template<class Type>
tmp
<
    GeometricField
    <
        typename outerProduct<vector,Type>::type, fvPatchField, volMesh
    >
>
reconstruct
(
    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
)
{
    typedef typename outerProduct<vector, Type>::type GradType;

    const fvMesh& mesh = ssf.mesh();

    surfaceVectorField faceVols
    (
        mesh.Sf()/(mesh.magSf()*mesh.nonOrthDeltaCoeffs())
    );
    
    const surfaceScalarField& weights = mesh.weights();
    
    volTensorField T(
        IOobject
            (
                "T",
                ssf.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
        mesh,
        dimensionedTensor("0", dimVolume, pTraits<tensor>::zero),
        "zeroGradient"
    );
    GeometricField <GradType, fvPatchField, volMesh> A(
        IOobject
            (
                "A",
                ssf.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
        mesh,
        dimensioned<GradType>("0", dimLength * ssf.dimensions(), pTraits<GradType>::zero),
        "zeroGradient"
    );
    const labelUList& owner = mesh.owner();
    const labelUList& neighbour = mesh.neighbour();
    forAll(faceVols, facei) {
        T[owner[facei]] += mesh.Sf()[facei] * faceVols[facei] * (1-weights[facei]);
        T[neighbour[facei]] += mesh.Sf()[facei] * faceVols[facei] * weights[facei];
        
        A[owner[facei]] += faceVols[facei] * ssf[facei]* (1-weights[facei]);
        A[neighbour[facei]] += faceVols[facei] * ssf[facei]* weights[facei];
    }
    forAll(faceVols.boundaryField(), patchi) {
        const vectorField& pSf = mesh.Sf().boundaryField()[patchi];
        const vectorField& pFaceVols = faceVols.boundaryField()[patchi];
        const scalarField& pWeights = weights.boundaryField()[patchi];
        const Field<Type>& pssf = ssf.boundaryField()[patchi];
        const labelUList& pFaceCells = mesh.boundary()[patchi].faceCells();
        
        forAll(faceVols.boundaryField()[patchi], facei) {
            tensor addT = pSf[facei] * pFaceVols[facei];
            GradType addA = pFaceVols[facei] * pssf[facei];
            if (faceVols.boundaryField()[patchi].coupled()) {
                addT *= (1-pWeights[facei]);
                addA *= (1-pWeights[facei]);
            }
            
            T[pFaceCells[facei]] += addT;
            A[pFaceCells[facei]] += addA;
        }
    }
    T.correctBoundaryConditions();
    A.correctBoundaryConditions();
            
    tmp<GeometricField<GradType, fvPatchField, volMesh> > treconField
    (
        new GeometricField<GradType, fvPatchField, volMesh>
        (
            IOobject
            (
                "volIntegrate("+ssf.name()+')',
                ssf.instance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            inv(T) & A,
            zeroGradientFvPatchField<GradType>::typeName
        )
    );

    treconField().correctBoundaryConditions();

    return treconField;
}

template<class Type>
tmp
<
    GeometricField
    <
        typename outerProduct<vector, Type>::type, fvPatchField, volMesh
    >
>
reconstruct
(
    const tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >& tssf
)
{
    typedef typename outerProduct<vector, Type>::type GradType;
    tmp<GeometricField<GradType, fvPatchField, volMesh> > tvf
    (
        fvc::reconstruct(tssf())
    );
    tssf.clear();
    return tvf;
}


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

} // End namespace fvc

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

} // End namespace Foam

// ************************************************************************* //
fvcReconstruct.C (5,389 bytes)   

user528

2014-02-18 15:21

  ~0002849

By the way, I have suddenly realized that faceVols is indeed (as the name suggest) something like the volume of tetrahedron formed by the face and cell center. But in such a case it is clear that it is different for different sides of mesh, and the correction I applied seems to be correct.

user4

2014-03-07 10:08

  ~0002936

Fixed in 2.3.0

Issue History

Date Modified Username Field Change
2014-02-17 15:35 user528 New Issue
2014-02-17 15:35 user528 File Added: ReconstructTest.zip
2014-02-18 15:02 user528 Note Added: 0002846
2014-02-18 15:06 user528 File Added: fvcReconstruct.C
2014-02-18 15:21 user528 Note Added: 0002849
2014-03-07 10:08 user4 Note Added: 0002936
2014-03-07 10:08 user4 Status new => resolved
2014-03-07 10:08 user4 Fixed in Version => 2.3.x
2014-03-07 10:08 user4 Resolution open => fixed
2014-03-07 10:08 user4 Assigned To => user4