- **OpenFOAM Running, Solving & CFD**
(*http://www.cfd-online.com/Forums/openfoam-solving/*)

- - **A doubt on solve**
(*http://www.cfd-online.com/Forums/openfoam-solving/58831-doubt-solve.html*)

Hello
I am trying to write Hello
I am trying to write a few simple codes to understand OpenFoam. I can't seem to find details about solve(), in all it's overloaded forms in the documentation. So i looked at scalarTransportFoam and icoFoam. Here are 2 snippets of code taken from OpenFoam-1.4.1/applications/solvers/incompressible OpenFoam-1.4.1/applications/solvers/basic. 1)solve ( fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(DT, T) ); Does this mean solve takes an input of type fvscalarMatrix, and somehow modifies value of T. So does that mean T has solve as it's friend? 2)Also in icoFoam i come across the line fvVectorMatrix UEqn ( fvm::ddt(U) + fvm::div(phi, U) - fvm::laplacian(nu, U) ); solve(UEqn == -fvc::grad(p)); The == really confuses me, since i can't find the relevant overloaded form for this operator. Also what is the return type of the expression Ueqn ==-fvc::grad(p) ? I am assuming solve magically changes values of U etc when it is invoked. I would appreciate it if someone would clarify this or at least point me to the right place to look/ documentation, as the programmersguide doesn't seem to touch this, and i have been staring at these lines for a couple of hours. Regards Srinath |

Hi Srinath
the operators +,Hi Srinath
the operators +,-,==, are defined in the class of fvMatrix located at /src/finiteVolume/fvMatrices/fvMatrix/ It seems that fvMatrix contains two main private parts 1: the coeficient matrix "A" of all cells 2: the source field "d" all the operations (div,ddt,etc) on the fields will generate a coeficient matrix or source field depending on the explicit(fvc::...) or implicit(fvm::...) feature. for fvm::div(..) +fvm::ddt(..)+fvm::... all the coeficient matrix for each term added together. for fvc::div(..) +fvc::... all the source field for each term added together for fvm::.. +fvc::.. Just add the field generated from the sencond term into the source field in the first term, leave the coeficient unchanged. the operator == has the lowest PRI, and thus it's operation was carried out in the final step. It's function is substract the term on the right from the term on the left handside. like this (coeficient matrix) - (coeficient matrix) (source field)- (source field) Once all the operations are done, we can solve the equations, because coeficient matrix and source term are assembled into fvmatrix. use global function solve(fvmatrix) or local function fvmatrix.solve() to solve these equations. Hope this help. Su Junwei |

Hi Su Junwei
Thanks, that cHi Su Junwei
Thanks, that clears a lot of things. I have one further doubt, How does the solution U, get changed when it is not passed to solve directly? I thought maybe if U(Vector of unknowns), was also a member of fvMatrix, it could be acted upon by solve. But i don't think this assumption of mine is correct. If we say solve(Ueqn == fvc:: ), Ueqn==fvc:: will generate a temporary variable of type fvMatrix, with a coefficientMatrix, source term and (vector of unknowns?). This vector of unknowns if altered by solve, will go out of scope. So this can't be right? I would appreciate it if you could clear this doubt, because i find it quite magical as to how solve changes the solution, when the solution vector is not passed to it.(At least i don't see how it is passed to it) Regards Srinath |

I'm no C++ expert, but I thinkI'm no C++ expert, but I think this[1] might help.
[1] http://powerlab.fsb.hr/ped/kturbo/Op...apers/Foam.pdf |

Hi Srinath
It should be notHi Srinath
It should be noted that in all the operations of fvm::... the volVectorField U were passed as a reference to the private member GeometricField<type,>& psi_; in the fvMatrix So, Once equation solved the psi_ will be updated, i.e. the field U will be changed. Of course,following code may also happen solve(fvm::div(U1) + fvm::ddt(U2)) U1 and U2 are different fields To avoid this, in all the operator function about fvmatrix, global function checkMethod() is activated ensuring &U1.psi_==&U2.psi_, i.e. U1 and U2 are the same variable Best,Su Junwei |

Thanks su junwei and Srinath
Thanks su junwei and Srinath
Now, the whole formulation doesn't seem so magical. It is in fact extremely elegant. Do the solvers(specifically the compressible/incompressible ones) have documentation on the formulation employed? I was looking at icofoam.C, and the following lines confused me. Could you tell me what they do? volScalarField rUA = 1.0/UEqn.A(); U = rUA*UEqn.H(); phi = (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, U, phi); adjustPhi(phi, U, p); Are we inverting the matrix formed by implicit(fvm) discretizatin in the first step? If so why is it of type volScalarField? If we are solving for U once more, in step 2, why not just say solve(Ueqn == ...) like in the previous step? The adjustPhi function seems to be used a lot in various solvers. Is it a fluxlimiter? Thanks srinath |

Hi srinath
This concerns wiHi srinath
This concerns with the PISO loop, please refer phd thesis of Hrvoje Jasak Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows |

I think I can help with a few I think I can help with a few of your queries:
UEqn.A(); is the matrix diagonal only. UEqn.H() is the right hand (constant + off-diagonal) side of the matrix. U = rUA*UEqn.H(); does not solve for U, it produces an intermediate variable U* (stored in U to save space) that excludes the effect of pressure gradient. The flux produced by this intermediate variable: phi* = (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, U, phi); is used to calculate the pressure via the continuity equation. adjustPhi is used to adjust the balance of fluxes in cases where no pressure boundary is present. See src/finiteVolume/lnInclude/adjustPhi.C |

Thanks Eugene
I looked at HThanks Eugene
I looked at Hrv's paper foam.pdf, his thesis and Issa's original paper on PISO. I understand and agree with most of what you said, However i have the following doubts.. In Issa's paper, the predictor step does include the effect of pressure gradient(from the previous timestep in the predictor). That is how, it seems to be implemented in icoFoam, as Ueqn is set, before entering the PISO section. What does the term fvc::ddtPhiCorr... do? Also are there references to the compressible flow solvers. I am reading Turkels review paper on densityBasedsolvers, but an exact reference would be helpful. Thanks Srinath |

Hello
My previous post is Hello
My previous post is not entirely correct. I found this good writeup on the net for PISO as implemented in OpenFoam. www.tfd.chalmers.se/~hani/kurser/OF_phD_2007/rhiechow.pdf Things are quite clear after reading this. Only queries are What does peqn.flux() do? Does this get grad(p)? What does fvc::ddtPhiCorr achieve?--Is this to correct mass errors due to interpolating cell centred values onto the face? Thanks Srinath |

deferred-correction=>What does fvc::ddtPhiCorr achieve?--Is this to correct mass errors due to interpolating cell centred values onto the face?
fvc::ddtPhiCor seems like a deferred-correction approach, see Peric's book (3rd rev.), Page122 |

All times are GMT -4. The time now is 03:26. |