[nimasam]

hi openFoamer

in compressibleInterFoam and in pEqn.H there are following syntax:

 dgdt =
                (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
               *(pEqnComp & p);

can any body tell me whats the meaning of "pEqnComp & p" ?

pEqnComp : fvScalarMatrix
p : volScalarField

so whats "&" here ?

[marupio]

I think it's the inner product.

[nimasam]

im confused how we can inner product fvScalarMatrix with volScalarField?
what is the output and how this calculation can be done? give me a mathematical concept please

[marupio]

fvMatrix.C implements operator& around line 2270. The result is a GeometricField. I'm not exactly sure what it does, but it looks like it affects the matrix. I don't have my Foaming computer in front of me at the moment.

[nimasam]

& operator is defined in fvMatrix.C like below:

template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::operator&
(
    const fvMatrix<Type>& M,
    const DimensionedField<Type, volMesh>& psi
)
{
    tmp<GeometricField<Type, fvPatchField, volMesh> > tMphi
    (
        new GeometricField<Type, fvPatchField, volMesh>
        (
            IOobject
            (
                "M&" + psi.name(),
                psi.instance(),
                psi.mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            psi.mesh(),
            M.dimensions()/dimVol,
            zeroGradientFvPatchScalarField::typeName
        )
    );
    GeometricField<Type, fvPatchField, volMesh>& Mphi = tMphi();

    // Loop over field components
    for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
    {
        scalarField psiCmpt = psi.field().component(cmpt);
        scalarField boundaryDiagCmpt(M.diag());
        M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
        Mphi.internalField().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
    }

    Mphi.internalField() += M.lduMatrix::H(psi.field()) + M.source();
    M.addBoundarySource(Mphi.internalField());

    Mphi.internalField() /= -psi.mesh().V();
    Mphi.correctBoundaryConditions();

    return tMphi;
}

1) could you please tell me what does this piece of code do?

it seems it loops over the field to make it diagonal predominate, but its confusing!

2)why it returns tMphi while it does all operation on Mphi? i feel it should be Mphi

[ngj]

add 2) tMphi is a tmp<Type>, whereas Mphi (of type Type) is a reference to tMphi. Hence all changes to Mphi is directly a change to tMphi.

tMphi is returned, as the tmp class is included in OpenFOAM to optimise the memory handling. Since volFieldsType> can potentially be very large, the return of a tmp<volField<Type> > speeds up the run-time.

Kind regards,

Niels

[marupio]

It's probably the inner product. When it involves a fvMatrix, it apparently does something to the Matrix's diagonal. I'm not going to try to decode it because I agree: it is confusing.

[marupio]

Hi Nima, sorry for being dismissive in the last response. If you want to go into detail about what the & operator does, I've added some comments in blue that might help you out. It looks like it doesn't actually affect the matrix. It looks like it is the inner product between the field associated with the matirx, with the field p you mention. Since a matrix is involved, it has to account for boundary conditions, and I think that's where it gets confusing, especially when it comes to boundary patches that have out-of-core multiplication (cyclic, processor, etc). Also see below for links to useful pages that might help further explain things.

template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
Foam::operator&
(
    const fvMatrix<Type>& M,
    const DimensionedField<Type, volMesh>& psi
)
{
    // Here we create the return object
    tmp<GeometricField<Type, fvPatchField, volMesh> > tMphi
    (
        new GeometricField<Type, fvPatchField, volMesh>
        (
            IOobject
            (
                "M&" + psi.name(),
                psi.instance(),
                psi.mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            psi.mesh(),
            M.dimensions()/dimVol,
            zeroGradientFvPatchScalarField::typeName
        )
    );
    // This is for convenience, because it is awkward dealing directly with tmp<>
    // See link below.
     GeometricField<Type, fvPatchField, volMesh>& Mphi = tMphi();

    // Loop over field components
    for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
    {
        // psiCmpt is your "p"
        scalarField psiCmpt = psi.field().component(cmpt);
        // Here we're grabbing a copy of the matrix diagonal of your "pEqnComp"
        scalarField boundaryDiagCmpt(M.diag());
        // This function brings in out-of-core multiplication effects from the boundary
        // (e.g. cyclic patches) and puts it into our copy of the diagonal
        M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
        // Multiplying p by the modified diagonal, store result in output object
        Mphi.internalField().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
    }

    // H operator stuff.  Shorthand operation built into the matrix that I never
    // remember what it does... see link below
     Mphi.internalField() += M.lduMatrix::H(psi.field()) + M.source();
    M.addBoundarySource(Mphi.internalField());

    // Per unit volume - the matrix is integrated over the volume; a field is not
    Mphi.internalField() /= -psi.mesh().V();
    Mphi.correctBoundaryConditions();

    return tMphi;
}

For H operator stuff, see this page:
http://openfoamwiki.net/index.php/Op...ide/H_operator

For what tmp<> is all about, see here:
http://openfoamwiki.net/index.php/OpenFOAM_guide/tmp