Wrong fvm::div assembling
1 Attachment(s)
I posted this bug today in OpenFOAM's bug tracker
Quote:
Quote:
We detected this problem working with a mate, but my advisor warned me about the thing too. 
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 reenters the domain. Such a condition exists already in OpenFOAM (inletOutlet). P.S. I did not see the bugreport on OpenCFD mantis. Best, 
2 Attachment(s)
Alberto,
Quote:
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. 
OK. I found the bug report. Sorry, for some reason my browser was not refreshing the page.
I see the problem too. 
Maybe we should link Hentry to this post?

I think they follow this forum.

Similar problem was reported previously by another independent person:
http://www.cfdonline.com/Forums/ope...edeltax.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. 
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, 
1 Attachment(s)
Alberto, your observation is very clever
Quote:
Regards. 
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. 
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_{i1} is used for phi_{i1/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. 
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. 
Quote:
Quote:
Quote:
Regards. P.S. I'll comment you about the progress. 
Thanks eugene for the clarification.
One question. Does this problem affect the dev/ext release too? Quote:
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. 
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. 
I would like to ask some more comment on Mantis, but it seems the bug cannot be reopened once closed (or I am not familiar enough with that bugreporting tool). Any idea if that's possible?

Since I had some doubt this "problem" has been left there for so long, I wrote to Hrv to ask for explanations :D
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: Quote:
Quote:
Best, Alberto 
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. 
I just had a more indepth 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. 
Thanks, Alberto & Eugene, some comments,
Quote:
Code:
fvm.internalCoeffs()[patchI] = patchFlux*psf.valueInternalCoeffs(pw); Quote:
Quote:
Quote:
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 advectivediffusive 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. 
All times are GMT 4. The time now is 15:00. 