Line data    Source code

       1             : /**********************************************************************/
       2             : /*                     DO NOT MODIFY THIS HEADER                      */
       3             : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
       4             : /*                                                                    */
       5             : /*            Copyright 2017 Battelle Energy Alliance, LLC            */
       6             : /*                        ALL RIGHTS RESERVED                         */
       7             : /**********************************************************************/
       8             : #ifdef FFTW3_ENABLED
       9             : 
      10             : #include "FourierPowerSpectrum.h"
      11             : #include "FourierTransform.h"
      12             : #include "MyTRIMMesh.h"
      13             : 
      14             : #include "libmesh/utility.h"
      15             : 
      16             : #include <functional>
      17             : 
      18             : registerMooseObject("MagpieApp", FourierPowerSpectrum);
      19             : 
      20             : InputParameters
      21           8 : FourierPowerSpectrum::validParams()
      22             : {
      23           8 :   InputParameters params = GeneralVectorPostprocessor::validParams();
      24           8 :   params.addClassDescription("Compute the power spectrum of a given fourier transform. The "
      25             :                              "resulting frequency is in reciprocal length units.");
      26          16 :   params.addRequiredParam<UserObjectName>(
      27             :       "fourier_transform", "FourierTransform user object to compute the power spectrum of");
      28           8 :   return params;
      29           0 : }
      30             : 
      31           4 : FourierPowerSpectrum::FourierPowerSpectrum(const InputParameters & parameters)
      32             :   : GeneralVectorPostprocessor(parameters),
      33           4 :     _fourier_transform(getUserObject<FourierTransform>("fourier_transform")),
      34           4 :     _dim(_fourier_transform.dimension()),
      35           4 :     _grid(_fourier_transform.getGrid()),
      36           4 :     _box_size(_fourier_transform.getBoxSize()),
      37           4 :     _buffer(_fourier_transform.getBuffer()),
      38           4 :     _spectrum_size(0),
      39           4 :     _max_frequency(0),
      40           4 :     _power_spectrum(declareVector("power_spectrum")),
      41           8 :     _frequency(declareVector("frequency"))
      42             : {
      43             :   // estimate spectrum bins
      44          12 :   for (unsigned int i = 0; i < _dim; ++i)
      45             :   {
      46           8 :     _max_frequency += Utility::pow<2>(_grid[i] * 0.5 / _box_size(i));
      47           8 :     _spectrum_size += Utility::pow<2>(_grid[i]);
      48             :   }
      49             : 
      50           4 :   _max_frequency = std::sqrt(_max_frequency);
      51           4 :   _spectrum_size = std::sqrt(_spectrum_size) + 1;
      52             : 
      53             :   // resize power spectrum vector
      54           4 :   _power_spectrum.resize(_spectrum_size);
      55           4 :   _frequency.resize(_spectrum_size);
      56           4 :   _n_samples.resize(_spectrum_size);
      57             : 
      58             :   // set frequencies
      59         152 :   for (auto i = beginIndex(_frequency); i < _spectrum_size; ++i)
      60         148 :     _frequency[i] = (i * _max_frequency) / (_spectrum_size - 1);
      61           4 : }
      62             : 
      63             : void
      64          63 : FourierPowerSpectrum::computePowerSpectrum(std::vector<int> & c,
      65             :                                            std::vector<Real> & F,
      66             :                                            std::size_t i)
      67             : {
      68          63 :   if (i == 0)
      69             :   {
      70             :     Real sum_high = 0.0;
      71         120 :     for (unsigned int j = 1; j < _dim; ++j)
      72          60 :       sum_high += F[j];
      73             : 
      74        1860 :     for (c[0] = 0; c[0] < _grid[0]; ++c[0])
      75             :     {
      76             :       // do stuff at lowest dimension
      77             :       const Real sum =
      78        1800 :           sum_high + Utility::pow<2>((c[0] * 2 > _grid[0] ? _grid[0] - c[0] : c[0]) / _box_size(0));
      79             : 
      80             :       // recalculate index (row major)
      81             :       int index = 0;
      82        5400 :       for (unsigned int j = 0; j < _dim; ++j)
      83        3600 :         index = index * _grid[j] + c[j];
      84             : 
      85             :       // find bin and add to spectrum
      86        1800 :       const std::size_t bin = std::floor((std::sqrt(sum) * (_spectrum_size - 1)) / _max_frequency);
      87             :       mooseAssert(bin < _power_spectrum.size(), "Bin out of bounds in power spectrum");
      88        1800 :       _power_spectrum[bin] += Utility::pow<2>(_buffer[index]);
      89        1800 :       ++_n_samples[bin];
      90             :     }
      91             :   }
      92             :   else
      93             :     // iterate over lower dimensions
      94          63 :     for (c[i] = 0; c[i] < _grid[i]; ++c[i])
      95             :     {
      96          60 :       F[i] = Utility::pow<2>((c[i] * 2 > _grid[i] ? _grid[i] - c[i] : c[i]) / _box_size(i));
      97          60 :       computePowerSpectrum(c, F, i - 1);
      98             :     }
      99          63 : }
     100             : 
     101             : void
     102           4 : FourierPowerSpectrum::execute()
     103             : {
     104             :   // clear spectrum
     105           4 :   std::fill(_power_spectrum.begin(), _power_spectrum.end(), 0.0);
     106             : 
     107             :   // for now only compute spectrum on processor 0
     108           4 :   if (processor_id() == 0)
     109             :   {
     110             :     // clear sample count
     111             :     std::fill(_n_samples.begin(), _n_samples.end(), 0);
     112             : 
     113             :     // compute the power spectrum
     114           3 :     std::vector<int> c(_dim, 0);
     115           3 :     std::vector<Real> F(_dim);
     116           3 :     computePowerSpectrum(c, F, _dim - 1);
     117             : 
     118             :     // normalize by sample count
     119             :     mooseAssert(_power_spectrum.size() == _n_samples.size(),
     120             :                 "Size mismatch between spectrum and sample number vectors");
     121         114 :     for (auto i = beginIndex(_power_spectrum); i < _power_spectrum.size(); ++i)
     122         111 :       if (_n_samples[i] > 0)
     123         108 :         _power_spectrum[i] /= _n_samples[i];
     124             :   }
     125           4 : }
     126             : 
     127             : void
     128           4 : FourierPowerSpectrum::finalize()
     129             : {
     130             :   // broadcast result to all cores (lazily using gatherSum - all other ranks have zero vectors)
     131           4 :   gatherSum(_power_spectrum);
     132           4 : }
     133             : 
     134             : #endif // FFTW3_ENABLED