Refactor examples to utilize World object for improved context handling and streamline decomposition and spacing access

Author: ahojukka5

examples/01_hello_world.cpp

@@ -39,10 +39,12 @@ int main() {
 
   // Retrieve world properties
   cout << "World properties:" << endl;
-  cout << "Dimensions: " << world.get_Lx() << " x " << world.get_Ly() << " x " << world.get_Lz() << endl;
-  cout << "Origin: (" << world.get_x0() << ", " << world.get_y0() << ", " << world.get_z0() << ")" << endl;
-  cout << "Discretization: dx = " << world.get_dx() << ", dy = " << world.get_dy() << ", dz = " << world.get_dz()
-       << endl;
+  auto size = world.size();
+  cout << "Dimensions: " << size[0] << " x " << size[1] << " x " << size[2] << endl;
+  auto world_origin = world.origin();
+  cout << "Origin: (" << world_origin[0] << ", " << world_origin[1] << ", " << world_origin[2] << ")" << endl;
+  auto spacing = world.spacing();
+  cout << "Spacing: dx = " << spacing[0] << ", dy = " << spacing[1] << ", dz = " << spacing[2] << endl;
 
   // Simpler way:
   cout << world << endl;

examples/02_domain_decomposition.cpp

@@ -3,8 +3,8 @@
 
 #include <iostream>
 #include <mpi.h>
+#include <openpfc/core/decomposition.hpp>
 #include <openpfc/core/world.hpp>
-#include <openpfc/decomposition.hpp>
 
 using namespace std;
 using namespace pfc;

examples/03_parallel_fft.cpp

@@ -4,8 +4,8 @@
 #include <iostream>
 #include <vector>
 
+#include <openpfc/core/decomposition.hpp>
 #include <openpfc/core/world.hpp>
-#include <openpfc/decomposition.hpp>
 #include <openpfc/fft.hpp>
 
 using namespace std;
@@ -113,7 +113,8 @@ int main(int argc, char *argv[]) {
   // Construct world, decomposition and fft
   World world({8, 1, 1});
   Decomposition decomposition(world, comm);
-  FFT fft(decomposition, comm);
+  auto plan_options = heffte::default_options<heffte::backend::fftw>();
+  FFT fft(decomposition, comm, plan_options, world);
 
   // Create two vectors; in contains input data and results are stored to out
   vector<double> in(fft.size_inbox());

examples/04_diffusion_model.cpp

@@ -5,8 +5,8 @@
 #include <iostream>
 #include <limits>
 
+#include <openpfc/core/decomposition.hpp>
 #include <openpfc/core/world.hpp>
-#include <openpfc/decomposition.hpp>
 #include <openpfc/fft.hpp>
 #include <openpfc/model.hpp>
 
@@ -125,9 +125,11 @@ class Diffusion : public Model {
     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++) {
-          double x = w.x0 + i * w.dx;
-          double y = w.y0 + j * w.dy;
-          double z = w.z0 + k * w.dz;
+          auto origin = w.origin();
+          auto spacing = w.spacing();
+          double x = origin[0] + i * spacing[0];
+          double y = origin[1] + j * spacing[1];
+          double z = origin[2] + k * spacing[2];
           psi[idx] = exp(-(x * x + y * y + z * z) / (4.0 * D));
           idx += 1;
         }
@@ -141,16 +143,17 @@ class Diffusion : public Model {
     if (rank0) cout << "Prepare operators" << endl;
     idx = 0;
     double pi = std::atan(1.0) * 4.0;
-    double fx = 2.0 * pi / (w.dx * w.Lx);
-    double fy = 2.0 * pi / (w.dy * w.Ly);
-    double fz = 2.0 * pi / (w.dz * w.Lz);
+    auto spacing = w.spacing();
+    auto size = w.size();
+    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]);
     for (int k = o_low[2]; k <= o_high[2]; k++) {
       for (int j = o_low[1]; j <= o_high[1]; j++) {
         for (int i = o_low[0]; i <= o_high[0]; i++) {
-          // Laplacian operator -k^2
-          double ki = (i <= w.Lx / 2) ? i * fx : (i - w.Lx) * fx;
-          double kj = (j <= w.Ly / 2) ? j * fy : (j - w.Ly) * fy;
-          double kk = (k <= w.Lz / 2) ? k * fz : (k - w.Lz) * fz;
+          double ki = (i <= size[0] / 2) ? i * fx : (i - size[0]) * fx;
+          double kj = (j <= size[1] / 2) ? j * fy : (j - size[1]) * fy;
+          double kk = (k <= size[2] / 2) ? k * fz : (k - size[2]) * fz;
           double kLap = -(ki * ki + kj * kj + kk * kk);
           opL[idx++] = 1.0 / (1.0 - dt * kLap);
         }
@@ -227,8 +230,9 @@ void run() {
   // Construct world, decomposition, fft and model
   World world({Lx, Ly, Lz}, {x0, y0, z0}, {dx, dy, dz});
   Decomposition decomp(world);
-  FFT fft(decomp);
-  Diffusion model(fft);
+  auto plan_options = heffte::default_options<heffte::backend::fftw>();
+  FFT fft(decomp, MPI_COMM_WORLD, plan_options, world);
+  Diffusion model(world);
 
   // Define time
   double t = 0.0;

examples/05_simulator.cpp

@@ -5,8 +5,8 @@
 #include <iostream>
 #include <limits>
 #include <memory>
+#include <openpfc/core/decomposition.hpp>
 #include <openpfc/core/world.hpp>
-#include <openpfc/decomposition.hpp>
 #include <openpfc/fft.hpp>
 #include <openpfc/field_modifier.hpp>
 #include <openpfc/model.hpp>
@@ -59,9 +59,11 @@ class GaussianIC : public FieldModifier {
     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++) {
-          double x = w.x0 + i * w.dx;
-          double y = w.y0 + j * w.dy;
-          double z = w.z0 + k * w.dz;
+          auto origin = w.origin();
+          auto spacing = w.spacing();
+          double x = origin[0] + i * spacing[0];
+          double y = origin[1] + j * spacing[1];
+          double z = origin[2] + k * spacing[2];
           field[idx] = exp(-(x * x + y * y + z * z) / (4.0 * D));
           idx += 1;
         }
@@ -102,16 +104,17 @@ class Diffusion : public Model {
 
     if (is_rank0()) std::cout << "Prepare operators" << std::endl;
     size_t idx = 0;
-    double fx = 2.0 * PI / (w.dx * w.Lx);
-    double fy = 2.0 * PI / (w.dy * w.Ly);
-    double fz = 2.0 * PI / (w.dz * w.Lz);
+    auto spacing = w.spacing();
+    auto size = w.size();
+    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]);
     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++) {
-          // Laplacian operator -k^2
-          double ki = (i <= w.Lx / 2) ? i * fx : (i - w.Lx) * fx;
-          double kj = (j <= w.Ly / 2) ? j * fy : (j - w.Ly) * fy;
-          double kk = (k <= w.Lz / 2) ? k * fz : (k - w.Lz) * fz;
+          double ki = (i <= size[0] / 2) ? i * fx : (i - size[0]) * fx;
+          double kj = (j <= size[1] / 2) ? j * fy : (j - size[1]) * fy;
+          double kk = (k <= size[2] / 2) ? k * fz : (k - size[2]) * fz;
           double kLap = -(ki * ki + kj * kj + kk * kk);
           opL[idx++] = 1.0 / (1.0 - dt * kLap);
         }
@@ -180,8 +183,9 @@ void run() {
   World world(dimensions, origin, discretization);
 
   Decomposition decomp(world);
-  FFT fft(decomp);
-  Diffusion model(fft);
+  auto plan_options = heffte::default_options<heffte::backend::fftw>();
+  FFT fft(decomp, MPI_COMM_WORLD, plan_options, world);
+  Diffusion model(world);
 
   // Define time
   double t0 = 0.0;

examples/08_discrete_fields.cpp

@@ -2,8 +2,8 @@
 // SPDX-License-Identifier: AGPL-3.0-or-later
 
 #include <iostream>
+#include <openpfc/core/decomposition.hpp>
 #include <openpfc/core/world.hpp>
-#include <openpfc/decomposition.hpp>
 #include <openpfc/discrete_field.hpp>
 #include <openpfc/utils.hpp>
 
@@ -48,14 +48,21 @@ int main() {
   Decomposition decomp3(world, 2, 4);
   Decomposition decomp4(world, 3, 4);
   std::cout << decomp1 << std::endl;
-  DiscreteField<double, 3> field1(decomp1.get_inbox_size(), decomp1.get_inbox_offset(), world.get_origin(),
-                                  world.get_discretization());
-  DiscreteField<double, 3> field2(decomp2.get_inbox_size(), decomp2.get_inbox_offset(), world.get_origin(),
-                                  world.get_discretization());
-  DiscreteField<double, 3> field3(decomp3.get_inbox_size(), decomp3.get_inbox_offset(), world.get_origin(),
-                                  world.get_discretization());
-  DiscreteField<double, 3> field4(decomp4.get_inbox_size(), decomp4.get_inbox_offset(), world.get_origin(),
-                                  world.get_discretization());
+  const std::array<int, 3> &inbox_size1 = decomp1.get_inbox_size();
+  const std::array<int, 3> &inbox_offset1 = decomp1.get_inbox_offset();
+  DiscreteField<double, 3> field1(inbox_size1, inbox_offset1, world.origin(), world.spacing());
+
+  const std::array<int, 3> &inbox_size2 = decomp2.get_inbox_size();
+  const std::array<int, 3> &inbox_offset2 = decomp2.get_inbox_offset();
+  DiscreteField<double, 3> field2(inbox_size2, inbox_offset2, world.origin(), world.spacing());
+
+  const std::array<int, 3> &inbox_size3 = decomp3.get_inbox_size();
+  const std::array<int, 3> &inbox_offset3 = decomp3.get_inbox_offset();
+  DiscreteField<double, 3> field3(inbox_size3, inbox_offset3, world.origin(), world.spacing());
+
+  const std::array<int, 3> &inbox_size4 = decomp4.get_inbox_size();
+  const std::array<int, 3> &inbox_offset4 = decomp4.get_inbox_offset();
+  DiscreteField<double, 3> field4(inbox_size4, inbox_offset4, world.origin(), world.spacing());
   std::cout << field1 << std::endl;
   std::cout << field2 << std::endl;
   std::cout << field3 << std::endl;

examples/09_parallel_fft_high_level.cpp

@@ -35,8 +35,9 @@ int main(int argc, char *argv[]) {
   // Array<double, 3> output2(decomp);
   // std::cout << output2 << std::endl; // this is {4, 3, 2}
 
-  // Create FFT object and perform parallel FFT
-  FFT fft(decomp);
+  auto plan_options = heffte::default_options<heffte::backend::fftw>();
+  FFT fft(decomp, MPI_COMM_WORLD, plan_options, world);
+
   fft.forward(input, output);
 
   // Display results

examples/10_ui_register_ic.cpp

@@ -40,6 +40,15 @@ void from_json(const json &params, MyIC &ic) {
 // Define model
 class MyModel : public Model {
 public:
+  /**
+   * @brief Constructs a MyModel instance with the given World object.
+   *
+   * @param world The World object to initialize the model.
+   */
+  explicit MyModel(const World &world) : Model(world) {
+    // Additional initialization if needed
+  }
+
   void initialize(double) override { std::cout << "initialize()" << std::endl; }
   void step(double) override { std::cout << "MyModel.step()" << std::endl; }
 };

examples/11_write_results.cpp

@@ -26,7 +26,7 @@ int main(int argc, char **argv) {
   writer.set_field_name("density");
   writer.set_domain(world.get_size(), field.get_size(), field.get_offset());
   writer.set_origin(world.get_origin());
-  writer.set_spacing(world.get_discretization());
+  writer.set_spacing(world.spacing());
   std::cout << "Writing results to file: " << writer.get_uri() << "\n";
   writer.initialize();
   writer.write(field.get_array().get_data());

examples/12_cahn_hilliard.cpp

@@ -19,6 +19,15 @@ class CahnHilliard : public Model {
   double D = 1.0;                              // Diffusion coefficient
 
 public:
+  /**
+   * @brief Constructs a CahnHilliard instance with the given World object.
+   *
+   * @param world The World object to initialize the model.
+   */
+  explicit CahnHilliard(const World &world) : Model(world) {
+    // Additional initialization if needed
+  }
+
   void initialize(double dt) override {
     FFT &fft = get_fft();
     const Decomposition &decomp = get_decomposition();
@@ -37,16 +46,17 @@ class CahnHilliard : public Model {
     std::array<int, 3> o_high = decomp.outbox.high;
     size_t idx = 0;
     double pi = std::atan(1.0) * 4.0;
-    double fx = 2.0 * pi / (w.dx * w.Lx);
-    double fy = 2.0 * pi / (w.dy * w.Ly);
-    double fz = 2.0 * pi / (w.dz * w.Lz);
+    auto spacing = w.spacing();
+    auto size = w.size();
+    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]);
     for (int k = o_low[2]; k <= o_high[2]; k++) {
       for (int j = o_low[1]; j <= o_high[1]; j++) {
         for (int i = o_low[0]; i <= o_high[0]; i++) {
-          // Laplacian operator -k^2
-          double ki = (i <= w.Lx / 2) ? i * fx : (i - w.Lx) * fx;
-          double kj = (j <= w.Ly / 2) ? j * fy : (j - w.Ly) * fy;
-          double kk = (k <= w.Lz / 2) ? k * fz : (k - w.Lz) * fz;
+          double ki = (i <= size[0] / 2) ? i * fx : (i - size[0]) * fx;
+          double kj = (j <= size[1] / 2) ? j * fy : (j - size[1]) * fy;
+          double kk = (k <= size[2] / 2) ? k * fz : (k - size[2]) * fz;
           double kLap = -(ki * ki + kj * kj + kk * kk);
           double L = kLap * (-D - D * gamma * kLap);
           opL[idx] = std::exp(L * dt);
@@ -104,8 +114,9 @@ int main(int argc, char **argv) {
   // Construct world, decomposition, fft and model
   World world({Lx, Ly, Lz}, {x0, y0, z0}, {dx, dy, dz});
   Decomposition decomp(world);
-  FFT fft(decomp);
-  CahnHilliard model;
+  auto plan_options = heffte::default_options<heffte::backend::fftw>();
+  FFT fft(decomp, MPI_COMM_WORLD, plan_options, world);
+  CahnHilliard model(world);
   model.set_fft(fft);
 
   // Define time
@@ -129,9 +140,9 @@ int main(int argc, char **argv) {
   // set uri as format cahn_hilliard_%04i.vti, where %04i is replaced by file_count
   writer.set_uri(sprintf("cahn_hilliard_%04i.vti", file_count));
   writer.set_field_name("concentration");
-  writer.set_domain(world.get_size(), decomp.get_inbox_size(), decomp.get_inbox_offset());
-  writer.set_origin(world.get_origin());
-  writer.set_spacing(world.get_discretization());
+  writer.set_domain(world.size(), decomp.get_inbox_size(), decomp.get_inbox_offset());
+  writer.set_origin(world.origin());
+  writer.set_spacing(world.spacing());
   writer.initialize();
   writer.write(field);