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

Possible bug (limitation?) in setValues

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

Reply
 
LinkBack Thread Tools Display Modes
Old   March 25, 2010, 14:37
Default 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
alberto will become famous soon enoughalberto will become famous soon enough
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
alberto is offline   Reply With Quote

Old   March 26, 2010, 04:16
Default
  #2
Senior Member
 
Join Date: Mar 2009
Posts: 854
Rep Power: 13
henry is on a distinguished road
setValues currently does not reset the diagonal, it uses the value the matrix provides to set the value.

H
henry is offline   Reply With Quote

Old   March 26, 2010, 05:29
Default
  #3
Senior Member
 
Hrvoje Jasak
Join Date: Mar 2009
Location: London, England
Posts: 1,763
Rep Power: 21
hjasak will become famous soon enough
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
hjasak is offline   Reply With Quote

Old   March 26, 2010, 12:04
Default
  #4
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
Thanks to both of you for your replies.

Quote:
Originally Posted by hjasak View Post
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
alberto is offline   Reply With Quote

Old   March 26, 2010, 20:45
Default
  #5
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 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
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
Serious bug in LES interface fs82 OpenFOAM Bugs 21 November 16, 2009 09:15
surfaceToPatch bug? bruce OpenFOAM Bugs 4 November 12, 2009 09:23
Please report this bug egp OpenFOAM Installation 5 December 8, 2006 13:56
Bug reports Mattijs Janssens (Mattijs) OpenFOAM 0 January 10, 2005 11:05
Forum y2k Bug Jonas Larsson Main CFD Forum 1 January 5, 2000 11:22


All times are GMT -4. The time now is 17:42.