source: sasmodels/sasmodels/kernel_iq.cl @ 7b7da6b

core_shell_microgelscostrafo411magnetic_modelrelease_v0.94release_v0.95ticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since 7b7da6b was 7b7da6b, checked in by Paul Kienzle <pkienzle@…>, 8 years ago

remove unused variable warnings for pd_values/pd_weights
  • Property mode set to 100644
File size: 10.9 KB
Line 
1
2/*
3    ##########################################################
4    #                                                        #
5    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
6    #   !!                                              !!   #
7    #   !!  KEEP THIS CODE CONSISTENT WITH KERNELPY.PY  !!   #
8    #   !!                                              !!   #
9    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   #
10    #                                                        #
11    ##########################################################
12*/
13
14#ifndef _PAR_BLOCK_ // protected block so we can include this code twice.
15#define _PAR_BLOCK_
16
17typedef struct {
18#if MAX_PD > 0
19    int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable
20    int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector
21    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
22    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
23#endif // MAX_PD > 0
24    int32_t pd_prod;            // total number of voxels in hypercube
25    int32_t pd_sum;             // total length of the weights vector
26    int32_t num_active;         // number of non-trivial pd loops
27    int32_t theta_par;          // id of spherical correction variable
28} ProblemDetails;
29
30typedef struct {
31    PARAMETER_TABLE;
32} ParameterBlock;
33#endif // _PAR_BLOCK_
34
35
36#ifdef MAGNETIC
37
38// Return value restricted between low and high
39static double clip(double value, double low, double high)
40{
41  return (value < low ? low : (value > high ? high : value));
42}
43
44// Compute spin cross sections given in_spin and out_spin
45// To convert spin cross sections to sld b:
46//     uu * (sld - m_sigma_x);
47//     dd * (sld + m_sigma_x);
48//     ud * (m_sigma_y + 1j*m_sigma_z);
49//     du * (m_sigma_y - 1j*m_sigma_z);
50static void spins(double in_spin, double out_spin,
51    double *uu, double *dd, double *ud, double *du)
52{
53  in_spin = clip(in_spin, 0.0, 1.0);
54  out_spin = clip(out_spin, 0.0, 1.0);
55  *uu = sqrt(sqrt(in_spin * out_spin));
56  *dd = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin)));
57  *ud = sqrt(sqrt(in_spin * (1.0-out_spin)));
58  *du = sqrt(sqrt((1.0-in_spin) * out_spin));
59}
60
61#endif // MAGNETIC
62
63kernel
64void KERNEL_NAME(
65    int32_t nq,                 // number of q values
66    const int32_t pd_start,     // where we are in the polydispersity loop
67    const int32_t pd_stop,      // where we are stopping in the polydispersity loop
68    global const ProblemDetails *details,
69    global const double *values,
70    global const double *q, // nq q values, with padding to boundary
71    global double *result,  // nq+1 return values, again with padding
72    const double cutoff     // cutoff in the polydispersity weight product
73    )
74{
75
76  // who we are and what element we are working with
77  const int q_index = get_global_id(0);
78  if (q_index >= nq) return;
79
80  // Storage for the current parameter values.  These will be updated as we
81  // walk the polydispersity cube.  local_values will be aliased to pvec.
82  ParameterBlock local_values;
83  double *pvec = (double *)&local_values;
84
85  // Fill in the initial variables
86  for (int i=0; i < NPARS; i++) {
87    pvec[i] = values[2+i];
88//if (q_index==0) printf("p%d = %g\n",i, pvec[i]);
89  }
90
91#ifdef MAGNETIC
92  // Location of the sld parameters in the parameter pvec.
93  // These parameters are updated with the effective sld due to magnetism.
94  const int32_t slds[] = { MAGNETIC_PARS };
95
96  const double up_frac_i = values[NPARS+2];
97  const double up_frac_f = values[NPARS+3];
98  const double up_angle = values[NPARS+4];
99  #define MX(_k) (values[NPARS+5+3*_k])
100  #define MY(_k) (values[NPARS+6+3*_k])
101  #define MZ(_k) (values[NPARS+7+3*_k])
102
103  // TODO: could precompute these outside of the kernel.
104  // Interpret polarization cross section.
105  double uu, dd, ud, du;
106  double cos_mspin, sin_mspin;
107  spins(up_frac_i, up_frac_f, &uu, &dd, &ud, &du);
108  SINCOS(-up_angle*M_PI_180, sin_mspin, cos_mspin);
109#endif // MAGNETIC
110
111  double pd_norm, this_result;
112  if (pd_start == 0) {
113    pd_norm = this_result = 0.0;
114  } else {
115    pd_norm = result[nq];
116    this_result = result[q_index];
117  }
118//if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, this_result);
119
120#if MAX_PD>0
121  global const double *pd_value = values + NUM_VALUES + 2;
122  global const double *pd_weight = pd_value + details->pd_sum;
123#endif
124
125  // Jump into the middle of the polydispersity loop
126#if MAX_PD>4
127  int n4=details->pd_length[4];
128  int i4=(pd_start/details->pd_stride[4])%n4;
129  const int p4=details->pd_par[4];
130  global const double *v4 = pd_value + details->pd_offset[4];
131  global const double *w4 = pd_weight + details->pd_offset[4];
132#endif
133#if MAX_PD>3
134  int n3=details->pd_length[3];
135  int i3=(pd_start/details->pd_stride[3])%n3;
136  const int p3=details->pd_par[3];
137  global const double *v3 = pd_value + details->pd_offset[3];
138  global const double *w3 = pd_weight + details->pd_offset[3];
139//if (q_index==0) printf("offset %d: %d %d\n", 3, details->pd_offset[3], NUM_VALUES);
140#endif
141#if MAX_PD>2
142  int n2=details->pd_length[2];
143  int i2=(pd_start/details->pd_stride[2])%n2;
144  const int p2=details->pd_par[2];
145  global const double *v2 = pd_value + details->pd_offset[2];
146  global const double *w2 = pd_weight + details->pd_offset[2];
147#endif
148#if MAX_PD>1
149  int n1=details->pd_length[1];
150  int i1=(pd_start/details->pd_stride[1])%n1;
151  const int p1=details->pd_par[1];
152  global const double *v1 = pd_value + details->pd_offset[1];
153  global const double *w1 = pd_weight + details->pd_offset[1];
154#endif
155#if MAX_PD>0
156  int n0=details->pd_length[0];
157  int i0=(pd_start/details->pd_stride[0])%n0;
158  const int p0=details->pd_par[0];
159  global const double *v0 = pd_value + details->pd_offset[0];
160  global const double *w0 = pd_weight + details->pd_offset[0];
161#endif
162
163
164  double spherical_correction=1.0;
165  const int theta_par = details->theta_par;
166#if MAX_PD>0
167  const int fast_theta = (theta_par == p0);
168  const int slow_theta = (theta_par >= 0 && !fast_theta);
169#else
170  const int slow_theta = (theta_par >= 0);
171#endif
172
173  int step = pd_start;
174
175
176#if MAX_PD>4
177  const double weight5 = 1.0;
178  while (i4 < n4) {
179    pvec[p4] = v4[i4];
180    double weight4 = w4[i4] * weight5;
181//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, pvec[p4], weight4);
182#elif MAX_PD>3
183    const double weight4 = 1.0;
184#endif
185#if MAX_PD>3
186  while (i3 < n3) {
187    pvec[p3] = v3[i3];
188    double weight3 = w3[i3] * weight4;
189//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, pvec[p3], weight3);
190#elif MAX_PD>2
191    const double weight3 = 1.0;
192#endif
193#if MAX_PD>2
194  while (i2 < n2) {
195    pvec[p2] = v2[i2];
196    double weight2 = w2[i2] * weight3;
197//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, pvec[p2], weight2);
198#elif MAX_PD>1
199    const double weight2 = 1.0;
200#endif
201#if MAX_PD>1
202  while (i1 < n1) {
203    pvec[p1] = v1[i1];
204    double weight1 = w1[i1] * weight2;
205//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, pvec[p1], weight1);
206#elif MAX_PD>0
207    const double weight1 = 1.0;
208#endif
209    if (slow_theta) { // Theta is not in inner loop
210      spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6);
211    }
212#if MAX_PD>0
213  while(i0 < n0) {
214    pvec[p0] = v0[i0];
215    double weight0 = w0[i0] * weight1;
216//if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, pvec[p0], weight0);
217    if (fast_theta) { // Theta is in inner loop
218      spherical_correction = fmax(fabs(cos(M_PI_180*pvec[p0])), 1.e-6);
219    }
220#else
221    const double weight0 = 1.0;
222#endif
223
224//if (q_index == 0) {printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NPARS; i++) printf("p%d=%g ",i, pvec[i]); printf("\n"); }
225//if (q_index == 0) printf("sphcor: %g\n", spherical_correction);
226
227    #ifdef INVALID
228    if (!INVALID(local_values))
229    #endif
230    {
231      // Accumulate I(q)
232      // Note: weight==0 must always be excluded
233      if (weight0 > cutoff) {
234        // spherical correction has some nasty effects when theta is +90 or -90
235        // where it becomes zero.
236        const double weight = weight0 * spherical_correction;
237        pd_norm += weight * CALL_VOLUME(local_values);
238
239#ifdef MAGNETIC
240          const double qx = q[2*q_index];
241          const double qy = q[2*q_index+1];
242          const double qsq = qx*qx + qy*qy;
243
244          // Constant across orientation, polydispersity for given qx, qy
245          double px, py, pz;
246          if (qsq > 1.e-16) {
247            px = (qy*cos_mspin + qx*sin_mspin)/qsq;
248            py = (qy*sin_mspin - qx*cos_mspin)/qsq;
249            pz = 1.0;
250          } else {
251            px = py = pz = 0.0;
252          }
253
254          double scattering = 0.0;
255          if (uu > 1.e-8) {
256            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
257                const double perp = (qy*MX(sk) - qx*MY(sk));
258                pvec[slds[sk]] = (values[slds[sk]+2] - perp*px)*uu;
259            }
260            scattering += CALL_IQ(q, q_index, local_values);
261          }
262
263          if (dd > 1.e-8){
264            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
265                const double perp = (qy*MX(sk) - qx*MY(sk));
266                pvec[slds[sk]] = (values[slds[sk]+2] + perp*px)*dd;
267            }
268            scattering += CALL_IQ(q, q_index, local_values);
269          }
270          if (ud > 1.e-8){
271            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
272                const double perp = (qy*MX(sk) - qx*MY(sk));
273                pvec[slds[sk]] = perp*py*ud;
274            }
275            scattering += CALL_IQ(q, q_index, local_values);
276            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
277                pvec[slds[sk]] = MZ(sk)*pz*ud;
278            }
279            scattering += CALL_IQ(q, q_index, local_values);
280          }
281          if (du > 1.e-8) {
282            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
283                const double perp = (qy*MX(sk) - qx*MY(sk));
284                pvec[slds[sk]] = perp*py*du;
285            }
286            scattering += CALL_IQ(q, q_index, local_values);
287            for (int sk=0; sk<NUM_MAGNETIC; sk++) {
288                pvec[slds[sk]] = -MZ(sk)*pz*du;
289            }
290            scattering += CALL_IQ(q, q_index, local_values);
291          }
292#else  // !MAGNETIC
293        const double scattering = CALL_IQ(q, q_index, local_values);
294#endif // !MAGNETIC
295        this_result += weight * scattering;
296      }
297    }
298    ++step;
299#if MAX_PD>0
300    if (step >= pd_stop) break;
301    ++i0;
302  }
303  i0 = 0;
304#endif
305#if MAX_PD>1
306    if (step >= pd_stop) break;
307    ++i1;
308  }
309  i1 = 0;
310#endif
311#if MAX_PD>2
312    if (step >= pd_stop) break;
313    ++i2;
314  }
315  i2 = 0;
316#endif
317#if MAX_PD>3
318    if (step >= pd_stop) break;
319    ++i3;
320  }
321  i3 = 0;
322#endif
323#if MAX_PD>4
324    if (step >= pd_stop) break;
325    ++i4;
326  }
327  i4 = 0;
328#endif
329
330//if (q_index==0) printf("res: %g/%g\n", this_result, pd_norm);
331  // Remember the current result and the updated norm.
332  result[q_index] = this_result;
333  if (q_index == 0) result[nq] = pd_norm;
334}
Note: See TracBrowser for help on using the repository browser.