Implement fft::create

ahojukka5 · ahojukka5 · commit fc5ebc2c7be8 · 2025-05-10T16:03:29.000+03:00
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -103,6 +103,7 @@ add_library(openpfc
     src/openpfc/core/box3d.cpp
     src/openpfc/core/decomposition.cpp
     src/openpfc/factory/decomposition_factory.cpp
+    src/openpfc/fft.cpp
     # Add more .cpp files as you go
 )
 
diff --git a/apps/aluminumNew/aluminumTest.cpp b/apps/aluminumNew/aluminumTest.cpp
@@ -44,10 +44,9 @@ using namespace Catch::Matchers;
 TEST_CASE("Aluminum functionality", "[Aluminum]") {
   SECTION("Step model and calculate norm of the result") {
     MPI_Worker worker(0, nullptr);
-    World world = world::create({32, 32, 32});
-    Decomposition decomp = make_decomposition(world);
-    auto plan_options = heffte::default_options<heffte::backend::fftw>();
-    FFT fft(decomp, MPI_COMM_WORLD, plan_options, world);
+    auto world = world::create({32, 32, 32});
+    auto decomp = make_decomposition(world);
+    auto fft = fft::create(decomp);
 
     Aluminum aluminum(world);
     aluminum.set_n0(-0.0060);
diff --git a/examples/02_domain_decomposition.cpp b/examples/02_domain_decomposition.cpp
@@ -31,20 +31,20 @@ int main(int argc, char *argv[]) {
   // Domain decomposition can be done in a manual manner, just by giving the
   // size of the calculation domain, id number and total number of subdomains:
   int comm_rank = 0, comm_size = 2;
-  World world1 = world::create({32, 4, 4});
-  Decomposition decomp1 = make_decomposition(world1, comm_rank, comm_size);
+  auto world1 = world::create({32, 4, 4});
+  auto decomp1 = make_decomposition(world1, comm_rank, comm_size);
   // cout << decomp1 << endl;
 
   // In practice, we let MPI communicator to decide the number of subdomains.
   MPI_Init(&argc, &argv);
   MPI_Comm comm = MPI_COMM_WORLD;
   MPI_Comm_rank(comm, &comm_rank);
-  World world2 = world::create({32, 4, 4});
-  Decomposition decomp2 = make_decomposition(world2, comm);
-  // if (comm_rank == 0) cout << decomp2 << endl;
+  auto world2 = world::create({32, 4, 4});
+  auto decomp2 = make_decomposition(world2, comm);
+  if (comm_rank == 0) cout << decomp2 << endl;
 
   // By default, MPI_COMM_WORLD is used, so the above example can be simplified:
-  // cout << make_decomposition(world2) << endl;
+  cout << make_decomposition(world2) << endl;
 
   MPI_Finalize();
   return 0;
diff --git a/examples/03_parallel_fft.cpp b/examples/03_parallel_fft.cpp
@@ -126,10 +126,9 @@ int main(int argc, char *argv[]) {
   MPI_Comm_size(comm, &num_procs);
 
   // Construct world, decomposition and fft
-  World world = world::create({8, 1, 1});
-  Decomposition decomposition = make_decomposition(world, comm);
-  auto plan_options = heffte::default_options<heffte::backend::fftw>();
-  FFT fft(decomposition, comm, plan_options, world);
+  auto world = world::create({8, 1, 1});
+  auto decomposition = make_decomposition(world, comm);
+  auto fft = fft::create(decomposition);
 
   // Create two vectors; in contains input data and results are stored to out
   vector<double> in(fft.size_inbox());
diff --git a/examples/04_diffusion_model.cpp b/examples/04_diffusion_model.cpp
@@ -230,10 +230,9 @@ void run() {
   double z0 = -0.5 * Lz * dz;
 
   // Construct world, decomposition, fft and model
-  World world = world::create({Lx, Ly, Lz}, {x0, y0, z0}, {dx, dy, dz});
-  Decomposition decomp = make_decomposition(world);
-  auto plan_options = heffte::default_options<heffte::backend::fftw>();
-  FFT fft(decomp, MPI_COMM_WORLD, plan_options, world);
+  auto world = world::create({Lx, Ly, Lz}, {x0, y0, z0}, {dx, dy, dz});
+  auto decomp = make_decomposition(world);
+  auto fft = fft::create(decomp);
   Diffusion model(world);
 
   // Define time
diff --git a/examples/05_simulator.cpp b/examples/05_simulator.cpp
@@ -183,11 +183,10 @@ void run() {
   std::array<int, 3> dimensions = {L, L, L};
   std::array<double, 3> discretization = {h, h, h};
   std::array<double, 3> origin = {o, o, o};
-  World world = world::create(dimensions, origin, discretization);
+  auto world = world::create(dimensions, origin, discretization);
 
-  Decomposition decomp = make_decomposition(world);
-  auto plan_options = heffte::default_options<heffte::backend::fftw>();
-  FFT fft(decomp, MPI_COMM_WORLD, plan_options, world);
+  auto decomp = make_decomposition(world);
+  auto fft = fft::create(decomp);
   Diffusion model(world);
 
   // Define time
diff --git a/examples/09_parallel_fft_high_level.cpp b/examples/09_parallel_fft_high_level.cpp
@@ -11,8 +11,8 @@ int main(int argc, char *argv[]) {
 
   // Create MPI session, World and Decomposition
   MPI_Worker worker(argc, argv);
-  World world = world::create({4, 3, 2});
-  Decomposition decomp = make_decomposition(world);
+  auto world = world::create({4, 3, 2});
+  auto decomp = make_decomposition(world);
 
   // Create input field
   // DiscreteField<double, 3> input(decomp);
@@ -30,8 +30,7 @@ int main(int argc, char *argv[]) {
   // Fill input field with random numbers
   apply(input, [&](auto, auto, auto) { return dis(gen); });
 
-  auto plan_options = heffte::default_options<heffte::backend::fftw>();
-  FFT fft(decomp, MPI_COMM_WORLD, plan_options, world);
+  auto fft = fft::create(decomp);
 
   // 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
diff --git a/examples/12_cahn_hilliard.cpp b/examples/12_cahn_hilliard.cpp
@@ -113,10 +113,9 @@ int main(int argc, char **argv) {
   double z0 = 0.0;
 
   // Construct world, decomposition, fft and model
-  World world = world::create({Lx, Ly, Lz}, {x0, y0, z0}, {dx, dy, dz});
-  Decomposition decomp = make_decomposition(world);
-  auto plan_options = heffte::default_options<heffte::backend::fftw>();
-  FFT fft(decomp, MPI_COMM_WORLD, plan_options, world);
+  auto world = world::create({Lx, Ly, Lz}, {x0, y0, z0}, {dx, dy, dz});
+  auto decomposition = make_decomposition(world);
+  auto fft = fft::create(decomposition);
   CahnHilliard model(world);
   model.set_fft(fft);
 
diff --git a/examples/diffusion_model_with_custom_initial_condition.cpp b/examples/diffusion_model_with_custom_initial_condition.cpp
@@ -92,9 +92,8 @@ void run() {
   Time time(tspan, saveat);
 
   MPI_Comm comm = MPI_COMM_WORLD;
-  Decomposition decomposition = make_decomposition(world, comm);
-  auto plan_options = heffte::default_options<heffte::backend::fftw>();
-  FFT fft(decomposition, comm, plan_options, world);
+  auto decomposition = make_decomposition(world, comm);
+  auto fft = fft::create(decomposition);
   Diffusion model(world);
   Simulator simulator(model, time);
 
diff --git a/include/openpfc/boundary_conditions/fixed_bc.hpp b/include/openpfc/boundary_conditions/fixed_bc.hpp
@@ -1,13 +1,15 @@
 // SPDX-FileCopyrightText: 2025 VTT Technical Research Centre of Finland Ltd
 // SPDX-License-Identifier: AGPL-3.0-or-later
 
-#ifndef PFC_BOUNDARY_CONDITIONS_FIXED_BC_HPP
-#define PFC_BOUNDARY_CONDITIONS_FIXED_BC_HPP
+#pragma once
 
-#include "../field_modifier.hpp"
+#include "openpfc/core/types.hpp"
+#include "openpfc/field_modifier.hpp"
 
 namespace pfc {
 
+using pfc::types::Int3;
+
 class FixedBC : public FieldModifier {
 
 private:
@@ -56,5 +58,3 @@ class FixedBC : public FieldModifier {
 };
 
 } // namespace pfc
-
-#endif // PFC_BOUNDARY_CONDITIONS_FIXED_BC_HPP
diff --git a/include/openpfc/core/decomposition.hpp b/include/openpfc/core/decomposition.hpp
@@ -11,6 +11,10 @@
 #include <vector>
 
 namespace pfc {
+
+/**
+ * @brief Namespace for decomposition-related classes and functions.
+ */
 namespace decomposition {
 
 using heffte::box3d;
diff --git a/include/openpfc/core/types.hpp b/include/openpfc/core/types.hpp
@@ -4,6 +4,7 @@
 #pragma once
 
 #include <array>
+#include <stdexcept>
 
 namespace pfc {
 
@@ -131,4 +132,9 @@ struct Spacing3 {
 };
 
 } // namespace types
+
+using Int3 = types::Int3;
+using Real3 = types::Real3;
+using Bool3 = types::Bool3;
+
 } // namespace pfc
diff --git a/include/openpfc/fft.hpp b/include/openpfc/fft.hpp
@@ -1,8 +1,7 @@
 // SPDX-FileCopyrightText: 2025 VTT Technical Research Centre of Finland Ltd
 // SPDX-License-Identifier: AGPL-3.0-or-later
 
-#ifndef PFC_FFT_HPP
-#define PFC_FFT_HPP
+#pragma once
 
 #include "core/decomposition.hpp"
 #include "openpfc/backends/heffte_adapter.hpp" // Ensure this is included for the conversion operator
@@ -13,6 +12,41 @@
 #include <mpi.h>
 
 namespace pfc {
+namespace fft {
+
+namespace layout {
+
+using box3di = heffte::box3d<int>;
+using Decomposition = pfc::decomposition::Decomposition<pfc::csys::CartesianTag>;
+using pfc::types::Int3;
+
+/**
+ * @brief Structure to hold the layout of FFT data.
+ *
+ * This structure contains the decomposition object, the direction of
+ * real-to-complex symmetry, and the boxes for real and complex FFT data.
+ */
+struct FFTLayout {
+  const Decomposition m_decomposition;       ///< The Decomposition object.
+  const int m_r2c_direction = 0;             ///< Real-to-complex symmetry direction.
+  const std::vector<box3di> m_real_boxes;    ///< Real boxes for FFT.
+  const std::vector<box3di> m_complex_boxes; ///< Complex boxes for FFT.
+};
+
+/**
+ * @brief Creates an FFTLayout object based on the given decomposition and
+ * parameters.
+ *
+ * @param decomposition The Decomposition object defining the domain
+ * decomposition.
+ * @param r2c_direction The direction of real-to-complex symmetry.
+ * @param num_domains The number of domains for the FFT layout.
+ * @return An FFTLayout object containing the layout information.
+ */
+const FFTLayout create(const Decomposition &decomposition, int r2c_direction,
+                       int num_domains);
+
+} // namespace layout
 
 using heffte::box3d;
 using pfc::types::Bool3;
@@ -39,10 +73,8 @@ struct FFT {
   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. */
-  ComplexVector m_wrk;       /**< Workspace vector for FFT computations. */
-  double m_fft_time = 0.0;   /**< Recorded FFT computation time. */
-  const World &m_world;      /**< Reference to the World object. */
-  heffte::box3d<int> domain; /**< Domain converted from World object. */
+  ComplexVector m_wrk;     /**< Workspace vector for FFT computations. */
+  double m_fft_time = 0.0; /**< Recorded FFT computation time. */
 
   /**
    * @brief Constructs an FFT object with the given Decomposition and MPI
@@ -55,15 +87,11 @@ struct FFT {
    * @param world The World object providing the domain size information.
    */
   FFT(const Decomposition &decomposition, MPI_Comm comm,
-      heffte::plan_options plan_options, const World &world)
+      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())),
-        m_world(world), domain(to_heffte_box(world)){
-                            // Use to_heffte_box for conversion
-                            // Explicit conversion
-                        };
+        m_wrk(std::vector<std::complex<double>>(m_fft.size_workspace())){};
 
   /**
    * @brief Performs the forward FFT transformation.
@@ -147,6 +175,13 @@ inline const Int3 &get_outbox_offset(const FFT &fft) noexcept {
   return get_outbox(fft).low;
 }
 
-} // namespace pfc
+FFT create(const Decomposition &decomposition);
+FFT create(const Decomposition &decomposition, MPI_Comm comm,
+           heffte::plan_options options);
 
-#endif
+} // namespace fft
+
+using FFT = fft::FFT;                     ///< Type alias for FFT class.
+using FFTLayout = fft::layout::FFTLayout; ///< Type alias for FFTLayout class.
+
+} // namespace pfc
diff --git a/include/openpfc/ui.hpp b/include/openpfc/ui.hpp
@@ -651,10 +651,10 @@ template <class ConcreteModel> class App {
     World world(ui::from_json<World>(m_settings));
     std::cout << "World: " << world << std::endl;
 
-    Decomposition decomp = make_decomposition(world, m_comm);
+    auto decomp = make_decomposition(world, m_comm);
     auto plan_options =
         ui::from_json<heffte::plan_options>(m_settings["plan_options"]);
-    FFT fft(decomp, m_comm, plan_options, world);
+    auto fft = fft::create(decomp, m_comm, plan_options);
     Time time(ui::from_json<Time>(m_settings));
     ConcreteModel model(world);
     model.set_fft(fft);
diff --git a/src/openpfc/fft.cpp b/src/openpfc/fft.cpp
@@ -0,0 +1,50 @@
+// SPDX-FileCopyrightText: 2025 VTT Technical Research Centre of Finland Ltd
+// SPDX-License-Identifier: AGPL-3.0-or-later
+
+#include "openpfc/fft.hpp"
+
+namespace pfc {
+namespace fft {
+namespace layout {
+
+const FFTLayout create(const Decomposition &decomposition, int r2c_direction,
+                       int num_domains) {
+  if (num_domains <= 0) {
+    throw std::logic_error("Cannot construct domain decomposition: !(nprocs > 0)");
+  }
+  auto world = get_world(decomposition);
+  auto [N1r, N2r, N3r] = get_size(world); // real size
+  auto [N1c, N2c, N3c] = get_size(world); // complex size
+  if (r2c_direction == 0) {
+    N1c = N1c / 2 + 1;
+  } else if (r2c_direction == 1) {
+    N2c = N2c / 2 + 1;
+  } else if (r2c_direction == 2) {
+    N3c = N3c / 2 + 1;
+  } else {
+    throw std::logic_error("Invalid r2c_direction: " +
+                           std::to_string(r2c_direction));
+  }
+  box3di real_indices({0, 0, 0}, {N1r - 1, N2r - 1, N3r - 1});
+  box3di complex_indices({0, 0, 0}, {N1c - 1, N2c - 1, N3c - 1});
+  Int3 grid = heffte::proc_setup_min_surface(real_indices, num_domains);
+  std::vector<box3di> real_boxes = heffte::split_world(real_indices, grid);
+  std::vector<box3di> complex_boxes = heffte::split_world(complex_indices, grid);
+  return FFTLayout{decomposition, r2c_direction, real_boxes, complex_boxes};
+}
+
+} // namespace layout
+
+FFT create(const Decomposition &decomposition, MPI_Comm comm,
+           heffte::plan_options options) {
+  return FFT(decomposition, comm, options);
+}
+
+FFT create(const Decomposition &decomposition) {
+  MPI_Comm comm = MPI_COMM_WORLD;
+  auto options = heffte::default_options<heffte::backend::fftw>();
+  return create(decomposition, comm, options);
+}
+
+} // namespace fft
+} // namespace pfc
diff --git a/tests/test_fft.cpp b/tests/test_fft.cpp
diff --git a/tests/test_fft_setting.cpp b/tests/test_fft_setting.cpp
diff --git a/tests/test_field_modifier.cpp b/tests/test_field_modifier.cpp
diff --git a/tests/test_ic_constant.cpp b/tests/test_ic_constant.cpp
diff --git a/tests/test_model.cpp b/tests/test_model.cpp
diff --git a/tests/test_simulator.cpp b/tests/test_simulator.cpp

Original file line number	Diff line number	Diff line change
`@@ -103,6 +103,7 @@ add_library(openpfc`
`103`	`103`	`src/openpfc/core/box3d.cpp`
`104`	`104`	`src/openpfc/core/decomposition.cpp`
`105`	`105`	`src/openpfc/factory/decomposition_factory.cpp`
	`106`	`+ src/openpfc/fft.cpp`
`106`	`107`	`# Add more .cpp files as you go`
`107`	`108`	`)`
`108`	`109`