corotation_common.h

 
#ifndef VA_COROTATION_COMMON_H
#define VA_COROTATION_COMMON_H
 
#include "common_defs.h"
#include "linalg3.h"
#include "spherical_integrators.h"
 
/***************************************************************************/
/***************************************************************************/
VA_DEVICE_FUN void valib_resvortstrain2D(VA_REAL *SO,
                                         VA_REAL *residual_vorticity,
                                         VA_REAL *residual_strain,
                                         VA_REAL *principal_axis)
{
    VA_REAL S_D_11, S_D_12;
    VA_REAL vort;
    VA_REAL as_D;
    VA_REAL length;
 
    // omega
    vort  = (SO[1] - SO[3]) / 2.;
 
    // values of the S_D tensor - deviatoric symmetric part of nabla u
    S_D_11 = 0.5 * (SO[0] - SO[4]);
    S_D_12 = 0.5 * (SO[1] + SO[3]);
 
    // |s_D| - positive eigenvalue of deviatoric strain-rate tensor S_D
    // version from the paper
    // as_D = sqrt ( (SO[0] - SO[4])*(SO[0] - SO[4])
    //              +(SO[1] + SO[3])*(SO[1] + SO[3])) / 2.;
    // using intermediate values
    as_D = sqrt( S_D_12 * S_D_12 + S_D_11 * S_D_11 );
 
    // try to use the first row for determining the eigenvector
    length = sqrt(((S_D_11 - as_D)*(S_D_11 - as_D) ) + (S_D_12 * S_D_12));
    if (length >= (VA_REAL) 1.e-5) {
        // diagonal is not an eigenvalue
        principal_axis[0] = S_D_12 / length;
        principal_axis[1] = - (S_D_11 - as_D) / length;
    }
    else {
        // try to use the second row for determining the eigenvector
        length = sqrt(((-S_D_11 - as_D)*(-S_D_11 - as_D))
                      +( S_D_12 * S_D_12 ) );
        if (length >= (VA_REAL) 1.e-5) {
           principal_axis[0] = (S_D_11 + as_D) / length;
           principal_axis[1] =  S_D_12 / length;
        }
        else {
        // the matrix contains just zeros, 0 is a multiple eigenvalue and any
        // vector is an eigenvector
           principal_axis[0] = (VA_REAL) 1.;
           principal_axis[1] = (VA_REAL) 0.;
        }
    }
 
    // omega_res and s_res
    if (fabs (vort) >= as_D) {
        *residual_vorticity = valib_sign( fabs(vort) - as_D,vort );
        *residual_strain    = (VA_REAL) 0.;
    }
    else {
        *residual_vorticity = (VA_REAL) 0.;
        *residual_strain    = as_D - fabs(vort);
    }
}
 
/***************************************************************************/
/***************************************************************************/
VA_DEVICE_FUN void valib_integrand_resvort(VA_REAL *n, VA_REAL *weight,
                                           VA_DEVICE_ADDR VA_REAL *A,
                                           VA_DEVICE_ADDR VA_REAL *avgcorot)
{
#if defined(CUDA)
    // put Q into shared memory
    extern __shared__ VA_REAL s_data[];
    VA_REAL *Q = &s_data[0];
#elif defined(OPENCL)
    // put Q into shared memory
    __private VA_REAL Q[9];
#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 vort_res, s_res, principal_axis[3];
 
    VA_REAL sina, cosa;
    VA_REAL sinb, cosb;
 
    cosb = n[2];
    sinb = sqrt(1. - cosb*cosb); // this is safe, points are never 0
    sina = n[1] / sinb;
    cosa = n[0] / sinb;
 
    // generate matrix Q
#if defined(CUDA)
    if (threadIdx.x == 0) {
#elif defined(OPENCL)
    if (get_local_id(0) == 0) {
#endif
       valib_constructQ3(sina, cosa, sinb, cosb, 0., 1., Q);
#if defined(CUDA)
    }
    __syncthreads();
#elif defined(OPENCL)
    }
    barrier(CLK_LOCAL_MEM_FENCE);
#endif
 
    // SO = QAQ^T
    valib_tatt3(A, Q, SO);
 
    // get residual vorticity in 2D
    valib_resvortstrain2D(SO, &vort_res, &s_res, principal_axis);
    // previous function only works on first two entries, correct the third
    principal_axis[2] = 0.;
 
    // add contribution to the surface integral in spherical coordinates
    avgcorot[0] = avgcorot[0]  + 6. * ( vort_res * n[0] ) * *weight;
    avgcorot[1] = avgcorot[1]  + 6. * ( vort_res * n[1] ) * *weight;
    avgcorot[2] = avgcorot[2]  + 6. * ( vort_res * n[2] ) * *weight;
 
#if defined(CUDA)
    __syncthreads();
#elif defined(OPENCL)
    barrier(CLK_LOCAL_MEM_FENCE);
#endif
}
 
/***************************************************************************/
VA_DEVICE_FUN void valib_integrand_resstrain(VA_REAL *n, VA_REAL *weight,
                                             VA_DEVICE_ADDR VA_REAL *A,
                                             VA_DEVICE_ADDR VA_REAL *S_RAVG)
{
    int i;
 
#if defined(CUDA)
    // put Q into shared memory
    extern __shared__ VA_REAL s_data[];
    VA_REAL *Q = &s_data[0];
#elif defined(OPENCL)
    // put Q into shared memory
    __private VA_REAL Q[9];
#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 localResStrain[9]; // tensor of residual strain returned
                               // to original coordinate system
 
    VA_REAL vort_res, s_res, principal_axis[3];
 
    VA_REAL sina, cosa;
    VA_REAL sinb, cosb;
 
    cosb = n[2];
    sinb = sqrt(1. - cosb*cosb); // this is safe, points are never 0
    sina = n[1] / sinb;
    cosa = n[0] / sinb;
 
    // generate matrix Q
#if defined(CUDA)
    if (threadIdx.x == 0) {
#elif defined(OPENCL)
    if (get_local_id(0) == 0) {
#endif
       valib_constructQ3(sina, cosa, sinb, cosb, 0., 1., Q);
#if defined(CUDA)
    }
    __syncthreads();
#elif defined(OPENCL)
    }
    barrier(CLK_LOCAL_MEM_FENCE);
#endif
 
    // SO = QAQ^T
    valib_tatt3(A, Q, SO);
 
    // get residual strain in 2D
    valib_resvortstrain2D(SO, &vort_res, &s_res, principal_axis);
    // previous function only works on first two entries
    principal_axis[2] = 0.;
 
    // now use SO to transform residual strain back
    for (i = 0; i < 9; i++) {
       SO[i] = (VA_REAL) 0.;
    }
    valib_rank1_update3(s_res, principal_axis, principal_axis, SO);
    // do the same for the second eigenvalue
    s_res = -s_res;
    principal_axis[2] =  principal_axis[0]; // to tmp
    principal_axis[0] = -principal_axis[1]; // orthogonal vector to [v1, v2]
                                            // in 2D is [-v2, v1]
    principal_axis[1] =  principal_axis[2]; // from tmp
    principal_axis[2] =  (VA_REAL) 0.;      // correct tmp
    valib_rank1_update3(s_res, principal_axis, principal_axis, SO);
 
    // transpose Q to multiply with it back
#if defined(CUDA)
    __syncthreads();
    if (threadIdx.x == 0) {
#elif defined(OPENCL)
    if (get_local_id(0) == 0) {
#endif
       valib_transpose3(Q);
#if defined(CUDA)
    }
    __syncthreads();
#elif defined(OPENCL)
    }
    barrier(CLK_LOCAL_MEM_FENCE);
#endif
 
    // result = Q^T*SO*Q
    valib_tatt3(SO, Q, localResStrain);
 
    // add contribution to the surface integral in spherical coordinates
    for (i = 0; i < 9; i++) {
       S_RAVG[i] = S_RAVG[i] + (VA_REAL) 2.5 * localResStrain[i] * *weight;
    }
 
#if defined(CUDA)
    __syncthreads();
#elif defined(OPENCL)
    barrier(CLK_LOCAL_MEM_FENCE);
#endif
}
 
/***************************************************************************/
VA_DEVICE_FUN void valib_integrand_surface(VA_REAL *n, VA_REAL *weight,
                                           VA_DEVICE_ADDR VA_REAL *A,
                                           VA_DEVICE_ADDR VA_REAL *result)
{
    result[0] = result[0] + *weight;
}
 
/***************************************************************************/
#endif // VA_COROTATION_COMMON_H