CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM

why in solve(UEqn == -fvc::grad(p)) UEqn.source() is not modified and UEqn.H() chang

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree4Likes
  • 4 Post By mAlletto

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   January 1, 2019, 11:26
Default why in solve(UEqn == -fvc::grad(p)) UEqn.source() is not modified and UEqn.H() chang
  #1
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Posts: 184
Rep Power: 2
mAlletto is on a distinguished road
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



Code:

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:


Code:
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 is offline   Reply With Quote

Old   January 1, 2019, 11:59
Default
  #2
Senior Member
 
Michael Alletto
Join Date: Jun 2018
Posts: 184
Rep Power: 2
mAlletto is on a distinguished road
Ok maybe I can answer my question by myself:


When the matrix



Code:
    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:


Code:
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:


Code:
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
hjasak, snak, rob3rt 0ng and 1 others like this.
mAlletto is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 06:33.