[Bernhard]

Recently, I have been using the components utility of foamCalc on a vector field (which itself is a gradient of a scalar). I am only interested in one of the components of this vector (or gradient) field. I calculate this component using two methods

1. Directly in the solver (i.e. by writing field.component(0))
2. Using foamCalc components field. Which, essentially, should do the same thing, but read the grad/vector-field from disk.

The cell values are identical for both approaches, however, looking at the value of the component at a patch, then both approach give me a different result.

The component value returned by the solver is the correct one, as it matches the analytical solution for a test case. The value obtained by calculating the components after reading the vector-field from disk are different. In the latter case, foamCalc seems to take into account that the (default?) zero-gradient condition, and thus returns me the cell value. This is expected behavior, and probably correct based on the input vector-field read from disk. However, it should not make a difference whether or not the field is re-read or not (in my understanding). Thus, is there information not written, when writing the gradient field?

I've attached a case and solver (modified laplacianFoam solver), with which the issue can be easily reproduced (for 2.1.0 and 2.2.0 at least). The tarball contains a "Allrun" script, which compiles, run and compares automatically.

Am I overlooking something? How can I get the foamCalc-route provide me with the correct solution? What information is lost when writing just the gradient-field?

[ngj]

Hi Bernhard,

I just gave your question a quick look. The issue is due to the file

OpenFOAM/OpenFOAM-2.2.0/src/postProcessing/foamCalcFunctions/field/components/writeComponentFields.C

which is called from

OpenFOAM/OpenFOAM-2.2.0/src/postProcessing/foamCalcFunctions/field/components/components.C

i.e. the method you invoke with foamCalc. The above mentioned file only contains the following method:

template <class Type>
void Foam::calcTypes::components::writeComponentFields
(
    const IOobject& header,
    const fvMesh& mesh,
    bool& processed
)
{
    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;

    if (header.headerClassName() == fieldType::typeName)
    {
        Info<< "    Reading " << header.name() << endl;
        fieldType field(header, mesh);

        for (direction i=0; i<Type::nComponents; i++)
        {
            Info<< "    Calculating " << header.name()
                << Type::componentNames[i] << endl;

            volScalarField componentField
            (
                IOobject
                (
                    header.name() + word(Type::componentNames[i]),
                    mesh.time().timeName(),
                    mesh,
                    IOobject::NO_READ
                ),
                field.component(i)
            );
            componentField.write();
        }

        processed = true;
    }
}

It is worth to notice that no information on the boundary condition are transfered from the field to the component output, i.e. the outputted volScalarField will assume that a zeroGradient is valid for the boundary.

Essentially, the solution of your problem will be to rewrite this particular function, such that the boundary types are included. This will subsequently require a loop, where the boundary fields are populated with the correct values.

Kind regards,

Niels

[Bernhard]

I see the point you make. However, then I don't understand how the write.H in the solver I attached does transfer the boundary condition. This is not explicitly done, and the written "grad(T)"-file has already zeroGradient boundary conditions, so is there anything left to transfer here? Or should grad(T) be written differently? I don't see a difference between the write.H file and the foamCalc components functions. (except for AUTO_WRITE versus component.write(), but I checked, this is not causing the difference))

[ngj]

Hi Bernhard,

Oh, I see your point. I will try to get some time during the weekend for a closer investigation.

Kind regards,

Niels

[ngj]

Hi Bernhard,

I have had a bit of fun, have been quite confused and now I think I have a waterproof explanation:

First of all:
Forget my entire post above, because it is not related to that!

Explanation:
The only file, which we have to consider is:

OpenFOAM/OpenFOAM-2.2.0/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C

The gradient is computed in two steps via this method (I added the comments in the source):

template<class Type>
Foam::tmp
<
    Foam::GeometricField
    <
        typename Foam::outerProduct<Foam::vector, Type>::type,
        Foam::fvPatchField,
        Foam::volMesh
    >
>
Foam::fv::gaussGrad<Type>::calcGrad
(
    const GeometricField<Type, fvPatchField, volMesh>& vsf,
    const word& name
) const
{
    typedef typename outerProduct<vector, Type>::type GradType;

    // STEP 1:
    tmp<GeometricField<GradType, fvPatchField, volMesh> > tgGrad
    (
        gradf(tinterpScheme_().interpolate(vsf), name)
    );
    GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad();

    // STEP 2:
    correctBoundaryConditions(vsf, gGrad);

    return tgGrad;
}

Step 1: This step computes the internal gradient and the boundary condition is set as zeroGradient.

Step 2: This step adds the non-zeroGradient component to the boundary field, i.e. it is temporary non-zeroGradient!

In your code, where you have the line

volVectorField gradT(fvc::grad(T));

// write all component fields

this is successful, because the boundary condition of gradT is the non-zeroGradient value. This is where my previous post is wrong, because the boundary values are carried over in the construction of the volScalarField. A quick look in the source code confirmed this.

If, however, you replace your code by this:

volVectorField gradT(fvc::grad(T));
gradT.correctBoundaryConditions();

// write all component fields

it results in a re-evaluation of the boundary field, which is set to zeroGradient; which is correct considering the type, which is assigned during the construction in

OpenFOAM/OpenFOAM-2.2.0/src/finiteVolume/finiteVolume/gradSchemes/gaussGrad/gaussGrad.C

Upon writing the field gradT, the writer is only concerned with the fact that it is a zeroGradient field, i.e. this is how it is written.

I tried to see, whether I could get a fixedGradient in the construction instead, but I was not immediately successful. This does not, however, mean that it is not possible.

Kind regards,

Niels

[Bernhard]

Aah, I had some suspicion that somewhere a correctBoundaryCondition had to come into play. Thanks for digging into this issue. I will need some time to consider my options to circumvent this feature for my applications

If I come up with any useful implementation, I will obviously share it here.