Possible bug (limitation?) in setValues

 User Name Remember Me Password
 Register Blogs Members List Search Today's Posts Mark Forums Read

 LinkBack Thread Tools Display Modes
 March 25, 2010, 14:37 Possible bug (limitation?) in setValues #1 Senior Member   Alberto Passalacqua Join Date: Mar 2009 Location: Ames, Iowa, United States Posts: 1,894 Rep Power: 26 Hello, in trying to use setValues to exclude cells were the momentum equation is not defined (rho ~ 0) from the matrix, I found a potential bug/limitation (if I didn't do anything wrong in my code ). Let's assume the momentum equation is in the compressible form, and I want to exclude cells were rho < SMALL. I loop to identify those cells: Code: ```labelHashSet doNotSolveMomentum; forAll(rho, cellI) { if (rho[cellI] <= 1.0e-6) { doNotSolveMomentum.insert(cellI); } }``` Then I create the list of cell labels and the field of values I want to fix: Code: ```labelList noMomentumCellLabels = doNotSolveMomentum.toc(); vectorField fixedU (noMomentumCellLabels.size()); forAll(fixedU, cellI) { fixedU[cellI] = vector(0.0, 0.0, 0.0); } UEqn.setValues(noMomentumCellLabels, fixedU);``` Now, when I check the central coefficient, it is zero in the cells where the condition is satisfied with the code: Code: ```myA = UEqn.A(); forAll(rho,i) { if (rho[i] < 1.0e-6) { Info << "rho = " << rho[i] << " A = " << myA[i] << endl; } }``` I obtain a list of zeros for A, which causes a division by zero when computing U = H/A. Is this the expected behaviour of setValues? Shouldn't A be set to 1 in those cells, so that the equation actually degenerates in U = myValue? 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

 March 26, 2010, 04:16 #2 Senior Member   Join Date: Mar 2009 Posts: 854 Rep Power: 13 setValues currently does not reset the diagonal, it uses the value the matrix provides to set the value. H

 March 26, 2010, 05:29 #3 Senior Member   Hrvoje Jasak Join Date: Mar 2009 Location: London, England Posts: 1,758 Rep Power: 21 Hi Alberto, You cannot just go in and change the diagonal value to 1 because the order of magnitude of other terms in the matrix is unknown: other diagonals can be anything betwees, say 1e-10 and 1e10, depending what you're solving. If you give the diagonal of a wrong order of magnitude, it will change matrix normalization and that makes the solver behave in unexpected ways. I had this trouble with the FEM solvers a long time ago and the solution was to leave the diagonal as is. How come you have a zero diagonal in the first place? Can you estimate the size of the diagonal coefficient, eg. cell volume/delta t? Hope this helps, Hrv __________________ Hrvoje Jasak Providing commercial FOAM/OpenFOAM and CFD Consulting: http://wikki.co.uk

March 26, 2010, 12:04
#4
Senior Member

Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,894
Rep Power: 26
Thanks to both of you for your replies.

Quote:
 Originally Posted by hjasak How come you have a zero diagonal in the first place? Can you estimate the size of the diagonal coefficient, eg. cell volume/delta t?
The diagonal is zero because in those cells the equation is undefined (0 = 0), due to rho being zero (or very small).

At the moment I am using a flag to set up the fvMatrix, doing something like

Code:
```fvVectorMatrix UEqn
(
solveMomentum*
(
fvm::ddt(...).
+ ...
) +
doNotSolveMomentum*
(
fvm::Sp(Uunits,U)
)
)```
but I was looking for a cleaner and possibly less intrusive solution.

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

 March 26, 2010, 20:45 #5 Senior Member   Alberto Passalacqua Join Date: Mar 2009 Location: Ames, Iowa, United States Posts: 1,894 Rep Power: 26 I forgot to answer to one point: in my case both the cell volume and the time step are usually constant (static mesh, fixed dt), so their ratio is known. 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

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post fs82 OpenFOAM Bugs 21 November 16, 2009 09:15 bruce OpenFOAM Bugs 4 November 12, 2009 09:23 egp OpenFOAM Installation 5 December 8, 2006 13:56 Mattijs Janssens (Mattijs) OpenFOAM 0 January 10, 2005 11:05 Jonas Larsson Main CFD Forum 1 January 5, 2000 11:22

All times are GMT -4. The time now is 15:02.

 Contact Us - CFD Online - Top