# 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

 January 1, 2019, 12:26 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 Location: Bremen Posts: 603 Rep Power: 14 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 Foam::tmp> Foam::operator== ( const fvMatrix& A, const DimensionedField& su ) { checkMethod(A, su, "=="); tmp> tC(new fvMatrix(A)); tC.ref().source() += su.mesh().V()*su.field(); return tC; } template Foam::tmp> Foam::operator== ( const fvMatrix& A, const tmp>& tsu ) { checkMethod(A, tsu(), "=="); tmp> tC(new fvMatrix(A)); tC.ref().source() += tsu().mesh().V()*tsu().field(); tsu.clear(); return tC; } template Foam::tmp> Foam::operator== ( const fvMatrix& A, const tmp>& tsu ) { checkMethod(A, tsu(), "=="); tmp> tC(new fvMatrix(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(A)) unchanged the function solve executes fvm.solve, clears the space and returns some statistics: Code: ```template Foam::SolverPerformance Foam::solve(const tmp>& tfvm) { SolverPerformance solverPerf = const_cast&>(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)); amuzeshi and Kummi like this.
 January 1, 2019, 12:59 #2 Senior Member   Michael Alletto Join Date: Jun 2018 Location: Bremen Posts: 603 Rep Power: 14 Ok maybe I can answer my question by myself: When the matrix Code: ``` tmp 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 Foam::fvMatrix::fvMatrix ( const GeometricField& 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 for field " << psi_.name() << endl; } // Initialise coupling coefficients forAll(psi.mesh().boundary(), patchi) { internalCoeffs_.set ( patchi, new Field ( psi.mesh().boundary()[patchi].size(), Zero ) ); boundaryCoeffs_.set ( patchi, new Field ( psi.mesh().boundary()[patchi].size(), Zero ) ); } // Update the boundary coefficients of psi without changing its event No. GeometricField& psiRef = const_cast&>(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 Foam::tmp> Foam::fvMatrix::H() const { tmp> tHphi ( new GeometricField ( IOobject ( "H("+psi_.name()+')', psi_.instance(), psi_.mesh(), IOobject::NO_READ, IOobject::NO_WRITE ), psi_.mesh(), dimensions_/dimVol, extrapolatedCalculatedFvPatchScalarField::typeName ) ); GeometricField& Hphi = tHphi.ref(); // Loop over field components for (direction cmpt=0; cmpt() ); for (direction cmpt=0; cmpt