
[Sponsors] 
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,904
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.0e6) { doNotSolveMomentum.insert(cellI); } } 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); Code:
myA = UEqn.A(); forAll(rho,i) { if (rho[i] < 1.0e6) { Info << "rho = " << rho[i] << " A = " << myA[i] << endl; } } 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 in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

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,769
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 1e10 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,904
Rep Power: 26 
Thanks to both of you for your replies.
Quote:
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) ) ) Best,
__________________
Alberto Passalacqua GeekoCFD  A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

March 26, 2010, 20:45 

#5 
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,904
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 in both physical and virtual formats. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

Thread Tools  
Display Modes  


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 