CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Bugs

Wrong fvm::div assembling

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree17Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   November 3, 2010, 16:15
Default
  #41
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,763
Rep Power: 21
hjasak will become famous soon enough
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
Attached Images
File Type: jpg lineX1.jpg (62.1 KB, 76 views)
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   November 3, 2010, 16:32
Default
  #42
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
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 Convection-diffusion in 1D : wrong solution for a large Delta x
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   November 3, 2010, 16:34
Default
  #43
New Member
 
Duncan Roy van der Heul
Join Date: May 2010
Posts: 14
Rep Power: 7
Duncan_vdH is on a distinguished road
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
Duncan_vdH is offline   Reply With Quote

Old   November 3, 2010, 17:26
Default
  #44
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,763
Rep Power: 21
hjasak will become famous soon enough
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?
Attached Files
File Type: gz 1D.tar.gz (1.8 KB, 26 views)
hua1015 likes this.
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   November 3, 2010, 17:51
Default
  #45
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
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.
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   November 3, 2010, 18:59
Default
  #46
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   November 3, 2010, 22:50
Default
  #47
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
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.
Attached Images
File Type: jpg mesh_Peclet.jpg (19.1 KB, 53 views)
File Type: jpg Versteeg_fig.jpg (5.5 KB, 46 views)
File Type: jpg Versteeg_problem.jpg (46.8 KB, 44 views)
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   November 5, 2010, 09:03
Default
  #48
Senior Member
 
Eugene de Villiers
Join Date: Mar 2009
Posts: 725
Rep Power: 12
eugene is on a distinguished road
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.
eugene is offline   Reply With Quote

Old   November 5, 2010, 10:25
Default
  #49
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,763
Rep Power: 21
hjasak will become famous soon enough
... 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
Attached Images
File Type: jpg upwind.jpg (43.9 KB, 36 views)
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   November 5, 2010, 10:26
Default
  #50
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,763
Rep Power: 21
hjasak will become famous soon enough
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!
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   November 5, 2010, 10:37
Default
  #51
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   November 5, 2010, 10:54
Default
  #52
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,763
Rep Power: 21
hjasak will become famous soon enough
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.
__________________
Hrvoje Jasak
Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk
hjasak is offline   Reply With Quote

Old   November 5, 2010, 10:56
Default
  #53
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
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,
dk1 likes this.
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Old   November 5, 2010, 11:37
Default
  #54
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 418
Rep Power: 15
santiagomarquezd will become famous soon enough
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.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D.
Post-doctoral Fellow
Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL
T.E.: 54-342-4511594 Ext. 1005
Güemes 3450 - (3000) Santa Fe
Santa Fe - Argentina
http://www.cimec.org.ar
santiagomarquezd is offline   Reply With Quote

Old   November 5, 2010, 12:56
Default
  #55
New Member
 
Duncan Roy van der Heul
Join Date: May 2010
Posts: 14
Rep Power: 7
Duncan_vdH is on a distinguished road
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!
Duncan_vdH is offline   Reply With Quote

Old   November 5, 2010, 13:33
Default
  #56
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
Quote:
Originally Posted by santiagomarquezd View Post
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.
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,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods

Last edited by alberto; November 5, 2010 at 14:52.
alberto is offline   Reply With Quote

Old   November 5, 2010, 16:06
Default
  #57
New Member
 
Duncan Roy van der Heul
Join Date: May 2010
Posts: 14
Rep Power: 7
Duncan_vdH is on a distinguished road
Quote:
Originally Posted by alberto View Post
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.
Duncan_vdH is offline   Reply With Quote

Old   November 5, 2010, 16:16
Default
  #58
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
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,
fumiya likes this.
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods

Last edited by alberto; November 5, 2010 at 16:20. Reason: Added clarification
alberto is offline   Reply With Quote

Old   November 5, 2010, 16:28
Default
  #59
New Member
 
Duncan Roy van der Heul
Join Date: May 2010
Posts: 14
Rep Power: 7
Duncan_vdH is on a distinguished road
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.

Last edited by Duncan_vdH; November 5, 2010 at 16:29. Reason: added 'first order' to upwind
Duncan_vdH is offline   Reply With Quote

Old   November 5, 2010, 17:57
Default
  #60
Senior Member
 
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
alberto will become famous soon enoughalberto will become famous soon enough
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...
Attached Images
File Type: jpg DT_0_001_Steady_Upwind_T0_1_T1_2_NoRelax.jpg (11.1 KB, 22 views)
File Type: jpg DT_0_001_Steady_Upwind_T0_1_T1_2_Relax_0_001.jpg (14.8 KB, 18 views)
File Type: jpg DT_0_001_Steady_Upwind_T0_1_T1_2_Relax_0_00001.jpg (13.0 KB, 20 views)
Attached Files
File Type: gz scalarTransportFoam.tar.gz (4.3 KB, 7 views)
File Type: gz 1D.tar.gz (6.6 KB, 5 views)
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as live DVD/USB, hard drive image and virtual image.
OpenQBMM - An open-source implementation of quadrature-based moment methods
alberto is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


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 axial-symmetric 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


All times are GMT -4. The time now is 22:04.