I am getting a bit sick of this discussion: WHAT ON EARTH ARE YOU TALKING ABOUT???

So, I took a 1-D convection-diffusion equation and solved it in 1-D, with U = 1 and DT = 0.01. The boundary conditions are T = 1 at inlet (left) and T = 2 at outlet (right). Aaaaand: the solution is attached, not a single problem, no single issue, nothing: it just works!

So: what are you talking about? and Why do I have to waste my time with this?

Hrv

Hrv, this is we were talking about:

Test case explanation and exact solution is in #1 (FOAM case attached), examples of output in #2, correct formulation by Prof. Versteeg in #37.

BTW, What are your mesh step (in order to calculate Pe number) and div schemes?

Regards.

P.S. Another previous independent post on the same issue http://www.cfd-online.com/Forums/ope...tml#post280021

Dear Hrvoje,

I am really happy that you do get the right solution to the problem and was really surprised to initially read here that it was not, as I told you this afternoon (really enjoyed the course by the way).

So apparantly, the original posters input (and also my 'total novice OpenFOAM user' understanding of it) does not correspond to your input. Would you care to share the files?

Kind regards,
Duncan

The length is 1, number of cells is 100, U = 1, DT = 0.01, deltaT = 0.001. It is 1-D. I have increased the time-step to 0.1 (10 steps, maxCo = 10) and it still works perfectly fine.

If you have a VERY LARGE diffusivity, and a very large time-step, you will lose the diagonal dominance in the last cell, and the iterative solver will blow up (eventually). In fact, I am using ILU_0 preconditioning, so not even that will happen. If you are having trouble, either use a direct solver (tridiagonal matrix in 1-D) or reduce the time-step to recover diagonal dominance.

If you have very small diffusivity (0.0001), you need to be careful to preserve the matrix because with zero diffusivity the correct answer is a "shock" at the outlet - in fact this is the first solution reported as "wrong" and in fact it is right.

Hrv

P.S. The case is dead-trivial to set up and it is attached. I really do not understand how and why can someone get this wrong.

P.P.S. Where did you all go to school?

The original had a relatively low diffusivity, which probably made Santiago think to an instability. Then the discussion deviated on something else, admittedly not correct.

P.S. Please let us know when you will stop being so arrogant. You have been asked by e-mail about this problem before, and the answer you gave (posted above) was not exactly "it works", so it is possible to confuse things.

Hrv, Alberto, problem was:

V*dphi/dx-k*d^2phi/dx^2=0, phi(0)=1, phi(L)=0

exact solution,

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

completely bounded between 0 and 1. Until I understand this problem has no shock solution while diffusivity is not zero, it has a "boundary layer" at outlet whichs is so thinner as the magnitude of Pe (Problem Pe not mesh Peclet). We checked the matrix assembled in FOAM and it is clear that using the downwind value (when you really are using upwinded scheme) to assemble last equation is exactly the opposite that proffesor Versteeg explains in his book (see reference in post #37). In other words in Hrv Thesis is clear that BC's are the boss even in upwind scheme, but professor Versteeg get rid off the BC in the last cell and see upstream (to the upstream value at cell centroid), then diagonal dominance is assured again. It is quite interesting because he explains the issue exactly with this problem (i knew about this reference some weeks after I started this thread).

I ever maintained the topic centered in this issue, and Eugene proposed a solution. I agree with Alberto, following the answer he posted, it seemed to be something like "we kwnow about the issue but is not really important".

I'll download the Hrv's case at home and will test it, I hope a really academic and fundamented explanation could be given.

Regards.

P.S. My academic data is below, our alma mater is Sergio R. Idelsohn.

Hrv, I tested you case, try this:

-ddt(T) steadyState
-div(phi,T) Gauss Upwind
-DT={0.1, 0.01, 0.001}, so Pe_mesh={0.1, 1, 10}

Solutions are in the figure. Due the steady state, no tricks can be done with the timestep (i.e. increase the diagonal dominance by V/dt contribution), so if div is not correctly assembled for the upwind, solution goes wrong for high Peclet. The only contribution to diagonal you have in the last cell is the one from diffusivity, which is small. If you use the upwind cell centroid value for the face in the outlet (right border), the contribution to diagonal is Sf_{right}*U.

The same is explained by Versteeg in his book, he uses F_B*phi_P not F_B*phi_B (as in Eq. 3.58 of your thesis), see attached figures.

Regards.

Lets not get all excited.

Hrv, I can second Santiago's latest analysis. This is a fundamental issue, but it is only apparent if you use upwind or an upwind-biased scheme (more than 50% upwind) in conjunction with a fixed value for the transported property at the outlet. Without upwind face interpolation of the convected property there is no issue at all. I have personally had problems with this in several practical applications, so it is not just an academic concern.

PS. I went to the same school as you. ;)

... and here's one with Upwind as well.

Would you care to offer an explanation why this works for me and does not for you - having in mind we went to the same school?

Hrv

To reproduce this result, take the 1-D case above and change the convection discretisation to Gauss upwind.

I made no other changes and... it works!

Hrv, could you maintain all your settings and change

-ddt(T) steadyState
-DT={0.01, 0.001}, so Pe_mesh={1, 10}

and post your results?

Regards.

Only if I add under-relaxation to recover diagonal dominance in the matrix that you have just thrown away with the removal of the ddt-term.

If I do not do that, I will no longer preserve the M-matrix properties, and the ITERATIVE linear equation solver would blow up, JUST LIKE THEY TAUGHT ME AT SCHOOL!!!

My other option is to use a direct solver, so even that may work with TDMA. the point is that there is nothing wrong with discretisation nor with linear solvers nor with boundary conditions and you guys have been talking about this for a month.

Santiago, the case with DT = 0.01 works fine.

The case with DT = 0.001 will give you problems (you have no points in the BL). Refine the grid around the outlet (I know I am changing the conditions :-P), and it will still works.

Best,

Hrv, I resembled the FOAM problem in octave with a mate and we obtained the same matrix as in FOAM, then solving with an exact solver issue appears too (values goes out of bounds). Diagonal dominance is lost because FOAM assembles the last equation using downstream info instead of really doing upwind.

This is quite interesting, if one use Prof. Versteeg formulation for the last cell (post #47) all the things are good and solution from his book is obtained (bounded solution). We coded this in Matlab and FOAM and things were good.

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

Alberto, if you change the dx you are changing the mesh Peclet, then is a completely different numerical problem. Upwind ought give you a bounded solution for all finite Pe.

Regards.

I now get the correct solutions:
With central discretisation, when the mesh-Peclet number <2,
to guarantee the operator is of positive type, like expected
With upwind discretisation, irrespective of the mesh-Peclet number, like expected (operator will be of positive type irrespective of the mesh-Peclet number).
In both cases using the steady-state solver. Great, I am happy with this result!

Yes, I am aware of it, but your problem does not originate from the numerical scheme, but from the solution of the linear system originated from the discrete equation.

Reducing the step just brings back the matrix to a better situation, exactly as under-relaxing.

Best,

I do not agree. When you reduce the meshwidth with the central discretisation such that Pe_h<2, your spatial discetisation is of positive type, and your solution will obey a discrete maximum principle, just as the exact solution obeys a maximum principle: the solution will be monotonic.
Underrelaxation will make the iterative process converge....to a discrete solution that is not monotonic and therefore not physical.

I think we are mixing things (sorry if I contributed to that).

We are talking about upwind schemes, not central schemes, which we all know have a stability condition on the Peclet number.

Now, if you consider the upwind scheme, you have no limitation on Pe, but you observe, as Hrv pointed out, a problem due to the numerical iterative solvers of the associated linear system (loss of diagonal dominance => the iterative solver is not guaranteed to converge or be stable), when Pe is large. If you reduce the mesh size where this problem happens, you lower Pe, so you keep the problem under control (which does not mean upwind schemes have a limitation on Pe, it is just a trick to recover in some way the diagonal dominance). That was my observation.

Best,

I don't know exactly how it is implemented in OF, but normally the first order upwind scheme for the scalar CD equation will lead to an M-matrix irrespective of the mesh peclet number.
The solution will always be monotonic. The thickness of the boundary layer will be overestimated if the mesh-peclet number is large, but the qualitative behavior will be right (monotonic).

I used Prof. Hjasak input files, put the div scheme to upwind and got the correct solution for the upwind scheme.

If you set DT = 0.001, with that case, without relaxing, refining the grid next to the outlet and/or using a small time stepping in the case of unsteady runs, you won't get the correct solution for the reason explained above (you lose diagonal dominance).

The pictures attached show you what I get with OF 1.7.x, using scalarTransportFoam with

T0 = 1
T1 = 2
DT = 0.001
Upwind

and with the same settings, but relaxing with URF 0.001 (still wrong solution), and 0.00001. The case and code are attached, since I had to modify scalarTransportFoam to include relaxation.

P.S. Plots are made in paraview, so the scales are what they are...