CFD Online Discussion Forums - Gauss linear gradient calculation discrepancy

Hello everyone,
I am trying to calculate in OF301 the gradient of a volScalarField without using fvc::grad (I set Gauss linear as scheme in the case study). Also, I have fully orthogonal hexahedrons.

So I try something like below:

Code:

// this is in setField

volVectorField gradT

    (

                IOobject

                (

                    "gradT",

                    runTime.timeName(),

                    fluidRegions[i],

                    IOobject::NO_READ,

                    IOobject::NO_WRITE

                ),

                fluidRegions[i],

                dimensionedVector("gradT", dimensionSet(0,-1,0,1,0,0,0), vector::zero)

    );

//



surfaceScalarField Tface = linearInterpolate(T);

    

    const fvMesh& mesh = Tface.mesh();

        const labelList& own = mesh.owner();

        const labelList& nei = mesh.neighbour();

        const vectorField& Sf = mesh.Sf();

    

    forAll(own, facei)

    {

                gradT[own[facei]] += Sf[facei]*Tface[facei]/mesh.V()[own[facei]];

                gradT[nei[facei]] -= Sf[facei]*Tface[facei]/mesh.V()[nei[facei]];

    }

    

    forAll(mesh.boundary(), patchi)

    {

         const labelList& pFaceCells =

             mesh.boundary()[patchi].faceCells();

         const vectorField& pSf = mesh.Sf().boundaryField()[patchi];

         const scalarField& pTface = Tface.boundaryField()[patchi];

 

         forAll(mesh.boundary()[patchi], facei)

         {

                gradT[pFaceCells[facei]] += pSf[facei]*pTface[facei]/mesh.V()[pFaceCells[facei]];

         }

    }

    

    gradT.correctBoundaryConditions();

    

    Info<< "\nDIFFERENCE\n" << max(mag(fvc::grad(T)-gradT)) << "\n" << endl;

which is more or less the same thing done by the member function gradf in gaussGrad.C.

However, if I try to compare the resulting gradient to the one calculated using the function fvc::grad, I see a considerable difference.

Am I missing something? Maybe the surfaceScalarField is not calculated through the linearInterpolate function?

Thank you very much
Cheers