[Duncan_vdH]

Can you tell me where exactly you lose the diagonal dominance?

The stencil for first order upwind scheme is (u, epsilon)>0

Internal points:

[-u -epsilon/h u+2*epsilon/h -epsilon/h]

First unknown:

[0 u+3*epsilon/h -epsilon/h]

using phi_0+phi_1=2*a

Last unknown

[-u-epsilon/h u+3*epsilon/h 0]

using phi_J+phi_(J+1)=2*b

Only using a central scheme, the scheme loses positivity when Pe_h>2

again, I am a OF novice, but I did get the right solution using Prof. Hjasak files modified for upwind,
but maybe my input is wrong.

[alberto]

In his case DT is large, that's why you get the right answer.

[Duncan_vdH]

I am sorry but no, I also tried the smaller values 0.001. In the upwind scheme you always get an M matrix, irrespective of the value of DT.

[alberto]

Do you use under-relaxation or use an unsteady solver (did you set steadyState in fvScheme?)?

The results I posted above are obtained exactly using Hrv case, with upwind scheme, obtaining a steady solution. The first picture is what I obtain without relaxing...

P.S. If you have a case with DT = 0.001, steady, without relaxation, I would like to see it, since I cannot really make it work in OF 1.7.x without relaxing it extremely heavily.

[santiagomarquezd]

Hi all, please let me post again (see post # 54) the spirit which I started this thread:

The main matter is, What one must do if had a downstream BC and is doing upwind? To use the BC value or to see upstream to the cell centroid to take a value for the boundary? If you follow Jasak you have to do one thing (the first one), if you follow Versteeg do the opposite. Which is conceptually correct?

Can anybody answer this fundamental question? Then let's talk about other related topics.

Regards.

[Duncan_vdH]

I think you should formulate your requirements, and then decide what to do:
The choice of how you incorporate the BC in the first order upwind scheme should:

• not screw up the positivity of the scheme->
  • lead to an M-matrix irrespective of the mesh peclet number
  • should not deteriorate the accuracy of the scheme-> be at least first order accurate in this case (automatically fulfilled for each consistent approximation).

So the choice of using the bc only to approximate the diffusive part of the flux fulfills both requirements.Any other choice you can up with that does not fulfill the requirements is wrong in my opinion.
A reasoning that 'upwind' is a scheme that should use information from 'upwind' is in my opinion not that valid. You are solving an elliptic equation, that does not have characteristics like a hyperbolic equation does. 'Upwind' is a first order scheme that has positivity irrespective of the mesh peclet number would be my definition.
Just my opinion. Unfortunately I cant post my results before monday, when I am back in my office.

[santiagomarquezd]

Duncan, respect of:

So the choice of using the bc only to approximate the diffusive part of the flux fulfills both requirements

this is I'm doing, like in Versteeg book, for the last cell in diffusive term, use the value given by the BC, because a centered scheme is used (CD), the upwind is used only for divergence term discretization (which has problems with centered schemes for Pe>1).

I think things would became much more clear if all people in the thread takes only 10 minutes, a pencil and a paper and solve the 5 cells problem proposed in Versteeg book. Then do the same using the info in Hrv. thesis. Systems obtained will be the same, except for the last equation. Assembling only the divergence term:

Versteeg

S_f*U*(-phi_4+phi_5)=0

Hrv

-S_f*U*phi_4+0*phi_5=S_f*U*BC

the issue is not centered in solving but in assembling.

Regards.

[eugene]

It seems clear (as mud) to me that there is no right or wrong here, just two different ways of approximating the same thing. I am not aware of any physical systems that are even conceptually similar to a fixed value outlet as is being discussed here, so I am not even sure "right" and "wrong" have any meaning in the sense that one approach can be said to approximate reality better than the other. Are the Versteeg boundaries better than the OpenFOAM fixed value boundaries? For most practical applications, probably not. If you want to compare artificial systems, then the two implementations can give different results. But what does this mean given that the systems being compared are artificial? The argument becomes something like "OpenFOAM behaves incorrectly when I try to model non-physical systems". It all gets a bit philosophical if you ask me.

The OpenFOAM fixedValue boundaries are what they are and do what they say on the tin. This happens not to be identical to what is detailed in some text books, but it is an approach that works just fine most of the time. If you want or need Versteeg type boundaries, they can be implemented with a modicum of effort, but they should not replace the existing boundaries. Both can co-exist in the same framework. Continuing to debate the correctness of the different approaches is not going to solve any problems.

For those who are interested, the simplest way to implement Versteeg-like fixed value outlets, is to make a new boundary derived from fixedValue with valueInternalCoeffs and valueBoundaryCoeffs replaced by those used in zeroGradient boundaries. Of course, these boundaries will not be convection scheme aware and should only be used with pure upwind and on outlets where all flux is going out of the domain. If I implement anything more general at some point in the future, I will leave a notice here.

[santiagomarquezd]

Eugene, I think things are not so polarized, Versteeg way to assemble BC's is the same as FOAM, except for the case of downstream BC and upwind scheme. For this particular case I've found that FOAM solution is wrong for high Pe.

I only want an explanation about why to use downwind info if I'm doing upwind. Some people argues that if you have a BC you must use it, but the fact is that this lead to an unbounded solution.

There is a related topic to this one. If one would use node centered FVM, like in Patankar's book, it allows to put the unknown just on the boundary, using a "half cell". In this formulation Hrv's and Verteeg's upwind schemes performs equal and resembles a Finite Difference solution.

In cell centered FVM, BC' aren't imposed in real unknown positions but in faces of volumes, so from my point of view the strong sense of a "BC" like in FDM or FEM cannot directly be applied, particularly if the way of reconstructing the face value leads to an unstable scheme.

Anyway... Imagine you are teaching and you want to explain why FOAM assembles the downstream BC with upwind in the way what is done. What would you say? (It's very probable that I have to give this explanation in the next semester...)

In Versteeg way I have the book and the facts supporting me, at least.

Regards.

P.S. You said that this issue appears in industrial cases too, didn't you?

[alberto]

Originally Posted by santiagomarquezd

I only want an explanation about why to use downwind info if I'm doing upwind. Some people argues that if you have a BC you must use it, but the fact is that this lead to an unbounded solution.

Well, as I told you, actually this is not the leading way of thinking. In most of the cases the BC is imposed weakly exactly to avoid numerical difficulties.

Anyway, I did what you said (took paper, pencil... and dig out some old code, which I should have done before ), and tried the two implementations.
Well, results are indeed quite interesting, since with Versteeg implementation, I can easily obtain a solution with Pe' = 100, whatever linear solver I use (I tried Gauss-Seidel and TDMA), while with the "OpenFOAM-like" implementation, a Pe' of 10 the correct solution is not obtained without altering the conditions in some way.

Santiago: do the result I posted above (see plots / codes in one of my previous post) match what you found with OF (you have a private message btw ) ?

[eugene]

Santiago, if you look in the code, Versteeg's way and the OpenFOAM way are not equivalent. They might produce equivalent results most of the time, but there is no consideration of the convection scheme in any non-coupled boundaries in FOAM. This is a fundamentally different approach - the FOAM way says the boundary conditions take precedence over the discretisation scheme, Versteeg's way says it is the other way around. In FOAM, the value on the boundary is thus the same for all operators. In Versteeg, the boundary simultaneously has multiple values. You could easily argue that Versteeg is the one who is confused!

Here is what I would say in a classroom: OpenFOAM does not implement the Versteeg fixed value boundary. It has an alternative approach. In OpenFOAM you are not using upwind or downwind info, in fact you are not considering the interpolation scheme (which is by definition approximate) at all. Rather you are using the more accurate exact information, i.e. the specified fixed value on the boundary for all purposes. What is the justification in Versteeg for using upwind interpolation to obtain a less accurate value at the boundary than the one already supplied? Numerical stability? From my perspective this falls into the realms of expediency and certainly is no more "correct" than the alternative, but neither is it "wrong".

I think we must try to make a distinction between 1) Not matching accepted dogma on a subject and 2) Containing a fundamental flaw - i.e. absolute correctness. FOAM does not match the generally accepted way of doing things and thus does not produce results in line with Versteeg. However, assuming that Versteeg and other established sources must be correct and everything that differs from them must be wrong is not necessarily a supportable position. Having the book on your side does not carry that much weight in my book. If someone can unequivocally prove to me that there is only one way to assemble the div operator in boundary cells, then I will gladly put my hands up, but so far all I have seen are "he said/she said" arguments. All instances where "correct" results have been obtained so far is when comparisons have been made to other numerical results using the same implementation. This is a circular argument and doesn't prove anything. Neither do I think you can argue Versteeg is correct because his solutions remain bounded - who says boundedness is an arbiter of correctness?

Don't get me wrong, I think we should implement Versteeg-like boundaries and that they are appropriate for a good range of problems, where they represent a more convenient and stable approximation than the alternative. Consider for example a porous jump outlet boundary with infinite resistance in the tangential direction (i.e. Utangential = (0 0 0) ). For everyday operating conditions, this outlet will not behave stably in conjunction with upwind differencing using a single value of U on the boundary (irrespective of whether it is provided by a gradient or a fixed value). So you might say "Ah, this proves Versteeg is correct because this real physical system is modelled more accurately by his approach than the FOAM approach." But it does not. Applying the equivalent of infinite resistance over zero space is in no way physical. In reality there is no such thing as a porous jump, it is only an approximation we use to make modelling some systems more convenient. I think it is the same for Versteeg - his div assembly is an approximation we should use to make modelling some systems easier. Saying that it is "correct" for this reason is like saying upwind differencing is correct and central differencing is wrong.

That's just how I see it at the moment.

[santiagomarquezd]

Originally Posted by eugene

Santiago, if you look in the code, Versteeg's way and the OpenFOAM way are not equivalent. They might produce equivalent results most of the time, but there is no consideration of the convection scheme in any non-coupled boundaries in FOAM. This is a fundamentally different approach - the FOAM way says the boundary conditions take precedence over the discretisation scheme, Versteeg's way says it is the other way around. In FOAM, the value on the boundary is thus the same for all operators. In Versteeg, the boundary simultaneously has multiple values. You could easily argue that Versteeg is the one who is confused!

Aha, I see the point, I forgot that in deed the difference is not only in upwind but in all schemes with more of 50% upwind, as you explained before.

Here is what I would say in a classroom: OpenFOAM does not implement the Versteeg fixed value boundary. It has an alternative approach. In OpenFOAM you are not using upwind or downwind info, in fact you are not considering the interpolation scheme (which is by definition approximate) at all. Rather you are using the more accurate exact information, i.e. the specified fixed value on the boundary for all purposes. What is the justification in Versteeg for using upwind interpolation to obtain a less accurate value at the boundary than the one already supplied? Numerical stability? From my perspective this falls into the realms of expediency and certainly is no more "correct" than the alternative, but neither is it "wrong".

I think we must try to make a distinction between 1) Not matching accepted dogma on a subject and 2) Containing a fundamental flaw - i.e. absolute correctness. FOAM does not match the generally accepted way of doing things and thus does not produce results in line with Versteeg. However, assuming that Versteeg and other established sources must be correct and everything that differs from them must be wrong is not necessarily a supportable position. Having the book on your side does not carry that much weight in my book. If someone can unequivocally prove to me that there is only one way to assemble the div operator in boundary cells, then I will gladly put my hands up, but so far all I have seen are "he said/she said" arguments. All instances where "correct" results have been obtained so far is when comparisons have been made to other numerical results using the same implementation. This is a circular argument and doesn't prove anything. Neither do I think you can argue Versteeg is correct because his solutions remain bounded - who says boundedness is an arbiter of correctness?

Boundedness is related to stability, when we use FOAM to do upwind in this problem we obtain an upwinded scheme but in the last cell, where a sort of central difference is done, because the use of upwind and downwind info. This is the reason of instability like in a pure central difference scheme. FOAM instabilities for upwind appears just after of Pe=1, exactly like in CD. It sounds like a 'corrupted' upwind difference, which not reach the objective of being stable to all Pe range, like real upwind do. Obviously, one can say 'Doing upwind is doing up-wind', if it improves stability or not is other discussion... my central concern is that FOAM says to do upwind, but finally it does something more likely to central differencing in the last cell.

Don't get me wrong, I think we should implement Versteeg-like boundaries and that they are appropriate for a good range of problems, where they represent a more convenient and stable approximation than the alternative. Consider for example a porous jump outlet boundary with infinite resistance in the tangential direction (i.e. Utangential = (0 0 0) ). For everyday operating conditions, this outlet will not behave stably in conjunction with upwind differencing using a single value of U on the boundary (irrespective of whether it is provided by a gradient or a fixed value). So you might say "Ah, this proves Versteeg is correct because this real physical system is modelled more accurately by his approach than the FOAM approach." But it does not. Applying the equivalent of infinite resistance over zero space is in no way physical. In reality there is no such thing as a porous jump, it is only an approximation we use to make modelling some systems more convenient. I think it is the same for Versteeg - his div assembly is an approximation we should use to make modelling some systems easier. Saying that it is "correct" for this reason is like saying upwind differencing is correct and central differencing is wrong.

That's just how I see it at the moment.

Aha, the point is that FOAM upwind is the only one upwind scheme that behaves as central difference in some cells, which is something I had never seen before, it makes me to think that FOAM don't do upwind but another different thing, let's call it 'FOAMUpwind'.

I think the base problem is what I posted before:

"There is a related topic to this one. If one would use node centered FVM, like in Patankar's book, it allows to put the unknown just on the boundary, using a "half cell". In this formulation Hrv's and Verteeg's upwind schemes performs equal and resembles a Finite Difference solution.

In cell centered FVM, BCs aren't imposed in real unknown positions but in faces of volumes, so from my point of view the strong sense of a "BC" like in FDM or FEM cannot directly be applied, particularly if the way of reconstructing the face value leads to an unstable scheme."

Regards.

[eugene]

Well I think we are pretty much in agreement then. FOAM imposes boundaries in the strong sense, most other approaches do so weakly to ensure stability for upwind and Pe > 1. This doesn't mean the upwind scheme is different or that it does central differencing in the last cell (although the effect might be equivalent), it just means boundary conditions are applied in the strong sense. This applies to all schemes and all boundary conditions.

I don't see this as right or wrong, just different and as long as the difference is understood, absolutely fine. The problem comes in when you assume the boundaries are applied weakly and then get numerical instabilities because they are applied in a strong sense. So perhaps some information along these lines should be put in a public place and weak alternative boundaries developed for applications that require them.

For interest sake, are the results for strong and weak boundaries identical for Pe < 1?

[eugene]

Just a quick update on a practical application that benefits from the approach discussed here. I have a boundary similar to Fluent's porous jump outlet. The normal velocity is zero gradient, but the tangential velocity should be zero. I used a fixedGradient boundary with the normal gradient set to zero and the tangential gradient set such that the velocity at the outlet would be zero. Convection scheme for U is upwind. Re_inlet ~ 1000, laminar.

I ran two cases, one using upwind for the gauss convection scheme outlet face value and the other using the standard fixedGradient specification. The difference between the two is simply setting valueBoundaryCoeffs to zero to use the upwind approach. The normal approach produces large discontinuities and unphysical behaviour near the outlet, while the modified formulation produces a smooth and much more plausible solution. (I would upload some pictures, but embarrassingly, I don't know how!)

For this kind of boundary the discretisation linked face formulation is thus of critical importance.

[alberto]

Originally Posted by eugene

I would upload some pictures, but embarrassingly, I don't know how!

You have to upload them somewhere and use the yellow button with the mountains in the toolbar, or attach them as files (clip icon).

Best,

[eugene]

Thanks Alberto. Here are the images in question (I hope!).

The one with the high velocity streak along the bottom (outlet) and left wall is the normal convection assembly formulation. The upwind formulation is just much better. Now admittedly this situation was set up to favour the upwind formulation, but it does illustrate the point.

[santiagomarquezd]

Originally Posted by eugene

For interest sake, are the results for strong and weak boundaries identical for Pe < 1?

Sorry for the delay, I had prepared some of this pics almost a month ago but I hadn't enough time to upload them and to write the message. They show a broad range of Pe number showing that solution differs from one scheme to another in all range. For Pe > 1 it is clear that actual FOAM BC scheme gives an unbounded solution whereas the other solution remains bounded. For Pe=1 and lower, solution from the weak scheme is more diffusive than analytic exactly how it should be (really ever is more diffusive, but in this range is more clear), on the other hand FOAM scheme give us less diffusive values (because the using of downstream info to assemble the scheme in the last cell). There is another pic with the analytical solution for Pe=1.

I will prepare more pics solving the problem for upwind scheme using finite difference (node based scheme) to compare.

Regards.

[jorkolino]

As I am solving a SS problem with scalar transport and BC fixedValue at outlet, I can confirm that problem with divergence exists and that changing just the BC to zerogradient fixes things perfectly. The equation solved:

dimensionedScalar LMRAsource
    (
    "LMRAsource", dimensionSet(0, 0, -1, 0, 0, 0, 0), -1.0
    ); 

    fvScalarMatrix LMRAEqn
    (
        fvm::div(phi, LMRA)
      - fvm::laplacian(kappaEff, LMRA)
    ==
    LMRAsource
    );

    LMRAEqn.solve();
    LMRAEqn.relax();

The Solution scheme:

LMRA
    {    
    solver          PBiCG;
        preconditioner  DILU;
    //smoother    DILU;
        tolerance       1e-05;
        relTol          0.1;
    }

My question is, how can I check the coefficients of the resulting matrix, as you have managed to do above? A command help will be greatly appreciated.

It is an issue indeed, because I only use stable cshemes: upwind for div and Gauss linear corrected for Laplace. And I do believe that the matrix coefficients are not composed correct for the near outlet cell(s). It is these cells that receive fantastic values in the process of itteration and flip flop the values in all other cells. This happens even for strong relaxation 0.1. But above all, just changing the BC to fixedgradient makes things sing.

[santiagomarquezd]

George, maybe changing the outlet BC is proper for the description of your particular physical problem, but in the case of the problem of this thread, using another outlet BC is solving a complete different problem. The solution comes by changing the way of BC are applied from 'strong' BCs to 'weak' BCs.

Respect of the inspection of matrix coefficients you have to go with a debugger (gdb for example) into the solve method until reach the point where BC's are imposed and system is finally solved. In case of scalarTransportFoam this point is after line 141 of fvScalarMatrix.C.
Then in order to examine the coefficients:

p this->diag().v_[start]@end
p this->lower().v_[start]@end
p this->upper().v_[start]@end

where start and end are the starting and ending indices of the piece of vector for diagonal, lower diagonal and upper diagonal that you want to inspect. There are many other ways to do this but it is simple and powerful.

Regards.

[jorkolino]

Originally Posted by eugene

It seems clear (as mud) to me that there is no right or wrong here, just two different ways of approximating the same thing. I am not aware of any physical systems that are even conceptually similar to a fixed value outlet as is being discussed here, so I am not even sure "right" and "wrong" have any meaning in the sense that one approach can be said to approximate reality better than the other. Are the Versteeg boundaries better than the OpenFOAM fixed value boundaries? For most practical applications, probably not. If you want to compare artificial systems, then the two implementations can give different results. But what does this mean given that the systems being compared are artificial? The argument becomes something like "OpenFOAM behaves incorrectly when I try to model non-physical systems". It all gets a bit philosophical if you ask me.
....

Hi,

I have to disagree that fixed values at outlets are to be rendered purely artificial or philosophical. I need to solve a parameter called "residual age of air" in a room, which has neatly defined theretically proven equation (standard convection, diffusion and fixed source terms), but also a very sound physical sense. It shows how much time is left for air molecules until they leave the domain. If you are familiar with the age of air parameter, this one is kind if opposite to it. And of course, due to the pure definition of this parameter, it has to be zero everywhere where it leaves the room. This is precisely the BC I need: zero value at all outlets. The problem I get with OF: unrealistic fluctuations in this and only this parameter, regardles of scheme and solution method. The fluctuations are so wild in the cells adjacent to outlet, that all domain values are ripped off for that parameter.

Another example I would give for using fixed outlet BC is the pressure BC. It is (almost) always fixed at outlets. There is no problem in solving the pressure though, because of the SIMPLE and others which create very different final matrix, and that is for pressure correction, not pressure itself.

I am still on a quest of composing stable matrix with zero outlet BC. Will post if anything comes around.