Refactor: FFT use fft::create everywhere

ahojukka5 · ahojukka5 · commit 207bf391b5ed · 2025-05-10T19:06:19.000+03:00
diff --git a/apps/aluminumNew/Aluminum.hpp b/apps/aluminumNew/Aluminum.hpp
@@ -162,26 +162,18 @@ class Aluminum : public Model {
   }
 
   void prepare_operators(double dt) {
-    const Decomposition &decomp = get_decomposition();
-    const FFT &fft = get_fft();
-    World w = decomposition::get_world(decomp);
-    auto spacing = get_spacing(w);
-    auto size = get_size(w);
-    auto dx = spacing[0];
-    auto dy = spacing[1];
-    auto dz = spacing[2];
-    auto Lx = size[0];
-    auto Ly = size[1];
-    auto Lz = size[2];
-
-    std::array<int, 3> low = get_outbox(fft).low;
-    std::array<int, 3> high = get_outbox(fft).high;
+    auto &fft = get_fft();
+    auto world = get_world();
+    auto [dx, dy, dz] = get_spacing(world);
+    auto [Lx, Ly, Lz] = get_size(world);
+    auto low = get_outbox(fft).low;
+    auto high = get_outbox(fft).high;
 
     int idx = 0;
-    const double pi = std::atan(1.0) * 4.0;
-    const double fx = 2.0 * pi / (dx * Lx);
-    const double fy = 2.0 * pi / (dy * Ly);
-    const double fz = 2.0 * pi / (dz * Lz);
+    double pi = std::atan(1.0) * 4.0;
+    double fx = 2.0 * pi / (dx * Lx);
+    double fy = 2.0 * pi / (dy * Ly);
+    double fz = 2.0 * pi / (dz * Lz);
 
     for (int k = low[2]; k <= high[2]; k++) {
       for (int j = low[1]; j <= high[1]; j++) {
diff --git a/apps/tungsten.cpp b/apps/tungsten.cpp
@@ -102,20 +102,13 @@ class Tungsten : public Model {
   }
 
   void prepare_operators(double dt) {
-    const FFT &fft = get_fft();
-    const Decomposition &decomp = get_decomposition();
-    World w = decomposition::get_world(decomp);
-    auto spacing = get_spacing(w);
-    auto size = get_size(w);
-    auto dx = spacing[0];
-    auto dy = spacing[1];
-    auto dz = spacing[2];
-    auto Lx = size[0];
-    auto Ly = size[1];
-    auto Lz = size[2];
-
-    std::array<int, 3> low = get_outbox(fft).low;
-    std::array<int, 3> high = get_outbox(fft).high;
+    auto &fft = get_fft();
+    auto &world = get_world();
+    auto [dx, dy, dz] = get_spacing(world);
+    auto [Lx, Ly, Lz] = get_size(world);
+
+    auto low = get_outbox(fft).low;
+    auto high = get_outbox(fft).high;
 
     int idx = 0;
     const double pi = std::atan(1.0) * 4.0;
diff --git a/examples/02_domain_decomposition.cpp b/examples/02_domain_decomposition.cpp
@@ -33,7 +33,7 @@ int main(int argc, char *argv[]) {
   int comm_rank = 0, comm_size = 2;
   auto world1 = world::create({32, 4, 4});
   auto decomp1 = make_decomposition(world1, comm_rank, comm_size);
-  // cout << decomp1 << endl;
+  cout << decomp1 << endl;
 
   // In practice, we let MPI communicator to decide the number of subdomains.
   MPI_Init(&argc, &argv);
diff --git a/examples/04_diffusion_model.cpp b/examples/04_diffusion_model.cpp
@@ -83,9 +83,8 @@ class Diffusion : public Model {
     if (rank0) cout << "Allocate space" << endl;
 
     // Get references to world, fft and domain decomposition
-    const Decomposition &decomp = get_decomposition();
-    const World &w = decomposition::get_world(decomp);
-    FFT &fft = get_fft();
+    auto &world = get_world();
+    auto &fft = get_fft();
 
     // Allocate space for the main variable and it's fourier transform
     psi.resize(fft.size_inbox());
@@ -104,30 +103,33 @@ class Diffusion : public Model {
     World is defining the global dimensions of the problem as well as origin and
     chosen discretization parameters.
     */
-    if (rank0) cout << "World: " << w << endl;
+    if (rank0) cout << "World: " << world << endl;
 
     /*
     Upper and lower limits for this particular MPI rank, in both inbox and
-    outbox, are given by domain decomposition object
+    outbox, are given by fft object
     */
-    Int3 i_low = get_inbox(fft).low;
-    Int3 i_high = get_inbox(fft).high;
-    Int3 o_low = get_outbox(fft).low;
-    Int3 o_high = get_outbox(fft).high;
+    auto i_low = get_inbox(fft).low;
+    auto i_high = get_inbox(fft).high;
+    auto o_low = get_outbox(fft).low;
+    auto o_high = get_outbox(fft).high;
 
     /*
     Typically initial conditions are constructed elsewhere. However, to keep
     things as simple as possible, the initial condition can be also constructed
     here.
     */
     if (rank0) cout << "Create initial condition" << endl;
+
+    auto size = get_size(world);
+    auto origin = get_origin(world);
+    auto spacing = get_spacing(world);
+
     int idx = 0;
     double D = 1.0;
     for (int k = i_low[2]; k <= i_high[2]; k++) {
       for (int j = i_low[1]; j <= i_high[1]; j++) {
         for (int i = i_low[0]; i <= i_high[0]; i++) {
-          auto origin = get_origin(w);
-          auto spacing = get_spacing(w);
           double x = origin[0] + i * spacing[0];
           double y = origin[1] + j * spacing[1];
           double z = origin[2] + k * spacing[2];
@@ -144,8 +146,6 @@ class Diffusion : public Model {
     if (rank0) cout << "Prepare operators" << endl;
     idx = 0;
     double pi = std::atan(1.0) * 4.0;
-    auto spacing = get_spacing(w);
-    auto size = get_size(w);
     double fx = 2.0 * pi / (spacing[0] * size[0]);
     double fy = 2.0 * pi / (spacing[1] * size[1]);
     double fz = 2.0 * pi / (spacing[2] * size[2]);
diff --git a/examples/05_simulator.cpp b/examples/05_simulator.cpp
@@ -98,9 +98,8 @@ class Diffusion : public Model {
   }
 
   void prepare_operators(double dt) {
-    const Decomposition &decomp = get_decomposition();
-    const World &w = decomposition::get_world(decomp);
-    FFT &fft = get_fft();
+    auto &w = get_world();
+    auto &fft = get_fft();
     std::array<int, 3> low = get_outbox(fft).low;
     std::array<int, 3> high = get_outbox(fft).high;
 
diff --git a/examples/09_parallel_fft_high_level.cpp b/examples/09_parallel_fft_high_level.cpp
@@ -35,7 +35,7 @@ int main(int argc, char *argv[]) {
   // Create output array to store FFT results. If requested array is of type T =
   // complex<double>, then array will be constructed using complex indices so
   // that it matches the Fourier-space, i.e. first dimension is floor(Lx/2) + 1.
-  Array<complex<double>, 3> output(get_outbox_size(fft));
+  Array<complex<double>, 3> output(get_outbox(fft).size);
 
   std::cout << "input: " << input << std::endl;   // this is {4, 3, 2}
   std::cout << "output: " << output << std::endl; // this is {3, 3, 2}
diff --git a/examples/12_cahn_hilliard.cpp b/examples/12_cahn_hilliard.cpp
@@ -141,7 +141,7 @@ int main(int argc, char **argv) {
   // file_count
   writer.set_uri(sprintf("cahn_hilliard_%04i.vti", file_count));
   writer.set_field_name("concentration");
-  writer.set_domain(get_size(world), get_inbox_size(fft), get_inbox_offset(fft));
+  writer.set_domain(get_size(world), get_inbox(fft).size, get_inbox(fft).low);
   writer.set_origin(get_origin(world));
   writer.set_spacing(get_spacing(world));
   writer.initialize();
diff --git a/examples/diffusion_model_with_custom_initial_condition.cpp b/examples/diffusion_model_with_custom_initial_condition.cpp
@@ -49,22 +49,23 @@ class Gaussian : public FieldModifier {
     if (m.rank0) {
       cout << "Applying custom initial condition at time " << t << endl;
     }
-    const World &w = m.get_world();
-    const Decomposition &d = m.get_decomposition();
-    Field &f = m.get_field();
-    auto low = get_inbox(d).low;
-    auto high = get_inbox(d).high;
+    auto &world = m.get_world();
+    auto &field = m.get_field();
+    auto &fft = m.get_fft();
+    auto origin = get_origin(world);
+    auto spacing = get_spacing(world);
+
+    auto low = get_inbox(fft).low;
+    auto high = get_inbox(fft).high;
     long int idx = 0;
 
     for (int k = low[2]; k <= high[2]; k++) {
       for (int j = low[1]; j <= high[1]; j++) {
         for (int i = low[0]; i <= high[0]; i++) {
-          auto origin = get_origin(w);
-          auto spacing = get_spacing(w);
           double x = origin[0] + i * spacing[0];
           double y = origin[1] + j * spacing[1];
           double z = origin[2] + k * spacing[2];
-          f[idx] = exp(-(x * x + y * y + z * z) / (4.0 * m_D));
+          field[idx] = exp(-(x * x + y * y + z * z) / (4.0 * m_D));
           idx += 1;
         }
       }
diff --git a/include/openpfc/array.hpp b/include/openpfc/array.hpp
@@ -34,7 +34,7 @@ template <typename T, size_t D> class Array {
    * @param fft
    */
   Array(const FFT &fft, std::true_type)
-      : index(get_outbox_size(fft), get_outbox_offset(fft)) {}
+      : index(get_outbox(fft).size, get_outbox(fft).low) {}
 
   /**
    * @brief Construct a new Array object from FFT, using inbox, if
@@ -43,7 +43,7 @@ template <typename T, size_t D> class Array {
    * @param decomp
    */
   Array(const FFT &fft, std::false_type)
-      : index(get_inbox_size(fft), get_inbox_offset(fft)) {}
+      : index(get_inbox(fft).size, get_inbox(fft).low) {}
 
   // Custom type trait to check if a type is complex
   template <typename U> struct is_complex : std::false_type {};
diff --git a/include/openpfc/fft.hpp b/include/openpfc/fft.hpp
@@ -46,33 +46,34 @@ struct FFTLayout {
 const FFTLayout create(const Decomposition &decomposition, int r2c_direction,
                        int num_domains);
 
+inline auto get_real_box(const FFTLayout &layout, int i) {
+  return layout.m_real_boxes.at(i);
+}
+
+inline auto get_complex_box(const FFTLayout &layout, int i) {
+  return layout.m_complex_boxes.at(i);
+}
+
 } // namespace layout
 
-using heffte::box3d;
-using pfc::types::Bool3;
 using pfc::types::Int3;
 using pfc::types::Real3;
-using Box3D = heffte::box3d<int>; ///< Type alias for 3D integer box.
 
-inline heffte::fft3d_r2c<heffte::backend::fftw>
-make_fft(const Box3D &inbox, const Box3D &outbox, const int &r2c_direction,
-         MPI_Comm comm, heffte::plan_options plan_options) {
-  return heffte::fft3d_r2c<heffte::backend::fftw>(inbox, outbox, r2c_direction, comm,
-                                                  plan_options);
-}
+using Decomposition = pfc::decomposition::Decomposition<pfc::csys::CartesianTag>;
+using ComplexVector = std::vector<std::complex<double>>;
+using fft_r2c = heffte::fft3d_r2c<heffte::backend::fftw>;
+using box3di = heffte::box3d<int>; ///< Type alias for 3D integer box.
 
 /**
  * @brief FFT class for performing forward and backward Fast Fourier
  * Transformations.
  */
 struct FFT {
 
-  using ComplexVector = std::vector<std::complex<double>>;
+  // const Decomposition m_decomposition; /**< The Decomposition object. */
+  // const box3di m_inbox, m_outbox;      /**< Local inbox and outbox boxes. */
 
-  const Decomposition m_decomposition;        /**< The Decomposition object. */
-  const heffte::box3d<int> m_inbox, m_outbox; /**< Local inbox and outbox boxes. */
-  const int m_r2c_direction; /**< Real-to-complex symmetry direction. */
-  const heffte::fft3d_r2c<heffte::backend::fftw> m_fft; /**< HeFFTe FFT object. */
+  const fft_r2c m_fft;     /**< HeFFTe FFT object. */
   ComplexVector m_wrk;     /**< Workspace vector for FFT computations. */
   double m_fft_time = 0.0; /**< Recorded FFT computation time. */
 
@@ -86,21 +87,15 @@ struct FFT {
    * @param plan_options Optional plan options for configuring the FFT behavior.
    * @param world The World object providing the domain size information.
    */
-  FFT(const Decomposition &decomposition, MPI_Comm comm,
-      heffte::plan_options plan_options)
-      : m_decomposition(decomposition), m_inbox(decomposition.m_inbox),
-        m_outbox(decomposition.m_outbox), m_r2c_direction(0),
-        m_fft(make_fft(m_inbox, m_outbox, m_r2c_direction, comm, plan_options)),
-        m_wrk(std::vector<std::complex<double>>(m_fft.size_workspace())){};
+  FFT(fft_r2c fft) : m_fft(std::move(fft)), m_wrk(m_fft.size_workspace()) {}
 
   /**
    * @brief Performs the forward FFT transformation.
    *
    * @param in Input vector of real values.
    * @param out Output vector of complex values.
    */
-  void forward(const std::vector<double> &in,
-               std::vector<std::complex<double>> &out) {
+  void forward(const std::vector<double> &in, ComplexVector &out) {
     m_fft_time -= MPI_Wtime();
     m_fft.forward(in.data(), out.data(), m_wrk.data());
     m_fft_time += MPI_Wtime();
@@ -112,8 +107,7 @@ struct FFT {
    * @param in Input vector of complex values.
    * @param out Output vector of real values.
    */
-  void backward(const std::vector<std::complex<double>> &in,
-                std::vector<double> &out) {
+  void backward(const ComplexVector &in, std::vector<double> &out) {
     m_fft_time -= MPI_Wtime();
     m_fft.backward(in.data(), out.data(), m_wrk.data(), heffte::scale::full);
     m_fft_time += MPI_Wtime();
@@ -136,7 +130,7 @@ struct FFT {
    *
    * @return Reference to the Decomposition object.
    */
-  const Decomposition &get_decomposition() { return m_decomposition; }
+  // const Decomposition &get_decomposition() { return m_decomposition; }
 
   /**
    * @brief Returns the size of the inbox used for FFT computations.
@@ -160,24 +154,36 @@ struct FFT {
   size_t size_workspace() const { return m_fft.size_workspace(); }
 };
 
-inline const Box3D &get_inbox(const FFT &fft) noexcept { return fft.m_inbox; }
-inline const Box3D &get_outbox(const FFT &fft) noexcept { return fft.m_outbox; }
-inline const Int3 &get_inbox_size(const FFT &fft) noexcept {
-  return get_inbox(fft).size;
-}
-inline const Int3 &get_inbox_offset(const FFT &fft) noexcept {
-  return get_inbox(fft).low;
-}
-inline const Int3 &get_outbox_size(const FFT &fft) noexcept {
-  return get_outbox(fft).size;
+inline const auto &get_fft_object(const FFT &fft) noexcept { return fft.m_fft; }
+
+inline const auto get_inbox(const FFT &fft) noexcept {
+  return get_fft_object(fft).inbox();
 }
-inline const Int3 &get_outbox_offset(const FFT &fft) noexcept {
-  return get_outbox(fft).low;
+
+inline const auto get_outbox(const FFT &fft) noexcept {
+  return get_fft_object(fft).outbox();
 }
 
-FFT create(const Decomposition &decomposition);
+/**
+ * @brief Creates an FFT object based on the given decomposition and MPI
+ * communicator.
+ *
+ * @param decomposition The Decomposition object defining the domain
+ * decomposition.
+ * @param comm The MPI communicator for parallel computations.
+ * @param options Optional plan options for configuring the FFT behavior.
+ * @return An FFT object containing the FFT configuration and data.
+ */
 FFT create(const Decomposition &decomposition, MPI_Comm comm,
            heffte::plan_options options);
+/**
+ * @brief Creates an FFT object based on the given decomposition.
+ *
+ * @param decomposition The Decomposition object defining the domain
+ * decomposition.
+ * @return An FFT object containing the FFT configuration and data.
+ */
+FFT create(const Decomposition &decomposition);
 
 } // namespace fft
 
diff --git a/include/openpfc/initial_conditions/file_reader.hpp b/include/openpfc/initial_conditions/file_reader.hpp
@@ -25,11 +25,10 @@ class FileReader : public FieldModifier {
 
   void apply(Model &m, double) override {
     const FFT &fft = m.get_fft();
-    const Decomposition &decomp = m.get_decomposition();
-    const auto &world = decomposition::get_world(decomp);
-    auto world_size = get_size(world);
-    auto inbox_size = get_inbox(fft).size;
-    auto inbox_offset = get_inbox(fft).low;
+    const auto &world = m.get_world();
+    const auto world_size = get_size(world);
+    const auto inbox_size = get_inbox(fft).size;
+    const auto inbox_offset = get_inbox(fft).low;
 
     Field &f = m.get_real_field(get_field_name());
     std::cout << "Reading initial condition from file" << get_filename()
diff --git a/include/openpfc/initial_conditions/random_seeds.hpp b/include/openpfc/initial_conditions/random_seeds.hpp
diff --git a/include/openpfc/model.hpp b/include/openpfc/model.hpp
diff --git a/include/openpfc/simulator.hpp b/include/openpfc/simulator.hpp
diff --git a/src/openpfc/fft.cpp b/src/openpfc/fft.cpp
