Line data Source code
1 : // This file contains an implementation of a vectorized cosine, which
2 : // is based in part on the implementations in the library "SLEEF" by
3 : // Naoki Shibata. SLEEF was used under the Boost Software License,
4 : // Version 1.0. The original source file contained the following
5 : // copyright notice:
6 : //
7 : // // Copyright Naoki Shibata 2010 - 2018.
8 : // // Distributed under the Boost Software License, Version 1.0.
9 : // // (See accompanying file LICENSE.txt or copy at
10 : // // http://www.boost.org/LICENSE_1_0.txt)
11 : //
12 : // SLEEF was used under the following license, which is not necessarily the
13 : // license that applies to this file:
14 : //
15 : // Boost Software License - Version 1.0 - August 17th, 2003
16 : //
17 : // Permission is hereby granted, free of charge, to any person or
18 : // organization obtaining a copy of the software and accompanying
19 : // documentation covered by this license (the "Software") to use,
20 : // reproduce, display, distribute, execute, and transmit the Software,
21 : // and to prepare derivative works of the Software, and to permit
22 : // third-parties to whom the Software is furnished to do so, all subject
23 : // to the following:
24 : //
25 : // The copyright notices in the Software and this entire statement,
26 : // including the above license grant, this restriction and the following
27 : // disclaimer, must be included in all copies of the Software, in whole
28 : // or in part, and all derivative works of the Software, unless such
29 : // copies or derivative works are solely in the form of
30 : // machine-executable object code generated by a source language
31 : // processor.
32 : //
33 : // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
34 : // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
35 : // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND
36 : // NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR ANYONE
37 : // DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR OTHER
38 : // LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT
39 : // OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
40 : // THE SOFTWARE.
41 :
42 : #include "crpropa/magneticField/turbulentField/PlaneWaveTurbulence.h"
43 : #include "crpropa/GridTools.h"
44 : #include "crpropa/Random.h"
45 : #include "crpropa/Units.h"
46 :
47 : #include "kiss/logger.h"
48 :
49 : #include <iostream>
50 :
51 : #if defined(FAST_WAVES)
52 : #if defined(__SSE__) && defined(__SSE2__) && defined(__SSE3__) && defined(__SSE4_1__) && defined(__SSE4_2__) && defined(__AVX__)
53 : #define ENABLE_FAST_WAVES
54 : #else
55 : #error "FAST_WAVES is enabled, but it appears that not all required SIMD extensions are enabled in the compiler. Without these extensions, the FAST_WAVES optimization cannot be used. Please make sure that the SIMD_EXTENSIONS option in cmake matches the capabilities of your target CPU, or (if your target CPU does not support the required extensions), disable the FAST_WAVES flag in cmake."
56 : #endif
57 : #endif
58 :
59 :
60 : #ifdef ENABLE_FAST_WAVES
61 : #include <immintrin.h>
62 : #endif
63 :
64 : namespace crpropa {
65 : #ifdef ENABLE_FAST_WAVES
66 : // see
67 : // https://stackoverflow.com/questions/49941645/get-sum-of-values-stored-in-m256d-with-sse-avx
68 : double hsum_double_avx(__m256d v) {
69 : __m128d vlow = _mm256_castpd256_pd128(v);
70 : __m128d vhigh = _mm256_extractf128_pd(v, 1); // high 128
71 : vlow = _mm_add_pd(vlow, vhigh); // reduce down to 128
72 :
73 : __m128d high64 = _mm_unpackhi_pd(vlow, vlow);
74 : return _mm_cvtsd_f64(_mm_add_sd(vlow, high64)); // reduce to scalar
75 : }
76 : #endif // defined(ENABLE_FAST_WAVES)
77 :
78 2 : PlaneWaveTurbulence::PlaneWaveTurbulence(const TurbulenceSpectrum &spectrum,
79 2 : int Nm, int seed)
80 2 : : TurbulentField(spectrum), Nm(Nm) {
81 :
82 : #ifdef ENABLE_FAST_WAVES
83 : KISS_LOG_INFO << "PlaneWaveTurbulence: Using SIMD TD13 implementation"
84 : << std::endl;
85 :
86 : // There used to be a cpuid check here, to see if the cpu running
87 : // this code would support SIMD (SSE + AVX). However, the library providing
88 : // the relevant function is no longer being used, and doing this manually
89 : // might be a bit too much work.
90 : #endif
91 :
92 2 : if (Nm <= 1) {
93 0 : throw std::runtime_error(
94 : "PlaneWaveTurbulence: Nm <= 1. Specify at least two wavemodes in "
95 0 : "order to generate the k distribution properly.");
96 : }
97 :
98 2 : Random random;
99 2 : if (seed != 0)
100 2 : random.seed(seed);
101 :
102 2 : double kmax = 2 * M_PI / spectrum.getLmin();
103 2 : double kmin = 2 * M_PI / spectrum.getLmax();
104 :
105 4 : xi = std::vector<Vector3d>(Nm, Vector3d(0.));
106 4 : kappa = std::vector<Vector3d>(Nm, Vector3d(0.));
107 2 : phi = std::vector<double>(Nm, 0.);
108 2 : costheta = std::vector<double>(Nm, 0.);
109 2 : beta = std::vector<double>(Nm, 0.);
110 2 : Ak = std::vector<double>(Nm, 0.);
111 2 : k = std::vector<double>(Nm, 0.);
112 :
113 2 : double delta = log10(kmax / kmin);
114 22 : for (int i = 0; i < Nm; i++) {
115 20 : k[i] = pow(10, log10(kmin) + ((double)i) / ((double)(Nm - 1)) * delta);
116 : }
117 :
118 : // * compute Ak *
119 :
120 : double delta_k0 =
121 2 : (k[1] - k[0]) / k[1]; // multiply this by k[i] to get delta_k[i]
122 : // Note: this is probably unnecessary since it's just a factor
123 : // and will get normalized out anyways. It's not like this is
124 : // performance-sensitive code though, and I don't want to change
125 : // anything numerical now.
126 :
127 : // For this loop, the Ak array actually contains Gk*delta_k (ie
128 : // non-normalized Ak^2). Normalization happens in a second loop,
129 : // once the total is known.
130 : double Ak2_sum = 0; // sum of Ak^2 over all k
131 22 : for (int i = 0; i < Nm; i++) {
132 20 : double k = this->k[i];
133 20 : double kHat = k * spectrum.getLbendover();
134 20 : double Gk = spectrum.energySpectrum(k) * (1 + kHat * kHat); // correct different implementation in TD 13 (eq. 5, missing + 1 in the denuminators exponent)
135 20 : Ak[i] = Gk * delta_k0 * k;
136 20 : Ak2_sum += Ak[i];
137 :
138 : // phi, costheta, and sintheta are for drawing vectors with
139 : // uniform distribution on the unit sphere.
140 : // This is similar to Random::randVector(): their t is our phi,
141 : // z is costheta, and r is sintheta. Our kappa is equivalent to
142 : // the return value of randVector(); however, TD13 then reuse
143 : // these values to generate a random vector perpendicular to kappa.
144 20 : double phi = random.randUniform(-M_PI, M_PI);
145 20 : double costheta = random.randUniform(-1., 1.);
146 20 : double sintheta = sqrt(1 - costheta * costheta);
147 :
148 20 : double alpha = random.randUniform(0, 2 * M_PI);
149 20 : double beta = random.randUniform(0, 2 * M_PI);
150 :
151 : Vector3d kappa =
152 20 : Vector3d(sintheta * cos(phi), sintheta * sin(phi), costheta);
153 :
154 : // NOTE: all other variable names match the ones from the TD13 paper.
155 : // However, our xi is actually their psi, and their xi isn't used at
156 : // all. (Though both can be used for the polarization vector, according
157 : // to the paper.) The reason for this discrepancy is that this code
158 : // used to be based on the original GJ99 paper, which provided only a
159 : // xi vector, and this xi happens to be almost the same as TD13's psi.
160 : Vector3d xi =
161 20 : Vector3d(costheta * cos(phi) * cos(alpha) + sin(phi) * sin(alpha),
162 20 : costheta * sin(phi) * cos(alpha) - cos(phi) * sin(alpha),
163 20 : -sintheta * cos(alpha));
164 :
165 : this->xi[i] = xi;
166 : this->kappa[i] = kappa;
167 20 : this->phi[i] = phi;
168 20 : this->costheta[i] = costheta;
169 20 : this->beta[i] = beta;
170 : }
171 :
172 : // Only in this loop are the actual Ak computed and stored.
173 : // This two-step process is necessary in order to normalize the values
174 : // properly.
175 22 : for (int i = 0; i < Nm; i++) {
176 20 : Ak[i] = sqrt(2 * Ak[i] / Ak2_sum) * spectrum.getBrms();
177 : }
178 :
179 : #ifdef ENABLE_FAST_WAVES
180 : // * copy data into AVX-compatible arrays *
181 : //
182 : // AVX requires all data to be aligned to 256 bit, or 32 bytes, which is the
183 : // same as 4 double precision floating point numbers. Since support for
184 : // alignments this big seems to be somewhat tentative in C++ allocators,
185 : // we're aligning them manually by allocating a normal double array, and
186 : // then computing the offset to the first value with the correct alignment.
187 : // This is a little bit of work, so instead of doing it separately for each
188 : // of the individual data arrays, we're doing it once for one big array that
189 : // all of the component arrays get packed into.
190 : //
191 : // The other thing to keep in mind is that AVX always reads in units of 256
192 : // bits, or 4 doubles. This means that our number of wavemodes must be
193 : // divisible by 4. If it isn't, we simply pad it out with zeros. Since the
194 : // final step of the computation of each wavemode is multiplication by the
195 : // amplitude, which will be set to 0, these padding wavemodes won't affect
196 : // the result.
197 :
198 : avx_Nm = ((Nm + 4 - 1) / 4) * 4; // round up to next larger multiple of 4:
199 : // align is 256 = 4 * sizeof(double) bit
200 : avx_data = std::vector<double>(itotal * avx_Nm + 3, 0.);
201 :
202 : // get the first 256-bit aligned element
203 : size_t size = avx_data.size() * sizeof(double);
204 : void *pointer = avx_data.data();
205 : align_offset =
206 : (double *)std::align(32, 32, pointer, size) - avx_data.data();
207 :
208 : // copy into the AVX arrays
209 : for (int i = 0; i < Nm; i++) {
210 : avx_data[i + align_offset + avx_Nm * iAxi0] = Ak[i] * xi[i].x;
211 : avx_data[i + align_offset + avx_Nm * iAxi1] = Ak[i] * xi[i].y;
212 : avx_data[i + align_offset + avx_Nm * iAxi2] = Ak[i] * xi[i].z;
213 :
214 : // The cosine implementation computes cos(pi*x), so we'll divide out the
215 : // pi here.
216 : avx_data[i + align_offset + avx_Nm * ikkappa0] =
217 : k[i] / M_PI * kappa[i].x;
218 : avx_data[i + align_offset + avx_Nm * ikkappa1] =
219 : k[i] / M_PI * kappa[i].y;
220 : avx_data[i + align_offset + avx_Nm * ikkappa2] =
221 : k[i] / M_PI * kappa[i].z;
222 :
223 : // We also need to divide beta by pi, since that goes into the argument
224 : // of the cosine as well.
225 : avx_data[i + align_offset + avx_Nm * ibeta] = beta[i] / M_PI;
226 : }
227 : #endif // ENABLE_FAST_WAVES
228 2 : }
229 :
230 45 : Vector3d PlaneWaveTurbulence::getField(const Vector3d &pos) const {
231 :
232 : #ifndef ENABLE_FAST_WAVES
233 : Vector3d B(0.);
234 495 : for (int i = 0; i < Nm; i++) {
235 450 : double z_ = pos.dot(kappa[i]);
236 450 : B += xi[i] * Ak[i] * cos(k[i] * z_ + beta[i]);
237 : }
238 45 : return B;
239 :
240 : #else // ENABLE_FAST_WAVES
241 :
242 : // Initialize accumulators
243 : //
244 : // There is one accumulator per component of the result vector.
245 : // Note that each accumulator contains four numbers. At the end of
246 : // the loop, each of these numbers will contain the sum of every
247 : // fourth wavemode, starting at a different offset. In the end, each
248 : // of the accumulator's numbers are added together (using
249 : // hsum_double_avx), resulting in the total sum for that component.
250 :
251 : __m256d acc0 = _mm256_setzero_pd();
252 : __m256d acc1 = _mm256_setzero_pd();
253 : __m256d acc2 = _mm256_setzero_pd();
254 :
255 : // broadcast position into AVX registers
256 : __m256d pos0 = _mm256_set1_pd(pos.x);
257 : __m256d pos1 = _mm256_set1_pd(pos.y);
258 : __m256d pos2 = _mm256_set1_pd(pos.z);
259 :
260 : for (int i = 0; i < avx_Nm; i += 4) {
261 :
262 : // Load data from memory into AVX registers:
263 : // - the three components of the vector A * xi
264 : __m256d Axi0 =
265 : _mm256_load_pd(avx_data.data() + i + align_offset + avx_Nm * iAxi0);
266 : __m256d Axi1 =
267 : _mm256_load_pd(avx_data.data() + i + align_offset + avx_Nm * iAxi1);
268 : __m256d Axi2 =
269 : _mm256_load_pd(avx_data.data() + i + align_offset + avx_Nm * iAxi2);
270 :
271 : // - the three components of the vector k * kappa
272 : __m256d kkappa0 = _mm256_load_pd(avx_data.data() + i + align_offset +
273 : avx_Nm * ikkappa0);
274 : __m256d kkappa1 = _mm256_load_pd(avx_data.data() + i + align_offset +
275 : avx_Nm * ikkappa1);
276 : __m256d kkappa2 = _mm256_load_pd(avx_data.data() + i + align_offset +
277 : avx_Nm * ikkappa2);
278 :
279 : // - the phase beta.
280 : __m256d beta =
281 : _mm256_load_pd(avx_data.data() + i + align_offset + avx_Nm * ibeta);
282 :
283 : // Then, do the computation.
284 :
285 : // This is the scalar product between k*kappa and pos:
286 : __m256d z = _mm256_add_pd(_mm256_mul_pd(pos0, kkappa0),
287 : _mm256_add_pd(_mm256_mul_pd(pos1, kkappa1),
288 : _mm256_mul_pd(pos2, kkappa2)));
289 :
290 : // Here, the phase is added on. This is the argument of the cosine.
291 : __m256d cos_arg = _mm256_add_pd(z, beta);
292 :
293 : // ********
294 : // * Computing the cosine
295 : // * Part 1: Argument reduction
296 : //
297 : // To understand the computation of the cosine, first note that the
298 : // cosine is periodic and we thus only need to model its behavior
299 : // between 0 and 2*pi to be able compute the function anywhere. In
300 : // fact, by mirroring the function along the x and y axes, even the
301 : // range between 0 and pi/2 is sufficient for this purpose. In this
302 : // range, the cosine can be efficiently evaluated with high precision
303 : // by using a polynomial approximation. Thus, to compute the cosine,
304 : // the input value is first reduced so that it lies within this range.
305 : // Then, the polynomial approximation is evaluated. Finally, if
306 : // necessary, the sign of the result is flipped (mirroring the function
307 : // along the x axis).
308 : //
309 : // The actual computation is slightly more involved. First, argument
310 : // reduction can be simplified drastically by computing cos(pi*x),
311 : // such that the values are reduced to the range [0, 0.5) instead of
312 : // [0, pi/2). Since the cosine is even (independent of the sign), we
313 : // can first reduce values to [-0.5, 0.5) – that is, a simple rounding
314 : // operation – and then neutralize the sign. In fact, precisely because
315 : // the cosine is even, all terms of the polynomial are powers of x^2,
316 : // so the value of x^2 (computed as x*x) forms the basis for the
317 : // polynomial approximation. If I understand things correctly, then (in
318 : // IEEE-754 floating point) x*x and (-x)*(-x) will always result in the
319 : // exact same value, which means that any error bound over [0, 0.5)
320 : // automatically applies to (-0.5, 0] as well.
321 :
322 : // First, compute round(x), and store it in q. If this value is odd,
323 : // we're looking at the negative half-wave of the cosine, and thus
324 : // will have to invert the sign of the result.
325 : __m256d q = _mm256_round_pd(
326 : cos_arg, (_MM_FROUND_TO_NEAREST_INT | _MM_FROUND_NO_EXC));
327 :
328 : // Since we're computing cos(pi*x), round(x) always yields the center of
329 : // a half-wave (where cos(pi*x) achieves an extremum). This point
330 : // logically corresponds to x=0. Therefore, we subtract this center from
331 : // the actual input argument to find the corresponding point on the
332 : // half-wave that is centered around zero.
333 : __m256d s = _mm256_sub_pd(cos_arg, q);
334 :
335 : // We now want to check whether q (the index of our half-wave) is even
336 : // or odd, since all of the odd-numbered half-waves are negative, so
337 : // we'll have to flip the final result. On an int, this is as simple as
338 : // checking the 0th bit. Idea: manipulate the double in such a way that
339 : // we can do this. So, we add 2^52, such that the last digit of the
340 : // mantissa is actually in the ones' position. Since q may be negative,
341 : // we'll also add 2^51 to make sure it's positive. Note that 2^51 is
342 : // even and thus leaves evenness invariant, which is the only thing we
343 : // care about here.
344 : //
345 : // This is based on the int extraction process described here:
346 : // https://stackoverflow.com/questions/41144668/how-to-efficiently-perform-double-int64-conversions-with-sse-avx/41223013
347 : //
348 : // We assume -2^51 <= q < 2^51 for this, which is unproblematic, as
349 : // double precision has decayed far enough at that point that the
350 : // usefulness of the cosine becomes limited.
351 : //
352 : // Explanation: The mantissa of a double-precision float has 52 bits
353 : // (excluding the implicit first bit, which is always one). If |q| >
354 : // 2^51, this implicit first bit has a place value of at least 2^51,
355 : // while the first stored bit of the mantissa has a place value of at
356 : // least 2^50. This means that the LSB of the mantissa has a place value
357 : // of at least 2^(-1), or 0.5. For a cos(pi*x), this corresponds to a
358 : // quarter of a cycle (pi/2), so at this point the precision of the
359 : // input argument is so low that going from one representable number to
360 : // the next causes the result to jump by +/-1.
361 :
362 : q = _mm256_add_pd(q, _mm256_set1_pd(0x0018000000000000));
363 :
364 : // Unfortunately, integer comparisons were only introduced in AVX2, so
365 : // we'll have to make do with a floating point comparison to check
366 : // whether the last bit is set. However, masking out all but the last
367 : // bit will result in a denormal float, which may either result in
368 : // performance problems or just be rounded down to zero, neither of
369 : // which is what we want here. To fix this, we'll mask in not only bit
370 : // 0, but also the exponent (and sign, but that doesn't matter) of q.
371 : // Luckily, the exponent of q is guaranteed to have the fixed value of
372 : // 1075 (corresponding to 2^52) after our addition.
373 :
374 : __m256d invert = _mm256_and_pd(
375 : q, _mm256_castsi256_pd(_mm256_set1_epi64x(0xfff0000000000001)));
376 :
377 : // If we did have a one in bit 0, our result will be equal to 2^52 + 1.
378 : invert = _mm256_cmp_pd(
379 : invert, _mm256_castsi256_pd(_mm256_set1_epi64x(0x4330000000000001)),
380 : _CMP_EQ_OQ);
381 :
382 : // Now we know whether to flip the sign of the result. However, remember
383 : // that we're working on multiple values at a time, so an if statement
384 : // won't be of much use here (plus it might perform badly). Instead,
385 : // we'll make use of the fact that the result of the comparison is all
386 : // ones if the comparison was true (i.e. q is odd and we need to flip
387 : // the result), and all zeroes otherwise. If we now mask out all bits
388 : // except the sign bit, we get something that, when xor'ed into our
389 : // final result, will flip the sign exactly when q is odd.
390 : invert = _mm256_and_pd(invert, _mm256_set1_pd(-0.0));
391 : // (Note that the binary representation of -0.0 is all 0 bits, except
392 : // for the sign bit, which is set to 1.)
393 :
394 : // TODO: clamp floats between 0 and 1? This would ensure that we never
395 : // see inf's, but maybe we want that, so that things dont just fail
396 : // silently...
397 :
398 : // * end of argument reduction
399 : // *******
400 :
401 : // ******
402 : // * Evaluate the cosine using a polynomial approximation for the zeroth
403 : // half-wave.
404 : // * The coefficients for this were generated using sleefs gencoef.c.
405 : // * These coefficients are probably far from optimal; however, they
406 : // should be sufficient for this case.
407 : s = _mm256_mul_pd(s, s);
408 :
409 : __m256d u = _mm256_set1_pd(+0.2211852080653743946e+0);
410 :
411 : u = _mm256_add_pd(_mm256_mul_pd(u, s),
412 : _mm256_set1_pd(-0.1332560668688523853e+1));
413 : u = _mm256_add_pd(_mm256_mul_pd(u, s),
414 : _mm256_set1_pd(+0.4058509506474178075e+1));
415 : u = _mm256_add_pd(_mm256_mul_pd(u, s),
416 : _mm256_set1_pd(-0.4934797516664651162e+1));
417 : u = _mm256_add_pd(_mm256_mul_pd(u, s), _mm256_set1_pd(1.));
418 :
419 : // Then, flip the sign of each double for which invert is not zero.
420 : // Since invert has only zero bits except for a possible one in bit 63,
421 : // we can xor it onto our result to selectively invert the 63rd (sign)
422 : // bit in each double where invert is set.
423 : u = _mm256_xor_pd(u, invert);
424 :
425 : // * end computation of cosine
426 : // **********
427 :
428 : // Finally, Ak*xi is multiplied on. Since this is a vector, the
429 : // multiplication needs to be done for each of the three
430 : // components, so it happens separately.
431 : acc0 = _mm256_add_pd(_mm256_mul_pd(u, Axi0), acc0);
432 : acc1 = _mm256_add_pd(_mm256_mul_pd(u, Axi1), acc1);
433 : acc2 = _mm256_add_pd(_mm256_mul_pd(u, Axi2), acc2);
434 : }
435 :
436 : return Vector3d(hsum_double_avx(acc0), hsum_double_avx(acc1),
437 : hsum_double_avx(acc2));
438 : #endif // ENABLE_FAST_WAVES
439 : }
440 :
441 : } // namespace crpropa
|
