Add FFTData (idaholab#401)

dschwen · Apr 29, 2020 · ee72a7f · ee72a7f
1 parent 4ff3845
commit ee72a7f
Show file tree

Hide file tree

Showing 11 changed files with 290 additions and 26 deletions.
diff --git a/include/problems/FFTProblem.h b/include/problems/FFTProblem.h
@@ -39,8 +39,6 @@ class FFTProblem : public FEProblem
 
   /// dummy system for the FFT variables
   AuxiliarySystem _fft_dummy_system;
-
-  unsigned int _fft_var_number;
 };
 
 template <typename T>

diff --git a/include/userobjects/FFTBufferBase.h b/include/userobjects/FFTBufferBase.h
@@ -18,8 +18,7 @@ class FFTBufferBase;
 #define usingFFTBufferBaseMembers                                                                  \
   using ElementUserObject::_perf_graph;                                                            \
   using FFTBufferBase<T>::_dim;                                                                    \
-  using FFTBufferBase<T>::_real_space_grid;                                                        \
-  using FFTBufferBase<T>::_reciprocal_space_grid;                                                  \
+  using FFTBufferBase<T>::_grid;                                                                   \
   using FFTBufferBase<T>::_real_space_data;                                                        \
   using FFTBufferBase<T>::_reciprocal_space_data;                                                  \
   using FFTBufferBase<T>::_real_space_data_start;                                                  \

diff --git a/include/utils/ComplexTypes.h b/include/utils/ComplexTypes.h
@@ -17,10 +17,10 @@
 #include "libmesh/vector_value.h"
 
 typedef std::complex<Real> Complex;
-typedef VectorValue<std::complex<Real>> ComplexVectorValue;
-typedef RankTwoTensorTempl<std::complex<Real>> ComplexRankTwoTensor;
-typedef RankThreeTensorTempl<std::complex<Real>> ComplexRankThreeTensor;
-typedef RankFourTensorTempl<std::complex<Real>> ComplexRankFourTensor;
+typedef VectorValue<Complex> ComplexVectorValue;
+typedef RankTwoTensorTempl<Complex> ComplexRankTwoTensor;
+typedef RankThreeTensorTempl<Complex> ComplexRankThreeTensor;
+typedef RankFourTensorTempl<Complex> ComplexRankFourTensor;
 
 // helper template to select the corresponding complex tensor type
 template <typename T>

diff --git a/include/utils/FFTData.h b/include/utils/FFTData.h
@@ -0,0 +1,65 @@
+/**********************************************************************/
+/*                     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 "ComplexTypes.h"
+
+// helper template to select the corresponding scalar type
+template <typename T>
+struct FFTScalarType
+{
+  using type = Real;
+};
+
+template <typename T>
+struct FFTScalarType<std::complex<T>>
+{
+  using type = Complex;
+};
+
+/**
+ * Helper class to hold the reciprocal space data
+ */
+template <typename T>
+class FFTData
+{
+  using ScalarT = typename FFTScalarType<T>::type;
+
+public:
+  FFTData(std::size_t size = 0) { resize(size); }
+  void resize(std::size_t size) { _buffer.resize(size); }
+
+  ///@{ data access by index
+  const T & operator[](std::size_t i) const { return _buffer[i]; }
+  T & operator[](std::size_t i) { return _buffer[i]; }
+  ///@}
+
+  ///@{ convenience math operators
+  FFTData<T> & operator+=(FFTData<T> const & rhs);
+  FFTData<T> & operator-=(FFTData<T> const & rhs);
+  FFTData<T> & operator*=(FFTData<ScalarT> const & rhs);
+  FFTData<T> & operator/=(FFTData<ScalarT> const & rhs);
+  FFTData<T> & operator*=(Real rhs);
+  FFTData<T> & operator/=(Real rhs);
+  FFTData<T> & operator=(FFTData<T> const & rhs);
+  ///@}
+
+  /// return the number of proper grid cells
+  std::size_t size() const { return _buffer.size(); }
+
+  /// get the addres of the first data element of the ith object in the buffer
+  void * start(std::size_t i);
+
+  /// get the number of transforms required for type T
+  std::size_t howMany() const;
+
+protected:
+  /// FFT data buffer
+  std::vector<T> _buffer;
+};
diff --git a/include/variables/MooseFFTVariable.h b/include/variables/MooseFFTVariable.h
@@ -53,7 +53,7 @@ class MooseFFTVariable : public MooseVariable
   {
     mooseError("Nodal FFT variables are not supported.");
   }
-  virtual void reinitNodesNeighbor(const std::vector<dof_id_type> & nodes)
+  virtual void reinitNodesNeighbor(const std::vector<dof_id_type> &)
   {
     mooseError("Nodal FFT variables are not supported.");
   }

diff --git a/src/executioners/SpectralExecutionerBase.C b/src/executioners/SpectralExecutionerBase.C
@@ -59,6 +59,7 @@ SpectralExecutionerBase::execute()
   // back and forth test
   auto & c_buffer = getFFTBuffer<Real>("c");
   auto c = c_buffer.realSpace();
+  auto c_tilde = c_buffer.reciprocalSpace();
   c_buffer.forward();
 
   auto & R_buffer = getFFTBuffer<RealVectorValue>("R");
@@ -114,24 +115,25 @@ SpectralExecutionerBase::kVectorMultiply(const FFTBufferBase<Real> & in_buffer,
   {
     case 1:
     {
+      const auto & ivec = in_buffer.kTable(0);
       const int ni = grid[0];
-      for (int i = 0; i < ni; ++i)
-      {
-        out[i](0) = in[i] * i;
-      }
+      for (int i = 0; i * 2 <= ni; ++i)
+        out[i](0) = in[i] * ivec[i];
       return;
     }
 
     case 2:
     {
       std::size_t index = 0;
+      const auto & ivec = in_buffer.kTable(0);
+      const auto & jvec = in_buffer.kTable(1);
       const int ni = grid[0];
       const int nj = grid[1];
       for (int i = 0; i < ni; ++i)
-        for (int j = 0; j < nj; ++j)
+        for (int j = 0; j * 2 <= nj; ++j)
         {
-          out[index](0) = in[index] * i;
-          out[index](1) = in[index] * j;
+          out[index](0) = in[index] * ivec[i];
+          out[index](1) = in[index] * jvec[j];
           index++;
         }
       return;
@@ -140,16 +142,19 @@ SpectralExecutionerBase::kVectorMultiply(const FFTBufferBase<Real> & in_buffer,
     case 3:
     {
       std::size_t index = 0;
+      const auto & ivec = in_buffer.kTable(0);
+      const auto & jvec = in_buffer.kTable(1);
+      const auto & kvec = in_buffer.kTable(2);
       const int ni = grid[0];
       const int nj = grid[1];
       const int nk = grid[2];
       for (int i = 0; i < ni; ++i)
         for (int j = 0; j < nj; ++j)
-          for (int k = 0; k < nk; ++k)
+          for (int k = 0; k * 2 <= nk; ++k)
           {
-            out[index](0) = in[index] * i;
-            out[index](1) = in[index] * j;
-            out[index](2) = in[index] * k;
+            out[index](0) = in[index] * ivec[i];
+            out[index](1) = in[index] * jvec[j];
+            out[index](2) = in[index] * kvec[k];
             index++;
           }
       return;

diff --git a/src/problems/FFTProblem.C b/src/problems/FFTProblem.C
@@ -85,7 +85,7 @@ FFTProblem::getVariable(THREAD_ID tid,
     auto params = MooseVariableBase::validParams();
     params.set<MooseEnum>("order") = "CONSTANT";
     params.set<MooseEnum>("family") = "MONOMIAL";
-    params.set<unsigned int>("_var_num") = _fft_var_number;
+    params.set<unsigned int>("_var_num") = fft_var_number;
     params.set<THREAD_ID>("tid") = tid;
     params.set<THREAD_ID>("_tid") = tid;
     params.set<Moose::VarKindType>("_var_kind") = Moose::VarKindType::VAR_AUXILIARY;

diff --git a/src/userobjects/FFTBufferBase.C b/src/userobjects/FFTBufferBase.C
@@ -45,6 +45,7 @@ FFTBufferBase<T>::FFTBufferBase(const InputParameters & parameters)
     _dim(_mesh.dimension()),
     _cell_volume(1.0),
     _moose_variable(coupledComponents("moose_variable")),
+    _k_table(_dim),
     _how_many(_real_space_data.howMany())
 {
   // make sure Real is double

diff --git a/src/utils/ComplexTypes.C b/src/utils/ComplexTypes.C
@@ -37,8 +37,11 @@ absoluteFuzzyEqual(const std::complex<Real> & var1,
 #include "RankThreeTensorImplementation.h"
 #include "RankFourTensorImplementation.h"
 
-template class VectorValue<std::complex<Real>>;
+namespace libMesh
+{
+template class VectorValue<Complex>;
+}
 
-template class RankTwoTensorTempl<std::complex<Real>>;
-template class RankThreeTensorTempl<std::complex<Real>>;
-template class RankFourTensorTempl<std::complex<Real>>;
+template class RankTwoTensorTempl<Complex>;
+template class RankThreeTensorTempl<Complex>;
+template class RankFourTensorTempl<Complex>;