Line data Source code
1 : #include "crpropa/magneticField/turbulentField/HelicalGridTurbulence.h"
2 : #include "crpropa/GridTools.h"
3 : #include "crpropa/Random.h"
4 :
5 : #ifdef CRPROPA_HAVE_FFTW3F
6 : #include "fftw3.h"
7 :
8 : namespace crpropa {
9 :
10 0 : HelicalGridTurbulence::HelicalGridTurbulence(const SimpleTurbulenceSpectrum &spectrum,
11 : const GridProperties &gridProp,
12 0 : double H, unsigned int seed)
13 0 : : SimpleGridTurbulence(spectrum, gridProp, seed), H(H) {
14 0 : initTurbulence(gridPtr, spectrum.getBrms(), spectrum.getLmin(),
15 0 : spectrum.getLmax(), -spectrum.getSindex() - 2, seed, H);
16 0 : }
17 :
18 0 : void HelicalGridTurbulence::initTurbulence(ref_ptr<Grid3f> grid, double Brms,
19 : double lMin, double lMax,
20 : double alpha, int seed, double H) {
21 :
22 0 : checkGridRequirements(grid, lMin, lMax);
23 :
24 : Vector3d spacing = grid->getSpacing();
25 : size_t n = grid->getNx(); // size of array
26 0 : size_t n2 = (size_t)floor(n / 2) +
27 0 : 1; // size array in z-direction in configuration space
28 :
29 : // arrays to hold the complex vector components of the B(k)-field
30 : fftwf_complex *Bkx, *Bky, *Bkz;
31 0 : Bkx = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
32 0 : Bky = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
33 0 : Bkz = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
34 :
35 0 : Random random;
36 0 : if (seed != 0)
37 0 : random.seed(seed); // use given seed
38 :
39 : // calculate the n possible discrete wave numbers
40 0 : double K[n];
41 0 : for (size_t i = 0; i < n; i++)
42 0 : K[i] = (double)i / n - i / (n / 2);
43 :
44 : // only used if there is a helicity
45 : double Bktot, Bkplus, Bkminus, thetaplus, thetaminus;
46 :
47 0 : double kMin = spacing.x / lMax;
48 0 : double kMax = spacing.x / lMin;
49 : Vector3f b; // real b-field vector
50 : Vector3f ek, e1, e2; // orthogonal base
51 : Vector3f n0(1, 1, 1); // arbitrary vector to construct orthogonal base
52 :
53 0 : for (size_t ix = 0; ix < n; ix++) {
54 0 : for (size_t iy = 0; iy < n; iy++) {
55 0 : for (size_t iz = 0; iz < n2; iz++) {
56 :
57 0 : size_t i = ix * n * n2 + iy * n2 + iz;
58 0 : ek.setXYZ(K[ix], K[iy], K[iz]);
59 0 : double k = ek.getR();
60 :
61 : // wave outside of turbulent range -> B(k) = 0
62 0 : if ((k < kMin) || (k > kMax)) {
63 0 : Bkx[i][0] = 0;
64 0 : Bkx[i][1] = 0;
65 0 : Bky[i][0] = 0;
66 0 : Bky[i][1] = 0;
67 0 : Bkz[i][0] = 0;
68 0 : Bkz[i][1] = 0;
69 : continue;
70 : }
71 :
72 : // construct an orthogonal base ek, e1, e2
73 : // (for helical fields together with the real transform the
74 : // following convention must be used: e1(-k) = e1(k), e2(-k) = -
75 : // e2(k)
76 0 : if (ek.getAngleTo(n0) < 1e-3) { // ek parallel to (1,1,1)
77 : e1.setXYZ(-1, 1, 0);
78 : e2.setXYZ(1, 1, -2);
79 : } else { // ek not parallel to (1,1,1)
80 : e1 = n0.cross(ek);
81 : e2 = ek.cross(e1);
82 : }
83 : e1 /= e1.getR();
84 : e2 /= e2.getR();
85 :
86 :
87 0 : double Bkprefactor = mu0 / (4 * M_PI * pow(k, 3));
88 0 : Bktot = fabs(random.randNorm() * pow(k, alpha / 2));
89 0 : Bkplus = Bkprefactor * sqrt((1 + H) / 2) * Bktot;
90 0 : Bkminus = Bkprefactor * sqrt((1 - H) / 2) * Bktot;
91 0 : thetaplus = 2 * M_PI * random.rand();
92 0 : thetaminus = 2 * M_PI * random.rand();
93 0 : double ctp = cos(thetaplus);
94 0 : double stp = sin(thetaplus);
95 0 : double ctm = cos(thetaminus);
96 0 : double stm = sin(thetaminus);
97 :
98 0 : Bkx[i][0] = ((Bkplus * ctp + Bkminus * ctm) * e1.x +
99 0 : (-Bkplus * stp + Bkminus * stm) * e2.x) /
100 : sqrt(2);
101 0 : Bkx[i][1] = ((Bkplus * stp + Bkminus * stm) * e1.x +
102 0 : (Bkplus * ctp - Bkminus * ctm) * e2.x) /
103 : sqrt(2);
104 0 : Bky[i][0] = ((Bkplus * ctp + Bkminus * ctm) * e1.y +
105 0 : (-Bkplus * stp + Bkminus * stm) * e2.y) /
106 : sqrt(2);
107 0 : Bky[i][1] = ((Bkplus * stp + Bkminus * stm) * e1.y +
108 0 : (Bkplus * ctp - Bkminus * ctm) * e2.y) /
109 : sqrt(2);
110 0 : Bkz[i][0] = ((Bkplus * ctp + Bkminus * ctm) * e1.z +
111 0 : (-Bkplus * stp + Bkminus * stm) * e2.z) /
112 : sqrt(2);
113 0 : Bkz[i][1] = ((Bkplus * stp + Bkminus * stm) * e1.z +
114 0 : (Bkplus * ctp - Bkminus * ctm) * e2.z) /
115 : sqrt(2);
116 :
117 : Vector3f BkRe(Bkx[i][0], Bky[i][0], Bkz[i][0]);
118 : Vector3f BkIm(Bkx[i][1], Bky[i][1], Bkz[i][1]);
119 : } // for iz
120 : } // for iy
121 : } // for ix
122 :
123 0 : executeInverseFFTInplace(grid, Bkx, Bky, Bkz);
124 :
125 0 : scaleGrid(grid, Brms / rmsFieldStrength(grid)); // normalize to Brms
126 0 : }
127 :
128 : } // namespace crpropa
129 :
130 : #endif // CRPROPA_HAVE_FFTW3F
|
