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

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 July 24, 2015, 10:50 fvc::laplacian(rAUf, p_rgh) versus fvm::laplacian(rAUf, p_rgh) #1 Senior Member   Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 14 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

 July 28, 2015, 03:54 #2 Senior Member   Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 14 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, 04:19 #3 Senior Member   Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 14 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, 10:14 #4 New Member   Jan Join Date: Oct 2014 Posts: 8 Rep Power: 11 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, 11:54 #5 Senior Member   Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 14 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

 July 31, 2015, 05:58 #6 Senior Member   Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 14 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::New(const fvMesh&, Istream&) : " "constructing laplacianScheme" << 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, 05:42 #7 Senior Member   Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 14 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 tmp > laplacianScheme::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 tmp > gaussLaplacianScheme::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 20, 2015, 16:31
#8
Retired Super Moderator

Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,975
Blog Entries: 45
Rep Power: 128
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() +=
}
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
__________________

 August 24, 2015, 11:39 #9 Senior Member   Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 14 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 24, 2015, 11:45
#10
Retired Super Moderator

Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,975
Blog Entries: 45
Rep Power: 128
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

 August 25, 2015, 12:38 #11 Senior Member   Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 14 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 tmp > gaussLaplacianScheme::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 26, 2015, 05:57
#12
Senior Member

Daniel Witte
Join Date: Nov 2011
Posts: 148
Rep Power: 14
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
 gaussLaplacianScheme.C (7.1 KB, 3 views) laplacianScheme.C (4.0 KB, 1 views)

 August 27, 2015, 09:03 #13 Senior Member   Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 14 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::fv::gaussLaplacianScheme::fvmLaplacian The one not used is: Code: template tmp > gaussLaplacianScheme::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

 August 30, 2015, 18:21 #14 Retired Super Moderator   Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 10,975 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 gaussLaplacianScheme where later in the file "gaussLaplacianScheme.H" is has these three methods: Code:  tmp > fvcLaplacian ( const GeometricField& ); tmp > fvmLaplacian ( const GeometricField&, const GeometricField& ); tmp > fvcLaplacian ( const GeometricField&, const GeometricField& ); 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::fv::gaussLaplacianScheme::fvmLaplacian \ ( \ const GeometricField& gamma, \ const GeometricField& vf \ ) \ ... template<> \ Foam::tmp >\ Foam::fv::gaussLaplacianScheme::fvcLaplacian \ ( \ const GeometricField& gamma, \ const GeometricField& 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, 12:01 #15 Senior Member   Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 14 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 void surfaceIntegrate ( Field& ivf, const GeometricField& ssf ) { const fvMesh& mesh = ssf.mesh(); const labelUList& owner = mesh.owner(); const labelUList& neighbour = mesh.neighbour(); const Field& 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& 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, 16:36 #16 Retired Super Moderator   Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 10,975 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 19, 2015, 12:20 #18 Retired Super Moderator   Bruno Santos Join Date: Mar 2009 Location: Lisbon, Portugal Posts: 10,975 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 12:22. Reason: see "edit:"

 September 21, 2015, 04:44 #19 Senior Member   Daniel Witte Join Date: Nov 2011 Posts: 148 Rep Power: 14 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.

September 21, 2015, 13:49
#20
Retired Super Moderator

Bruno Santos
Join Date: Mar 2009
Location: Lisbon, Portugal
Posts: 10,975
Blog Entries: 45
Rep Power: 128
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

 Tags laplacian operator

 Thread Tools Search this Thread Search this Thread: Advanced Search Display Modes Linear Mode

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

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

 Contact Us - CFD Online - Privacy Statement - Top