problems of `patchNeighbourField()` in jump condition

chengdi · January 5, 2018, 01:04

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;
}

hjasak · January 7, 2018, 10:19

> 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

chengdi · January 7, 2018, 10:46

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.

January 5, 2018, 01:04	problems of `patchNeighbourField()` in jump condition	#1
chengdi Member Di Cheng Join Date: May 2010 Location: Beijing, China Posts: 47 Rep Power: 16	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; }

January 7, 2018, 10:19		#2
hjasak Senior Member Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,905 Rep Power: 33	> 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 __________________ Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
[Commercial meshers] Boundary condition problems (OpenFOAM)	Milos	OpenFOAM Meshing & Mesh Conversion	13	October 13, 2016 19:58
Divergent with rhoSimpleFoam and the boundary condition problems	qjh888	OpenFOAM Running, Solving & CFD	0	May 17, 2016 20:31
boundary condition: problems with outlet at 0.99 atm when doing a trasient simulation	umair64	CFX	1	August 17, 2015 18:11
Subsea hole release boundary condition problems	dotapro	CFX	5	March 22, 2015 04:45
CFX problems with supersonic inlet condition - Inlet values in CFX-Post are wrong	jannnesss	CFX	5	February 25, 2011 16:24

January 7, 2018, 10:46		#3
chengdi Member Di Cheng Join Date: May 2010 Location: Beijing, China Posts: 47 Rep Power: 16	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.