fvc::laplacian(rAUf, p_rgh) versus fvm::laplacian(rAUf, p_rgh)

danny123 · July 24, 2015, 11:50

Dear all,

I have a problem to understand the way how OpenFoam dicretizes the pressure equation. There is already a lot of information on this forum, but always on the explicit formulation. E.g.

Code:

fvc::laplacian(rAUf, p_rgh)

is identical to:

Code:

fvc::surfaceIntegrate(rAUf*fvc::snGrad(p_rgh)*mag(mesh.Sf()))

I verified this and it is correct.

Also

Code:

fvc::surfaceIntegrate(rAUf*fvc::snGrad(p_rgh)*mag(mesh.Sf())-phiHbyA)

is identical to

Code:

fvc::laplacian(rAUf, p_rgh) - fvc::div(phiHbyA)

. So the sum rule of the divergence does apply indeed.

But if I use

fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
);

the equation is decitized into a linear matrix equation type

Code:

S - A * p_rgh = rA

whereas rA is the residual that is iterated to be small (A is the matrix). In my case I get an average number per cell of 10^-20, which is really small. However, if compare this by the explicit divergence of the calulated phi field, I get an error of around 10^-8.

Is there any explanation to this?

Regards,

Daniel

danny123 · July 28, 2015, 04:54

I think I have an idea where this may come from. OpenFoam uses a feature called "non-othogonal correction". It supposely takes care of this issue. In order to see, if it is doing what I need, I can add non-orthogonal corrector loops and hope this does change something.

It would be nice to have some measurement of effictiveness of non-othogonal correction. I would need the residuals of the matrix iteration and, divide them by the expicit error (div (phi) in this case and some the absolute of these local errors up, eventually divide it by the number of cells. This would give me some measure if the non-orthogonal correction does something.

Is there some method to get the matrix residuals. I do not see any of such object in the solverperformance declaration.

Regards,

Daniel

danny123 · July 29, 2015, 05:19

I checked the non orthogonal correction. Even this swithed off, I get the same inconsistency. So, I am kind of stuck here. Is there anybody who has an explanation? some higher order terms, which are neglegted in fvm discretization procedure?

Curreously: if I set use the explicit error and add it to a 2nd fvm discretization, the deviation remains.

I am kind of stuck here. Is there somebody who has an idea?

Regards,

Daniel

riddim · July 30, 2015, 11:14

Hey danny,

I'm very new here and as you will have recognized I myself am also a little bit lost

To be honest I don't really see what your question is ...

As I understand it, you want to know what the difference between the explicit and the implicit class is.
I am not sure if I got it right myself ... I am asking myself the same thing for some time now

but as the "Programmers Guide" says:
"vm and fvc contain static functions, representing differential
operators, e.g. ∇ 2 , ∇ • and ∂/∂t, that discretise geometricField < Type > s. The purpose of defining these functions within two classes, fvm and fvc, rather than one, is to distinguish:
• functions of fvm that calculate implicit derivatives of and return an fvMatrix < Type >
• some functions of fvc that calculate explicit derivatives and other explicit calcula-
tions, returning a geometricField < Type > ."
(P-36)

So I think the main difference is that derivatives are calculated in different ways. (see http://www.sosmath.com/calculus/diff/der05/der05.html )
Since the path to your solution is a different one this may (or will) cause a small difference in the numerical solution.

Though there is a small difference both solutions should be pretty similar and about equally "correct" (and roughly equally wrong)

I don't know much about the "non-orthogonal correction" loops but they seem to be important in some cases ... I didn't really get what they do yet ... (P-51)

Sorry - I know I'm not very helpful here ... but as I mentioned I'm not even sure what exactly your problem is ... and I don't know how you can sum up exact differences ...

danny123 · July 30, 2015, 12:54

It is getting even more confusing since I added a comment line in the gausslapacianscheme and it is not to be found in the output. I think I have to sleep this over.

Since CFD uses linear discretization exclusively, You can always represent a differential equation as A * x = b, A being a matrix, x the value vector and b the source term. So, there should not be no difference to the explicit calculation unless there is some multiplication.

E.g. in the convective term there is a quadratic term (rho phi U), but this is treated is a linear relationship too, since only U is treated implicetly.

By the way, fvm::lapalcian(A) does exists already in the namespace. It is treated the same way I told you, but consideres non-othogonal correction. You can switch it off though.

Regards,

Daniel

danny123 · July 31, 2015, 06:58

I put info statements into laplacianScheme.C. It compiles, but there is no output. Is there any reason to this? Or do I look into the wrong file? To test, I simply did this:

Code:

//    if (fv::debug)
    {
        Info<< "laplacianScheme<Type, GType>::New(const fvMesh&, Istream&) : "
               "constructing laplacianScheme<Type, GType>"
            << endl;
    }

This should basically cancel out the debug switch.

In the application solver, I call

Code:

        fvScalarMatrix pcorrEqn
        (
            fvm::laplacian(rAUf, pcorr) == fvc::div(phiHbyA)
        );

This works and I get an ouput, but the info statement above is not plotted.

Is there anybody who has an idea? The same happens, if I add info statements into gaussLaplacianScheme.C

Regards,

Daniel

danny123 · August 4, 2015, 06:42

Hello,

Now it becomes more and more tricky. I have 2 implicit laplacians in my code, one explicit laplacian. So far, so simple.

I put info statement into the code. It seems that

Code:

template<class Type, class GType>
tmp<fvMatrix<Type> >
laplacianScheme<Type, GType>::fvmLaplacian

is called located in laplacianScheme.C. It gets called only once. The explicit call within the same file does not get called

.

When the individual laplacians are built,

Code:

template<class Type, class GType>
tmp<fvMatrix<Type> >
gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected

gets called, located in gaussLaplacianScheme.C. It does not matter if I specify "Gauss linear corrected" or "Gauss linear uncorrected"

.
The explicit laplacian call also does not get called by the pcorrEq.flux() command.

Is there anybody who knows how the laplacian call logic works? It would be nice to have that of the snGrad logic too, since the explicit laplacian calls snGrade (even though maybe not, since the fvc laplacian does not get called).

Regards,

Daniel

wyldckat · August 20, 2015, 17:31

Greetings to all!

@Daniel: Sorry, but only today did I manage to look into this and I have to say that I'm only able to answer questions regarding code itself, I don't know enough how OpenFOAM discretizes the equations

Quote:

Originally Posted by danny123

I put info statements into laplacianScheme.C. It compiles, but there is no output. Is there any reason to this?

I don't know if you ever found this thread: http://www.cfd-online.com/Forums/ope...epository.html - but the explanation is simple enough: these classes you're modifying are template classes; what happens is that when they are included by the solver, these template classes are re-implemented into the solver itself and the ones on the binary build of "finiteVolume" library are not the ones used by the solver itself. This is why that you do not see any messages when you modify the code and build only the library "finiteVolume". You will have to build the solver as well, for the changes to propagate.

Quote:

Originally Posted by danny123

Is there some method to get the matrix residuals. I do not see any of such object in the solverperformance declaration.

You will have to hack directly the class for the matrix solver being used for an equation. For example, internal residuals for the GAMG solver are in this file: "src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverSolve.C"

Quote:

Originally Posted by danny123

Now it becomes more and more tricky. I have 2 implicit laplacians in my code, one explicit laplacian. So far, so simple.

[...]gets called, located in gaussLaplacianScheme.C. It does not matter if I specify "Gauss linear corrected" or "Gauss linear uncorrected"

.
The explicit laplacian call also does not get called by the pcorrEq.flux() command.

Is there anybody who knows how the laplacian call logic works? It would be nice to have that of the snGrad logic too, since the explicit laplacian calls snGrade (even though maybe not, since the fvc laplacian does not get called).

From what I can see in the source code, the uncorrected version needs to be called even when correction is applied, because the correction is added afterwards: https://github.com/OpenFOAM/OpenFOAM...nScheme.C#L184

Code:

    if (this->tsnGradScheme_().corrected())
    {
        tfaceFluxCorrection() +=
            SfGammaSn*this->tsnGradScheme_().correction(vf);
    }

As for the "pEqn.flux()", it's done in this method: https://github.com/OpenFOAM/OpenFOAM...vMatrix.C#L883 - it only manipulates/sets-up a matrix to be delivered, it does not do any discretization, i.e. it bases itself on whatever the result was from the call to "solve".

If you're trying to find your bearings in the source code, hopefully this description will put you in the right direction.

Best regards,
Bruno

danny123 · August 24, 2015, 12:39

Thanks Bruno,

This is very helpful. Only a little question, rebuilding the solver, you mean the application solver (interFoam bla. bla.) or the matrix solver, meaning rebuilding within src/OpenFOAM?

Regards,

Daniel

wyldckat · August 24, 2015, 12:45

Hi Daniel,

I think you're referring to this paragraph:

Quote:

Originally Posted by wyldckat

I don't know if you ever found this thread: http://www.cfd-online.com/Forums/ope...epository.html - but the explanation is simple enough: these classes you're modifying are template classes; what happens is that when they are included by the solver, these template classes are re-implemented into the solver itself and the ones on the binary build of "finiteVolume" library are not the ones used by the solver itself. This is why that you do not see any messages when you modify the code and build only the library "finiteVolume". You will have to build the solver as well, for the changes to propagate.

Yes, sorry about not having be clearer, but yes, I'm referring to the solver application, such as interFoam, simpleFoam, etc...

Best regards,
Bruno

danny123 · August 25, 2015, 13:38

Hello Bruno,

Well, I did what you said, but there is no difference. My application solver is located in run/myworks, it is home build one based on interDymFoam. Maybe this makes a difference... It should not though. I also tried wclean upfront of wmake, nothing different though.

Reading the code, I thought the laplacian is structured the following way:

fvm:

the header files are located in laplacianScheme.H, which then calls laplacianScheme.C for the definitions. There are a number of different definitions, basically depending whether Gamma is present or if Gamma is not (otherwise Gamma is defined as unity), and whether Gamma is taken at cell center, so a face interpolation is needed. They all converge into a face interpolated Gamma laplacian. This then should call (somehow?)

Code:

template<class Type, class GType>
tmp<fvMatrix<Type> >
gaussLaplacianScheme<Type, GType>::fvmLaplacian

This template class calls the uncorrected laplacian in gaussLaplacianScheme.C for the implicit (matrix) part and gammaSnGradCorr for the explicit correction, which is added to the source term.

After solving the flux function is called directly on the matrix equation, as you correctly pointed out. This is where I suspect a mistake in the code, but regardless, I need to be sure that the way how it is described here is correct. Adding a comment line in either fvmLaplacian or gammaSnGradCorr does not show up, in fvmLaplacianUncorrected it does. They all are template classes.

fvc:

very much the same way, only there is no matrix built, but simply the div operator is used and snGrade to determine the face fluxes. Curiously, this operator does not get called at all. I obviously put an fvc::Laplacian into the application solver code.

If you could help what else I could do, I would appreciate. I will try to compile directly in the native interFoam solver, just to see if the problem stems from there.

Regards,

Daniel

danny123 · August 26, 2015, 06:57

So, now I compiled interFoam new and ran the application. I attach the 2 modified laplacian files, where I put comment lines in.

The output is:

HTML Code:

...  fvmLaplacian is here 3... fvmLaplacian is here 1 fvmLaplacian is here 3 fvmLaplacian is here 3... fvmLaplacian is here 1 fvmLaplacian is here 3 fvmLaplacian is here 3... etc.

This is identical to my custom application solver. In interDymFoam (as in interFoam), there are 3 laplacian fvm calls, the 1st in pcorrEqn, the 2nd in divDevRhoReff located in kEpsilon.C and the 3rd one is in p_rghEq. The 2nd one uses a cell centred gamma, this is why you see "fvmLaplacian is here 1" upfront of "fvmLaplacian is here 3". For the 2 others face interpolation of rau to rauf (=gamma) is done before the call of laplacian.

In my opinion, looking into the code, for each "fvmLaplacian is here 3" mentioned above, the output should be:

HTML Code:

 fvmLaplacian is here 5 tgammaSnGradCorr is here fvmLaplacian is here 3

For each fvc laplacian, there should be:

HTML Code:

 fvcLaplacian is here 4

preceded by

HTML Code:

 fvcLaplacian is here 2

if face interpolation of gamma is required. I checked this too, "fvcLaplacian is here 2" shows up as expected, but no " fvcLaplacian is here 4" (There is no fvcLaplacian in interFoam and also no cell one for cell centered gamma, so I added the code into the application solver).

Looking back into the code, it is striking that only fvmLaplacianUncorrected shows up within gaussLaplacianScheme. C. If one looks into the header file, this is the only static template (code in gaussLaplacianScheme.H). Moreover, at the bottom there is following code

Code:

// Use macros to emulate partial-specialisation of the the Laplacian functions
// for scalar diffusivity gamma

#define defineFvmLaplacianScalarGamma(Type)                                 \
                                                                            \
template<>                                                                  \
tmp<fvMatrix<Type> > gaussLaplacianScheme<Type, scalar>::fvmLaplacian       \
(                                                                           \
    const GeometricField<scalar, fvsPatchField, surfaceMesh>&,              \
    const GeometricField<Type, fvPatchField, volMesh>&                      \
);                                                                          \
                                                                            \
template<>                                                                  \
tmp<GeometricField<Type, fvPatchField, volMesh> >                           \
gaussLaplacianScheme<Type, scalar>::fvcLaplacian                            \
(                                                                           \
    const GeometricField<scalar, fvsPatchField, surfaceMesh>&,              \
    const GeometricField<Type, fvPatchField, volMesh>&                      \
);

So, are there some macros, that get called and prevent the source code being active? Where are they, how do I get around them?

Regards,

Daniel

danny123 · August 27, 2015, 10:03

Hello Bruno and you other alls,

I found the solution, it is a tricky one, not that obvious and, honestly, I do not understand the reason why it is implemented this way. Maybe, Bruno, can you shed some light on this?

The fvmlaplacian as well as the fvclaplacian actually exists twice in the source code. The alternative code is located in gaussLaplacianSchemes.C. It is written after a #define statement, so everything is in pink and backslashes at the end of lines. It can be edited as any other code (and compiles too).

The declaration line the the fvm laplacian used is

Code:

template<>                                                                   \
Foam::tmp<Foam::fvMatrix<Foam::Type> >                                       \
Foam::fv::gaussLaplacianScheme<Foam::Type, Foam::scalar>::fvmLaplacian

The one not used is:

Code:

template<class Type, class GType>
tmp<fvMatrix<Type> >
gaussLaplacianScheme<Type, GType>::fvmLaplacian

So, I assume this difference determines, which one is selected (not sure).

I now can figure out the difference where the difference in fvm and fvc comes from and will report, if succesful.

Regards,

Daniel

wyldckat · August 30, 2015, 19:21

Hi Daniel,

The file "gaussLaplacianSchemes.C" has dedicated specializations of the main template. This is a C++ feature, where we can define a main template structure and then have a specific implementation have a slightly different thing done in some situations.

In this example, we have the main template class "gaussLaplacianScheme", with this template structure:

Code:

template<class Type, class GType> class gaussLaplacianScheme

where later in the file "gaussLaplacianScheme.H" is has these three methods:

Code:

        tmp<GeometricField<Type, fvPatchField, volMesh> > fvcLaplacian
        (
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        tmp<fvMatrix<Type> > fvmLaplacian
        (
            const GeometricField<GType, fvsPatchField, surfaceMesh>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

        tmp<GeometricField<Type, fvPatchField, volMesh> > fvcLaplacian
        (
            const GeometricField<GType, fvsPatchField, surfaceMesh>&,
            const GeometricField<Type, fvPatchField, volMesh>&
        );

It's hard to spot, but the names "Type" and "GType" are placeholders for any data structures, such as "scalar" and "vector".
Then in "gaussLaplacianSchemes.C" we have dedicated implementations for when "GType" is "scalar":

Code:

template<>                                                                   \
Foam::tmp<Foam::fvMatrix<Foam::Type> >                                       \
Foam::fv::gaussLaplacianScheme<Foam::Type, Foam::scalar>::fvmLaplacian       \
(                                                                            \
    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,         \
    const GeometricField<Type, fvPatchField, volMesh>& vf                    \
)                                                                            \

...

template<>                                                                   \
Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh> >\
Foam::fv::gaussLaplacianScheme<Foam::Type, Foam::scalar>::fvcLaplacian       \
(                                                                            \
    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,         \
    const GeometricField<Type, fvPatchField, volMesh>& vf                    \
)                                                                            \

These two methods replace the generic implementation that is defined in "gaussLaplacianScheme.C".

Comparing the source code, the idea is that the "gamma" field (generic name) that this specialization uses is something that is previously calculated... probably the "phi" field. While the generic template class has to always calculate the "gamma" field.

Best regards,
Bruno

danny123 · September 9, 2015, 13:01

Hello Bruno,

This definetly helped me to understand the code. I found the missmatch, but I do not know the motivation of it.

Basically the fvm::laplacian does a balance over the cell faces of the product of the gradient of some surface field and another surface field. That other surface field is called gamma, sometimes called diffusity. The gradient of the 1st surface field is derived from the cell values of the field.

In order to do a correct balance, you need some kind of flow over a face. Then you sum them up over each faces over a cell (and apply boundary conditions). The result of that sum has the same unit as the individual face flows. For pcorr.Eq in interFoam this is m3/s. You now may divide each balance value by some norming factor. But this would mean that you cannot add or substract additional terms. So, my guess is that the implementation is correct.

What is not correct is the explicit implementation. You can see this in fvcSurfaceIntegrate.C:

Code:

template<class Type>
void surfaceIntegrate
(
    Field<Type>& ivf,
    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
)
{
    const fvMesh& mesh = ssf.mesh();

    const labelUList& owner = mesh.owner();
    const labelUList& neighbour = mesh.neighbour();

    const Field<Type>& issf = ssf;

    forAll(owner, facei)
    {
        ivf[owner[facei]] += issf[facei];
        ivf[neighbour[facei]] -= issf[facei];
    }

    forAll(mesh.boundary(), patchi)
    {
        const labelUList& pFaceCells =
            mesh.boundary()[patchi].faceCells();

        const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];

        forAll(mesh.boundary()[patchi], facei)
        {
            ivf[pFaceCells[facei]] += pssf[facei];
        }
    }

    ivf /= mesh.Vsc();
}

The explicit Lapacian uses the divergence operator, which is identical to surfaceIntegrate. As one can see, the result is normed on the cell volume mesh.Vsc(). So, this means any explicit fvc:div or fvc:laplacian (which contains a div) is potentially wrong. We can correct the code by multiplying mesh.Vsc() or alter the code stated above.

I tried it out and indeed fvc::div(phi) has the unit 1/s.

I have not checked all the other implementations (fvm:div). Do you have an idea why the fvc::surfaceIntegrate is implemented such?

Regards,

Daniel

wyldckat · September 9, 2015, 17:36

Hi Daniel,

I won't be able to look into this in more detail any time soon, but I do suggest you check Hrvoje Jasak's PhD Thesis: http://powerlab.fsb.hr/ped/kturbo/Op...jeJasakPhD.pdf - hopefully this and other details are explained there.

Best regards,
Bruno

danny123 · September 18, 2015, 06:24

Hello,

For those who are interested in this topic, this is what I figured out:

First, and most importantly, the implementation is fundametally correct. On can debate on some minor points, see below.

When you do a balance on a control volume, or as it is named in OpenFoam, a cell, you have to establish an equation that has the dimension flow. Basically, the property of the cells changes because something flows inside the cell or outside therefrom. This corresponds to the integral form of the Narvier Stakes Equation. As an example, take:

$\nabla \cdot U = 0$

The integral form is:

$\int_{V} (\nabla \cdot U) dV \approx \sum_{faces} U{_f} S{_f} = 0$

The unit of this equation is m^3/s. U^f is the volumetric face flux and the product of face surface and flux is the volumetric face flow phi.

fvm operators and fvc both need to have the same dimension, in this case m^3/s. What is confusing is, that the for fvc::div, you get the dimension 1/s, since the result of the balance is divided by the cell volume. I can only guess why this is done so, maybe it is a convention or whatever (since nabla dot U also is of dimension 1/s). Now, if one adds a fvc::div term to the balance equation, the result is automatically remultiplied by the cell volume. So, you divide by the cell volume and remultiply and get back to the same result. If you look at the code in your application solver, the same code line fvc::div within a balance equation (matrix equation) and as regular calculation have different results. This is what did confuse me.

Some other points:

- I found small deviations because of the numerical precision since multiplying and redividing by a number does not lead to identical results. These errors are very small, but can be annoying.
- the non-orthogonal correction has been removed from the snGrad and therefore does not apply for the explicit laplacian. This term should be present in the snGradScheme.C.
- if non orthogonal correction is to be used, you need to add "Gauss linear uncorrected" for the laplacianSchemes in fvSchemes.C to apply the corretor for the fvm::div operator as well for the flux correction, and "corrected" for the snGradSchemes to apply for the fvc::div operator. Note that there is a fvc::div embedded in the fvm::div operator (the non orthogonal correction), so you get some kind of "correction of the correction". Note also that the correction will apply for all "stand alone" snGrad operators. I do not know the thinking behind this is since what needs to be corrected is the gradient and not the divergence operator. So, if you want gradient correction, lapalcians always need to be corrected too since the gradient is an integral part of the laplacian.

Hope these explanations help. Any comment is welcome.

Regards,

Daniel

wyldckat · September 19, 2015, 13:20

Hi Daniel,

Many thanks for the detailed explanation!
About your questions regarding how the correctors were implemented, have you read Hrvoje Jasak's PhD thesis "Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows"?
It has a lot of these kinds of details explained in there. I've only managed to do a very quick search right now... and I ended up getting lost

Nonetheless, subsections 1.2 and 3.5 seem to be relevant to your questions.

edit: Helpful thesis links here: http://openfoamwiki.net/index.php/IcoFoam

Best regards,
Bruno

danny123 · September 21, 2015, 05:44

Hello Bruno,

Yes, I read this thesis, at least the relevant parts that you mention. They are very helpful in understanding the code. Unfortunatly, the code that is currently implemented is slightly different than the original "over relaxed" method. You can switch to other non orthogonal correction methods by changing the comment "/" sign location in surfaceInterpolation.C.

There is one relevant equation that hit me, which is eq. 3.138. The author forgot to put the integral sign over volume on the left hand side of the equation, So, left hand side you get unit 1/s and right hand side m^3/s...

It is also important to understand why the balance is not done on a reduced by cell volume basis. The matrix in pEq would not be symmetric anymore, increasing computional cost for solving that matrix equation.

This being said, I still not understand why there is this double coding of the laplacian operator. The code in gaussLaplacianSchemes.C covers both scalar fields as vector scalar fields.

Regards,

Daniel

wyldckat · September 21, 2015, 14:49

Hi Daniel,

Quote:

Originally Posted by danny123

There is one relevant equation that hit me, which is eq. 3.138. The author forgot to put the integral sign over volume on the left hand side of the equation, So, left hand side you get unit 1/s and right hand side m^3/s...

Well, a 394 page document has a very high probability of having typos

If we compare with the definition in equation 3.15, the missing term on the left side is "Vp", which is the cell volume, which would fix the problem with the units. Since it's equal to zero either way, it's possible that this was the reason the typo wasn't spotted sooner.

Quote:

Originally Posted by danny123

It is also important to understand why the balance is not done on a reduced by cell volume basis. The matrix in pEq would not be symmetric anymore, increasing computional cost for solving that matrix equation.

Sorry, you lost me here

What do you mean by the following expression?

Quote:

why the balance is not done on a reduced by cell volume basis.

Quote:

Originally Posted by danny123

This being said, I still not understand why there is this double coding of the laplacian operator. The code in gaussLaplacianSchemes.C covers both scalar fields as vector scalar fields.

Are you referring to the second version that uses the "gamma" variable? Or am I missing something?

Best regards,
Bruno

July 24, 2015, 11:50	fvc::laplacian(rAUf, p_rgh) versus fvm::laplacian(rAUf, p_rgh)	#1
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	Dear all, I have a problem to understand the way how OpenFoam dicretizes the pressure equation. There is already a lot of information on this forum, but always on the explicit formulation. E.g. Code: fvc::laplacian(rAUf, p_rgh) is identical to: Code: fvc::surfaceIntegrate(rAUffvc::snGrad(p_rgh)mag(mesh.Sf())) I verified this and it is correct. Also Code: fvc::surfaceIntegrate(rAUffvc::snGrad(p_rgh)mag(mesh.Sf())-phiHbyA) is identical to Code: fvc::laplacian(rAUf, p_rgh) - fvc::div(phiHbyA) . So the sum rule of the divergence does apply indeed. But if I use fvScalarMatrix p_rghEqn ( fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA) ); the equation is decitized into a linear matrix equation type Code: S - A * p_rgh = rA whereas rA is the residual that is iterated to be small (A is the matrix). In my case I get an average number per cell of 10^-20, which is really small. However, if compare this by the explicit divergence of the calulated phi field, I get an error of around 10^-8. Is there any explanation to this? Regards, Daniel

July 31, 2015, 06:58		#6
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	I put info statements into laplacianScheme.C. It compiles, but there is no output. Is there any reason to this? Or do I look into the wrong file? To test, I simply did this: Code: // if (fv::debug) { Info<< "laplacianScheme<Type, GType>::New(const fvMesh&, Istream&) : " "constructing laplacianScheme<Type, GType>" << endl; } This should basically cancel out the debug switch. In the application solver, I call Code: fvScalarMatrix pcorrEqn ( fvm::laplacian(rAUf, pcorr) == fvc::div(phiHbyA) ); This works and I get an ouput, but the info statement above is not plotted. Is there anybody who has an idea? The same happens, if I add info statements into gaussLaplacianScheme.C Regards, Daniel

August 4, 2015, 06:42		#7
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	Hello, Now it becomes more and more tricky. I have 2 implicit laplacians in my code, one explicit laplacian. So far, so simple. I put info statement into the code. It seems that Code: template<class Type, class GType> tmp<fvMatrix<Type> > laplacianScheme<Type, GType>::fvmLaplacian is called located in laplacianScheme.C. It gets called only once. The explicit call within the same file does not get called . When the individual laplacians are built, Code: template<class Type, class GType> tmp<fvMatrix<Type> > gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected gets called, located in gaussLaplacianScheme.C. It does not matter if I specify "Gauss linear corrected" or "Gauss linear uncorrected" . The explicit laplacian call also does not get called by the pcorrEq.flux() command. Is there anybody who knows how the laplacian call logic works? It would be nice to have that of the snGrad logic too, since the explicit laplacian calls snGrade (even though maybe not, since the fvc laplacian does not get called). Regards, Daniel

August 25, 2015, 13:38		#11
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	Hello Bruno, Well, I did what you said, but there is no difference. My application solver is located in run/myworks, it is home build one based on interDymFoam. Maybe this makes a difference... It should not though. I also tried wclean upfront of wmake, nothing different though. Reading the code, I thought the laplacian is structured the following way: fvm: the header files are located in laplacianScheme.H, which then calls laplacianScheme.C for the definitions. There are a number of different definitions, basically depending whether Gamma is present or if Gamma is not (otherwise Gamma is defined as unity), and whether Gamma is taken at cell center, so a face interpolation is needed. They all converge into a face interpolated Gamma laplacian. This then should call (somehow?) Code: template<class Type, class GType> tmp<fvMatrix<Type> > gaussLaplacianScheme<Type, GType>::fvmLaplacian This template class calls the uncorrected laplacian in gaussLaplacianScheme.C for the implicit (matrix) part and gammaSnGradCorr for the explicit correction, which is added to the source term. After solving the flux function is called directly on the matrix equation, as you correctly pointed out. This is where I suspect a mistake in the code, but regardless, I need to be sure that the way how it is described here is correct. Adding a comment line in either fvmLaplacian or gammaSnGradCorr does not show up, in fvmLaplacianUncorrected it does. They all are template classes. fvc: very much the same way, only there is no matrix built, but simply the div operator is used and snGrade to determine the face fluxes. Curiously, this operator does not get called at all. I obviously put an fvc::Laplacian into the application solver code. If you could help what else I could do, I would appreciate. I will try to compile directly in the native interFoam solver, just to see if the problem stems from there. Regards, Daniel

August 27, 2015, 10:03		#13
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	Hello Bruno and you other alls, I found the solution, it is a tricky one, not that obvious and, honestly, I do not understand the reason why it is implemented this way. Maybe, Bruno, can you shed some light on this? The fvmlaplacian as well as the fvclaplacian actually exists twice in the source code. The alternative code is located in gaussLaplacianSchemes.C. It is written after a #define statement, so everything is in pink and backslashes at the end of lines. It can be edited as any other code (and compiles too). The declaration line the the fvm laplacian used is Code: template<> \ Foam::tmp<Foam::fvMatrix<Foam::Type> > \ Foam::fv::gaussLaplacianScheme<Foam::Type, Foam::scalar>::fvmLaplacian The one not used is: Code: template<class Type, class GType> tmp<fvMatrix<Type> > gaussLaplacianScheme<Type, GType>::fvmLaplacian So, I assume this difference determines, which one is selected (not sure). I now can figure out the difference where the difference in fvm and fvc comes from and will report, if succesful. Regards, Daniel

July 28, 2015, 04:54		#2
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	I think I have an idea where this may come from. OpenFoam uses a feature called "non-othogonal correction". It supposely takes care of this issue. In order to see, if it is doing what I need, I can add non-orthogonal corrector loops and hope this does change something. It would be nice to have some measurement of effictiveness of non-othogonal correction. I would need the residuals of the matrix iteration and, divide them by the expicit error (div (phi) in this case and some the absolute of these local errors up, eventually divide it by the number of cells. This would give me some measure if the non-orthogonal correction does something. Is there some method to get the matrix residuals. I do not see any of such object in the solverperformance declaration. Regards, Daniel

July 29, 2015, 05:19		#3
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	I checked the non orthogonal correction. Even this swithed off, I get the same inconsistency. So, I am kind of stuck here. Is there anybody who has an explanation? some higher order terms, which are neglegted in fvm discretization procedure? Curreously: if I set use the explicit error and add it to a 2nd fvm discretization, the deviation remains. I am kind of stuck here. Is there somebody who has an idea? Regards, Daniel

July 30, 2015, 11:14		#4
riddim New Member Jan Join Date: Oct 2014 Posts: 9 Rep Power: 12	Hey danny, I'm very new here and as you will have recognized I myself am also a little bit lost To be honest I don't really see what your question is ... As I understand it, you want to know what the difference between the explicit and the implicit class is. I am not sure if I got it right myself ... I am asking myself the same thing for some time now but as the "Programmers Guide" says: "vm and fvc contain static functions, representing differential operators, e.g. ∇ 2 , ∇ • and ∂/∂t, that discretise geometricField < Type > s. The purpose of defining these functions within two classes, fvm and fvc, rather than one, is to distinguish: • functions of fvm that calculate implicit derivatives of and return an fvMatrix < Type > • some functions of fvc that calculate explicit derivatives and other explicit calcula- tions, returning a geometricField < Type > ." (P-36) So I think the main difference is that derivatives are calculated in different ways. (see http://www.sosmath.com/calculus/diff/der05/der05.html ) Since the path to your solution is a different one this may (or will) cause a small difference in the numerical solution. Though there is a small difference both solutions should be pretty similar and about equally "correct" (and roughly equally wrong) I don't know much about the "non-orthogonal correction" loops but they seem to be important in some cases ... I didn't really get what they do yet ... (P-51) Sorry - I know I'm not very helpful here ... but as I mentioned I'm not even sure what exactly your problem is ... and I don't know how you can sum up exact differences ...

July 30, 2015, 12:54		#5
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	It is getting even more confusing since I added a comment line in the gausslapacianscheme and it is not to be found in the output. I think I have to sleep this over. Since CFD uses linear discretization exclusively, You can always represent a differential equation as A * x = b, A being a matrix, x the value vector and b the source term. So, there should not be no difference to the explicit calculation unless there is some multiplication. E.g. in the convective term there is a quadratic term (rho phi U), but this is treated is a linear relationship too, since only U is treated implicetly. By the way, fvm::lapalcian(A) does exists already in the namespace. It is treated the same way I told you, but consideres non-othogonal correction. You can switch it off though. Regards, Daniel

August 24, 2015, 12:39		#9
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	Thanks Bruno, This is very helpful. Only a little question, rebuilding the solver, you mean the application solver (interFoam bla. bla.) or the matrix solver, meaning rebuilding within src/OpenFOAM? Regards, Daniel

August 30, 2015, 19:21		#14
wyldckat Retired Super Moderator Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 10,981 Blog Entries: 45 Rep Power: 128	Hi Daniel, The file "gaussLaplacianSchemes.C" has dedicated specializations of the main template. This is a C++ feature, where we can define a main template structure and then have a specific implementation have a slightly different thing done in some situations. In this example, we have the main template class "gaussLaplacianScheme", with this template structure: Code: template<class Type, class GType> class gaussLaplacianScheme where later in the file "gaussLaplacianScheme.H" is has these three methods: Code: tmp<GeometricField<Type, fvPatchField, volMesh> > fvcLaplacian ( const GeometricField<Type, fvPatchField, volMesh>& ); tmp<fvMatrix<Type> > fvmLaplacian ( const GeometricField<GType, fvsPatchField, surfaceMesh>&, const GeometricField<Type, fvPatchField, volMesh>& ); tmp<GeometricField<Type, fvPatchField, volMesh> > fvcLaplacian ( const GeometricField<GType, fvsPatchField, surfaceMesh>&, const GeometricField<Type, fvPatchField, volMesh>& ); It's hard to spot, but the names "Type" and "GType" are placeholders for any data structures, such as "scalar" and "vector". Then in "gaussLaplacianSchemes.C" we have dedicated implementations for when "GType" is "scalar": Code: template<> \ Foam::tmp<Foam::fvMatrix<Foam::Type> > \ Foam::fv::gaussLaplacianScheme<Foam::Type, Foam::scalar>::fvmLaplacian \ ( \ const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma, \ const GeometricField<Type, fvPatchField, volMesh>& vf \ ) \ ... template<> \ Foam::tmp<Foam::GeometricField<Foam::Type, Foam::fvPatchField, Foam::volMesh> >\ Foam::fv::gaussLaplacianScheme<Foam::Type, Foam::scalar>::fvcLaplacian \ ( \ const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma, \ const GeometricField<Type, fvPatchField, volMesh>& vf \ ) \ These two methods replace the generic implementation that is defined in "gaussLaplacianScheme.C". Comparing the source code, the idea is that the "gamma" field (generic name) that this specialization uses is something that is previously calculated... probably the "phi" field. While the generic template class has to always calculate the "gamma" field. Best regards, Bruno __________________ OpenFOAM: FAQ \| Getting started Forum: How to get help, to post code/output and forum guide Read this before sending me PM

September 9, 2015, 13:01		#15
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	Hello Bruno, This definetly helped me to understand the code. I found the missmatch, but I do not know the motivation of it. Basically the fvm::laplacian does a balance over the cell faces of the product of the gradient of some surface field and another surface field. That other surface field is called gamma, sometimes called diffusity. The gradient of the 1st surface field is derived from the cell values of the field. In order to do a correct balance, you need some kind of flow over a face. Then you sum them up over each faces over a cell (and apply boundary conditions). The result of that sum has the same unit as the individual face flows. For pcorr.Eq in interFoam this is m3/s. You now may divide each balance value by some norming factor. But this would mean that you cannot add or substract additional terms. So, my guess is that the implementation is correct. What is not correct is the explicit implementation. You can see this in fvcSurfaceIntegrate.C: Code: template<class Type> void surfaceIntegrate ( Field<Type>& ivf, const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf ) { const fvMesh& mesh = ssf.mesh(); const labelUList& owner = mesh.owner(); const labelUList& neighbour = mesh.neighbour(); const Field<Type>& issf = ssf; forAll(owner, facei) { ivf[owner[facei]] += issf[facei]; ivf[neighbour[facei]] -= issf[facei]; } forAll(mesh.boundary(), patchi) { const labelUList& pFaceCells = mesh.boundary()[patchi].faceCells(); const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi]; forAll(mesh.boundary()[patchi], facei) { ivf[pFaceCells[facei]] += pssf[facei]; } } ivf /= mesh.Vsc(); } The explicit Lapacian uses the divergence operator, which is identical to surfaceIntegrate. As one can see, the result is normed on the cell volume mesh.Vsc(). So, this means any explicit fvc:div or fvc:laplacian (which contains a div) is potentially wrong. We can correct the code by multiplying mesh.Vsc() or alter the code stated above. I tried it out and indeed fvc::div(phi) has the unit 1/s. I have not checked all the other implementations (fvm:div). Do you have an idea why the fvc::surfaceIntegrate is implemented such? Regards, Daniel

September 9, 2015, 17:36		#16
wyldckat Retired Super Moderator Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 10,981 Blog Entries: 45 Rep Power: 128	Hi Daniel, I won't be able to look into this in more detail any time soon, but I do suggest you check Hrvoje Jasak's PhD Thesis: http://powerlab.fsb.hr/ped/kturbo/Op...jeJasakPhD.pdf - hopefully this and other details are explained there. Best regards, Bruno

September 18, 2015, 06:24		#17
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	Hello, For those who are interested in this topic, this is what I figured out: First, and most importantly, the implementation is fundametally correct. On can debate on some minor points, see below. When you do a balance on a control volume, or as it is named in OpenFoam, a cell, you have to establish an equation that has the dimension flow. Basically, the property of the cells changes because something flows inside the cell or outside therefrom. This corresponds to the integral form of the Narvier Stakes Equation. As an example, take: $\nabla \cdot U = 0$ The integral form is: $\int_{V} (\nabla \cdot U) dV \approx \sum_{faces} U{_f} S{_f} = 0$ The unit of this equation is m^3/s. U^f is the volumetric face flux and the product of face surface and flux is the volumetric face flow phi. fvm operators and fvc both need to have the same dimension, in this case m^3/s. What is confusing is, that the for fvc::div, you get the dimension 1/s, since the result of the balance is divided by the cell volume. I can only guess why this is done so, maybe it is a convention or whatever (since nabla dot U also is of dimension 1/s). Now, if one adds a fvc::div term to the balance equation, the result is automatically remultiplied by the cell volume. So, you divide by the cell volume and remultiply and get back to the same result. If you look at the code in your application solver, the same code line fvc::div within a balance equation (matrix equation) and as regular calculation have different results. This is what did confuse me. Some other points: - I found small deviations because of the numerical precision since multiplying and redividing by a number does not lead to identical results. These errors are very small, but can be annoying. - the non-orthogonal correction has been removed from the snGrad and therefore does not apply for the explicit laplacian. This term should be present in the snGradScheme.C. - if non orthogonal correction is to be used, you need to add "Gauss linear uncorrected" for the laplacianSchemes in fvSchemes.C to apply the corretor for the fvm::div operator as well for the flux correction, and "corrected" for the snGradSchemes to apply for the fvc::div operator. Note that there is a fvc::div embedded in the fvm::div operator (the non orthogonal correction), so you get some kind of "correction of the correction". Note also that the correction will apply for all "stand alone" snGrad operators. I do not know the thinking behind this is since what needs to be corrected is the gradient and not the divergence operator. So, if you want gradient correction, lapalcians always need to be corrected too since the gradient is an integral part of the laplacian. Hope these explanations help. Any comment is welcome. Regards, Daniel wyldckat likes this.

September 19, 2015, 13:20		#18
wyldckat Retired Super Moderator Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 10,981 Blog Entries: 45 Rep Power: 128	Hi Daniel, Many thanks for the detailed explanation! About your questions regarding how the correctors were implemented, have you read Hrvoje Jasak's PhD thesis "Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows"? It has a lot of these kinds of details explained in there. I've only managed to do a very quick search right now... and I ended up getting lost Nonetheless, subsections 1.2 and 3.5 seem to be relevant to your questions. edit: Helpful thesis links here: http://openfoamwiki.net/index.php/IcoFoam Best regards, Bruno Last edited by wyldckat; September 19, 2015 at 13:22. Reason: see "edit:"

September 21, 2015, 05:44		#19
danny123 Senior Member Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 15	Hello Bruno, Yes, I read this thesis, at least the relevant parts that you mention. They are very helpful in understanding the code. Unfortunatly, the code that is currently implemented is slightly different than the original "over relaxed" method. You can switch to other non orthogonal correction methods by changing the comment "/" sign location in surfaceInterpolation.C. There is one relevant equation that hit me, which is eq. 3.138. The author forgot to put the integral sign over volume on the left hand side of the equation, So, left hand side you get unit 1/s and right hand side m^3/s... It is also important to understand why the balance is not done on a reduced by cell volume basis. The matrix in pEq would not be symmetric anymore, increasing computional cost for solving that matrix equation. This being said, I still not understand why there is this double coding of the laplacian operator. The code in gaussLaplacianSchemes.C covers both scalar fields as vector scalar fields. Regards, Daniel amuzeshi likes this.