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 "FourierLengthScale.h"
      11             : #include "FourierTransform.h"
      12             : #include "MyTRIMMesh.h"
      13             : 
      14             : #include "libmesh/utility.h"
      15             : 
      16             : #include <functional>
      17             : 
      18             : registerMooseObject("MagpieApp", FourierLengthScale);
      19             : 
      20             : InputParameters
      21          16 : FourierLengthScale::validParams()
      22             : {
      23          16 :   InputParameters params = GeneralPostprocessor::validParams();
      24          16 :   params.addClassDescription(
      25             :       "Compute the average length scale form a given Fourier transform in length units.");
      26          32 :   params.addRequiredParam<UserObjectName>(
      27             :       "fourier_transform", "FourierTransform user object to compute the length scale of");
      28          16 :   return params;
      29           0 : }
      30             : 
      31           8 : FourierLengthScale::FourierLengthScale(const InputParameters & parameters)
      32             :   : GeneralPostprocessor(parameters),
      33           8 :     _fourier_transform(getUserObject<FourierTransform>("fourier_transform")),
      34           8 :     _dim(_fourier_transform.dimension()),
      35           8 :     _grid(_fourier_transform.getGrid()),
      36           8 :     _box_size(_fourier_transform.getBoxSize()),
      37           8 :     _buffer(_fourier_transform.getBuffer())
      38             : {
      39           8 : }
      40             : 
      41             : void
      42         171 : FourierLengthScale::computeLengthScale(std::vector<int> & c, std::vector<Real> & F, std::size_t i)
      43             : {
      44         171 :   if (i == _dim - 1)
      45             :   {
      46             :     Real sum_high = 0.0;
      47         330 :     for (unsigned int j = 1; j < _dim; ++j)
      48         165 :       sum_high += F[j];
      49             : 
      50             :     auto & cc = c[_dim - 1];
      51        6615 :     for (cc = 0; cc < _grid[_dim - 1]; ++cc)
      52             :     {
      53             :       // do stuff at highest dimension
      54        6450 :       F[_dim - 1] =
      55        6450 :           Utility::pow<2>((cc * 2 > _grid[_dim - 1] ? _grid[_dim - 1] - cc : cc) / _box_size(0));
      56        6450 :       Real sum = sum_high + F[_dim - 1];
      57             : 
      58             :       // find bin and add to spectrum
      59        6450 :       if (sum > 0)
      60             :       {
      61             :         // compute frequency
      62        6444 :         const Real frequency = std::sqrt(sum);
      63             : 
      64             :         // compute number of wavevectors contributing to this frequency (frequency^(_dim-1))
      65             :         Real normalization = 1.0;
      66       12888 :         for (unsigned int j = 1; j < _dim; ++j)
      67        6444 :           normalization *= frequency;
      68             : 
      69             :         // weight is the normalized intensity at this frequency
      70        6444 :         const Real weight = Utility::pow<2>(_buffer[_index]) / normalization;
      71             : 
      72             :         // compute weighted average of frequencies
      73        6444 :         _length_scale += frequency * weight;
      74        6444 :         _weight_sum += weight;
      75             :       }
      76        6450 :       _index++;
      77             :     }
      78             :   }
      79             :   else
      80             :     // iterate over lower dimensions
      81         171 :     for (c[i] = 0; c[i] < _grid[i]; ++c[i])
      82             :     {
      83         165 :       F[i] = Utility::pow<2>((c[i] * 2 > _grid[i] ? _grid[i] - c[i] : c[i]) / _box_size(i));
      84         165 :       computeLengthScale(c, F, i + 1);
      85             :     }
      86         171 : }
      87             : 
      88             : void
      89           8 : FourierLengthScale::execute()
      90             : {
      91             :   // clear data
      92           8 :   _length_scale = 0.0;
      93           8 :   _weight_sum = 0.0;
      94             : 
      95             :   // for now only compute length scale on processor 0
      96           8 :   if (processor_id() == 0)
      97             :   {
      98             :     // compute the length scale
      99           6 :     std::vector<int> c(_dim, 0);
     100           6 :     std::vector<Real> F(_dim);
     101           6 :     _index = 0;
     102           6 :     computeLengthScale(c, F, 0);
     103             : 
     104             :     // divide by weight and take reciprocal (frequency -> length)
     105           6 :     _length_scale = _weight_sum / _length_scale;
     106             :   }
     107           8 : }
     108             : 
     109             : void
     110           8 : FourierLengthScale::finalize()
     111             : {
     112             :   // broadcast result to all cores (lazily using gatherSum - all other ranks have zero length scale)
     113           8 :   gatherSum(_length_scale);
     114           8 : }
     115             : 
     116             : #endif // FFTW3_ENABLED