diff --git a/include/userobjects/FFTBufferBase.h b/include/userobjects/FFTBufferBase.h new file mode 100644 index 00000000..b328fe43 --- /dev/null +++ b/include/userobjects/FFTBufferBase.h @@ -0,0 +1,73 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ + +#pragma once + +#include "GeneralUserObject.h" + +template +class FFTBufferBase; + +/** + * Generic FFT interleaved data buffer base class + */ +template +class FFTBufferBase : public GeneralUserObject +{ +public: + static InputParameters validParams(); + + FFTBufferBase(const InputParameters & parameters); + + virtual void initialize() {} + virtual void execute() {} + virtual void finalize() {} + + // transforms + virtual void forward() = 0; + virtual void backward() = 0; + + // data access + const T & operator[](std::size_t i) const { return _buffer[i]; } + T & operator[](std::size_t i) { return _buffer[i]; } + +protected: + /// get the addres of the first data element of the ith object in the bufefr + Real * start(std::size_t i); + + /// get the number of transforms required for type T + std::size_t howMany() const; + + ///@{ mesh data + MooseMesh & _mesh; + unsigned int _dim; + ///@} + + /// grid size for FFT (needs to be signed for FFTW) + std::vector _grid; + + ///@{ simulation box extents + Point _min_corner; + Point _max_corner; + Point _box_size; + ///@} + + /// FFT grid cell volume + Real _cell_volume; + + ///@{ FFT data buffer + std::vector _buffer; + std::size_t _buffer_size; + ///@} + + /// pointer to the start of the data + Real * _start; + + /// stride in units of double size + std::ptrdiff_t _stride; +}; diff --git a/include/userobjects/FFTWBufferBase.h b/include/userobjects/FFTWBufferBase.h new file mode 100644 index 00000000..54fc7229 --- /dev/null +++ b/include/userobjects/FFTWBufferBase.h @@ -0,0 +1,39 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ +#ifdef FFTW3_ENABLED + +#pragma once + +#include "FFTBufferBase.h" + +#include "fftw3.h" + +template +class FFTWBufferBase; + +/** + * FFTW specific interleaved data buffer base class + */ +template +class FFTWBufferBase : public FFTBufferBase +{ +public: + FFTWBufferBase(const InputParameters & parameters); + ~FFTWBufferBase(); + + // transforms + void forward() override; + void backward() override; + +protected: + /// FFTW plans + fftw_plan _forward_plan; + fftw_plan _backward_plan; +}; + +#endif diff --git a/src/userobjects/FFTBufferBase.C b/src/userobjects/FFTBufferBase.C new file mode 100644 index 00000000..44d77ce6 --- /dev/null +++ b/src/userobjects/FFTBufferBase.C @@ -0,0 +1,176 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ +#include "FFTBufferBase.h" +#include "MooseTypes.h" +#include "MyTRIMMesh.h" +#include "RankTwoTensor.h" +#include "RankThreeTensor.h" +#include "RankFourTensor.h" + +#include + +// template <> +// InputParameters +// validParams() +// { +// InputParameters params = validParams(); +// params.addClassDescription(""); +// return params; +// } + +template +InputParameters +FFTBufferBase::validParams() +{ + InputParameters params = GeneralUserObject::validParams(); + params.addClassDescription("Fourier transform data buffer object."); + params.addRangeCheckedParam>( + "grid", + "grid > 0", + "Number of grid cells in each dimension to compute " + "the FFT on (can be omitted when using MyTRIMMesh)"); + params.addParam( + "max_h_level", 0, "Further grid refinement to apply when using MyTRIMMesh with adaptivity"); + return params; +} + +template +FFTBufferBase::FFTBufferBase(const InputParameters & parameters) + : GeneralUserObject(parameters), + _mesh(_subproblem.mesh()), + _dim(_mesh.dimension()), + _cell_volume(1.0), + _buffer_size(1) + +{ + // make sure Real is double + static_assert( + std::is_same::value, + "Libmesh must be build with double precision floating point numbers as the 'Real' type."); + + // FFT needs a regular grid of values to operate on + if (isParamValid("grid")) + { + // if valid use the user-supplied sampling grid + _grid = getParam>("grid"); + if (_grid.size() != _dim) + paramError("grid", "Number of entries must match the mesh dimension ", _dim); + } + else + { + auto * mytrim_mesh = dynamic_cast(&_mesh); + if (!mytrim_mesh) + mooseError("Must specify 'grid' parameter if not using MyTRIMMesh"); + + // querying the mesh to set grid + for (unsigned int i = 0; i < _dim; ++i) + _grid.push_back(mytrim_mesh->getCellCountInDimension(i)); + } + + // refine grid by max_h_level powers of two + auto max_h_level = getParam("max_h_level"); + for (auto & count : _grid) + count = count << max_h_level; + + // get mesh extents and calculate space required and estimate spectrum bins + for (unsigned int i = 0; i < _dim; ++i) + { + _min_corner(i) = _mesh.getMinInDimension(i); + _max_corner(i) = _mesh.getMaxInDimension(i); + _box_size(i) = _max_corner(i) - _min_corner(i); + _cell_volume *= _box_size(i) / _grid[i]; + + // last direction needs to be padded for in-place transforms + _buffer_size *= (i == _dim - 1) ? ((_grid[i] >> 1) + 1) << 1 : _grid[i]; + } + _buffer.resize(_buffer_size); + + // compute stride and start pointer + auto istart = start(0); + std::ptrdiff_t istride = reinterpret_cast(start(1)) - reinterpret_cast(istart); + if (istride % sizeof(Real) != 0) + mooseError("Invalid data alignment"); + istride /= sizeof(Real); +} + +template <> +Real * +FFTBufferBase::start(std::size_t i) +{ + return &_buffer[i]; +} + +template <> +Real * +FFTBufferBase::start(std::size_t i) +{ + return &_buffer[i](0); +} + +template <> +Real * +FFTBufferBase::start(std::size_t i) +{ + return &_buffer[i](0, 0); +} + +template <> +Real * +FFTBufferBase::start(std::size_t i) +{ + return &_buffer[i](0, 0, 0); +} + +template <> +Real * +FFTBufferBase::start(std::size_t i) +{ + return &_buffer[i](0, 0, 0, 0); +} + +template <> +std::size_t +FFTBufferBase::howMany() const +{ + return 1; +} + +template <> +std::size_t +FFTBufferBase::howMany() const +{ + return LIBMESH_DIM; +} + +template <> +std::size_t +FFTBufferBase::howMany() const +{ + return LIBMESH_DIM * LIBMESH_DIM; +} + +template <> +std::size_t +FFTBufferBase::howMany() const +{ + return LIBMESH_DIM * LIBMESH_DIM * LIBMESH_DIM; +} + +template <> +std::size_t +FFTBufferBase::howMany() const +{ + return LIBMESH_DIM * LIBMESH_DIM * LIBMESH_DIM * LIBMESH_DIM; +} + +// explicit instantiation +template class FFTBufferBase; +template class FFTBufferBase; +template class FFTBufferBase; +template class FFTBufferBase; +template class FFTBufferBase; diff --git a/src/userobjects/FFTWBufferBase.C b/src/userobjects/FFTWBufferBase.C new file mode 100644 index 00000000..87344bfd --- /dev/null +++ b/src/userobjects/FFTWBufferBase.C @@ -0,0 +1,60 @@ +/**********************************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */ +/* */ +/* Copyright 2017 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/**********************************************************************/ +#ifdef FFTW3_ENABLED + +#include "FFTWBufferBase.h" +#include "MyTRIMMesh.h" +#include "MooseTypes.h" +#include "RankTwoTensor.h" +#include "RankThreeTensor.h" +#include "RankFourTensor.h" + +template +FFTWBufferBase::FFTWBufferBase(const InputParameters & parameters) + : FFTBufferBase(parameters), _forward_plan(nullptr), _backward_plan(nullptr) +{ + // create plans +} + +template +FFTWBufferBase::~FFTWBufferBase() +{ + // destroy FFTW plans + fftw_destroy_plan(_forward_plan); + fftw_destroy_plan(_backward_plan); +} + +template +void +FFTWBufferBase::forward() +{ + // execute plan + fftw_execute(_forward_plan); +} + +template +void +FFTWBufferBase::backward() +{ + // execute plan + fftw_execute(_backward_plan); +} + +// explicit instantiation and registration +#define FFTWBufferInstance(T) \ + template class FFTWBufferBase; \ + using T##FFTWBuffer = FFTWBufferBase; \ + registerMooseObject("MagpieApp", T##FFTWBuffer) + +FFTWBufferInstance(Real); +FFTWBufferInstance(RealVectorValue); +FFTWBufferInstance(RankTwoTensor); +FFTWBufferInstance(RankThreeTensor); +FFTWBufferInstance(RankFourTensor); + +#endif