diff --git a/include/userobjects/FFTBufferBase.h b/include/userobjects/FFTBufferBase.h index c7e33aab..e7db1c80 100644 --- a/include/userobjects/FFTBufferBase.h +++ b/include/userobjects/FFTBufferBase.h @@ -8,7 +8,9 @@ #pragma once +#include "ComplexTypes.h" #include "ElementUserObject.h" +#include "FFTData.h" template class FFTBufferBase; @@ -16,10 +18,14 @@ class FFTBufferBase; #define usingFFTBufferBaseMembers \ using ElementUserObject::_perf_graph; \ using FFTBufferBase::_dim; \ - using FFTBufferBase::_grid; \ - using FFTBufferBase::_buffer; \ - using FFTBufferBase::_start; \ - using FFTBufferBase::_stride; \ + using FFTBufferBase::_real_space_grid; \ + using FFTBufferBase::_reciprocal_space_grid; \ + using FFTBufferBase::_real_space_data; \ + using FFTBufferBase::_reciprocal_space_data; \ + using FFTBufferBase::_real_space_data_start; \ + using FFTBufferBase::_reciprocal_space_data_start; \ + using FFTBufferBase::_real_space_data_stride; \ + using FFTBufferBase::_reciprocal_space_data_stride; \ using FFTBufferBase::_how_many /** @@ -28,6 +34,8 @@ class FFTBufferBase; template class FFTBufferBase : public ElementUserObject { + using ComplexT = typename ComplexType::type; + public: static InputParameters validParams(); @@ -43,31 +51,20 @@ class FFTBufferBase : public ElementUserObject virtual void backward() = 0; ///@} - ///@{ data access by index - const T & operator[](std::size_t i) const { return _buffer[i]; } - T & operator[](std::size_t i) { return _buffer[i]; } - ///@} + ///@{ buffer access + FFTData & realSpace() { return _real_space_data; } + FFTData & reciprocalSpace() { return _reciprocal_space_data; } + const FFTData & realSpace() const { return _real_space_data; } + const FFTData & reciprocalSpace() const { return _reciprocal_space_data; } - ///@{ data access by location + ///@{ real space data access by location const T & operator()(const Point & p) const; T & operator()(const Point & p); ///@} - ///@{ convenience math operators - FFTBufferBase & operator+=(FFTBufferBase const & rhs); - FFTBufferBase & operator-=(FFTBufferBase const & rhs); - FFTBufferBase & operator*=(FFTBufferBase const & rhs); - FFTBufferBase & operator/=(FFTBufferBase const & rhs); - FFTBufferBase & operator*=(Real rhs); - FFTBufferBase & operator/=(Real rhs); - ///@} - /// return the number of grid cells along each dimension without padding const std::vector & grid() const { return _grid; } - /// return the number of proper grid cells without padding - const std::size_t & size() const { return _grid_size; } - /// return the buffer dimension const unsigned int & dim() const { return _dim; } @@ -79,12 +76,6 @@ class FFTBufferBase : public ElementUserObject } 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; @@ -103,15 +94,17 @@ class FFTBufferBase : public ElementUserObject Real _cell_volume; ///@{ FFT data buffer and unpadded number of grid cells - std::vector _buffer; - std::size_t _grid_size; + FFTData _real_space_data; + FFTData _reciprocal_space_data; ///@} /// pointer to the start of the data - Real * _start; + Real * _real_space_data_start; + Complex * _reciprocal_space_data_start; /// stride in units of double size - std::ptrdiff_t _stride; + std::ptrdiff_t _real_space_data_stride; + std::ptrdiff_t _reciprocal_space_data_stride; /// optional moose sister variabe (to obtain IC from) std::vector _moose_variable; diff --git a/src/executioners/SpectralExecutionerBase.C b/src/executioners/SpectralExecutionerBase.C index e7b0bdbb..09067b3e 100644 --- a/src/executioners/SpectralExecutionerBase.C +++ b/src/executioners/SpectralExecutionerBase.C @@ -57,16 +57,22 @@ SpectralExecutionerBase::execute() _fe_problem.advanceState(); // back and forth test - auto & c = getFFTBuffer("c"); - c.forward(); - auto & R = getFFTBuffer("R"); - R.forward(); + auto & c_buffer = getFFTBuffer("c"); + auto c = c_buffer.realSpace(); + c_buffer.forward(); + + auto & R_buffer = getFFTBuffer("R"); + auto R = R_buffer.realSpace(); + R_buffer.forward(); // gradient test - auto & u = getFFTBuffer("u"); - u.forward(); - auto & grad_u = getFFTBuffer("grad_u"); - kVectorMultiply(u, grad_u); + auto & u_buffer = getFFTBuffer("u"); + auto u = u_buffer.realSpace(); + u_buffer.forward(); + + auto & grad_u_buffer = getFFTBuffer("grad_u"); + auto grad_u = grad_u_buffer.realSpace(); + kVectorMultiply(u_buffer, grad_u_buffer); _time_step = 1; _fe_problem.execute(EXEC_FINAL); @@ -75,14 +81,15 @@ SpectralExecutionerBase::execute() _fe_problem.advanceState(); // back and forth test - c.backward(); - R.backward(); + c_buffer.backward(); + R_buffer.backward(); R /= 10000.0; c /= 10000.0; // gradient test - u.backward(); - grad_u.backward(); + u_buffer.backward(); + grad_u_buffer.backward(); + u /= 10000.0; grad_u /= 100.0; @@ -93,14 +100,17 @@ SpectralExecutionerBase::execute() } void -SpectralExecutionerBase::kVectorMultiply(const FFTBufferBase & in, - FFTBufferBase & out) const +SpectralExecutionerBase::kVectorMultiply(const FFTBufferBase & in_buffer, + FFTBufferBase & out_buffer) const { + mooseAssert(in_buffer.dim() == out_buffer.dim(), "Buffer dimensions must be equal"); + + const FFTData & in = in_buffer.reciprocalSpace(); + FFTData & out = out_buffer.reciprocalSpace(); mooseAssert(in.size() == out.size(), "Buffer sizes must be equal"); - mooseAssert(in.dim() == out.dim(), "Buffer dimensions must be equal"); - const auto & grid = in.grid(); - switch (in.dim()) + const auto & grid = in_buffer.grid(); + switch (in_buffer.dim()) { case 1: { @@ -108,8 +118,6 @@ SpectralExecutionerBase::kVectorMultiply(const FFTBufferBase & in, for (int i = 0; i < ni; ++i) { out[i](0) = in[i] * i; - out[i](1) = 0.0; - out[i](2) = 0.0; } return; } @@ -124,7 +132,6 @@ SpectralExecutionerBase::kVectorMultiply(const FFTBufferBase & in, { out[index](0) = in[index] * i; out[index](1) = in[index] * j; - out[index](2) = 0.0; index++; } return; diff --git a/src/userobjects/FFTBufferBase.C b/src/userobjects/FFTBufferBase.C index af984c19..5b37241f 100644 --- a/src/userobjects/FFTBufferBase.C +++ b/src/userobjects/FFTBufferBase.C @@ -44,16 +44,17 @@ FFTBufferBase::FFTBufferBase(const InputParameters & parameters) _mesh(_subproblem.mesh()), _dim(_mesh.dimension()), _cell_volume(1.0), - _grid_size(1), _moose_variable(coupledComponents("moose_variable")), - _k_table(_dim), - _how_many(howMany()) + _how_many(_real_space_data.howMany()) { // 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."); + static_assert(sizeof(std::complex) == 2 * sizeof(Real), + "Complex numbers based on 'Real' should have the size of two 'Real' types."); + // FFT needs a regular grid of values to operate on if (isParamValid("grid")) { @@ -99,7 +100,8 @@ FFTBufferBase::FFTBufferBase(const InputParameters & parameters) } // get mesh extents and calculate space required and estimate spectrum bins - std::size_t buffer_size = 1; + std::size_t real_space_buffer_size = 1; + std::size_t reciprocal_space_buffer_size = 1; for (unsigned int i = 0; i < _dim; ++i) { _min_corner(i) = _mesh.getMinInDimension(i); @@ -107,25 +109,35 @@ FFTBufferBase::FFTBufferBase(const InputParameters & parameters) _box_size(i) = _max_corner(i) - _min_corner(i); _cell_volume *= _box_size(i) / _grid[i]; - // unpadded size (number of physical grid cells) - _grid_size *= _grid[i]; - - // last direction needs to be padded for in-place transforms - buffer_size *= (i == _dim - 1) ? ((_grid[i] >> 1) + 1) << 1 : _grid[i]; + // unumber of physical grid cells + real_space_buffer_size *= _grid[i]; + // last direction needs to be roughly cut in half + reciprocal_space_buffer_size *= (i == _dim - 1) ? ((_grid[i] >> 1) + 1) : _grid[i]; // precompute kvector components for current direction _k_table[i].resize(_grid[i]); for (int j = 0; j < _grid[i]; ++j) _k_table[i][j] = (j * 2 > _grid[i] ? Real(_grid[i] - j) : Real(j)) / _box_size(i); } - _buffer.resize(buffer_size); + _real_space_data.resize(real_space_buffer_size); + _reciprocal_space_data.resize(reciprocal_space_buffer_size); // compute stride and start pointer - _start = reinterpret_cast(start(0)); - _stride = reinterpret_cast(start(1)) - reinterpret_cast(_start); - if (_stride % sizeof(Real) != 0) + _real_space_data_start = reinterpret_cast(_real_space_data.start(0)); + _reciprocal_space_data_start = reinterpret_cast(_reciprocal_space_data.start(0)); + + _real_space_data_stride = reinterpret_cast(_real_space_data.start(1)) - + reinterpret_cast(_real_space_data_start); + _reciprocal_space_data_stride = reinterpret_cast(_reciprocal_space_data.start(1)) - + reinterpret_cast(_reciprocal_space_data_start); + + if (_real_space_data_stride % sizeof(Real) != 0) mooseError("Invalid data alignment"); - _stride /= sizeof(Real); + _real_space_data_stride /= sizeof(Real); + + if (_reciprocal_space_data_stride % sizeof(Complex) != 0) + mooseError("Invalid data alignment"); + _reciprocal_space_data_stride /= sizeof(Complex); } template @@ -140,7 +152,7 @@ FFTBufferBase::execute() // copy solution value from the variables into the buffer for (unsigned i = 0; i < _moose_variable.size(); ++i) - _start[a * _how_many + i] = (*_moose_variable[i])[0]; + _real_space_data_start[a * _how_many + i] = (*_moose_variable[i])[0]; } template @@ -150,7 +162,7 @@ FFTBufferBase::operator()(const Point & p) const std::size_t a = 0; for (unsigned int i = 0; i < _dim; ++i) a = a * _grid[i] + std::floor(((p(i) - _min_corner(i)) * _grid[i]) / _box_size(i)); - return _buffer[a]; + return _real_space_data[a]; } template @@ -160,135 +172,10 @@ FFTBufferBase::operator()(const Point & p) std::size_t a = 0; for (unsigned int i = 0; i < _dim; ++i) a = a * _grid[i] + std::floor(((p(i) - _min_corner(i)) * _grid[i]) / _box_size(i)); - return _buffer[a]; -} - -template -FFTBufferBase & -FFTBufferBase::operator+=(FFTBufferBase const & rhs) -{ - for (std::size_t i = 0; i < _grid_size; ++i) - _buffer[i] += rhs[i]; - return *this; -} - -template -FFTBufferBase & -FFTBufferBase::operator-=(FFTBufferBase const & rhs) -{ - for (std::size_t i = 0; i < _grid_size; ++i) - _buffer[i] -= rhs[i]; - return *this; -} - -template -FFTBufferBase & -FFTBufferBase::operator*=(FFTBufferBase const & rhs) -{ - for (std::size_t i = 0; i < _grid_size; ++i) - _buffer[i] *= rhs[i]; - return *this; -} - -template -FFTBufferBase & -FFTBufferBase::operator/=(FFTBufferBase const & rhs) -{ - for (std::size_t i = 0; i < _grid_size; ++i) - _buffer[i] /= rhs[i]; - return *this; -} - -template -FFTBufferBase & -FFTBufferBase::operator*=(Real rhs) -{ - for (std::size_t i = 0; i < _grid_size; ++i) - _buffer[i] *= rhs; - return *this; -} - -template -FFTBufferBase & -FFTBufferBase::operator/=(Real rhs) -{ - const Real reciprocal = 1 / rhs; - for (std::size_t i = 0; i < _grid_size; ++i) - _buffer[i] *= reciprocal; - return *this; -} - -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; + return _real_space_data[a]; } -// explicit instantiation +// explicit instantiations template class FFTBufferBase; template class FFTBufferBase; template class FFTBufferBase; diff --git a/src/userobjects/FFTWBufferBase.C b/src/userobjects/FFTWBufferBase.C index 420e419e..59ff60ce 100644 --- a/src/userobjects/FFTWBufferBase.C +++ b/src/userobjects/FFTWBufferBase.C @@ -24,35 +24,33 @@ FFTWBufferBase::FFTWBufferBase(const InputParameters & parameters) { TIME_SECTION(_perf_plan); - std::vector forward_kind(_dim, FFTW_R2HC); - _forward_plan = fftw_plan_many_r2r(_dim, - _grid.data(), - _how_many, - _start, - _grid.data(), - _stride, - 1, - _start, - _grid.data(), - _stride, - 1, - forward_kind.data(), - FFTW_ESTIMATE); + _forward_plan = + fftw_plan_many_dft_r2c(_dim, + _grid.data(), + _how_many, + _real_space_data_start, + nullptr, + _real_space_data_stride, + 1, + reinterpret_cast(_reciprocal_space_data_start), + nullptr, + _reciprocal_space_data_stride, + 1, + FFTW_ESTIMATE); - std::vector backward_kind(_dim, FFTW_HC2R); - _backward_plan = fftw_plan_many_r2r(_dim, - _grid.data(), - _how_many, - _start, - _grid.data(), - _stride, - 1, - _start, - _grid.data(), - _stride, - 1, - backward_kind.data(), - FFTW_ESTIMATE); + _backward_plan = + fftw_plan_many_dft_c2r(_dim, + _grid.data(), + _how_many, + reinterpret_cast(_reciprocal_space_data_start), + nullptr, + _reciprocal_space_data_stride, + 1, + _real_space_data_start, + nullptr, + _real_space_data_stride, + 1, + FFTW_ESTIMATE); } }