[santiagomarquezd]

I posted this bug today in OpenFOAM's bug tracker

When fvm::div term is assembled in steady advection diffusion equation values in faces from BC are used even though the face cell is downstream, so if you use Gauss upwind scheme this value is not correct, because value at face must be taken from upstream cell center and not from BC.

Solve (case attached) the problem (1D) V*dphi/dx-k*d^2phi/dx^2=0, phi(0)=1, phi(L)=0. Exact solution must be:

phi(x)=(-e^(V/k*L)+e^(V/k*x))/(1-e^(V/k*L))

In case of high Peclet numbers, solution blows up, for example use V=1, L=2, k=0.001, dx=0.1, Pe=100 then solution goes to 17.0882 in x=1.95. Correct solution must be bounded between boundary conditions, i.e. 0<phi<1.
We tested this problem in 1.6.x and 1.7.x

case was closed and comment from Henry was:

It depends on the boundary conditions you choose.

I think is quite obvious it depends on the boundary conditions, but particularly for this boundary conditions FOAM does wrong in assembling this term, so it might be corrected in order to manage this kind of BCs that, BTW, are very simple.

We detected this problem working with a mate, but my advisor warned me about the thing too.

[alberto]

Hello,

you are saying that if a fixedValue condition is specified, and the condition is basically an outlet with a fixed value for a given property, that value is enforced?

If so, isn't it the expected behaviour with that BC, which all what it does is exactly to set the value on the face to the specified one?
If I understand your problem correctly, what you are asking is a mixed condition, which is zeroGradient when the flow is exiting, and fixedValue when there is backmixing, so the flow re-enters the domain. Such a condition exists already in OpenFOAM (inletOutlet).

P.S. I did not see the bugreport on OpenCFD mantis.

Best,

[santiagomarquezd]

Alberto,

Originally Posted by alberto

Hello,

you are saying that if a fixedValue condition is specified, and the condition is basically an outlet with a fixed value for a given property, that value is enforced?

If so, isn't it the expected behaviour with that BC, which all what it does is exactly to set the value on the face to the specified one?
If I understand your problem correctly, what you are asking is a mixed condition, which is zeroGradient when the flow is exiting, and fixedValue when there is backmixing, so the flow re-enters the domain. Such a condition exists already in OpenFOAM (inletOutlet).

P.S. I did not see the bugreport on OpenCFD mantis.

Best,

problem is simply steady advection-diffusion so you can set two BCs. You have an advective velocitiy V, a fixed value at inlet and a fixed value at outlet. Due to diffusion both fixed values can be satisfied. If Peclet number is below 1, solutions tends to be a line between 1 (fixed value at inlet) and 0 (fixed value at outlet). If Pe >> 1, fixed value at inlet is transported to outlet creating a 'boundary layer', this boudary layer is thinner when Pe is greater. Velocity is not being calculated, is a data from the problem (see attached case).

Correct solution have to be like advDiff.jpg (bounded between 0 and 1), but running the attached case you obtain something like test_1.7.png. We checked the resulting matrix in FOAM against the same matrix assembled using the theory and the problem is that FOAM is using the downstream BC to assemble the div operator, if you use Gauss Upwind doing this is wrong. In FOAM the problem is correctly assembled, but the last cell.

Problem is so simple and it's mathematically well posed... bug is actually in the bug tracker (# 44) but Henry closed it as I commented before in this post.

Regards.

[alberto]

OK. I found the bug report. Sorry, for some reason my browser was not refreshing the page.

I see the problem too.

[bastil]

Maybe we should link Hentry to this post?

[alberto]

I think they follow this forum.

[santiagomarquezd]

Similar problem was reported previously by another independent person:

http://www.cfd-online.com/Forums/ope...e-delta-x.html

Could you Alberto and Bastil check the attached model, specially in BCs and run it? It's small and run with scalarTransportFoam.

Thx in advance.

[alberto]

Hi,

you said the problem happens with upwind scheme. Is it limited to it? Did you try to use the linear scheme with an appropriate grid to satisfy the stability criterion?

Best,

[santiagomarquezd]

Alberto, your observation is very clever

Originally Posted by alberto

Hi,

you said the problem happens with upwind scheme. Is it limited to it? Did you try to use the linear scheme with an appropriate grid to satisfy the stability criterion?

Best,

I tried going down in Pe changing the diffusion and results were OK with central difference scheme (gauss linear), i.e. linear solution between inlet BC and outlet BC. I tried your suggestion too, Peclet number was 1, to do so I changed the grid spacing as you said and with central difference results were OK (see attached figure). Problems appears using upwind.

Regards.

[eugene]

I have have run into this issue a few times myself without knowing what the source of the problem was (other than "knowing" that advecting something into a fixed value outlet is bad). Some of my workarounds have been rather ... convoluted.

So you are saying the solution (and correct approach) is to ignore the specified boundary face value of the convected property in the div operator if that face is down stream of the cell centre and replace it with the upwind derived value?

I.e. for an outlet

fvm::div(phi, T)

T should use upwind interpolation get the face value of T instead of the specified boundary condition. However, for the other terms like diffusion that also depend on the boundary face value, the normal boundary value of T is used. Is this a correct summation?

This would be like having a separate zero gradient boundary condition for the div term part and a fixed value BC for the rest of the equation terms. I guess the magic would have to happen in the convection scheme code, i.e. gaussConvectionScheme.C.

No idea how to do it elegantly right now, but I'm pretty sure it could be hacked in.

[santiagomarquezd]

Eugene, your approach is interesting, few days ago we (the group in which I'm working) were wondering how to circumvent the issue, and it is related in fact with gaussConvection scheme as you said. My advisor is implementing a kind of FOAM emulator in Matlab and he asks to us for a lot of details of the code.
In Eq. 3.21, from Hrv, thesis it is clear that for upwind it is necessary to set phi_f as the upwind node value, nevertheless this criteria is not taken in account in Eq. 3.58 (when he discusses how to set BCs).
If term for node i is being assembled and phi_{i-1} is used for phi_{i-1/2} (upstream face) and phi_BC for phi_{i+1/2} (downstream face), finally no contribution for node i is done by fvm::div operator. i.e. advective contribution for K_{ii} is null, which we observed in FOAM using a debugger. This is like using some sort of "central difference" which is, as all of us know, unstable for Peclet number > 1. We implemented the problem by hand using Matlab, and if one take the resulting matrix and cancel the contribution of proper upwind sheme in K_{ii} it's possible to obtain the same wrong values.
Problem has a second derivative, then it accepts two fixed value BCs, but if upwind is used, downstream BC cannot be used to assemble fvm::div.
We prefer that developers fix the problem, I think they can do it elegantly like nobody more, but my attempt to report the bug in the bugtracker was useless. Maybe opening the issue to all the people can help the community or at least, if we are wrong, we could have a little bit more info that: "It depends on the boundary conditions you choose."

Regards.

[eugene]

I think the light has gone on. Its not a problem with convection in particular, although this is where you are most likely to encounter the issue, because it is the only operator which commonly uses a flux biased interpolation. It is a fundamental problem with the definition of the boundary conditions in FOAM. You just don't see it in the diffusion term because diffusion generally uses linear interpolation (not upwind). The great thing about it being a general problem, is that it can be fixed universally within the current framework. The bad thing is that a lot of OpenFOAM results so far would have been contaminated by this error (fortunately most to a very small degree).

The first question we have to ask, is "What is a boundary condition in FV?"
Second, "How does the value of the boundary condition interact with the boundary face value?" From what you are saying, the boundary condition is never the boss of the interpolation scheme, rather the reverse is true.

In all the cases I have checked, only the internal face values of a variable phi are affected by the interpolation scheme. This should be extended so that the boundary face values are similarly influenced.

If you look at the construction of the matrix boundary in gaussConvectionScheme, you will see that the internal and patch coefficients can be influenced by two things:
1. The boundary condition
2. The weights from the interpolation scheme.

However, when you check the code for fixedValue and zero gradient boundary conditions, you see firstly that the weights (and thus the interpolation scheme) are completely ignored. Secondly, if you look at the the interpolation schemes like linear upwind, you see that corrections are not applied to boundary values at all.

So, to correct the behaviour of uncorrected schemes like upwind, all you would have to do is change the boundary condition functions valueInternalCoeffs and valueBoundaryCoeffs such that they include the correct effects of the weights. If you want to fix corrected schemes like linear upwind, then you also need to fix the correction function in the interpolations schemes in question.

I have not yet though about what the exact form of these valueCoeff functions should be, but if you assume that for linear interpolation the boundary condition is exact, then some investigation of the boundary interpolation weights for the different schemes should tell you how to modify the functions to give the right results (at least for uncorrected schemes). This of course assumes that the boundary weights are consistent with their interior counterparts. I have to get back to work, but let me know if you make any progress, I might be able to provide more practical help at a later stage.

On the issue of getting the "developers" to fix this: I would say Henry is no more obligated to do something about this issue than me or anyone else, especially if he does not consider it a problem. As long as you are not his client there is no responsibility on his part either. That is the nature of open source: sometimes you have to get your fixes where you can find them and sometimes (if no one agrees with your assessment) you have to make them yourself.

[santiagomarquezd]

Originally Posted by eugene

I think the light has gone on. Its not a problem with convection in particular, although this is where you are most likely to encounter the issue, because it is the only operator which commonly uses a flux biased interpolation. It is a fundamental problem with the definition of the boundary conditions in FOAM. You just don't see it in the diffusion term because diffusion generally uses linear interpolation (not upwind). The great thing about it being a general problem, is that it can be fixed universally within the current framework. The bad thing is that a lot of OpenFOAM results so far would have been contaminated by this error (fortunately most to a very small degree).

I agree, it's more a conceptual error rather than an implementation's one. We use FOAM in teaching (among other things), and it is interesting that it is not possible to solve such simple problem (1D adv-diff).

The first question we have to ask, is "What is a boundary condition in FV?"
Second, "How does the value of the boundary condition interact with the boundary face value?" From what you are saying, the boundary condition is never the boss of the interpolation scheme, rather the reverse is true.

In all the cases I have checked, only the internal face values of a variable phi are affected by the interpolation scheme. This should be extended so that the boundary face values are similarly influenced.

If you look at the construction of the matrix boundary in gaussConvectionScheme, you will see that the internal and patch coefficients can be influenced by two things:
1. The boundary condition
2. The weights from the interpolation scheme.

However, when you check the code for fixedValue and zero gradient boundary conditions, you see firstly that the weights (and thus the interpolation scheme) are completely ignored. Secondly, if you look at the the interpolation schemes like linear upwind, you see that corrections are not applied to boundary values at all.

So, to correct the behaviour of uncorrected schemes like upwind, all you would have to do is change the boundary condition functions valueInternalCoeffs and valueBoundaryCoeffs such that they include the correct effects of the weights. If you want to fix corrected schemes like linear upwind, then you also need to fix the correction function in the interpolations schemes in question.

I have not yet though about what the exact form of these valueCoeff functions should be, but if you assume that for linear interpolation the boundary condition is exact, then some investigation of the boundary interpolation weights for the different schemes should tell you how to modify the functions to give the right results (at least for uncorrected schemes). This of course assumes that the boundary weights are consistent with their interior counterparts. I have to get back to work, but let me know if you make any progress, I might be able to provide more practical help at a later stage.

Ok, we'll consider your suggestions in order to implement a bugfix.

On the issue of getting the "developers" to fix this: I would say Henry is no more obligated to do something about this issue than me or anyone else, especially if he does not consider it a problem. As long as you are not his client there is no responsibility on his part either. That is the nature of open source: sometimes you have to get your fixes where you can find them and sometimes (if no one agrees with your assessment) you have to make them yourself.

Yes I know that. Even when you are paying for the software sometimes bug reports are not considered. The spirit of the comment was: "We think if bug is solved by developers it might benefit more users. If it is not possible, let know the community about it and maybe the users could give us some feedback". I think we have received a lot of, then, let's code!.

Regards.

P.S. I'll comment you about the progress.

[alberto]

Thanks eugene for the clarification.

One question. Does this problem affect the -dev/-ext release too?

Originally Posted by eugene

On the issue of getting the "developers" to fix this: I would say Henry is no more obligated to do something about this issue than me or anyone else, especially if he does not consider it a problem. As long as you are not his client there is no responsibility on his part either. That is the nature of open source: sometimes you have to get your fixes where you can find them and sometimes (if no one agrees with your assessment) you have to make them yourself.

It might be the nature of open source, but OpenFOAM is also a CFD code used for scientific purposes, and as such you do not expect a bug closed without any valid motivation.

Conceptual errors like this one throw more than one shadow on the quality of the code and its validation, and are detrimental to its reputation.

Developers do not have to fix this, but they cannot just close a bug ignoring a problem which is clear, and affects "textbook" cases, which should have been run as test. Claiming it is not a problem just contributes to make things worse.

[santiagomarquezd]

I'm downloading the -dev/-ext versión, it is taking a very long while... I'll post the news once download and compiling are finished. Meanwhile, I'm studying Eugene suggestions and checking them with gdb.

Best.

[alberto]

I would like to ask some more comment on Mantis, but it seems the bug cannot be re-opened once closed (or I am not familiar enough with that bug-reporting tool). Any idea if that's possible?

[alberto]

Since I had some doubt this "problem" has been left there for so long, I wrote to Hrv to ask for explanations

He replied that indeed it is not a bug in the fvm::div operator, but it depends on the boundary condition used, as already said by Henry on Mantis.

The quote of his answer follows:

The boundary condition is there to specify HOW to determine the boundary value and it will decide this by its type. The divergence term will use any face values that are served to it and is not allowed to interfere with the operation of the boundary condition class.

and

There are very good (stability, accuracy, boundedness) reasons WHY you are not allowed to use fixed value outlet boundaries; if you wish to have them, they require special treatment (in fact, they internally convert the "value" into a fixed gradient condition and adjust accordingly. I think we've got examples like this, eg. in advective conditions, but this is quite complex and serious.

So, basically you have to write an appropriate boundary condition for your case. In your textbook example, pay attention to how the boundary is treated when you discretize it, so that the problem is solved correctly.

Best,
Alberto

[eugene]

This was my point. It isn't 100% clear to me that there is a right or wrong in this case. I wrote a long reply about it on Saturday evening, but then lost it due to some browser hiccup. Grrr.

There is however a fundamental difference between what Hrv's reply and what Santiago is proposing. Santiago's outlet boundary would be stable, bounded and accurate, yet it would still be a fixed value. The conceptual difference comes from how one prioritises the interpolation scheme relative to the boundary value. Traditionally in FOAM, the boundary conditions are completely insensitive to the interpolation scheme. Thus, a fixed value is fixed on the boundary faces irrespective of the scheme used to calculate face values in the rest of the domain. Santiago is saying this is not the way it is supposed to be - interpolation schemes always take precedence over boundary conditions (as much as they can anyway).

This means that if the face interpolation scheme is pure upwind the resulting face values at an outlet are completely independent of the boundary condition for that particular operation. In fact, they are not only independent, the cannot know anything about what has been specified at the boundary! If a different operation uses a different interpolation scheme then the boundary face value changes too. The important thing to note is that the boundary value of the face can have different values for different operators in the same equation.

To my knowledge this has never been tried in FOAM and it may be a small yet pervasive and fundamental flaw. The reason why I don't make more noise about it is because I am not sure its not just a matter of opinion. There might be no "right" or "wrong" in this case, just schools of thought.

Santiago's approach also leaves me with a few uncertainties. If upwind imposes the requirement to use the cell centred value at the outlet face, what does it do at the inlet? Should we extrapolate to get the inlet face value assuming all fixed value specifications are linear interpolations or should we just "limit" the inlet value to the linear equivalent? What will this do to stability? I'm note sure, but these are interesting questions, the solution of which would solve some long standing issues and I for one would love to see some numerical experiments along these lines. I am thus waiting in some anticipation for Santiago's feedback.

[eugene]

I just had a more in-depth look at how the matrix is constructed from the boundary conditions and the values of the weighting field (for upwind at least).

Some things to note:
1. gaussConvection uses values on the boundary, gaussLaplacian uses gradients
2. The boundary condition interface supports passing weights to the value functions, but not the gradient functions. In practice this wont have much of an influence since the interpolation schemes for the Laplacian are generally linear.
3. The upwind weighting on the outlet face for convection is 1. I.e. on an internal face the value would have been set to the face owner value. Unfortunately, the weight for all boundary faces using linear interpolation is also 1.

This has to be changed. For this approach to work, the entire system has to be reconfigured with linear weight = 0.5 on the boundary (fvPatch.C). If a boundary weight of 0.5 = linear and 1 = upwind, the rest should just fall into place. Of course you would have to modify all boundaries with the valueInternalCoeffs and valueBoundaryCoeffs interfaces to reflect this change (fortunately there are only 10 or so of them). You also have to make sure the behaviour of inlets and other boundaries are correct as well.

[santiagomarquezd]

Thanks, Alberto & Eugene, some comments,

He replied that indeed it is not a bug in the fvm::div operator, but it depends on the boundary condition used, as already said by Henry on Mantis.

I was thinking a lot in Henry's answer and in writing to Hvr. about the matter, thanks for do it. Looking at the code is true that fvm::div is assembled in a generic way and does not have nothing to do with BCs. BCs influence is applied at the solving stage via addBoundaryDiag and addBoundarySource methods. Particularly addBoundaryDiag usses internalCoeffs. These values are calculated in fvm::div method using valueInternalCoeffs, this method is owned by BC class, in the case of fixedValue it returns zero, so that, we have:

fvm.internalCoeffs()[patchI] = patchFlux*psf.valueInternalCoeffs(pw);

internalCoeffs are zero and no contribution is done to the diagonal in case of downstream patch. In upstream patch (i.e. inlet) BC is assembled properly, because correct diagonal coefficient is already assembled by fvm::div.

The divergence term will use any face values that are served to it and is not allowed to interfere with the operation of the boundary condition class.

Following this spirit, again, the issue should be solved in BC class. Also problem maybe must named "Wrong fixedValue BC assignment" or something like this.

There are very good (stability, accuracy, boundedness) reasons WHY you are not allowed to use fixed value outlet boundaries; if you wish to have them, they require special treatment (in fact, they internally convert the "value" into a fixed gradient condition and adjust accordingly. I think we've got examples like this, eg. in advective conditions, but this is quite complex and serious.

Taking this comment and Henry's one, one can say "Ok problem is not in fvm::div, is in the BC, but at the same time correct BC is not implemented...". Now, what is the path? Enrich existing fixedValue BC to cope this type of problems of generate something like convectiveOutletFixedValue BC? Following Hrv. problem exists but it is usually circumvented due this kind of BC are not usual (and desirable) in CFD, nevertheless, it appears in textbook problems from numerical analysis.

Santiago's approach also leaves me with a few uncertainties. If upwind imposes the requirement to use the cell centred value at the outlet face, what does it do at the inlet? Should we extrapolate to get the inlet face value assuming all fixed value specifications are linear interpolations or should we just "limit" the inlet value to the linear equivalent? What will this do to stability? I'm note sure, but these are interesting questions, the solution of which would solve some long standing issues and I for one would love to see some numerical experiments along these lines. I am thus waiting in some anticipation for Santiago's feedback.

In FOAM matrices are assembled in a face basis, but if assembling is viewed from cell point of view things are clear. For first cell, whose usptream face is inlet, values for T * U & Sf (where T is the advected scalar) are zero for faces perpendicular to U (supposing 1D geometry and U going from inlet to outlet alingned with geometry). For inlet value is T_BC * U & Sf, for downstream face, T_i * U & Sf (where T_i is upstream cell centroid T value). On the other hand at the outlet values are: for upstream face T_(n-1) * U & Sf, and for downstream one T_n * U & Sf in stead of BC * U & Sf, because the upwind scheme.

If the problem is only advective, the last is wrong but because setting a BC at inlet is wrong (purely hyperbolic problem, only BC at inlet). If problem is advective-diffusive solution is correct. Nevertheless question is open, Have the BCs to be the boss in each term discretization or is it enough with having them in account in any part of the complete discretization (laplancian term in this case)?

Now I have a preliminar code of the BC, but the new valueInternalCoeffs method needs at least one more parameter to know what to do as Eugene said. I have to send the code to Nisi to fix some implementation issues that I cannont manage at all (he is in informatics, not me...).

Regards.