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

floating point exception error for eigenvalue calculation in applications/test/tensor

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

Reply
 
LinkBack Thread Tools Display Modes
Old   November 16, 2016, 11:11
Default floating point exception error for eigenvalue calculation in applications/test/tensor
  #1
New Member
 
SSSSS
Join Date: Jun 2011
Posts: 14
Rep Power: 7
doctorWho is on a distinguished road
I am using OpenFOAM 4.x and found that applications/test/tensor make a floating point exception error.

in the below code, vector e = eigenValues(t6) makes the error
Code:
    tensor t6(1,0,-4,0,5,4,-4,4,3);
    //tensor t6(1,2,0,2,5,0,0,0,0);

    Info<< "tensor " << t6 << endl;
    vector e = eigenValues(t6);
    Info<< "eigenvalues " << e << endl;
problematic part is in the member function eigenValues(const tensor& t) of src/OpenFOAM/primitives/Tensor/tensor/tensor.C

in the code below, Q makes zero for the test matrix (1 0 -4 0 5 4 -4 4 3) which has eigen values (-3, 3, 9). So when checking mag(PPP/QQ -1), it leads to floating point exception error before going to the third conditional PPP > QQ.

The same test code works for OpenFOAM 2.3.0 but not work for OpenFOAM 4.x (also 3.0.x). I don't know much about the OpenFOAM code but it looks like a bug after version 3.0.x. Maybe Henry Weller can fix this behavior.



Code:
    // non-diagonal matrix
    else
    {
        // Coefficients of the characteristic polynmial
        // x^3 + a*x^2 + b*x + c = 0
        scalar a =
           - t.xx() - t.yy() - t.zz();

        scalar b =
            t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
          - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz();

        scalar c =
          - t.xx()*t.yy()*t.zz()
          - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx()
          + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx();

        // Auxillary variables
        scalar aBy3 = a/3;

        scalar P = (a*a - 3*b)/9; // == -p_wikipedia/3
        scalar PPP = P*P*P;

        scalar Q = (2*a*a*a - 9*a*b + 27*c)/54; // == q_wikipedia/2
        scalar QQ = Q*Q;

        // Three identical roots
        if (mag(P) < SMALL && mag(Q) < SMALL)
        {
            return vector(- aBy3, - aBy3, - aBy3);
        }

        // Two identical roots and one distinct root
        else if (mag(PPP/QQ - 1) < SMALL)
        {
            scalar sqrtP = sqrt(P);
            scalar signQ = sign(Q);

            i = ii = signQ*sqrtP - aBy3;
            iii = - 2*signQ*sqrtP - aBy3;
        }
doctorWho 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
A floating point exception - SEM Model yansheng STAR-CCM+ 1 April 4, 2016 04:57
fluent divergence for no reason sufjanst FLUENT 2 March 23, 2016 17:08
Floating point exception from twoPhaseEulerFoam openfoammaofnepo OpenFOAM Running, Solving & CFD 1 March 19, 2016 14:56
icoFoam floating point exception (8) leizhao512 OpenFOAM Running, Solving & CFD 6 September 10, 2015 15:33
Floating point exception (core dumped) for GAMG solver yuhou1989 OpenFOAM Running, Solving & CFD 2 March 24, 2015 20:28


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