CFD Online Discussion Forums - problems of `patchNeighbourField()` in jump condition

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - problems of `patchNeighbourField()` in jump condition (https://www.cfd-online.com/Forums/openfoam-programming-development/197433-problems-patchneighbourfield-jump-condition.html)

problems of `patchNeighbourField()` in jump condition

In the implementation of `patchNeighbourField()` at line 94 of file jumpCyclicFvPatchField.C

It add or subtract the jump to the value in neighbor cells.

In `fvMatrix`, which represent $A\cdot x=s, x_0 = \psi$ , so if the psi_ is zero, the `Amul()` result, which represents $A\cdot x$ must be zero.

However, if the problem contains `jumpCyclicFvPatchField` boundary condition, which is a coupled interface, and `Amul()` will represent $A\cdot x - B^CT(x')$ , which contains coupled interface's contribution. Now, if the psi_ is zero, `Amul()` is no longer zero. Which violates the linear algebra meaning of `Amul()`.

I think there is a limitation on the design of fvMatrix's treatment of coupled BC, which cannot include boundary source term independent of the values of both sides.

Another problem is: mathematically, jumpCyclic BC must enforce the jump condition between patch faces. However, in the implementation, it seems that the jump condition are enforced between the inner cell and neighbour cells, which is obviously incorrect.

Code:

template<class Type>

Foam::tmp<Foam::Field<Type>>

Foam::jumpCyclicFvPatchField<Type>::patchNeighbourField() const

{

    const Field<Type>& iField = this->primitiveField();

    const labelUList& nbrFaceCells =

        this->cyclicPatch().neighbFvPatch().faceCells();



    tmp<Field<Type>> tpnf(new Field<Type>(this->size()));

    Field<Type>& pnf = tpnf.ref();



    Field<Type> jf(this->jump());

    if (!this->cyclicPatch().owner())

    {

        jf *= -1.0;

    }



    if (this->doTransform())

    {

        forAll(*this, facei)

        {

            pnf[facei] = transform

            (

                this->forwardT()[0], iField[nbrFaceCells[facei]]

            ) - jf[facei];

        }

    }

    else

    {

        forAll(*this, facei)

        {

            pnf[facei] = iField[nbrFaceCells[facei]] - jf[facei];

        }

    }



    return tpnf;

}

> mathematically, jumpCyclic BC must enforce the jump condition between patch faces

Nonsense: the two faces are actually one. Therefore the jump is the jump between two cell centres.

> Which violates the linear algebra meaning of `Amul()`.

Another nonsense. What is the "mathematical definition" of the Amul function??? The code is exactly as it should be and it works fine. I use it every day.

Hrv

I think the jump condition (in mathematics) must be like this for 1D problem in interval [0,1]:
$p(x=0) = p(x=1)+\delta p$
However, the implementation is :
$p(x=0+\Delta x/2) = p(x=1-\Delta x/2)+\delta p$

Which is dependent on grid spacing $\Delta x$ . It may impair the spatial order of accuracy of this BC.

about Amul(), by the document and its usage in Krylov type solvers. It must be the definition of $A\cdot x$ . which means the output must be zero vector if the x is a zero vector, which is what I mean by the term "linear algebra meaning"

However, in the implementation of ldu::Amul(), it invokes the interfaces' "updateInterfaceMatrix()" method, which essentially is adding the product of boundaryCoeffs and neighbor side cell values (transformed if neccesary).

However, according to the implementation of jumpCyclic method "updateInterfaceMatrix()", if the interface is the jumpCyclic condition, The cell value transferred in is not zero if the state vector is zero. So I say it violates the "linear algebra meaning" of Amul().

Maybe the implementation of "Amul()" is compatible with other codes because most solvers actually works.