
[Sponsors] 
June 2, 2008, 13:50 
Hello
I am trying to write

#1 
Member
srinath
Join Date: Mar 2009
Location: Champaign, USA
Posts: 90
Rep Power: 10 
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 OpenFoam1.4.1/applications/solvers/incompressible OpenFoam1.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 

June 2, 2008, 21:45 
Hi Srinath
the operators +,

#2 
Senior Member

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 

June 3, 2008, 00:02 
Hi Su Junwei
Thanks, that c

#3 
Member
srinath
Join Date: Mar 2009
Location: Champaign, USA
Posts: 90
Rep Power: 10 
Hi 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 

June 3, 2008, 00:18 
I'm no C++ expert, but I think

#4 
Senior Member
Srinath Madhavan (a.k.a pUl)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 703
Rep Power: 14 
I'm no C++ expert, but I think this[1] might help.
[1] http://powerlab.fsb.hr/ped/kturbo/Op...apers/Foam.pdf 

June 3, 2008, 02:55 
Hi Srinath
It should be not

#5 
Senior Member

Hi 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 

June 3, 2008, 23:41 
Thanks su junwei and Srinath

#6 
Member
srinath
Join Date: Mar 2009
Location: Champaign, USA
Posts: 90
Rep Power: 10 
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 

June 4, 2008, 04:01 
Hi srinath
This concerns wi

#7 
Senior Member

Hi 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 

June 4, 2008, 05:29 
I think I can help with a few

#8 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 14 
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 + offdiagonal) 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 

June 5, 2008, 08:22 
Thanks Eugene
I looked at H

#9 
Member
srinath
Join Date: Mar 2009
Location: Champaign, USA
Posts: 90
Rep Power: 10 
Thanks 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 

June 10, 2008, 05:57 
Hello
My previous post is

#10 
Member
srinath
Join Date: Mar 2009
Location: Champaign, USA
Posts: 90
Rep Power: 10 
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 

April 12, 2009, 04:26 
deferredcorrection

#11 
New Member
Jiang Lijun
Join Date: Mar 2009
Posts: 14
Rep Power: 10 
=>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 deferredcorrection approach, see Peric's book (3rd rev.), Page122 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Doubt  New User  FLUENT  2  July 11, 2008 23:50 
doubt  chandu  FLUENT  2  August 23, 2005 05:17 
Doubt  Mukund  FLUENT  4  August 4, 2004 11:28 
another VOF doubt  some1  Main CFD Forum  6  April 16, 2002 06:52 
Doubt on VOF  some1  Main CFD Forum  12  April 6, 2002 14:18 