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