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

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

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   July 24, 2015, 11:50
Default fvc::laplacian(rAUf, p_rgh) versus fvm::laplacian(rAUf, p_rgh)
  #1
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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 is offline   Reply With Quote

Old   July 28, 2015, 04:54
Default
  #2
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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 is offline   Reply With Quote

Old   July 29, 2015, 05:19
Default
  #3
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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
danny123 is offline   Reply With Quote

Old   July 30, 2015, 11:14
Default
  #4
New Member
 
Jan
Join Date: Oct 2014
Posts: 8
Rep Power: 10
riddim is on a distinguished road
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 ...
riddim is offline   Reply With Quote

Old   July 30, 2015, 12:54
Default
  #5
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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 is offline   Reply With Quote

Old   July 31, 2015, 06:58
Default
  #6
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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 is offline   Reply With Quote

Old   August 4, 2015, 06:42
Default
  #7
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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
danny123 is offline   Reply With Quote

Old   August 20, 2015, 17:31
Default
  #8
Retired Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,969
Blog Entries: 45
Rep Power: 126
wyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to all
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 View Post
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 View Post
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 View Post
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
__________________
wyldckat is offline   Reply With Quote

Old   August 24, 2015, 12:39
Default
  #9
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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
danny123 is offline   Reply With Quote

Old   August 24, 2015, 12:45
Default
  #10
Retired Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,969
Blog Entries: 45
Rep Power: 126
wyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to all
Hi Daniel,

I think you're referring to this paragraph:
Quote:
Originally Posted by wyldckat View Post
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
wyldckat is offline   Reply With Quote

Old   August 25, 2015, 13:38
Default
  #11
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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 is offline   Reply With Quote

Old   August 26, 2015, 06:57
Default
  #12
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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
Attached Files
File Type: c gaussLaplacianScheme.C (7.1 KB, 3 views)
File Type: c laplacianScheme.C (4.0 KB, 1 views)
danny123 is offline   Reply With Quote

Old   August 27, 2015, 10:03
Default
  #13
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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
danny123 is offline   Reply With Quote

Old   August 30, 2015, 19:21
Default
  #14
Retired Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,969
Blog Entries: 45
Rep Power: 126
wyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to all
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
__________________
wyldckat is offline   Reply With Quote

Old   September 9, 2015, 13:01
Default
  #15
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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
danny123 is offline   Reply With Quote

Old   September 9, 2015, 17:36
Default
  #16
Retired Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,969
Blog Entries: 45
Rep Power: 126
wyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to all
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
wyldckat is offline   Reply With Quote

Old   September 18, 2015, 06:24
Default
  #17
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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.
danny123 is offline   Reply With Quote

Old   September 19, 2015, 13:20
Default
  #18
Retired Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,969
Blog Entries: 45
Rep Power: 126
wyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to all
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:"
wyldckat is offline   Reply With Quote

Old   September 21, 2015, 05:44
Default
  #19
Senior Member
 
Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 13
danny123 is on a distinguished road
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.
danny123 is offline   Reply With Quote

Old   September 21, 2015, 14:49
Default
  #20
Retired Super Moderator
 
Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,969
Blog Entries: 45
Rep Power: 126
wyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to allwyldckat is a name known to all
Hi Daniel,

Quote:
Originally Posted by danny123 View Post
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 View Post
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 View Post
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
wyldckat is offline   Reply With Quote

Reply

Tags
laplacian operator

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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



All times are GMT -4. The time now is 04:25.