[srinath]

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

[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

[srinath]

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

[msrinath80]

I'm no C++ expert, but I think this[1] might help.

[1] http://powerlab.fsb.hr/ped/kturbo/Op...apers/Foam.pdf

[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

[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

[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

[eugene]

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

[srinath]

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]

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

[L.J.Jiang]

=>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