valib
Vortex Analysis LIBrary
spherical_integrators.h
1 
5 #ifndef VA_SPHERICAL_INTEGRATORS_H
6 #define VA_SPHERICAL_INTEGRATORS_H
7 
8 #include "common_defs.h"
9 
10 /***************************************************************************/
15 /***************************************************************************/
20 VA_DEVICE_FUN void valib_integrator_fibonacci(
21  VA_DEVICE_ADDR void integrand(VA_REAL *n, VA_REAL *weight,
22  VA_DEVICE_ADDR VA_REAL *A,
23  VA_DEVICE_ADDR VA_REAL *result),
24  VA_DEVICE_ADDR VA_REAL *A, VA_DEVICE_ADDR VA_REAL *result)
25 {
26  int i,j;
27 
28  VA_REAL n[3];
29  VA_REAL weight;
30 
31  // Which two numbers from the Fibonacci sequence should be used for
32  // integration? A higher number corresponds to more integration points
33  // and higher accuracy of integration. Value 15 can be recommended
34  // for sigle precision, corresponding to
35  // (2 * (15-2)th-Fibonacci-number)-2, i.e. 464, integration points.
36  const int VA_FIBONACCI_LEVEL = 15;
37 
38  // find two Fibonacci numbers
39  int precedingFibonacciNumber = 0;
40  int aFibonacciNumber = 1;
41  for (i = 3; i < VA_FIBONACCI_LEVEL; i++ ) {
42  int newNumber = aFibonacciNumber + precedingFibonacciNumber;
43  precedingFibonacciNumber = aFibonacciNumber;
44  aFibonacciNumber = newNumber;
45  }
46 
47  for (j = 1; j < aFibonacciNumber; j++) {
48  for (i = 2*j-1; i < 2*j+1; i++) {
49  VA_REAL zj_star = 2.* j / (VA_REAL) aFibonacciNumber - 1.;
50  VA_REAL phij_star = VA_PI * j
51  * (VA_REAL) precedingFibonacciNumber
52  / (VA_REAL) aFibonacciNumber;
53 
54  // normal components
55  n[2] = zj_star + sin( VA_PI * zj_star ) / VA_PI;
56  n[0] = pow((VA_REAL) -1.,(VA_REAL) 2.+i-2.*j) * cos(phij_star)
57  * sqrt( 1. - n[2]*n[2] );
58  n[1] = pow((VA_REAL) -1.,(VA_REAL) 2.+i-2.*j) * sin(phij_star)
59  * sqrt( 1. - n[2]*n[2] );
60 
61  // weight
62  weight = (1+cos(VA_PI * zj_star)) / (2.*(VA_REAL) aFibonacciNumber);
63 
64  // evaluate integrand at the given point
65  integrand(n, &weight, A, result);
66  }
67  }
68 }
69 
70 /***************************************************************************/
74 VA_DEVICE_FUN void valib_integrator_bazant(
75  VA_DEVICE_ADDR void integrand(VA_REAL *n, VA_REAL *weight,
76  VA_DEVICE_ADDR VA_REAL *A,
77  VA_DEVICE_ADDR VA_REAL *result),
78  VA_DEVICE_ADDR VA_REAL *A, VA_DEVICE_ADDR VA_REAL *result)
79 {
80  #include "bazant_oh_2x61_integration_points.h"
81 
82  int iround;
83  int ipoint;
84 
85  VA_REAL n[3];
86  VA_REAL weight;
87 
88  for (iround = 0; iround < 2; iround++) {
89  for (ipoint = 0; ipoint < 61; ipoint++) {
90  // normal components
91  n[0] = gpoints[ipoint][0];
92  n[1] = gpoints[ipoint][1];
93  n[2] = gpoints[ipoint][2];
94 
95  if (iround == 1) {
96  // use opposite sign
97  n[0] = -n[0];
98  n[1] = -n[1];
99  n[2] = -n[2];
100  }
101 
102  // weight
103  weight = gpoints[ipoint][3];
104 
105  // evaluate integrand at the given point
106  integrand(n, &weight, A, result);
107  }
108  }
109 }
110 
111 /***************************************************************************/
116 VA_DEVICE_FUN void valib_integrator_rectangle(
117  VA_DEVICE_ADDR void integrand(VA_REAL *n, VA_REAL *weight,
118  VA_DEVICE_ADDR VA_REAL *A,
119  VA_DEVICE_ADDR VA_REAL *result),
120  VA_DEVICE_ADDR VA_REAL *A, VA_DEVICE_ADDR VA_REAL *result, int numIntervals)
121 {
122  VA_REAL n[3];
123  VA_REAL weight;
124 
125  VA_REAL phi, theta; // spherical coordinates on unit sphere
126  VA_REAL sinP, cosP; // sin(phi), cos(phi)
127  VA_REAL sinT, cosT; // sin(theta), cos(theta)
128 
129  // compute stepSize in radians
130  VA_REAL stepSize = VA_PI / (VA_REAL) numIntervals;
131 
132  VA_REAL corrFact = stepSize * stepSize / ( (VA_REAL) 4. * VA_PI );
133 
134  // Perform numerical integration in spherical coordinates
135  int iphi, itheta;
136  for (itheta = 0; itheta < numIntervals; itheta++) {
137  theta = itheta*stepSize + stepSize / (VA_REAL) 2.; // add half step to
138  // get the middle of
139  // the interval
140 #if defined(CUDA)
141  sincos(theta, &sinT, &cosT);
142 #elif defined(OPENCL)
143  sinT = native_sin(theta);
144  cosT = native_cos(theta);
145 #else
146  sinT = sin(theta);
147  cosT = cos(theta);
148 #endif
149 
150  weight = sinT * corrFact;
151 
152  for (iphi = 0; iphi < 2*numIntervals; iphi++) {
153  phi = iphi*stepSize + stepSize / (VA_REAL) 2.; // add half step to
154  // get the middle of
155  // the interval
156 #if defined(CUDA)
157  sincos(phi, &sinP, &cosP);
158 #elif defined(OPENCL)
159  sinP = native_sin(phi);
160  cosP = native_cos(phi);
161 #else
162  sinP = sin(phi);
163  cosP = cos(phi);
164 #endif
165 
166  n[0] = sinT*cosP;
167  n[1] = sinT*sinP;
168  n[2] = cosT;
169 
170  // evaluate integrand at the given point
171  integrand(n, &weight, A, result);
172  }
173  }
174 }
175 
176 /***************************************************************************/
180 #endif // VA_SPHERICAL_INTEGRATORS_H
valib_integrator_bazant
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)
Definition: spherical_integrators.h:74
valib_integrator_fibonacci
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)
Definition: spherical_integrators.h:20
valib_integrator_rectangle
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)
Definition: spherical_integrators.h:116