[pod]

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:

// 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