rortex.h

Go to the documentation of this file.

 
#ifndef VA_RORTEX_H
#define VA_RORTEX_H
 
#include "common_defs.h"
#include "linalg3.h"
 
/***************************************************************************/
VA_DEVICE_FUN void valib_rortex(VA_REAL *A, VA_REAL *rortex)
{
    int debug = 0;
 
    // get vector of vorticity
    VA_REAL vorticity[3];
    vorticity[0] = A[2+1*3] - A[1+2*3];
    vorticity[1] = A[0+2*3] - A[2+0*3];
    vorticity[2] = A[1+0*3] - A[0+1*3];
 
 
    if (debug) {
        printf("the vorticity vector is: [ %lf, %lf, %lf ]\n",
               vorticity[0], vorticity[1], vorticity[2]);
    }
 
    // cast the eigenvalue problem to finding roots of the cubic characteristic
    // equation
    // ax^3 + bx^2 + cx + d = 0
    //  x^3 - trace(A) x^2 + II_A x - det(A) = 0
    VA_REAL a, b, c, d;
 
    a = 1.;
 
    valib_trace3(A, &b);
    b = -b;
 
    valib_second_invariant3(A, &c);
 
    valib_determinant3(A, &d);
    d = -d;
 
    if (debug) {
       printf("  %s:\n", "A");
       int i, j;
       for (i = 0; i < 3; i++) {
          printf("   [ ");
          for (j = 0; j < 3; j++) {
             printf(" %13.6f ", A[j*3+i]);
          }
          printf(" ]\n");
       }
       printf("a, b, c, d:  %lf, %lf, %lf, %lf\n", a, b, c, d);
    }
 
    // find the three complex roots
    int info;
    VA_COMPLEX x1, x2, x3;
    valib_solve_cubic_equation(a, b, c, d, &x1, &x2, &x3, &info);
 
    VA_REAL lambda_r, lambda_ci;
 
    if (info == -1) {
        // only real eigenvalues exist, Rortex is zero
        rortex[0] = 0.;
        rortex[1] = 0.;
        rortex[2] = 0.;
 
        return;
    }
    else if (info == 0) {
        // search for eigenvalues was successful
        lambda_r  = creal(x1);
        lambda_ci = cimag(x2);
    }
    else {
        // there was another error in the search for eigenvalues
        printf("The solve of the cubic equation has failed with info: %d \n", info);
        rortex[0] = 0.;
        rortex[1] = 0.;
        rortex[2] = 0.;
        return;
    }
 
    // now find the corresponding eigenvector as the solution of
    // (A - lambda_r*I)x = 0
 
    // form matrix B = A - lambda_r*I
    VA_REAL B[9];
    for (int i = 0; i < 9; i++)
        B[i] = A[i];
    for (int i = 0; i < 3; i++)
        B[i+i*3] -= lambda_r;
 
    valib_find_nullspace3(B, rortex, &info);
    if (info < 0) {
        printf("Something is wrong, find_nullspace returned with info %d.\n", info);
        rortex[0] = 0.;
        rortex[1] = 0.;
        rortex[2] = 0.;
        return;
    }
 
    // print the eigenvector
    if (debug) {
        printf("the real eigenvector is: [ %lf, %lf, %lf ]\n",
               rortex[0], rortex[1], rortex[2]);
    }
 
    // ensure condition (14) by reverting the eigenvector if needed
    VA_REAL vr = valib_scalar_prod3(vorticity, rortex);
    if (vr < 0) {
        for (int i = 0; i < 3; i++)
            rortex[i] = -rortex[i];
    }
 
    vr = valib_scalar_prod3(vorticity, rortex);
    if (vr < 0) {
        printf("Something is wrong, the scalar product of vorticity and rortex should be positive, %lf.\n", vr);
        rortex[0] = 0.;
        rortex[1] = 0.;
        rortex[2] = 0.;
        return;
    }
 
    // compute the magnitude of Rortex by equation (13)
    VA_REAL diff = vr*vr - 4.*lambda_ci*lambda_ci;
    VA_REAL R;
    if (diff >= 0.)
       R = vr - sqrt(diff);
    else
       R = 0.;
 
    // scale the rortex
    for (int i = 0; i < 3; i++)
        rortex[i] *= R;
}
 
#endif // VA_RORTEX_H