# 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
 Thread Tools Search this Thread Show Printable Version Email this Page Search this Thread: Advanced Search Display Modes Linear Mode Switch to Hybrid Mode Switch to Threaded Mode 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 OffTrackbacks are Off Pingbacks are On Refbacks are On Forum Rules