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             :                     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