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             : 
       9             : #include "FFTBufferBase.h"
      10             : #include "MyTRIMMesh.h"
      11             : 
      12             : #include "MooseTypes.h"
      13             : #include "AuxiliarySystem.h"
      14             : #include "AllLocalDofIndicesThread.h"
      15             : #include "RankTwoTensor.h"
      16             : #include "RankThreeTensor.h"
      17             : #include "RankFourTensor.h"
      18             : 
      19             : #include <type_traits>
      20             : 
      21             : template <typename T>
      22             : InputParameters
      23          84 : FFTBufferBase<T>::validParams()
      24             : {
      25          84 :   InputParameters params = ElementUserObject::validParams();
      26          84 :   params.addClassDescription("Fourier transform data buffer object.");
      27         168 :   params.addRangeCheckedParam<std::vector<int>>(
      28             :       "grid",
      29             :       "grid > 0",
      30             :       "Number of grid cells in each dimension to compute "
      31             :       "the FFT on (can be omitted when using MyTRIMMesh)");
      32         168 :   params.addCoupledVar("moose_variable",
      33             :                        "Optional AuxVariable to take the initial condition of the buffer from.");
      34             : 
      35             :   // make sure we run the object on initial to apply the initial conditions
      36          84 :   params.set<ExecFlagEnum>("execute_on") = EXEC_INITIAL;
      37             : 
      38          84 :   return params;
      39           0 : }
      40             : 
      41             : template <typename T>
      42          42 : FFTBufferBase<T>::FFTBufferBase(const InputParameters & parameters)
      43             :   : ElementUserObject(parameters),
      44          84 :     _mesh(_subproblem.mesh()),
      45          42 :     _dim(_mesh.dimension()),
      46          42 :     _cell_volume(1.0),
      47          42 :     _moose_variable(coupledComponents("moose_variable")),
      48          42 :     _k_table(_dim),
      49         126 :     _how_many(_real_space_data.howMany())
      50             : {
      51             :   // make sure Real is double
      52             :   static_assert(
      53             :       std::is_same<Real, double>::value,
      54             :       "Libmesh must be build with double precision floating point numbers as the 'Real' type.");
      55             : 
      56             :   static_assert(sizeof(std::complex<Real>) == 2 * sizeof(Real),
      57             :                 "Complex numbers based on 'Real' should have the size of two 'Real' types.");
      58             : 
      59             :   // FFT needs a regular grid of values to operate on
      60          84 :   if (isParamValid("grid"))
      61             :   {
      62             :     // cannot specify a corresponding MOOSE Variable when running mesh-less
      63           0 :     if (!_moose_variable.empty())
      64           0 :       paramError("moose_variable",
      65             :                  "You cannot specify a corresponding MOOSE Variable when running mesh-less, i.e. "
      66             :                  "using the 'grid' parameter.");
      67             : 
      68             :     // if valid use the user-supplied sampling grid
      69           0 :     _grid = getParam<std::vector<int>>("grid");
      70           0 :     if (_grid.size() != _dim)
      71           0 :       paramError("grid", "Number of entries must match the mesh dimension ", _dim);
      72             :   }
      73             :   else
      74             :   {
      75          42 :     auto * mytrim_mesh = dynamic_cast<MyTRIMMesh *>(&_mesh);
      76          42 :     if (!mytrim_mesh)
      77           0 :       mooseError("Must specify 'grid' parameter if not using MyTRIMMesh");
      78             : 
      79             :     // querying the mesh to set grid
      80         144 :     for (unsigned int i = 0; i < _dim; ++i)
      81         102 :       _grid.push_back(mytrim_mesh->getCellCountInDimension(i));
      82             :   }
      83             : 
      84             :   // check optional coupled MOOSE variable (only for Real valued buffers)
      85          42 :   if (!_moose_variable.empty())
      86             :   {
      87          21 :     if (_moose_variable.size() != _how_many)
      88           0 :       paramError("moose_variable", "Variable needs to have ", _how_many, " components.");
      89             : 
      90          48 :     for (unsigned i = 0; i < _moose_variable.size(); ++i)
      91             :     {
      92             :       // check
      93          27 :       auto var = getVar("moose_variable", i);
      94          27 :       if (var->order() != 0)
      95           0 :         paramError("moose_variable", "Variable needs to have CONSTANT order.");
      96          27 :       if (var->isNodal())
      97           0 :         paramError("moose_variable", "Variable must be elemental.");
      98             : 
      99          27 :       _moose_variable[i] = &coupledValue("moose_variable", i);
     100             :     }
     101             :   }
     102             : 
     103             :   // get mesh extents and calculate space required and estimate spectrum bins
     104             :   std::size_t real_space_buffer_size = 1;
     105             :   std::size_t reciprocal_space_buffer_size = 1;
     106         144 :   for (unsigned int i = 0; i < _dim; ++i)
     107             :   {
     108         102 :     _min_corner(i) = _mesh.getMinInDimension(i);
     109         102 :     _max_corner(i) = _mesh.getMaxInDimension(i);
     110         102 :     _box_size(i) = _max_corner(i) - _min_corner(i);
     111         102 :     _cell_volume *= _box_size(i) / _grid[i];
     112             : 
     113             :     // unumber of physical grid cells
     114         102 :     real_space_buffer_size *= _grid[i];
     115             :     // last direction needs to be roughly cut in half
     116         102 :     reciprocal_space_buffer_size *= (i == _dim - 1) ? ((_grid[i] >> 1) + 1) : _grid[i];
     117             : 
     118             :     // precompute kvector components for current direction
     119         102 :     _k_table[i].resize(_grid[i]);
     120        4770 :     for (int j = 0; j < _grid[i]; ++j)
     121        4668 :       _k_table[i][j] = 2.0 * libMesh::pi *
     122        4668 :                        ((j < (_grid[i] >> 1) + 1) ? Real(j) : (Real(j) - _grid[i])) / _box_size(i);
     123             : 
     124         102 :     if (i == _dim - 1)
     125          42 :       _k_table[i].resize((_grid[i] >> 1) + 1);
     126             :   }
     127             :   _real_space_data.resize(real_space_buffer_size);
     128             :   _reciprocal_space_data.resize(reciprocal_space_buffer_size);
     129             : 
     130             :   // compute stride and start pointer
     131          42 :   _real_space_data_start = reinterpret_cast<Real *>(_real_space_data.start(0));
     132          42 :   _reciprocal_space_data_start = reinterpret_cast<Complex *>(_reciprocal_space_data.start(0));
     133             : 
     134          42 :   _real_space_data_stride = reinterpret_cast<char *>(_real_space_data.start(1)) -
     135          42 :                             reinterpret_cast<char *>(_real_space_data_start);
     136          42 :   _reciprocal_space_data_stride = reinterpret_cast<char *>(_reciprocal_space_data.start(1)) -
     137          42 :                                   reinterpret_cast<char *>(_reciprocal_space_data_start);
     138             : 
     139          42 :   if (_real_space_data_stride % sizeof(Real) != 0)
     140           0 :     mooseError("Invalid data alignment");
     141          42 :   _real_space_data_stride /= sizeof(Real);
     142             : 
     143          42 :   if (_reciprocal_space_data_stride % sizeof(Complex) != 0)
     144           0 :     mooseError("Invalid data alignment");
     145          42 :   _reciprocal_space_data_stride /= sizeof(Complex);
     146          42 : }
     147             : 
     148             : template <typename T>
     149             : void
     150      682224 : FFTBufferBase<T>::execute()
     151             : {
     152             :   // get  grid / buffer location
     153      682224 :   Point centroid = _current_elem->vertex_average();
     154             :   std::size_t a = 0;
     155     2636496 :   for (unsigned int i = 0; i < _dim; ++i)
     156     1954272 :     a = a * _grid[i] + std::floor(((centroid(i) - _min_corner(i)) * _grid[i]) / _box_size(i));
     157             : 
     158             :   // copy solution value from the variables into the buffer
     159      872928 :   for (unsigned i = 0; i < _moose_variable.size(); ++i)
     160      190704 :     _real_space_data_start[a * _how_many + i] = (*_moose_variable[i])[0];
     161      682224 : }
     162             : 
     163             : template <typename T>
     164             : void
     165          15 : FFTBufferBase<T>::forward()
     166             : {
     167          15 :   forwardRaw();
     168          15 : }
     169             : 
     170             : template <typename T>
     171             : void
     172          78 : FFTBufferBase<T>::backward()
     173             : {
     174          78 :   backwardRaw();
     175          78 :   _real_space_data *= backwardScale();
     176          78 : }
     177             : 
     178             : template <typename T>
     179             : const T &
     180     3267120 : FFTBufferBase<T>::operator()(const Point & p) const
     181             : {
     182             :   std::size_t a = 0;
     183    12750480 :   for (unsigned int i = 0; i < _dim; ++i)
     184     9483360 :     a = a * _grid[i] + std::floor(((p(i) - _min_corner(i)) * _grid[i]) / _box_size(i));
     185     3267120 :   return _real_space_data[a];
     186             : }
     187             : 
     188             : template <typename T>
     189             : T &
     190           0 : FFTBufferBase<T>::operator()(const Point & p)
     191             : {
     192             :   std::size_t a = 0;
     193           0 :   for (unsigned int i = 0; i < _dim; ++i)
     194           0 :     a = a * _grid[i] + std::floor(((p(i) - _min_corner(i)) * _grid[i]) / _box_size(i));
     195           0 :   return _real_space_data[a];
     196             : }
     197             : 
     198             : // explicit instantiations
     199             : template class FFTBufferBase<Real>;
     200             : template class FFTBufferBase<RealVectorValue>;
     201             : template class FFTBufferBase<RankTwoTensor>;
     202             : template class FFTBufferBase<RankThreeTensor>;
     203             : template class FFTBufferBase<RankFourTensor>;