# Bug in eigenValue(const tensor& t)

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

 February 21, 2012, 06:27 Bug in eigenValue(const tensor& t) #1 New Member   Evgeny Shchukin Join Date: Feb 2012 Posts: 1 Rep Power: 0 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]. fumiya likes this.

 February 21, 2012, 06:35 #2 Senior Member     Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,264 Rep Power: 23 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.

 Tags eigenvalues

 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 daZigeiner OpenFOAM Bugs 2 May 18, 2015 15:54 fs82 OpenFOAM Bugs 21 November 16, 2009 09:15 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 23:23.