spherical_integrators.h

Go to the documentation of this file.

 
#ifndef VA_SPHERICAL_INTEGRATORS_H
#define VA_SPHERICAL_INTEGRATORS_H
 
#include "common_defs.h"
 
/***************************************************************************/
/***************************************************************************/
VA_DEVICE_FUN void valib_integrator_fibonacci(
    VA_DEVICE_ADDR void integrand(VA_REAL *n, VA_REAL *weight,
                                  VA_DEVICE_ADDR VA_REAL *A,
                                  VA_DEVICE_ADDR VA_REAL *result),
    VA_DEVICE_ADDR VA_REAL *A, VA_DEVICE_ADDR VA_REAL *result)
{
    int i,j;
 
    VA_REAL n[3];
    VA_REAL weight;
 
    // Which two numbers from the Fibonacci sequence should be used for
    // integration? A higher number corresponds to more integration points
    // and higher accuracy of integration. Value 15 can be recommended
    // for sigle precision, corresponding to
    // (2 * (15-2)th-Fibonacci-number)-2, i.e. 464, integration points.
    const int VA_FIBONACCI_LEVEL = 15;
 
    // find two Fibonacci numbers
    int precedingFibonacciNumber = 0;
    int aFibonacciNumber         = 1;
    for (i = 3; i < VA_FIBONACCI_LEVEL; i++ ) {
        int newNumber = aFibonacciNumber + precedingFibonacciNumber;
        precedingFibonacciNumber = aFibonacciNumber;
        aFibonacciNumber         = newNumber;
    }
 
    for (j = 1; j < aFibonacciNumber; j++) {
        for (i = 2*j-1; i < 2*j+1; i++) {
            VA_REAL zj_star   = 2.* j / (VA_REAL) aFibonacciNumber - 1.;
            VA_REAL phij_star = VA_PI * j
                              * (VA_REAL) precedingFibonacciNumber
                              / (VA_REAL) aFibonacciNumber;
 
            // normal components
            n[2] = zj_star + sin( VA_PI * zj_star ) / VA_PI;
            n[0] = pow((VA_REAL) -1.,(VA_REAL) 2.+i-2.*j) * cos(phij_star)
                 * sqrt( 1. - n[2]*n[2] );
            n[1] = pow((VA_REAL) -1.,(VA_REAL) 2.+i-2.*j) * sin(phij_star)
                 * sqrt( 1. - n[2]*n[2] );
 
            // weight
            weight = (1+cos(VA_PI * zj_star)) / (2.*(VA_REAL) aFibonacciNumber);
 
            // evaluate integrand at the given point
            integrand(n, &weight, A, result);
        }
    }
}
 
/***************************************************************************/
VA_DEVICE_FUN void valib_integrator_bazant(
    VA_DEVICE_ADDR void integrand(VA_REAL *n, VA_REAL *weight,
                                  VA_DEVICE_ADDR VA_REAL *A,
                                  VA_DEVICE_ADDR VA_REAL *result),
    VA_DEVICE_ADDR VA_REAL *A, VA_DEVICE_ADDR VA_REAL *result)
{
    #include "bazant_oh_2x61_integration_points.h"
 
    int iround;
    int ipoint;
 
    VA_REAL n[3];
    VA_REAL weight;
 
    for (iround = 0; iround < 2; iround++) {
        for (ipoint = 0; ipoint < 61; ipoint++) {
            // normal components
            n[0] = gpoints[ipoint][0];
            n[1] = gpoints[ipoint][1];
            n[2] = gpoints[ipoint][2];
 
            if (iround == 1) {
                // use opposite sign
                n[0] = -n[0];
                n[1] = -n[1];
                n[2] = -n[2];
            }
 
            // weight
            weight = gpoints[ipoint][3];
 
            // evaluate integrand at the given point
            integrand(n, &weight, A, result);
        }
    }
}
 
/***************************************************************************/
VA_DEVICE_FUN void valib_integrator_rectangle(
    VA_DEVICE_ADDR void integrand(VA_REAL *n, VA_REAL *weight,
                                  VA_DEVICE_ADDR VA_REAL *A,
                                  VA_DEVICE_ADDR VA_REAL *result),
    VA_DEVICE_ADDR VA_REAL *A, VA_DEVICE_ADDR VA_REAL *result, int numIntervals)
{
    VA_REAL n[3];
    VA_REAL weight;
 
    VA_REAL phi, theta;   // spherical coordinates on unit sphere
    VA_REAL sinP, cosP;   // sin(phi), cos(phi)
    VA_REAL sinT, cosT;   // sin(theta), cos(theta)
 
    // compute stepSize in radians
    VA_REAL stepSize = VA_PI / (VA_REAL) numIntervals;
 
    VA_REAL corrFact = stepSize * stepSize / ( (VA_REAL) 4. * VA_PI );
 
    // Perform numerical integration in spherical coordinates
    int iphi, itheta;
    for (itheta = 0; itheta < numIntervals; itheta++) {
        theta = itheta*stepSize + stepSize / (VA_REAL) 2.; // add half step to
                                                           // get the middle of
                                                           // the interval
#if defined(CUDA)
        sincos(theta, &sinT, &cosT);
#elif defined(OPENCL)
        sinT = native_sin(theta);
        cosT = native_cos(theta);
#else
        sinT = sin(theta);
        cosT = cos(theta);
#endif
 
        weight = sinT * corrFact;
 
        for (iphi = 0; iphi < 2*numIntervals; iphi++) {
            phi = iphi*stepSize + stepSize / (VA_REAL) 2.; // add half step to
                                                               // get the middle of
                                                               // the interval
#if defined(CUDA)
            sincos(phi, &sinP, &cosP);
#elif defined(OPENCL)
            sinP = native_sin(phi);
            cosP = native_cos(phi);
#else
            sinP = sin(phi);
            cosP = cos(phi);
#endif
 
            n[0] = sinT*cosP;
            n[1] = sinT*sinP;
            n[2] = cosT;
 
            // evaluate integrand at the given point
            integrand(n, &weight, A, result);
        }
    }
}
 
/***************************************************************************/
#endif // VA_SPHERICAL_INTEGRATORS_H