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

A doubt on solve

Register Blogs Community New Posts Updated Threads Search

Like Tree2Likes
  • 1 Post By su_junwei
  • 1 Post By eugene

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 2, 2008, 14:50
Default Hello I am trying to write
  #1
Member
 
srinath
Join Date: Mar 2009
Location: Champaign, USA
Posts: 91
Rep Power: 17
srinath is on a distinguished road
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
srinath is offline   Reply With Quote

Old   June 2, 2008, 22:45
Default Hi Srinath the operators +,
  #2
Senior Member
 
su_junwei's Avatar
 
su junwei
Join Date: Mar 2009
Location: Xi'an China
Posts: 151
Rep Power: 20
su_junwei is on a distinguished road
Send a message via MSN to su_junwei
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
HuZexi likes this.
su_junwei is offline   Reply With Quote

Old   June 3, 2008, 01:02
Default Hi Su Junwei Thanks, that c
  #3
Member
 
srinath
Join Date: Mar 2009
Location: Champaign, USA
Posts: 91
Rep Power: 17
srinath is on a distinguished road
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
srinath is offline   Reply With Quote

Old   June 3, 2008, 01:18
Default 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: 21
msrinath80 is on a distinguished road
I'm no C++ expert, but I think this[1] might help.

[1] http://powerlab.fsb.hr/ped/kturbo/Op...apers/Foam.pdf
msrinath80 is offline   Reply With Quote

Old   June 3, 2008, 03:55
Default Hi Srinath It should be not
  #5
Senior Member
 
su_junwei's Avatar
 
su junwei
Join Date: Mar 2009
Location: Xi'an China
Posts: 151
Rep Power: 20
su_junwei is on a distinguished road
Send a message via MSN to su_junwei
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
su_junwei is offline   Reply With Quote

Old   June 4, 2008, 00:41
Default Thanks su junwei and Srinath
  #6
Member
 
srinath
Join Date: Mar 2009
Location: Champaign, USA
Posts: 91
Rep Power: 17
srinath is on a distinguished road
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
srinath is offline   Reply With Quote

Old   June 4, 2008, 05:01
Default Hi srinath This concerns wi
  #7
Senior Member
 
su_junwei's Avatar
 
su junwei
Join Date: Mar 2009
Location: Xi'an China
Posts: 151
Rep Power: 20
su_junwei is on a distinguished road
Send a message via MSN to su_junwei
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
su_junwei is offline   Reply With Quote

Old   June 4, 2008, 06:29
Default I think I can help with a few
  #8
Senior Member
 
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 21
eugene is on a distinguished road
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
openfoam_aero likes this.
eugene is offline   Reply With Quote

Old   June 5, 2008, 09:22
Default Thanks Eugene I looked at H
  #9
Member
 
srinath
Join Date: Mar 2009
Location: Champaign, USA
Posts: 91
Rep Power: 17
srinath is on a distinguished road
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
srinath is offline   Reply With Quote

Old   June 10, 2008, 06:57
Default Hello My previous post is
  #10
Member
 
srinath
Join Date: Mar 2009
Location: Champaign, USA
Posts: 91
Rep Power: 17
srinath is on a distinguished road
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
srinath is offline   Reply With Quote

Old   April 12, 2009, 05:26
Default deferred-correction
  #11
New Member
 
Jiang Lijun
Join Date: Mar 2009
Posts: 14
Rep Power: 17
L.J.Jiang is on a distinguished road
=>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
L.J.Jiang is offline   Reply With Quote

Reply


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 Off
Pingbacks are On
Refbacks are On


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


All times are GMT -4. The time now is 11:22.