[alexandross]

Hello everyone,

I am a new OpenFOAM user, but experienced finite volume programmer. Recently we have conducted a study that shows that the Gauss gradient used in OpenFOAM (and the cell-based Gauss gradient of Fluent) is inconsistent on unstructured meshes, i.e. with grid refinement it converges to a skewness-dependent operator that is different than the actual gradient. Furthermore, we have found, by solving a Poisson equation with an analytic solution, that use of this gradient scheme results in zeroth-order accuracy for the finite volume method as a whole. This means that no matter how fine a grid you use, the accuracy is comparable to that on a very coarse mesh. Some of these tests were conducted with OpenFOAM and Gmsh or Netgen grids. Our study can be found here:

https://hal.archives-ouvertes.fr/hal-01532882v2

(An expanded version that includes also an analysis of the least-squares gradient can be found here:

https://www.researchgate.net/publica...volume_methods

I know that this is probably surprising to many people given the popularity and long history of this scheme, and its reputation for being second-order accurate. However, the analysis is straightforward and the numerical tests are simple and easily repeatable. The problem goes away if one uses the least-squares gradient, which is always at least first-order accurate.

So, I would appreciate your feedback on this issue. Has anyone had a similar experience? Is it documented somewhere else? Since "Gauss linear" is the default scheme in OpenFOAM, and the least-squares scheme is mentioned under "Other schemes that are rarely used" in the cfd-direct users guide, it is puzzling to me how this issue has gone unnoticed so far given the popularity of OpenFOAM.

Of course, the Gauss gradient is not restricted to OpenFOAM but is used more widely in the finite volume realm (I already metioned Fluent for example) but since it is the default choice in OpenFOAM and we used the latter in our numerical tests, I think it is a good idea to post this here. This is a serious issue and it is important that awareness is raised and people switch to "leastSquares" if they use unstructured meshes.

[guin]

Hello Alexandros,
Thank you for this thorough analysis on the accuracy of the Gauss linear discretization method. I think that this is a very important topic, since it may be potentially affecting simulation results of almost everyone here. Would you mind to post your findings at the bug reporting system in order to inform the development team? They can provide you more details of the actual implementation and hopefully improve it if necessary
https://bugs.openfoam.org/

Sent from my HUAWEI RIO-L01 using CFD Online Forum mobile app

[alexandross]

Thanks Rodrigo,

I followed your advice and submitted a bug report (ID 0002596). Although strictly speaking this is not a code bug, informing the developers seems to be the right thing to do.

[arjun]

Originally Posted by alexandross

Thanks Rodrigo,

I followed your advice and submitted a bug report (ID 0002596). Although strictly speaking this is not a code bug, informing the developers seems to be the right thing to do.

This is not a bug. There is no better way where we could have grid independent gradients. The only way could be least square but that too depends on where the neighbour cell centers are located.

I am now trying to get another algorithm for gradient computation in FVUS/wildkatze solver that shall have better gradients and will also be doing similar tests to validate. Still the stability of this new algorithm remains a question. (other than the fact that it is going to be more expensive too). For this reason it would be provided as user model so that user can activate it when he wishes for higher accurate gradients and for the variables he wants to use.

This is certainly not a bug, because they can not fix it and code works as intended.

[alexandross]

Hi Arjun,

Thanks for the reply. As I said, this is not a bug in the code. However, there is a "bug" in the algorithm itself, as the Gauss gradient does not work as intended on unstructured grids. With grid refinement, it converges to a skewness-dependent operator different than the actual gradient. This makes the finite volume method unable to converge with grid refinement to the exact solution of the differential equation solved (at least in the case of the Poisson equation that we tried). Grid convergence is a property expected of any numerical method for solving differential equations, so this is a serious problem.

Originally Posted by arjun

There is no better way where we could have grid independent gradients.

It's hard to agree with this statement when the Gauss gradient is not "grid-independent" whereas many other schemes are; in our paper we show that the least-squares gradient and two alternative gradient schemes based on the Gauss (divergence) theorem are at least first-order accurate on any type of grid, and provide references to other works in the literature that also describe consistent gradient schemes.

I had a look at your code FVUS/wildkatze and it looks like a serious effort. Please let me know of the details once you finalise the new gradient scheme that you are developing.

Best regards,
Alexandros

[arjun]

Originally Posted by alexandross

Hi Arjun,


It's hard to agree with this statement when the Gauss gradient is not "grid-independent" whereas many other schemes are; in our paper we show that the least-squares gradient and two alternative gradient schemes based on the Gauss (divergence) theorem are at least first-order accurate on any type of grid, and provide references to other works in the literature that also describe consistent gradient schemes.

I had a look at your code FVUS/wildkatze and it looks like a serious effort. Please let me know of the details once you finalise the new gradient scheme that you are developing.

Best regards,
Alexandros

The reason is that in case of Green gauss you need face value which is usually computed based on linear interpolation (or say volume weighted ie in starccm+), this means this value is a value that you get somewhere on the line connecting two cell centers. When there is skew the face center does not lie on this line and thus there is now way linear interpolation could give you correct value. So this error that is introduced is skew dependent.

In case of least square the gradient is computed from cell centered and skew does not play direct role.

In FVUS, the algorithm then will be two steps, one find that linear interpolated value on the cell to cell line, then from there using gradient it will find the face center value. So in second step we remove the error introduced due to skew.

I am not the only one to think of it and trying it, in literature there are other people who have done this. The difference is in the method how the gradient at the point on line connected cells is computed. my method is more expensive and does not use cell center gradients here, that previous methods used.

I have once implemented it in starccm and compared against standard method and did find it to be much more accurate, now it will be brought to FVUS/wildkatze for the case people want more accurate solution on expense of 15 to 20% more efforts.

Having said this, how the limiting etc work is also an issue so don't know what will happen to convergence in complicated industrial cases. Need to have the algorithm to see how it works.

Originally Posted by alexandross

Hi Arjun,

I had a look at your code FVUS/wildkatze and it looks like a serious effort.

Yes, it is quite a lot of code and lots of interesting things (Summary here FVUS (Finite Volume Utilities and Solvers) / Wildkatze (solver) Updates
There are two major things that would come next (or least will be tried), one this gradients improvement (efforts to fight skew) and expert drivers to control solver. I wish that user shall not provide any solver setting, solver shall do this part. There is framework for it and briefly tried with flow model, but now efforts will be made to extend this to other models. Once this framework is working as intended , a good amount of cases solver will be able to run and converge the solution. Another 6 months will be very interesting.

[arjun]

Here is one of the attempts for example

http://onlinelibrary.wiley.com/doi/1....1189/abstract

[alexandross]

Thanks for the information Arjun,

Originally Posted by arjun

The reason is that in case of Green gauss you need face value which is usually computed based on linear interpolation (or say volume weighted ie in starccm+), this means this value is a value that you get somewhere on the line connecting two cell centers. When there is skew the face center does not lie on this line and thus there is now way linear interpolation could give you correct value. So this error that is introduced is skew dependent.

Yes, this is precisely the reason why the Gauss gradient fails on unstructured meshes.

Originally Posted by arjun

Here is one of the attempts for example

http://onlinelibrary.wiley.com/doi/1....1189/abstract

I am aware of this publication, it is ref. [18] of our extended paper. In fact this technique dates back at least to 1997, as it is implemented in the CAFFAC code that accompanies the famous book by Ferziger & Peric. The problem is that whereas each "corrector" step improves the accuracy, the order of accuracy remains zero, i.e. as you refine the grid the discretisation error drops only up to a certain limit, but further refinement does not reduce the error any more. Performing more corrector steps lowers this limit, but does not remove it. This is proved in p. 9-10 of our paper (I am referring to the extended version, the one on ResearchGate), and demonstrated in figs. 17 and 19. Theoretically, in the limit of infinite number of corrector steps, the gradient becomes first-order accurate. In eq. (22) of the paper we propose an iterative procedure to achieve this limit without the extra cost of the corrector steps.

There are also other gradient schemes that are derived from the Gauss theorem that you can try. For example, you can calculate the value of the variable at face vertices using second-order accurate interpolation from neighbour cells, and then use linear interpolation between these vertex values to obtain a second-order accurate value for the face centre, which would make the Gauss gradient at least first-order accurate on any grid. I think that this is the "vertex-based" Gauss gradient of Fluent. Also, you could apply the Gauss theorem to the auxiliary cell of Fig. 6 of our paper, which also results in a first-order accurate gradient. Have you thought about these possibilities?

[arjun]

Originally Posted by alexandross

Thanks for the information Arjun,

I am aware of this publication, it is ref. [18] of our extended paper. In fact this technique dates back at least to 1997, as it is implemented in the CAFFAC code that accompanies the famous book by Ferziger & Peric. The problem is that whereas each "corrector" step improves the accuracy, the order of accuracy remains zero, i.e. as you refine the grid the discretisation error drops only up to a certain limit, but further refinement does not reduce the error any more. Performing more corrector steps lowers this limit, but does not remove it. This is proved in p. 9-10 of our paper (I am referring to the extended version, the one on ResearchGate), and demonstrated in figs. 17 and 19. Theoretically, in the limit of infinite number of corrector steps, the gradient becomes first-order accurate. In eq. (22) of the paper we propose an iterative procedure to achieve this limit without the extra cost of the corrector steps.

Yes indeed the correctors don't really help. Thats because the reconstruction gradient still is skew dependent.

Originally Posted by alexandross

There are also other gradient schemes that are derived from the Gauss theorem that you can try. For example, you can calculate the value of the variable at face vertices using second-order accurate interpolation from neighbour cells, and then use linear interpolation between these vertex values to obtain a second-order accurate value for the face centre, which would make the Gauss gradient at least first-order accurate on any grid. I think that this is the "vertex-based" Gauss gradient of Fluent. Also, you could apply the Gauss theorem to the auxiliary cell of Fig. 6 of our paper, which also results in a first-order accurate gradient. Have you thought about these possibilities?

That fig 6 is an interesting approach, but I wonder how do you use this in 3D.

About the node based green gauss fluent mentions that they have comparable accuracy to least square method so interpolating onto nodes really does not make much sense.

In FVUS even if user selects least square method, the final gradient is still green gauss bases. What changes is the reconstruction gradient method.

For FVUS, the aim is to provide more accurate method to the user and thus i am not against least square as recon gradient. The problem is that least square too is control volume dependent (albeit much less than pure green gauss). So there are two approaches that would take place. One of them already was used in iNavier , in this case we add few more points to least sqr gradient stencil. So we look for skew and when skew is high we try to look for another point that could be better position for least sqr gradient.
Second approach is to compute the gradient at the cell to cell line to get face value.

After I wrote last post , i have added the framework in FVUS that provides me bigger stencil in parallel solver. This means I can now in position to soon apply the new gradient algorithm. It shall be soon added. (I will release it with next release in another 2 weeks time).

[alexandross]

Originally Posted by arjun

That fig 6 is an interesting approach, but I wonder how do you use this in 3D.

I haven't, actually! My code is 2D only. It would be more involved, depending on how you have organised your data structures. The auxiliary cell would be of polyhedral shape, whose faces are triangles and its vertices are the centroids of the neighbour cells of the original cell.

Originally Posted by arjun

About the node based green gauss fluent mentions that they have comparable accuracy to least square method so interpolating onto nodes really does not make much sense.

But the least-squares is always at least first-order accurate; do you hope to achieve more than that?

Originally Posted by arjun

In FVUS even if user selects least square method, the final gradient is still green gauss bases. What changes is the reconstruction gradient method.

OK, but the final gradient will still be first-order accurate on unstructured meshes, won't it? So what is its advantage over the least-squares method?

Originally Posted by arjun

For FVUS, the aim is to provide more accurate method to the user and thus i am not against least square as recon gradient. The problem is that least square too is control volume dependent (albeit much less than pure green gauss). So there are two approaches that would take place. One of them already was used in iNavier , in this case we add few more points to least sqr gradient stencil. So we look for skew and when skew is high we try to look for another point that could be better position for least sqr gradient.
Second approach is to compute the gradient at the cell to cell line to get face value.

I guess you mean that the least-squares gradient looses accuracy on skewed meshes, but not order of accuracy (it is at least first-order, independent of skewness). Is there somewhere I can find more information on this?

Originally Posted by arjun

After I wrote last post , i have added the framework in FVUS that provides me bigger stencil in parallel solver. This means I can now in position to soon apply the new gradient algorithm. It shall be soon added. (I will release it with next release in another 2 weeks time).

OK, keep us informed.

[arjun]

Originally Posted by alexandross

OK, but the final gradient will still be first-order accurate on unstructured meshes, won't it? So what is its advantage over the least-squares method?

I guess you mean that the least-squares gradient looses accuracy on skewed meshes, but not order of accuracy (it is at least first-order, independent of skewness). Is there somewhere I can find more information on this?

An example of problem is what happened at thin prism cells. You can see that there green gauss is better than least square gradients.

This is why StarCCM for example provides hybrid gradient computation when you chose least square method. Here in prism cells green gauss method is used even when user selects least square method.

If least square could maintain it's accuracy irrespective of cell shape you would not need to switch to green gauss in prism cells.

I don't have the papers on my tips but this thing is being echoed by papers and demonstrated.

This is the reason in FVUS the final gradient is always green gauss because in the cells the recon gradient fails (even using least square) we fall to green gauss method. When recon gradient is good the green gauss is good too.

[arjun]

i accidentally bumped here so posting it, but it kind of echo what i was writing.

https://bugs.openfoam.org/view.php?id=141

[alexandross]

Originally Posted by arjun

i accidentally bumped here so posting it, but it kind of echo what i was writing.

https://bugs.openfoam.org/view.php?id=141

Thanks. It seems like an implementation issue, as is mentioned in the discussion, because in exact precision arithmetic the least squares gradient is always at least first-order. Using the normal equations to solve the least squares system often results in very ill-conditioned systems. It is better to use a QR decomposition. It also depends on whether single- or double-precision was used. The volume weighting may also play a role; if you have very thin cells whose volume tends to zero, compared to other neighbour cells, this may create a problem. I've only used distance weighting.

[arjun]

Originally Posted by alexandross

Thanks. It seems like an implementation issue, as is mentioned in the discussion, because in exact precision arithmetic the least squares gradient is always at least first-order. Using the normal equations to solve the least squares system often results in very ill-conditioned systems. It is better to use a QR decomposition. It also depends on whether single- or double-precision was used. The volume weighting may also play a role; if you have very thin cells whose volume tends to zero, compared to other neighbour cells, this may create a problem. I've only used distance weighting.

This was closing off comment from Henry

"Yes, this is not a programming bug but a limitation of the current leastSquares implementation which accumulates round-off error on extreme meshes."

So I guess that by implementation issue he means the currently implemented version.

Anyway I have implemented the newer algorithm in FVUS. It will be benchmarked and compared like you did but for quick testing, I had a case with bad cells in it. It is about drag prediction of the car. The second order scheme so far did over predict but the bounded central scheme produced much better results.

With new gradient algorithm second order upwind could catch on with bounded central scheme results.

But as of today, the new gradient algorithm is available. It also opens up for third order solver in FVUS because now extended cell stencil is available in parallel solver too. (This was needed for this gradient).

[arjun]

Okay so I did little bit of benchmarking on tetrahedral mesh case. A simple bend pipe.
Using exact pressure from https://pdfs.semanticscholar.org/434...00ac52c18e.pdf

This is how total error sum looks like

Green Gauss
Weighted least squre
New Green Gauss algorithm.


dpdx dpdy dpdz
4.23014291071706e+03 4.09934782748109e+03 4.06669780949349e+03
1.36558073007596e+03 1.40226594356697e+03 1.32879637216518e+03
6.97935064452681e+02 7.03712372769450e+02 6.88826526239326e+02

You can see this beats least square algorithm by good margin.

[alexandross]

Hello everyone,

Since in a previous post I informed you that I submitted a bug report to bugs.openfoam.org on this issue, following the well advised recommendation of Rodrigo, I think it is fitting to inform you about the outcome of this bug report submission (bug ID 0002596). I apologise for not doing so earlier, but I had a high workload.

The reply said that "Gauss linear gradient is very stable and robust and a good starting point for any simulation and we recommend it for general purpose use", and that "leastSquares can be more accurate than Gauss linear but is a lot less stable and could not be recommended for all simulations and cannot be the default". The report was on for a couple of days, and then its status was changed to "closed", its priority was changed to "normal", and its category was changed from "bug" to "feature".

Therefore, the zeroth-order accuracy issue is not denied, but it is belittled. In my opinion, the zeroth-order accuracy of the Gauss gradient is a major issue but because of its widespread use and associated embarrassment of admitting it openly many people will try to keep it unrevealed, until there is a gradual clandestine transition from the Gauss scheme to a consistent scheme. We had such an experience in trying to publish the arXiv paper I mentioned in my original post. Originally, we submitted it to the Journal of Computational Physics and it was rejected because the reviewers didn't think that the gradient inconsistency would pass on to the finite volume method as a whole (and to give them some credit, it is true that the preliminary version of the paper that we submitted did not include the present section 6, where the gradient is used for the FVM solution of a Poisson equation). Then we added section 6 and submitted it to Computers & Fluids. After more than 6 months of review and repeated emails to the journal, we were finally sent the reviewers' comments. The comments of the first reviewer were reasonable, but the second reviewer asked us to remove all the sections dealing with the Gauss gradient and leave only those on the least-squares one, because the inconsistency issue is "already known" and "redundancy should not rule the research and only originality must be sought" – he/she referred us to the conference paper by Sozer et al (2014) (now ref. [17] of our paper), which we had missed in our literature review. However, since our analysis goes much deeper we refused to remove the sections on the Gauss gradient. We then submitted the revised manuscript, and confidentially warned the editor that the 2nd reviewer is probably someone who has invested in the Gauss gradient and has an interest in keeping the inconsistency issue hidden. The result is that we received a "reject" decision, based on the follow-up comments of only the 2nd reviewer (we were not even sent the follow-up comments of the first reviewer) who now changed his mind completely: he said that the paper of Sozer et al is only a conference paper, not peer-reviewed, and therefore not credible, while he deemed our results also untrustworthy. He stated that it is known from published papers (none of which he cited) that the Gauss gradient is second-order accurate on unstructured meshes.

So, my advice is to conduct your own tests on simple problems where the accuracy can be readily checked instead of just relying on the popularity of the software. It seems that there are still unresolved fundamental issues with finite volume methods on unstructured grids.

[arjun]

Originally Posted by alexandross

The reply said that "Gauss linear gradient is very stable and robust and a good starting point for any simulation and we recommend it for general purpose use", and that "leastSquares can be more accurate than Gauss linear but is a lot less stable and could not be recommended for all simulations and cannot be the default". The report was on for a couple of days, and then its status was changed to "closed", its priority was changed to "normal", and its category was changed from "bug" to "feature".

My 2 cents based on experience so far.

I doubt the stability claim of Green Gauss compared to Least Square. It seems to originates in the near wall prisms, where Green Gauss becomes more accurate and hence more stable.

While it is true for prism, the same thing applies for other part of the mesh where when skew is high least square out performs green gauss and are more stable there.

Over all the behaviour is mixed, one can create a test case where in the interior skew is very bad while keeping the prism of good quality (or have no prisms at all) and demonstrate that least square comes out to be more stable.

(I have fact have such case, i can find it put it only to try with open foam for comparison purposes and it will test the claim too).

Originally Posted by alexandross

Hello everyone,



Therefore, the zeroth-order accuracy issue is not denied, but it is belittled. In my opinion, the zeroth-order accuracy of the Gauss gradient is a major issue but because of its widespread use and associated embarrassment of admitting it openly many people will try to keep it unrevealed, until there is a gradual clandestine transition from the Gauss scheme to a consistent scheme. We had such an experience in trying to publish the arXiv paper I mentioned in my original post. Originally, we submitted it to the Journal of Computational Physics and it was rejected because the reviewers didn't think that the gradient inconsistency would pass on to the finite volume method as a whole (and to give them some credit, it is true that the preliminary version of the paper that we submitted did not include the present section 6, where the gradient is used for the FVM solution of a Poisson equation). Then we added section 6 and submitted it to Computers & Fluids. After more than 6 months of review and repeated emails to the journal, we were finally sent the reviewers' comments. The comments of the first reviewer were reasonable, but the second reviewer asked us to remove all the sections dealing with the Gauss gradient and leave only those on the least-squares one, because the inconsistency issue is "already known" and "redundancy should not rule the research and only originality must be sought" – he/she referred us to the conference paper by Sozer et al (2014) (now ref. [17] of our paper), which we had missed in our literature review. However, since our analysis goes much deeper we refused to remove the sections on the Gauss gradient. We then submitted the revised manuscript, and confidentially warned the editor that the 2nd reviewer is probably someone who has invested in the Gauss gradient and has an interest in keeping the inconsistency issue hidden. The result is that we received a "reject" decision, based on the follow-up comments of only the 2nd reviewer (we were not even sent the follow-up comments of the first reviewer) who now changed his mind completely: he said that the paper of Sozer et al is only a conference paper, not peer-reviewed, and therefore not credible, while he deemed our results also untrustworthy. He stated that it is known from published papers (none of which he cited) that the Gauss gradient is second-order accurate on unstructured meshes.

I agree with you here so next comment is just a side addition.

For these reasons as I mentioned earlier least square is only used for recon in FVUS. The final algorithm for gradient is green gauss with using least square recon gradients.

In the next versions of FVUS, the least square might take second stage and newer Green Gauss will take central role and this will be for the reasons of accuracy, contrary to what is been told there.

[akidess]

You left out that the bug report is to be reopened once you contribute your code. I agree with Henry that this is not a bug but a feature. There is no default setting in OpenFOAM - if you do not specify a gradient scheme in fvSchemes, the code will exit. It is the general philosophy of OpenFOAM to leave such decisions to the user. The user needs to be aware of the limitations of the selected schemes (thank you for your contribution in that regard), and make appropriate choices.

Originally Posted by alexandross

Hello everyone,

Since in a previous post I informed you that I submitted a bug report to bugs.openfoam.org on this issue, following the well advised recommendation of Rodrigo, I think it is fitting to inform you about the outcome of this bug report submission (bug ID 0002596). I apologise for not doing so earlier, but I had a high workload.

The reply said that "Gauss linear gradient is very stable and robust and a good starting point for any simulation and we recommend it for general purpose use", and that "leastSquares can be more accurate than Gauss linear but is a lot less stable and could not be recommended for all simulations and cannot be the default". The report was on for a couple of days, and then its status was changed to "closed", its priority was changed to "normal", and its category was changed from "bug" to "feature".

Therefore, the zeroth-order accuracy issue is not denied, but it is belittled. In my opinion, the zeroth-order accuracy of the Gauss gradient is a major issue but because of its widespread use and associated embarrassment of admitting it openly many people will try to keep it unrevealed, until there is a gradual clandestine transition from the Gauss scheme to a consistent scheme. We had such an experience in trying to publish the arXiv paper I mentioned in my original post. Originally, we submitted it to the Journal of Computational Physics and it was rejected because the reviewers didn't think that the gradient inconsistency would pass on to the finite volume method as a whole (and to give them some credit, it is true that the preliminary version of the paper that we submitted did not include the present section 6, where the gradient is used for the FVM solution of a Poisson equation). Then we added section 6 and submitted it to Computers & Fluids. After more than 6 months of review and repeated emails to the journal, we were finally sent the reviewers' comments. The comments of the first reviewer were reasonable, but the second reviewer asked us to remove all the sections dealing with the Gauss gradient and leave only those on the least-squares one, because the inconsistency issue is "already known" and "redundancy should not rule the research and only originality must be sought" – he/she referred us to the conference paper by Sozer et al (2014) (now ref. [17] of our paper), which we had missed in our literature review. However, since our analysis goes much deeper we refused to remove the sections on the Gauss gradient. We then submitted the revised manuscript, and confidentially warned the editor that the 2nd reviewer is probably someone who has invested in the Gauss gradient and has an interest in keeping the inconsistency issue hidden. The result is that we received a "reject" decision, based on the follow-up comments of only the 2nd reviewer (we were not even sent the follow-up comments of the first reviewer) who now changed his mind completely: he said that the paper of Sozer et al is only a conference paper, not peer-reviewed, and therefore not credible, while he deemed our results also untrustworthy. He stated that it is known from published papers (none of which he cited) that the Gauss gradient is second-order accurate on unstructured meshes.

So, my advice is to conduct your own tests on simple problems where the accuracy can be readily checked instead of just relying on the popularity of the software. It seems that there are still unresolved fundamental issues with finite volume methods on unstructured grids.

[alexandross]

Originally Posted by akidess

You left out that the bug report is to be reopened once you contribute your code. I agree with Henry that this is not a bug but a feature. There is no default setting in OpenFOAM - if you do not specify a gradient scheme in fvSchemes, the code will exit. It is the general philosophy of OpenFOAM to leave such decisions to the user. The user needs to be aware of the limitations of the selected schemes (thank you for your contribution in that regard), and make appropriate choices.

My code is written in Fortran. I currently do not have the experience or the time to implement it in OpenFOAM, otherwise I would gladly contribute it. But reopening the bug report in this case does not mean that the accuracy issue of the Gauss gradient will be reconsidered. It merely means that the method I contribute will be considered as an additional gradient option, while the Gauss gradient will remain the default option. The Gauss linear scheme is clearly favoured in the users guide,

https://cfd.direct/openfoam/user-guide/fvschemes/

where it says right at the top of section 4.4.2 that

"The default discretisation scheme that is primarily used for gradient terms is:
default Gauss linear;
"
whereas the leastSquares scheme is described as "rarely used". The Gauss gradient is pre-selected in most of the example files distributed with OpenFOAM on which the users base their test cases.

I agree that this is not a bug in the typical sense, i.e. a bug in the code. I wrote this in this forum, and also in the bug report. What I suggested in the bug report is that the users are made aware of the inconsistency issue so that they can choose whether to use it or not, depending on their case, or so that they know what accuracy to expect, but my suggestion has not been followed yet as far as I know. Most users do not have the skills to assess the limitations of the scheme (probably even experts may be surprised by this finding) and just rely on the suggestions of the OpenFOAM (or other software) developers. Furthermore, there are people who try with great effort to conceal this issue, as I explained with regards to our experience in trying to publish the paper. So, how is the average user supposed to be well enough informed in order to make a judicious choice about what gradient scheme to use?

So, I think I did my part by raising this issue in this forum, and I appreciate your acknowledgement of that.

[guin]

Hi Alexandros,
Thanks again for your insightful analysis.

Originally Posted by alexandross

My code is written in Fortran. I currently do not have the experience or the time to implement it in OpenFOAM, otherwise I would gladly contribute it.\[...\]

To be honest, I do not have the time either. But there may be anyone partially idle in the community, who can get this development done. I'd suggest to start by sharing your code, even if it is currently in FORTRAN form. This can be done either right here or by means of any of the main fork projects:

https://bugs.openfoam.org/view.php?id=2596 -> Your bug report at the OpenFOAM Foundation. Probably they will not be able to accomplish this action without the corresponding funding, because it is mostly just a 3-guys team at CFD-Direct + a few individual volunteers. Anyway, this can be used as support for future actions e.g. reopening the report or creating a new one.

http://openfoam.com/contact/ -> ESI-OpenCFD team is responsible for the OpenFOAM-plus versions. They have more resources (bigger dedicated team) and they seem to have more active interaction with the user community when it comes to integrate new models into their OpenFOAM releases.

https://sourceforge.net/p/foam-extend/_members/ -> You can contact any of the few current developers of the FOAM-Extend project. It is a fully community-driven project and it would seem a right place to share any new model. But the development is quite slow compared to the other two, which comes from the fact that everyone here contributes in their free time.

Just a question: If I understood correctly, your proposal was only developed in 2D. Do you know how to proceed in order to extend it to 3D?

Originally Posted by alexandross

\[...\]So, I think I did my part by raising this issue in this forum, and I appreciate your acknowledgement of that.

Indeed we do appreciate your effort.