triple_decomposition.h

Go to the documentation of this file.

 
#ifndef VA_TRIPLE_DECOMPOSITION_H
#define VA_TRIPLE_DECOMPOSITION_H
 
#include "common_defs.h"
#include "linalg3.h"
 
/***************************************************************************/
/***************************************************************************/
VA_DEVICE_FUN void valib_get_brf(VA_DEVICE_ADDR VA_REAL *A, int *num_intervals,
                                 VA_REAL *alpha_brf,
                                 VA_REAL *beta_brf,
                                 VA_REAL *gamma_brf,
                                 VA_REAL *goal_brf)
{
    int ialpha, ibeta, igamma;
 
#if defined(CUDA)
    // put Q into shared memory
    extern __shared__ VA_REAL s_data[];
    VA_REAL *Q  = &s_data[0];
#else
    VA_REAL Q[9];  // orthogonal matrices of transformations
#endif
    VA_REAL SO[9]; // symmetric (upper triangle) and antisymmetric (lower
                   // triangle) parts of velocity gradient
 
    VA_REAL alpha, beta, gamma; // rotation angles
    VA_REAL sina, cosa;         // sin(alpha), cos(alpha)
    VA_REAL sinb, cosb;         // sin(beta),  cos(beta)
    VA_REAL sing, cosg;         // sin(gamma), cos(gamma)
    VA_REAL goal_function;
 
    // compute stepSize in radians
    VA_REAL stepSize = VA_PI / (VA_REAL) *num_intervals;
 
    // initialize optimal angles
    *alpha_brf = 0.;
    *beta_brf  = 0.;
    *gamma_brf  = 0.;
 
    // initialize maxima
    *goal_brf = 0.;
 
 
    // perform search over all combinations of spherical coordinates
    for (ialpha = 0; ialpha < *num_intervals + 1; ialpha++) {
        alpha = ialpha*stepSize; // alpha in [0,pi]
#if defined(CUDA)
        sincos(alpha, &sina, &cosa);
#else
        sina = sin(alpha);
        cosa = cos(alpha);
#endif
        for (ibeta = 0; ibeta < *num_intervals + 1; ibeta++) {
            beta = ibeta*stepSize; // beta in [0,pi]
#if defined(CUDA)
            sincos(beta, &sinb, &cosb);
#else
            sinb = sin(beta);
            cosb = cos(beta);
#endif
            for (igamma = 0; igamma < *num_intervals / 2 + 1; igamma++) {
                gamma = igamma*stepSize; // gamma in [0,pi/2]
#if defined(CUDA)
                sincos(gamma, &sing, &cosg);
#else
                sing = sin(gamma);
                cosg = cos(gamma);
#endif
 
                // generate matrix Q
#if defined(CUDA)
                __syncthreads();
                if (threadIdx.x == 0) {
#endif
                    valib_constructQ3(sina, cosa, sinb, cosb, sing, cosg, Q);
#if defined(CUDA)
                }
                __syncthreads();
#endif
                // SO = QAQ^T
                valib_tatt3(A, Q, SO);
 
                // get residual vorticity in 2D
                valib_sym_antisym3(SO);
 
                // evaluate the goal function
                goal_function = fabs(SO[1]*SO[3])
                              + fabs(SO[2]*SO[6])
                              + fabs(SO[5]*SO[7]);
 
                // update maxima
                if (goal_function > *goal_brf) {
                    *alpha_brf = alpha;
                    *beta_brf  = beta;
                    *gamma_brf  = gamma;
 
                    *goal_brf = goal_function;
                }
            }
        }
    }
}
 
/***************************************************************************/
VA_DEVICE_FUN void valib_triple_decomposition_base(
    VA_DEVICE_ADDR VA_REAL *A,
    int *num_intervals,
    VA_DEVICE_ADDR VA_REAL *residual_vorticity_strain,
    VA_DEVICE_ADDR VA_REAL *shear_vorticity_strain,
    VA_DEVICE_ADDR VA_REAL *alpha,
    VA_DEVICE_ADDR VA_REAL *beta,
    VA_DEVICE_ADDR VA_REAL *gamma
)
{
    VA_REAL Q[9];                        // orthogonal matrices of
                                         // transformations
    VA_REAL QAQT[9];                     // transformed gradient QAQ^T
    VA_REAL alpha_brf, beta_brf, gamma_brf; // angles of basic reference
                                            // frame (BRF)
    VA_REAL goal_brf;                    // value of the shear scalar in BRF
    VA_REAL sina, cosa;                  // sin(alpha), cos(alpha)
    VA_REAL sinb, cosb;                  // sin(beta),  cos(beta)
    VA_REAL sing, cosg;                  // sin(gamma),  cos(gamma)
 
    int i;
    VA_REAL aux;
 
    // find BRF angles
    valib_get_brf(A, num_intervals, &alpha_brf, &beta_brf, &gamma_brf, &goal_brf);
 
    // Fill output variables with the results
    *alpha = alpha_brf;
    *beta  = beta_brf;
    *gamma = gamma_brf;
 
#if defined(CUDA)
    sincos(alpha_brf, &sina, &cosa);
    sincos(beta_brf,  &sinb, &cosb);
    sincos(gamma_brf, &sing, &cosg);
#else
    sina = sin(alpha_brf);
    cosa = cos(alpha_brf);
    sinb = sin(beta_brf);
    cosb = cos(beta_brf);
    sing = sin(gamma_brf);
    cosg = cos(gamma_brf);
#endif
 
    // generate matrix Q (this call is not shared for CUDA because this time,
    // each thread has different Q)
    valib_constructQ3(sina, cosa, sinb, cosb, sing, cosg, Q);
 
    // QAQ^T
    valib_tatt3(A, Q, QAQT);
 
    // Construct the residual tensor.
    // diagonal entries of the residual tensor are just copied
    for (i = 0; i < 3; i++) {
        residual_vorticity_strain[i*3+i] = QAQT[i*3+i];
    }
    // Find the residual vorticity and shear inside SO - filter out shear stress
    aux   = valib_min(fabs(QAQT[3]), fabs(QAQT[1]));
    residual_vorticity_strain[3] = valib_sign(aux, QAQT[3]);
    residual_vorticity_strain[1] = valib_sign(aux, QAQT[1]);
 
    aux   = valib_min(fabs(QAQT[6]), fabs(QAQT[2]));
    residual_vorticity_strain[6] = valib_sign(aux, QAQT[6]);
    residual_vorticity_strain[2] = valib_sign(aux, QAQT[2]);
 
    aux   = valib_min(fabs(QAQT[7]), fabs(QAQT[5]));
    residual_vorticity_strain[7] = valib_sign(aux, QAQT[7]);
    residual_vorticity_strain[5] = valib_sign(aux, QAQT[5]);
 
    // Construct the shear tensor as QAQT - residual tensor
    for (i = 0; i < 9; i++) {
        shear_vorticity_strain[i] = QAQT[i] - residual_vorticity_strain[i];
    }
 
    // Return to the original coordinate frame by applying the
    // inverse transformation
    valib_ttat3(residual_vorticity_strain, Q, QAQT);
    valib_mat_copy3(QAQT, residual_vorticity_strain);
    valib_ttat3(shear_vorticity_strain, Q, QAQT);
    valib_mat_copy3(QAQT, shear_vorticity_strain);
 
    // Find symmetric and antisymmetric parts of the residual tensor
    valib_sym_antisym3(residual_vorticity_strain);
    // Find symmetric and antisymmetric parts of the shear tensor
    valib_sym_antisym3(shear_vorticity_strain);
}
/***************************************************************************/
/***************************************************************************/
VA_DEVICE_FUN void valib_triple_decomposition(
    VA_DEVICE_ADDR VA_REAL *A,
    int *num_intervals,
    VA_DEVICE_ADDR VA_REAL *residual_vorticity,
    VA_DEVICE_ADDR VA_REAL *residual_strain,
    VA_DEVICE_ADDR VA_REAL *shear)
{
    VA_REAL residual_vorticity_strain[9];
    VA_REAL shear_vorticity_strain[9];
    VA_REAL alpha_brf, beta_brf, gamma_brf; // angles of basic reference
                                            // frame (BRF)
 
    valib_triple_decomposition_base(A, num_intervals, 
                                    residual_vorticity_strain,
                                    shear_vorticity_strain,
                                    &alpha_brf, &beta_brf, &gamma_brf);
    
    // ||Omega_RES||_F
    // read from the lower triangle of the residual_vorticity_strain array
    *residual_vorticity = sqrt(
         2.*valib_square(residual_vorticity_strain[1])
       + 2.*valib_square(residual_vorticity_strain[2])   
       + 2.*valib_square(residual_vorticity_strain[5])
    );
 
    // ||S_RES||_F
    // read from the upper triangle of the residual_strain array
    *residual_strain = sqrt(
         2.*valib_square(residual_vorticity_strain[3])
       + 2.*valib_square(residual_vorticity_strain[6])
       + 2.*valib_square(residual_vorticity_strain[7])
       +    valib_square(residual_vorticity_strain[0])
       +    valib_square(residual_vorticity_strain[4])
       +    valib_square(residual_vorticity_strain[8])
    );
 
    *shear = sqrt(
      // diagonal terms
        valib_square(shear_vorticity_strain[0])
      + valib_square(shear_vorticity_strain[4])
      + valib_square(shear_vorticity_strain[8])
      // upper square
      + valib_square(shear_vorticity_strain[3]-shear_vorticity_strain[1])
      + valib_square(shear_vorticity_strain[6]-shear_vorticity_strain[2])
      + valib_square(shear_vorticity_strain[7]-shear_vorticity_strain[5])
      // lower square
      + valib_square(shear_vorticity_strain[3]+shear_vorticity_strain[1])
      + valib_square(shear_vorticity_strain[6]+shear_vorticity_strain[2])
      + valib_square(shear_vorticity_strain[7]+shear_vorticity_strain[5])
    );
}
 
/***************************************************************************/
VA_DEVICE_FUN void valib_residual_vorticity(
    VA_DEVICE_ADDR VA_REAL *A,
    int *num_intervals,
    VA_DEVICE_ADDR VA_REAL *residual_vorticity,
    VA_DEVICE_ADDR VA_REAL *residual_strain,
    VA_DEVICE_ADDR VA_REAL *shear)
{
    VA_REAL residual_vorticity_strain[9];
    VA_REAL shear_vorticity_strain[9];
    VA_REAL alpha_brf, beta_brf, gamma_brf; // angles of basic reference
                                            // frame (BRF)
 
    valib_triple_decomposition_base(A, num_intervals, 
                                    residual_vorticity_strain,
                                    shear_vorticity_strain,
                                    &alpha_brf, &beta_brf, &gamma_brf);
    
    // Get omega_RES from the antisymmetric tensor.
    residual_vorticity[0] =  2. * residual_vorticity_strain[5];
    residual_vorticity[1] = -2. * residual_vorticity_strain[2];  
    residual_vorticity[2] =  2. * residual_vorticity_strain[1];
 
    // ||S_RES||_F
    // read from the upper triangle of the residual_strain array
    if (residual_strain != NULL) {
        *residual_strain = sqrt(
             2.*valib_square(residual_vorticity_strain[3])
           + 2.*valib_square(residual_vorticity_strain[6])
           + 2.*valib_square(residual_vorticity_strain[7])
           +    valib_square(residual_vorticity_strain[0])
           +    valib_square(residual_vorticity_strain[4])
           +    valib_square(residual_vorticity_strain[8])
        );
    }
 
    if (shear != NULL) {
        *shear = sqrt(
          // diagonal terms
            valib_square(shear_vorticity_strain[0])
          + valib_square(shear_vorticity_strain[4])
          + valib_square(shear_vorticity_strain[8])
          // upper square
          + valib_square(shear_vorticity_strain[3]-shear_vorticity_strain[1])
          + valib_square(shear_vorticity_strain[6]-shear_vorticity_strain[2])
          + valib_square(shear_vorticity_strain[7]-shear_vorticity_strain[5])
          // lower square
          + valib_square(shear_vorticity_strain[3]+shear_vorticity_strain[1])
          + valib_square(shear_vorticity_strain[6]+shear_vorticity_strain[2])
          + valib_square(shear_vorticity_strain[7]+shear_vorticity_strain[5])
        );
    }
}
 
/***************************************************************************/
VA_DEVICE_FUN void valib_triple_decomposition_4norms(
    VA_DEVICE_ADDR VA_REAL *A, int *num_intervals,
    VA_DEVICE_ADDR VA_REAL *residual_vorticity,
    VA_DEVICE_ADDR VA_REAL *residual_strain,
    VA_DEVICE_ADDR VA_REAL *shear_vorticity,
    VA_DEVICE_ADDR VA_REAL *shear_strain)
{
    VA_REAL residual_vorticity_strain[9];
    VA_REAL shear_vorticity_strain[9];
    VA_REAL alpha_brf, beta_brf, gamma_brf; // angles of basic reference
                                            // frame (BRF)
 
    valib_triple_decomposition_base(A, num_intervals, 
                                    residual_vorticity_strain,
                                    shear_vorticity_strain,
                                    &alpha_brf, &beta_brf, &gamma_brf);
    
    // ||Omega_RES||_F
    // read from the lower triangle of the residual_vorticity_strain array
    *residual_vorticity = sqrt(
         2.*valib_square(residual_vorticity_strain[1])
       + 2.*valib_square(residual_vorticity_strain[2])   
       + 2.*valib_square(residual_vorticity_strain[5])
    );
 
    // ||S_RES||_F
    // read from the upper triangle of the residual_strain array
    *residual_strain = sqrt(
         2.*valib_square(residual_vorticity_strain[3])
       + 2.*valib_square(residual_vorticity_strain[6])
       + 2.*valib_square(residual_vorticity_strain[7])
       +    valib_square(residual_vorticity_strain[0])
       +    valib_square(residual_vorticity_strain[4])
       +    valib_square(residual_vorticity_strain[8])
    );
 
    *shear_vorticity = sqrt(
      // lower square
         2.*valib_square(shear_vorticity_strain[1])
       + 2.*valib_square(shear_vorticity_strain[2])
       + 2.*valib_square(shear_vorticity_strain[5])
    );
 
    *shear_strain = sqrt(
         2.*valib_square(shear_vorticity_strain[3])
       + 2.*valib_square(shear_vorticity_strain[6])
       + 2.*valib_square(shear_vorticity_strain[7])
       +    valib_square(shear_vorticity_strain[0])
       +    valib_square(shear_vorticity_strain[4])
       +    valib_square(shear_vorticity_strain[8])
    );
}
 
/***************************************************************************/
VA_DEVICE_FUN void valib_triple_decomposition_with_angles(
    VA_DEVICE_ADDR VA_REAL *A,
    int *num_intervals,
    VA_DEVICE_ADDR VA_REAL *residual_vorticity,
    VA_DEVICE_ADDR VA_REAL *residual_strain,
    VA_DEVICE_ADDR VA_REAL *shear,
    VA_DEVICE_ADDR VA_REAL *alpha,
    VA_DEVICE_ADDR VA_REAL *beta,
    VA_DEVICE_ADDR VA_REAL *gamma
    )
{
    VA_REAL residual_vorticity_strain[9];
    VA_REAL shear_vorticity_strain[9];
    VA_REAL alpha_brf, beta_brf, gamma_brf; // angles of basic reference
                                            // frame (BRF)
 
    valib_triple_decomposition_base(A, num_intervals, 
                                    residual_vorticity_strain,
                                    shear_vorticity_strain,
                                    &alpha_brf, &beta_brf, &gamma_brf);
 
    // return the angles
    *alpha = alpha_brf;
    *beta  = beta_brf;
    *gamma = gamma_brf;
    
    // ||Omega_RES||_F
    // read from the lower triangle of the residual_vorticity_strain array
    *residual_vorticity = sqrt(
         2.*valib_square(residual_vorticity_strain[1])
       + 2.*valib_square(residual_vorticity_strain[2])   
       + 2.*valib_square(residual_vorticity_strain[5])
    );
 
    // ||S_RES||_F
    // read from the upper triangle of the residual_strain array
    *residual_strain = sqrt(
         2.*valib_square(residual_vorticity_strain[3])
       + 2.*valib_square(residual_vorticity_strain[6])
       + 2.*valib_square(residual_vorticity_strain[7])
       +    valib_square(residual_vorticity_strain[0])
       +    valib_square(residual_vorticity_strain[4])
       +    valib_square(residual_vorticity_strain[8])
    );
 
    *shear = sqrt(
      // diagonal terms
        valib_square(shear_vorticity_strain[0])
      + valib_square(shear_vorticity_strain[4])
      + valib_square(shear_vorticity_strain[8])
      // upper square
      + valib_square(shear_vorticity_strain[3]-shear_vorticity_strain[1])
      + valib_square(shear_vorticity_strain[6]-shear_vorticity_strain[2])
      + valib_square(shear_vorticity_strain[7]-shear_vorticity_strain[5])
      // lower square
      + valib_square(shear_vorticity_strain[3]+shear_vorticity_strain[1])
      + valib_square(shear_vorticity_strain[6]+shear_vorticity_strain[2])
      + valib_square(shear_vorticity_strain[7]+shear_vorticity_strain[5])
    );
}
 
#endif // VA_TRIPLE_DECOMPOSITION_H