[Mithrandirrr]

Dear all,

I need to compute fluxes over cell faces given a volVectorField. To fulfill the needs of my custom parallel code, the resulting surfaceScalarField requires consistent values across processor patches.

I would like to use a piece of code like this one:

 phi_ = fvc::interpolate(U_) & mesh_.Sf();

However, I'm running into problems since fvc::interpolate(U) does not yield consistent values on processor patches.

Hence my question: is there any direct way to interpolate a field from cell centers to face centers, with consistent values on processor patches?

To show the inconsistency that I mention, I attached a test case with a Cartesian mesh that has two grid cells and is decomposed onto two processors (i.e., one grid cell per processor). A velocity field U has value (0, 0, 0) on processor 0 and (1, 0, 0) on processor 1. As a result of fvc::interpolate(U), I would expected a value of (0.5, 0, 0) on the processor patch on both processors. However, the (stripped) output is as follows:

[1] phi = dimensions      [0 1 -2 0 0 0 0];
[1] 
[1] internalField   nonuniform 0();
[1] 
[1] boundaryField
[1] {
[1]     domainBoundary
[1]     {
[1]         type            calculated;
[1]         value           uniform (0 0 0);
[1]     }
[1]     procBoundary1to0
[1]     {
[1]         type            processor;
[1]         value           uniform (0.5 0 0);
[1]     }
[1] }
[0] phi = dimensions      [0 1 -2 0 0 0 0];
[0] 
[0] internalField   nonuniform 0();
[0] 
[0] boundaryField
[0] {
[0]     domainBoundary
[0]     {
[0]         type            calculated;
[0]         value           uniform (0 0 0);
[0]     }
[0]     procBoundary0to1
[0]     {
[0]         type            processor;
[0]         value           uniform (0 0 0);
[0]     }

The values on the patches procBoundary0to1 and procBoundary1to0 are not consistent. Please see the source file testProcessorPatchInterpolation.C for details. The test case can be compiled and run using

wmake
./Allrun

To circumvent this, I came up with a quick-and-dirty solution which is not satisfactory:

void
Foam::synchronizeProcessorPatchField
(
    volVectorField& vf
)
{
    const polyBoundaryMesh& patches = vf.mesh().boundaryMesh();

    FieldField<fvPatchField, vector> oldBoundaryField = vf.boundaryField();

    forAll(patches, patchI)
    {
        vf.boundaryField()[patchI].initEvaluate(Pstream::blocking);
    }
    
    forAll(patches, patchI)
    {
        const polyPatch& pp = patches[patchI];
        fvPatchField<vector>& pf = vf.boundaryField()[patchI];
        // The following call puts the values of the patch neighbour field
        // info pf.
        pf.evaluate(Pstream::blocking);

        if (isA<processorPolyPatch>(pp))
        {
            pf = 0.5*(pf.patchNeighbourField() + oldBoundaryField[patchI]);
        }
    }
}

This code is certainly a messy version of a feature which is present somewhere in the OpenFOAM source code. The use of 0.5 as a weighting factor is quite dirty, since it assumes linear interpolation on a Cartesian equidistant grid...

Any input is welcome!

Thanks,
Adrien

[Ranj]

Hi, has anyone ever found a solution this?

[olesen]

Originally Posted by Ranj

Hi, has anyone ever found a solution this?

You will need to check what is being done with the values later on. In general, the processor patches in OpenFOAM work a bit like a halo-value swap. That means that the the procBoundary0to1 contains values from the faces of processor1 and procBoundary1to0 contains values from the faces of processor0.


If, for your purposes, you need the averaged values for determining fluxes, you need to handle the coupled patches yourself.


For example,


https://develop.openfoam.com/Develop.../fluxSummary.C