Line data    Source code

       1             : #include "crpropa/magneticField/turbulentField/GridTurbulence.h"
       2             : #include "crpropa/GridTools.h"
       3             : #include "crpropa/Random.h"
       4             : 
       5             : #ifdef CRPROPA_HAVE_FFTW3F
       6             : 
       7             : namespace crpropa {
       8             : 
       9             : 
      10           6 : GridTurbulence::GridTurbulence(const TurbulenceSpectrum &spectrum,
      11             :                                const GridProperties &gridProp,
      12           6 :                                unsigned int seed)
      13           6 :     : TurbulentField(spectrum), seed(seed) {
      14           6 :         initGrid(gridProp);
      15           6 :         checkGridRequirements(gridPtr, spectrum.getLmin(), spectrum.getLmax());
      16           6 :         initTurbulence();
      17           6 : }
      18             : 
      19           6 : void GridTurbulence::initGrid(const GridProperties &p) {
      20           6 :         gridPtr = new Grid3f(p);
      21           6 : }
      22             : 
      23           4 : Vector3d GridTurbulence::getField(const Vector3d &pos) const {
      24           4 :         return gridPtr->interpolate(pos);
      25             : }
      26             : 
      27           1 : const ref_ptr<Grid3f> &GridTurbulence::getGrid() const { return gridPtr; }
      28             : 
      29           6 : void GridTurbulence::initTurbulence() {
      30             : 
      31             :         Vector3d spacing = gridPtr->getSpacing();
      32             :         size_t n = gridPtr->getNx(); // size of array
      33           6 :         size_t n2 = (size_t)floor(n / 2) +
      34             :                     1; // size array in z-direction in configuration space
      35             : 
      36             :         // arrays to hold the complex vector components of the B(k)-field
      37             :         fftwf_complex *Bkx, *Bky, *Bkz;
      38           6 :         Bkx = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
      39           6 :         Bky = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
      40           6 :         Bkz = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n2);
      41             : 
      42           6 :         Random random;
      43           6 :         if (seed != 0)
      44           4 :                 random.seed(seed); // use given seed
      45             : 
      46             :         // calculate the n possible discrete wave numbers
      47           6 :         double K[n];
      48         390 :         for (size_t i = 0; i < n; i++)
      49         384 :                 K[i] = ((double)i / n - i / (n / 2));
      50             : 
      51             :         // double kMin = 2*M_PI / lMax; // * 2 * spacing.x; // spacing.x / lMax;
      52             :         // double kMax = 2*M_PI / lMin; // * 2 * spacing.x; // spacing.x / lMin;
      53           6 :         double kMin = spacing.x / spectrum.getLmax();
      54           6 :         double kMax = spacing.x / spectrum.getLmin();
      55           6 :         auto lambda = 1 / spacing.x * 2 * M_PI;
      56             : 
      57             :         Vector3f n0(1, 1, 1); // arbitrary vector to construct orthogonal base
      58             : 
      59         390 :         for (size_t ix = 0; ix < n; ix++) {
      60       24960 :                 for (size_t iy = 0; iy < n; iy++) {
      61      835584 :                         for (size_t iz = 0; iz < n2; iz++) {
      62             :         
      63             :                                 Vector3f ek, e1, e2;  // orthogonal base
      64             : 
      65      811008 :                                 size_t i = ix * n * n2 + iy * n2 + iz;
      66      811008 :                                 ek.setXYZ(K[ix], K[iy], K[iz]);
      67      811008 :                                 double k = ek.getR();
      68             : 
      69             :                                 // wave outside of turbulent range -> B(k) = 0
      70      811008 :                                 if ((k < kMin) || (k > kMax)) {
      71      395939 :                                         Bkx[i][0] = 0;
      72      395939 :                                         Bkx[i][1] = 0;
      73      395939 :                                         Bky[i][0] = 0;
      74      395939 :                                         Bky[i][1] = 0;
      75      395939 :                                         Bkz[i][0] = 0;
      76      395939 :                                         Bkz[i][1] = 0;
      77             :                                         continue;
      78             :                                 }
      79             : 
      80             :                                 // construct an orthogonal base ek, e1, e2
      81      415069 :                                 if (ek.isParallelTo(n0, float(1e-3))) {
      82             :                                         // ek parallel to (1,1,1)
      83             :                                         e1.setXYZ(-1., 1., 0);
      84             :                                         e2.setXYZ(1., 1., -2.);
      85             :                                 } else {
      86             :                                         // ek not parallel to (1,1,1)
      87             :                                         e1 = n0.cross(ek);
      88             :                                         e2 = ek.cross(e1);
      89             :                                 }
      90             :                                 e1 /= e1.getR();
      91             :                                 e2 /= e2.getR();
      92             : 
      93             :                                 // random orientation perpendicular to k
      94      415069 :                                 double theta = 2 * M_PI * random.rand();
      95      415069 :                                 Vector3f b = e1 * std::cos(theta) + e2 * std::sin(theta); // real b-field vector
      96             : 
      97             :                                 // normal distributed amplitude with mean = 0
      98      415069 :                                 b *= std::sqrt(spectrum.energySpectrum(k*lambda));
      99             :                                 
     100             :                                 // uniform random phase
     101      415069 :                                 double phase = 2 * M_PI * random.rand();
     102      415069 :                                 double cosPhase = std::cos(phase); // real part
     103      415069 :                                 double sinPhase = std::sin(phase); // imaginary part
     104             : 
     105      415069 :                                 Bkx[i][0] = b.x * cosPhase;
     106      415069 :                                 Bkx[i][1] = b.x * sinPhase;
     107      415069 :                                 Bky[i][0] = b.y * cosPhase;
     108      415069 :                                 Bky[i][1] = b.y * sinPhase;
     109      415069 :                                 Bkz[i][0] = b.z * cosPhase;
     110      415069 :                                 Bkz[i][1] = b.z * sinPhase;
     111             :                         } // for iz
     112             :                 }     // for iy
     113             :         }         // for ix
     114             : 
     115           6 :         executeInverseFFTInplace(gridPtr, Bkx, Bky, Bkz);
     116             : 
     117           6 :         fftwf_free(Bkx);
     118           6 :         fftwf_free(Bky);
     119           6 :         fftwf_free(Bkz);
     120             : 
     121          24 :         scaleGrid(gridPtr, spectrum.getBrms() /
     122          12 :                                rmsFieldStrength(gridPtr)); // normalize to Brms
     123          12 : }
     124             : 
     125             : // Check the grid properties before the FFT procedure
     126          13 : void GridTurbulence::checkGridRequirements(ref_ptr<Grid3f> grid, double lMin,
     127             :                                            double lMax) {
     128             :         size_t Nx = grid->getNx();
     129             :         size_t Ny = grid->getNy();
     130             :         size_t Nz = grid->getNz();
     131             :         Vector3d spacing = grid->getSpacing();
     132             : 
     133          13 :         if ((Nx != Ny) or (Ny != Nz))
     134           0 :                 throw std::runtime_error("turbulentField: only cubic grid supported");
     135          13 :         if ((spacing.x != spacing.y) or (spacing.y != spacing.z))
     136           0 :                 throw std::runtime_error("turbulentField: only equal spacing suported");
     137          13 :         if (lMin < 2 * spacing.x)
     138           1 :                 throw std::runtime_error("turbulentField: lMin < 2 * spacing");
     139          12 :         if (lMax > Nx * spacing.x) // before was (lMax > Nx * spacing.x / 2), why?
     140           1 :                 throw std::runtime_error("turbulentField: lMax > size");
     141          11 :         if (lMax < lMin) 
     142           1 :                 throw std::runtime_error("lMax < lMin");
     143          10 : }
     144             : 
     145             : // Execute inverse discrete FFT in-place for a 3D grid, from complex to real
     146             : // space
     147          10 : void GridTurbulence::executeInverseFFTInplace(ref_ptr<Grid3f> grid,
     148             :                                               fftwf_complex *Bkx,
     149             :                                               fftwf_complex *Bky,
     150             :                                               fftwf_complex *Bkz) {
     151             : 
     152             :         size_t n = grid->getNx(); // size of array
     153          10 :         size_t n2 = (size_t)floor(n / 2) +
     154             :                     1; // size array in z-direction in configuration space
     155             : 
     156             :         // in-place, complex to real, inverse Fourier transformation on each
     157             :         // component note that the last elements of B(x) are unused now
     158             :         float *Bx = (float *)Bkx;
     159          10 :         fftwf_plan plan_x = fftwf_plan_dft_c2r_3d(n, n, n, Bkx, Bx, FFTW_ESTIMATE);
     160          10 :         fftwf_execute(plan_x);
     161          10 :         fftwf_destroy_plan(plan_x);
     162             : 
     163             :         float *By = (float *)Bky;
     164          10 :         fftwf_plan plan_y = fftwf_plan_dft_c2r_3d(n, n, n, Bky, By, FFTW_ESTIMATE);
     165          10 :         fftwf_execute(plan_y);
     166          10 :         fftwf_destroy_plan(plan_y);
     167             : 
     168             :         float *Bz = (float *)Bkz;
     169          10 :         fftwf_plan plan_z = fftwf_plan_dft_c2r_3d(n, n, n, Bkz, Bz, FFTW_ESTIMATE);
     170          10 :         fftwf_execute(plan_z);
     171          10 :         fftwf_destroy_plan(plan_z);
     172             : 
     173             :         // save to grid
     174         650 :         for (size_t ix = 0; ix < n; ix++) {
     175       41600 :                 for (size_t iy = 0; iy < n; iy++) {
     176     2662400 :                         for (size_t iz = 0; iz < n; iz++) {
     177     2621440 :                                 size_t i = ix * n * 2 * n2 + iy * 2 * n2 + iz;
     178             :                                 Vector3f &b = grid->get(ix, iy, iz);
     179     2621440 :                                 b.x = Bx[i];
     180     2621440 :                                 b.y = By[i];
     181     2621440 :                                 b.z = Bz[i];
     182             :                         }
     183             :                 }
     184             :         }
     185          10 : }
     186             : 
     187           0 : Vector3f GridTurbulence::getMeanFieldVector() const {
     188           0 :         return meanFieldVector(gridPtr);
     189             : }
     190             : 
     191           0 : double GridTurbulence::getMeanFieldStrength() const {
     192           0 :         return meanFieldStrength(gridPtr);
     193             : }
     194             : 
     195           0 : double GridTurbulence::getRmsFieldStrength() const {
     196           0 :         return rmsFieldStrength(gridPtr);
     197             : }
     198             : 
     199           0 : std::array<float, 3> GridTurbulence::getRmsFieldStrengthPerAxis() const {
     200           0 :         return rmsFieldStrengthPerAxis(gridPtr);
     201             : }
     202             :         
     203           0 : std::vector<std::pair<int, float>> GridTurbulence::getPowerSpectrum() const {
     204           0 :         return gridPowerSpectrum(gridPtr);
     205             : }
     206             : 
     207           0 : void GridTurbulence::dumpToFile(std::string filename) const {
     208           0 :         dumpGrid(gridPtr, filename);
     209           0 : }
     210             : 
     211             : } // namespace crpropa
     212             : 
     213             : #endif // CRPROPA_HAVE_FFTW3F