Line data    Source code

       1             : #include "crpropa/magneticField/turbulentField/SimpleGridTurbulence.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           3 : SimpleGridTurbulence::SimpleGridTurbulence(const SimpleTurbulenceSpectrum &spectrum,
      11             :                                            const GridProperties &gridProp,
      12           3 :                                            unsigned int seed)
      13           3 :     : GridTurbulence(spectrum, gridProp, seed) {
      14           6 :         initTurbulence(gridPtr, spectrum.getBrms(), spectrum.getLmin(),
      15           3 :                        spectrum.getLmax(), -spectrum.getSindex() - 2, seed);
      16           3 : }
      17             : 
      18           7 : void SimpleGridTurbulence::initTurbulence(ref_ptr<Grid3f> grid, double Brms,
      19             :                                           double lMin, double lMax,
      20             :                                           double alpha, int seed) {
      21             :         
      22             :         Vector3d spacing = grid->getSpacing();
      23             : 
      24          11 :         checkGridRequirements(grid, lMin, lMax);
      25             : 
      26             :         size_t n = grid->getNx(); // size of array
      27           4 :         size_t n2 = (size_t)floor(n / 2) +
      28             :                     1; // size array in z-direction in configuration space
      29             : 
      30             :         // arrays to hold the complex vector components of the B(k)-field
      31             :         fftwf_complex *Bkx, *Bky, *Bkz;
      32           4 :         Bkx = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
      33           4 :         Bky = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
      34           4 :         Bkz = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
      35             : 
      36           4 :         Random random;
      37           4 :         if (seed != 0)
      38           2 :                 random.seed(seed); // use given seed
      39             : 
      40             :         // calculate the n possible discrete wave numbers
      41           4 :         double K[n];
      42         260 :         for (size_t i = 0; i < n; i++)
      43         256 :                 K[i] = (double)i / n - i / (n / 2);
      44             : 
      45           4 :         double kMin = spacing.x / lMax;
      46           4 :         double kMax = spacing.x / lMin;
      47             :         Vector3f n0(1, 1, 1); // arbitrary vector to construct orthogonal base
      48             : 
      49         260 :         for (size_t ix = 0; ix < n; ix++) {
      50       16640 :                 for (size_t iy = 0; iy < n; iy++) {
      51      557056 :                         for (size_t iz = 0; iz < n2; iz++) {
      52             :         
      53             :                                 Vector3f ek, e1, e2;  // orthogonal base
      54             : 
      55      540672 :                                 size_t i = ix * n * n2 + iy * n2 + iz;
      56      540672 :                                 ek.setXYZ(K[ix], K[iy], K[iz]);
      57      540672 :                                 double k = ek.getR();
      58             : 
      59             :                                 // wave outside of turbulent range -> B(k) = 0
      60      540672 :                                 if ((k < kMin) || (k > kMax)) {
      61      264724 :                                         Bkx[i][0] = 0;
      62      264724 :                                         Bkx[i][1] = 0;
      63      264724 :                                         Bky[i][0] = 0;
      64      264724 :                                         Bky[i][1] = 0;
      65      264724 :                                         Bkz[i][0] = 0;
      66      264724 :                                         Bkz[i][1] = 0;
      67             :                                         continue;
      68             :                                 }
      69             : 
      70             :                                 // construct an orthogonal base ek, e1, e2
      71      275948 :                                 if (ek.isParallelTo(n0, float(1e-3))) {
      72             :                                         // ek parallel to (1,1,1)
      73             :                                         e1.setXYZ(-1., 1., 0);
      74             :                                         e2.setXYZ(1., 1., -2.);
      75             :                                 } else {
      76             :                                         // ek not parallel to (1,1,1)
      77             :                                         e1 = n0.cross(ek);
      78             :                                         e2 = ek.cross(e1);
      79             :                                 }
      80             :                                 e1 /= e1.getR();
      81             :                                 e2 /= e2.getR();
      82             : 
      83             :                                 // random orientation perpendicular to k
      84      275948 :                                 double theta = 2 * M_PI * random.rand();
      85      275948 :                                 Vector3f b = e1 * cos(theta) + e2 * sin(theta); // real b-field vector
      86             : 
      87             :                                 // normal distributed amplitude with mean = 0 and sigma =
      88             :                                 // k^alpha/2
      89      275948 :                                 b *= random.randNorm() * pow(k, alpha / 2);
      90             : 
      91             :                                 // uniform random phase
      92      275948 :                                 double phase = 2 * M_PI * random.rand();
      93      275948 :                                 double cosPhase = cos(phase); // real part
      94      275948 :                                 double sinPhase = sin(phase); // imaginary part
      95             : 
      96      275948 :                                 Bkx[i][0] = b.x * cosPhase;
      97      275948 :                                 Bkx[i][1] = b.x * sinPhase;
      98      275948 :                                 Bky[i][0] = b.y * cosPhase;
      99      275948 :                                 Bky[i][1] = b.y * sinPhase;
     100      275948 :                                 Bkz[i][0] = b.z * cosPhase;
     101      275948 :                                 Bkz[i][1] = b.z * sinPhase;
     102             :                         } // for iz
     103             :                 }     // for iy
     104             :         }         // for ix
     105             : 
     106           4 :         executeInverseFFTInplace(grid, Bkx, Bky, Bkz);
     107             : 
     108           4 :         fftwf_free(Bkx);
     109           4 :         fftwf_free(Bky);
     110           4 :         fftwf_free(Bkz);
     111             : 
     112          16 :         scaleGrid(grid, Brms / rmsFieldStrength(grid)); // normalize to Brms
     113          11 : }
     114             : 
     115             : } // namespace crpropa
     116             : 
     117             : #endif // CRPROPA_HAVE_FFTW3F