
[Sponsors] 
September 20, 2010, 17:25 
Wrong fvm::div assembling

#1  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

September 21, 2010, 02:43 

#2 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27 
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,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 21, 2010, 10:35 

#3  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

September 21, 2010, 13:27 

#4 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27 
OK. I found the bug report. Sorry, for some reason my browser was not refreshing the page.
I see the problem too.
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 21, 2010, 14:05 

#5 
Senior Member
BastiL
Join Date: Mar 2009
Posts: 488
Rep Power: 12 
Maybe we should link Hentry to this post?


September 21, 2010, 14:35 

#6 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27 
I think they follow this forum.
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 22, 2010, 11:42 

#7 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
Similar problem was reported previously by another independent person:
Convectiondiffusion in 1D : wrong solution for a large Delta x 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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

September 23, 2010, 01:13 

#8 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27 
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,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 23, 2010, 08:57 

#9  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
Alberto, your observation is very clever
Quote:
Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

September 23, 2010, 19:38 

#10 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 13 
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. 

September 23, 2010, 22:23 

#11 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

September 24, 2010, 06:18 

#12 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 13 
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. 

September 24, 2010, 09:42 

#13  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
Quote:
Quote:
Quote:
Regards. P.S. I'll comment you about the progress.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

September 24, 2010, 10:42 

#14  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27 
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.
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 24, 2010, 20:52 

#15 
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

September 26, 2010, 00:51 

#16 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27 
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?
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 26, 2010, 15:05 

#17  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,907
Rep Power: 27 
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: Quote:
Quote:
Best, Alberto
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

September 27, 2010, 06:55 

#18 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 13 
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. 

September 27, 2010, 08:12 

#19 
Senior Member
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 13 
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. 

September 27, 2010, 12:12 

#20  
Senior Member
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 430
Rep Power: 16 
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe  Argentina. http://www.cimec.org.ar 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Warning: Dynamic zone with wrong CG using 6DOF  Manoj Kumar  FLUENT  1  August 11, 2012 04:03 
BuoyantBoussinesqSimpleFoam and axialsymmetric results wrong mass flow  Thomas Baumann  OpenFOAM  6  December 21, 2009 11:31 
udf error  srihari  FLUENT  0  February 9, 2009 10:00 
Pressure contour seems wrong???  Harry Qiu  FLUENT  1  June 29, 2001 05:53 
What's wrong with my UDF?  olivia  FLUENT  1  June 23, 2001 17:06 