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 "FourierTransform.h"
      11             : #include "MyTRIMMesh.h"
      12             : 
      13             : #include "libmesh/utility.h"
      14             : 
      15             : #include <functional>
      16             : 
      17             : registerMooseObject("MagpieApp", FourierTransform);
      18             : 
      19             : InputParameters
      20          27 : FourierTransform::validParams()
      21             : {
      22          27 :   InputParameters params = ElementUserObject::validParams();
      23          27 :   params.addClassDescription("Compute the Fourier transform of a given variable field.");
      24          54 :   params.addCoupledVar("variable", "Variable field to compute the transform of");
      25          54 :   params.addRangeCheckedParam<std::vector<int>>(
      26             :       "grid",
      27             :       "grid > 0",
      28             :       "Number of grid cells in each dimension to compute "
      29             :       "the FFT on (can be omitted when using MyTRIMMesh)");
      30          54 :   params.addParam<unsigned int>(
      31          54 :       "max_h_level", 0, "Further grid refinement to apply when using MyTRIMMesh with adaptivity");
      32          27 :   return params;
      33           0 : }
      34             : 
      35          15 : FourierTransform::FourierTransform(const InputParameters & parameters)
      36             :   : ElementUserObject(parameters),
      37          15 :     _var(coupledValue("variable")),
      38          15 :     _dim(_mesh.dimension()),
      39          15 :     _cell_volume(1.0),
      40          15 :     _buffer_size(1),
      41          15 :     _perf_plan(registerTimedSection("fftw_plan_r2r", 2)),
      42          45 :     _perf_fft(registerTimedSection("fftw_execute", 2))
      43             : {
      44             :   // FFT needs a regular grid of values to operate on
      45          30 :   if (isParamValid("grid"))
      46             :   {
      47             :     // if valid use the user-supplied sampling grid
      48           0 :     _grid = getParam<std::vector<int>>("grid");
      49           0 :     if (_grid.size() != _dim)
      50           0 :       paramError("grid", "Number of entries must match the mesh dimension ", _dim);
      51             :   }
      52             :   else
      53             :   {
      54          15 :     auto * mytrim_mesh = dynamic_cast<MyTRIMMesh *>(&_mesh);
      55          15 :     if (!mytrim_mesh)
      56           0 :       mooseError("Must specify 'grid' parameter if not using MyTRIMMesh");
      57             : 
      58             :     // querying the mesh to set grid
      59          45 :     for (unsigned int i = 0; i < _dim; ++i)
      60          30 :       _grid.push_back(mytrim_mesh->getCellCountInDimension(i));
      61             :   }
      62             : 
      63             :   // refine grid by max_h_level powers of two
      64          30 :   auto max_h_level = getParam<unsigned int>("max_h_level");
      65          45 :   for (auto & count : _grid)
      66          30 :     count = count << max_h_level;
      67             : 
      68             :   // get mesh extents and calculate space required and estimate spectrum bins
      69          45 :   for (unsigned int i = 0; i < _dim; ++i)
      70             :   {
      71          30 :     _min_corner(i) = _mesh.getMinInDimension(i);
      72          30 :     _max_corner(i) = _mesh.getMaxInDimension(i);
      73          30 :     _box_size(i) = _max_corner(i) - _min_corner(i);
      74          30 :     _cell_volume *= _box_size(i) / _grid[i];
      75             : 
      76             :     // last direction needs to be padded for in-place transforms
      77          30 :     _buffer_size *= (i == _dim - 1) ? ((_grid[i] >> 1) + 1) << 1 : _grid[i];
      78             :   }
      79             : 
      80             :   // initialize FFTW buffer and plan
      81          15 :   _buffer.resize(_buffer_size);
      82          15 :   std::vector<fftw_r2r_kind> kind(_dim, FFTW_R2HC);
      83             :   {
      84          15 :     TIME_SECTION(_perf_plan);
      85          15 :     _plan = fftw_plan_r2r(
      86          15 :         _dim, _grid.data(), _buffer.data(), _buffer.data(), kind.data(), FFTW_ESTIMATE);
      87             :   }
      88          15 : }
      89             : 
      90          30 : FourierTransform::~FourierTransform()
      91             : {
      92             :   // destroy FFTW plan
      93          15 :   fftw_destroy_plan(_plan);
      94          30 : }
      95             : 
      96             : void
      97          15 : FourierTransform::initialize()
      98             : {
      99             :   // avoid assign here to make _sure_ the data address (which has already been passed to the plan)
     100             :   // stays the same (i.e. no reallocations)!
     101             :   std::fill(_buffer.begin(), _buffer.end(), 0.0);
     102          15 : }
     103             : 
     104             : void
     105        8250 : FourierTransform::execute()
     106             : {
     107             :   // loop over all quadrature points in the element
     108       41250 :   for (std::size_t qp = 0; qp < _qrule->n_points(); ++qp)
     109             :   {
     110             :     // get FFT grid cell index for current QP
     111             :     std::size_t index = 0;
     112       99000 :     for (unsigned int component = 0; component < _dim; ++component)
     113             :     {
     114             :       const unsigned int cell =
     115       66000 :           ((_q_point[qp](component) - _min_corner(component)) * _grid[component]) /
     116       66000 :           _box_size(component);
     117       66000 :       index = index * _grid[component] + cell;
     118             :     }
     119             : 
     120             :     // add qp contribution to FFT grid cell
     121       33000 :     _buffer[index] += _var[qp] * _JxW[qp] * _coord[qp] / _cell_volume;
     122             :   }
     123        8250 : }
     124             : 
     125             : void
     126          12 : FourierTransform::finalize()
     127             : {
     128             :   // get the sum of all buffers to all processes
     129          12 :   gatherSum(_buffer);
     130             : 
     131             :   // for now only run FFT on processor 0
     132          12 :   if (processor_id() != 0)
     133             :     return;
     134             : 
     135             :   // execute transform
     136             :   {
     137           9 :     TIME_SECTION(_perf_fft);
     138           9 :     fftw_execute(_plan);
     139             :   }
     140             : }
     141             : 
     142             : void
     143           3 : FourierTransform::threadJoin(const UserObject & y)
     144             : {
     145             :   const FourierTransform & vpp = static_cast<const FourierTransform &>(y);
     146             :   mooseAssert(vpp._buffer.size() == _buffer.size(), "Inconsistent buffer sizes across threads");
     147             : 
     148             :   // sum the two vectors using the binary functional std::plus
     149             :   std::transform(vpp._buffer.begin(),
     150             :                  vpp._buffer.end(),
     151             :                  _buffer.begin(),
     152             :                  _buffer.begin(),
     153             :                  std::plus<double>());
     154           3 : }
     155             : 
     156             : #endif