CFD Online Discussion Forums - Difference between UEqn().H() & UEqn().H1()?

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

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

- - Difference between UEqn().H() & UEqn().H1()? (https://www.cfd-online.com/Forums/openfoam-programming-development/141592-difference-between-ueqn-h-ueqn-h1.html)

Difference between UEqn().H() & UEqn().H1()?

Dear Foamers,

what is the difference between UEqn().H() & UEqn().H1() equations?

Hi,

I think (from some tests that i did) that:

H1() is simply the diagonal coeffients divided by the volume and thus it is equal to diag()/mesh.V(). Moreover, it does not contain the contribution of the BC cells. So that's why it is different from A().

H() ? i am not sure, i asked for it too. But, i guess it should be as the following:

lets say we want to solve for matrix system: K*psi=S. OpenFoam re-arrange the system as: (A+H')*psi=S with A: diagonal coeffs. H': non-diagonal coeffs., and S: source terms coeffs. Then one can write:
A*psi = (S - H'*psi). I think (not sure) that the H() is the one as:
H= (S - H'*psi).

FOAMERS, Correct me if i am mistaken..

Regards,

T.D.

Quote:

Hi
Hi,

I think (from some tests that i did) that:

H1() is simply the diagonal coeffients divided by the volume and thus it is equal to diag()/mesh.V(). Moreover, it does not contain the contribution of the BC cells. So that's why it is different from A().

H() ? i am not sure, i asked for it too. But, i guess it should be as the following:

lets say we want to solve for matrix system: K*psi=S. OpenFoam re-arrange the system as: (A+H')*psi=S with A: diagonal coeffs. H': non-diagonal coeffs., and S: source terms coeffs. Then one can write:
A*psi = (S - H'*psi). I think (not sure) that the H() is the one as:
H= (S - H'*psi).

FOAMERS, Correct me if i am mistaken..

Regards,

T.D.

Please, excuse for the late reply.

I think, UEqn.H1( ) instead depends on the non-diagonal terms. Or, in more general way it is equal to non-diagonal()/mesh.V().

What do you think?

-
Best Regards!

Hi,

So what you are saying is that according to my previous posts, it will be:
H'= H1() .

Foamers, any confirmation (proof) ?

Regards,

T.D.

Quote:

Originally Posted by Tushar@cfd (Post 514619)

Please, excuse for the late reply.

I think, UEqn.H1( ) instead depends on the non-diagonal terms. Or, in more general way it is equal to non-diagonal()/mesh.V().

What do you think?

-
Best Regards!

I am not sure either, but I think it is more related to non-diagonal terms.

Correct me if I am wrong?

-
Best Regards!

Quote:

Originally Posted by Tushar@cfd (Post 514625)

I am not sure either, but I think it is more related to non-diagonal terms.

Correct me if I am wrong?

-
Best Regards!

Hello guys,

Cause recently I wanna do a simplecFoam for fun. This may relate to this topic. I tried to search the code, after lots of nested class(right? bad for my english..) I got this:

Code:

H1_.internalField() = lduMatrix::H1();



Hphi.internalField() += lduMatrix::H(psi_.internalField()) + source_;//this is for H()

https://github.com/OpenFOAM/OpenFOAM...rix/fvMatrix.C

Code:

template<class Type>

Foam::tmp<Foam::Field<Type> > Foam::lduMatrix::H(const Field<Type>& psi) const

{

    tmp<Field<Type> > tHpsi

    (

        new Field<Type>(lduAddr().size(), pTraits<Type>::zero)

    );



    if (lowerPtr_ || upperPtr_)

    {

        Field<Type> & Hpsi = tHpsi();



        Type* __restrict__ HpsiPtr = Hpsi.begin();



        const Type* __restrict__ psiPtr = psi.begin();



        const label* __restrict__ uPtr = lduAddr().upperAddr().begin();

        const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();



        const scalar* __restrict__ lowerPtr = lower().begin();

        const scalar* __restrict__ upperPtr = upper().begin();



        register const label nFaces = upper().size();



        for (register label face=0; face<nFaces; face++)

        {

            HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];

            HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];

        }

    }



    return tHpsi;

}



Foam::tmp<Foam::scalarField > Foam::lduMatrix::H1() const

{

    tmp<scalarField > tH1

    (

        new scalarField(lduAddr().size(), 0.0)

    );



    if (lowerPtr_ || upperPtr_)

    {

        scalarField& H1_ = tH1();



        scalar* __restrict__ H1Ptr = H1_.begin();



        const label* __restrict__ uPtr = lduAddr().upperAddr().begin();

        const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();



        const scalar* __restrict__ lowerPtr = lower().begin();

        const scalar* __restrict__ upperPtr = upper().begin();



        register const label nFaces = upper().size();



        for (register label face=0; face<nFaces; face++)

        {

            H1Ptr[uPtr[face]] -= lowerPtr[face];

            H1Ptr[lPtr[face]] -= upperPtr[face];

        }

    }



    return tH1;

}

https://github.com/OpenFOAM/OpenFOAM...ixOperations.C
https://github.com/OpenFOAM/OpenFOAM...rixTemplates.C

For now too many codes Im a little dizzy. For my CFD level this is endless to dig deeper...

However, I think the main difference is that H() includes the source when U apply FVM discretation on UEqn. But H1() does not.

From this equations(No need to post all of it here cuz its chinese), H1() does not consider that source term in Ueqn.

http://www.cfd-online.com/Forums/att...1&d=1432987688

I think this is rite from T.D.

Code:

So what you are saying is that according to my previous posts, it will be:

H'= H1() .

Correct me if im wrong?:eek:

Best,:D

Hello!
UEqn().H1() should be sum of all non-diagonal coefficients of matrix not including multiplication by Un, so
H1 = sum{An}
whereas
H = sum{An*Un}+source

Here a good explanation of the usage and computation of H1 is given:

http://hobbyfoam.blogspot.com/2016/0...-openfoam.html