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

Bug in eigenValue(const tensor& t)

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

Reply
 
LinkBack Thread Tools Display Modes
Old   February 21, 2012, 06:27
Default Bug in eigenValue(const tensor& t)
  #1
New Member
 
Evgeny Shchukin
Join Date: Feb 2012
Posts: 1
Rep Power: 0
EvgenyS is on a distinguished road
I have found a subtle bug in the eigenValue(const tensor& t) function. It has the following code:

...
else
{
scalar Q = (a*a - 3*b)/9;
scalar R = (2*a*a*a - 9*a*b + 27*c)/54;

scalar R2 = sqr(R);
scalar Q3 = pow3(Q);

// Three different real roots
if (R2 < Q3)
{
scalar sqrtQ = sqrt(Q);
scalar theta = acos(R/(Q*sqrtQ));
...

It is assumed here that if R2 < Q3 then |R|<=Q*sqrtQ, so that the argument of acos never exceeds 1.0 by absolute value. But it is not always the case. Consider the following values:

R = 00111011001011011110100111010100100101111010111101 01101010000111
= (8419873611799175)/(680564733841876926926749214863536422912) ~ 1.23719e-23,

Q = 00111100110000110100010111001001101110000010100100 01101110001100 = (1356189309486819)/(2535301200456458802993406410752) ~= 5.34922e-16.

Here the first value is 64 bits representing the corresponding double value according to the IEEE 754 binary64 format and the second is the exact value of these 64 bits in the form of a rational number. The last value is the approximation of this rational number in scientific notation. We have

R2 = 00110110011010111111011001111101111110000101001100 01101001001101 = (7870845268728397)/(5142201741628768881734278695491720328071049580104 9370729644032),

Q3 = 00110110011010111111011001111101111110000101001100 01101001001110 = (3935422634364199)/(2571100870814384440867139347745860164035524790052 4685364822016),

so that

R2-Q3 = -(1)/(5142201741628768881734278695491720328071049580104 9370729644032) < 0,

and R2 < Q3. On the other hand, we have

Q*sqrtQ = 00111011001011011110100111010100100101111010111101 01101010000110 = (4209936805899587)/(340282366920938463463374607431768211456),

and

R - Q*sqrtQ = (1)/(680564733841876926926749214863536422912) > 0,

and thus R/(Q*sqrtQ) > 1.0. In fact, we have

R/(Q*sqrtQ) = 00111111111100000000000000000000000000000000000000 00000000000001 = (4503599627370497)/(4503599627370496) > 1.

Most probably, the simplest possible solution of this problem is to replace the call

acos(R/(Q*sqrtQ))

by the call

acos(sqrt(R2/Q3))

It is not clear that there are no new problems with this approach, but at least it avoids arguments outside the range [-1, 1].
EvgenyS is offline   Reply With Quote

Old   February 21, 2012, 06:35
Default
  #2
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 919
Rep Power: 17
akidess will become famous soon enough
Better post this on the bug tracker - http://openfoam.com/bugs/
__________________
*On twitter @akidTwit
*Spend as much time formulating your questions as you expect people to spend on their answer.
*Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology
akidess is offline   Reply With Quote

Reply

Tags
eigenvalues

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
Problem/ bug in mesh.cellCells() ?? daZigeiner OpenFOAM Bugs 2 May 18, 2015 15:54
Serious bug in LES interface fs82 OpenFOAM Bugs 21 November 16, 2009 09:15
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 14:31.