[mAlletto]

Hello Foamers,


I have a fundamental question understanding simpleFoam.


when solve(UEqn == -fvc::grad(p)); is called in simpleFoam the UEqn.source() not modified (i printed the UEqn.source() for a small cavity flow) but UEqn.H() is modified.


from the source code in fvMatrix.C this is not to see.



The operator == is defined as follows

template<class Type>
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
(
    const fvMatrix<Type>& A,
    const DimensionedField<Type, volMesh>& su
)
{
    checkMethod(A, su, "==");
    tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
    tC.ref().source() += su.mesh().V()*su.field();
    return tC;
}

template<class Type>
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
(
    const fvMatrix<Type>& A,
    const tmp<DimensionedField<Type, volMesh>>& tsu
)
{
    checkMethod(A, tsu(), "==");
    tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
    tC.ref().source() += tsu().mesh().V()*tsu().field();
    tsu.clear();
    return tC;
}

template<class Type>
Foam::tmp<Foam::fvMatrix<Type>> Foam::operator==
(
    const fvMatrix<Type>& A,
    const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu
)
{
    checkMethod(A, tsu(), "==");
    tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A));
    tC.ref().source() += tsu().mesh().V()*tsu().primitiveField();
    tsu.clear();
    return tC;
}

as far as I understood it returns a new fvMatrix with the source term modified leaving the
fvMatrix<Type>(A)) unchanged

the function solve executes fvm.solve, clears the space and returns some statistics:

template<class Type>
Foam::SolverPerformance<Type> Foam::solve(const tmp<fvMatrix<Type>>& tfvm)
{
    SolverPerformance<Type> solverPerf =
        const_cast<fvMatrix<Type>&>(tfvm()).solve();

    tfvm.clear();

    return solverPerf;
 }

From this analysis I would guess that the matrix UEqn is solved with the addition of the pressure gradient but the H operator should not be changed. This is obviously in contrast from the printing of H and source during the iteration.


What I also noticed is that the field U is also modified after calling solve(UEqn == -fvc::grad(p));

[mAlletto]

Ok maybe I can answer my question by myself:


When the matrix

    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
      + MRF.DDt(U)
      + turbulence->divDevReff(U)
     ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();

is constructed the solution vector psi_ is defined as a reference to U:

template<class Type>
Foam::fvMatrix<Type>::fvMatrix
(
    const GeometricField<Type, fvPatchField, volMesh>& psi,
    const dimensionSet& ds
)
:
    lduMatrix(psi.mesh()),
    psi_(psi),
    dimensions_(ds),
    source_(psi.size(), Zero),
    internalCoeffs_(psi.mesh().boundary().size()),
    boundaryCoeffs_(psi.mesh().boundary().size()),
    faceFluxCorrectionPtr_(nullptr)
{
    if (debug)
    {
        InfoInFunction
            << "Constructing fvMatrix<Type> for field " << psi_.name() << endl;
    }

    // Initialise coupling coefficients
    forAll(psi.mesh().boundary(), patchi)
    {
        internalCoeffs_.set
        (
            patchi,
            new Field<Type>
            (
                psi.mesh().boundary()[patchi].size(),
                Zero
            )
        );

        boundaryCoeffs_.set
        (
            patchi,
            new Field<Type>
            (
                psi.mesh().boundary()[patchi].size(),
                Zero
            )
        );
    }

    // Update the boundary coefficients of psi without changing its event No.
    GeometricField<Type, fvPatchField, volMesh>& psiRef =
       const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);

    label currentStatePsi = psiRef.eventNo();
    psiRef.boundaryFieldRef().updateCoeffs();
    psiRef.eventNo() = currentStatePsi;
}

after that the function solve(UEqn == -fvc::grad(p)); is called. It indeed solves a new matrix (the old one where to the already existing source the pressure gradient is added).


But when the matrix is solved, the solution vector psi_ is modified which is a reference to U. This is in line with the observation that U is modified after the solution of the matrix.


Hence also the solution vector psi_ of UEqn is also modified. But the source of UEqn remains unchanged.



since the H() is computed from the solution vector psi_ and the source:

template<class Type>
Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>>
Foam::fvMatrix<Type>::H() const
{
    tmp<GeometricField<Type, fvPatchField, volMesh>> tHphi
    (
        new GeometricField<Type, fvPatchField, volMesh>
        (
            IOobject
            (
                "H("+psi_.name()+')',
                psi_.instance(),
                psi_.mesh(),
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            psi_.mesh(),
            dimensions_/dimVol,
            extrapolatedCalculatedFvPatchScalarField::typeName
        )
    );
    GeometricField<Type, fvPatchField, volMesh>& Hphi = tHphi.ref();

    // Loop over field components
    for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
    {
        scalarField psiCmpt(psi_.primitiveField().component(cmpt));

        scalarField boundaryDiagCmpt(psi_.size(), 0.0);
        addBoundaryDiag(boundaryDiagCmpt, cmpt);
        boundaryDiagCmpt.negate();
        addCmptAvBoundaryDiag(boundaryDiagCmpt);

        Hphi.primitiveFieldRef().replace(cmpt, boundaryDiagCmpt*psiCmpt);
    }

    Hphi.primitiveFieldRef() += lduMatrix::H(psi_.primitiveField()) + source_;
    addBoundarySource(Hphi.primitiveFieldRef());

    Hphi.primitiveFieldRef() /= psi_.mesh().V();
    Hphi.correctBoundaryConditions();

    typename Type::labelType validComponents
    (
        psi_.mesh().template validComponents<Type>()
    );

    for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
    {
        if (validComponents[cmpt] == -1)
        {
            Hphi.replace
            (
                cmpt,
                dimensionedScalar("0", Hphi.dimensions(), 0.0)
            );
        }
    }

    return tHphi;
}

this explains why the H() field is changed and source not