allocate plans and memory

src/fluid_solver.cpp

@@ -21,44 +21,86 @@
 
 
 #include "fluid_solver.hpp"
-
-///*****************************************************************************/
-///* generic definitions, empty by design                                      */
-//
-//template <class rnumber>
-//fluid_solver<rnumber>::fluid_solver()
-//{}
-//template <class rnumber>
-//fluid_solver<rnumber>::~fluid_solver()
-//{}
-//template <class rnumber>
-//void fluid_solver<rnumber>::step()
-//{}
-///*****************************************************************************/
+#include "fftwf_tools.hpp"
 
 
 
 /*****************************************************************************/
 /* macro for specializations to numeric types compatible with FFTW           */
 
-#define FLUID_SOLVER_DEFINITIONS(X, R, C) \
+#define FLUID_SOLVER_DEFINITIONS(FFTW, R, C) \
  \
 template<> \
-fluid_solver<R>::fluid_solver() \
+fluid_solver<R>::fluid_solver( \
+        int nx, \
+        int ny, \
+        int nz) \
 { \
-    this->c2r_vorticity = new X(plan);\
-    this->r2c_vorticity = new X(plan);\
-    this->c2r_velocity  = new X(plan);\
-    this->r2c_velocity  = new X(plan);\
+    FFTW(get_descriptors_3D)(nz, ny, nx, &this->fr, &this->fc);\
+    this->cvorticity = FFTW(alloc_complex)(this->fc->local_size*3);\
+    this->cvelocity  = FFTW(alloc_complex)(this->fc->local_size*3);\
+    this->rvorticity = FFTW(alloc_real)(this->fc->local_size*6);\
+    this->rvelocity  = FFTW(alloc_real)(this->fc->local_size*6);\
+ \
+    this->c2r_vorticity = new FFTW(plan);\
+    this->r2c_vorticity = new FFTW(plan);\
+    this->c2r_velocity  = new FFTW(plan);\
+    this->r2c_velocity  = new FFTW(plan);\
+ \
+    ptrdiff_t sizes[] = {this->fr->sizes[0], \
+                         this->fr->sizes[1], \
+                         this->fr->sizes[2]};\
+ \
+    *(FFTW(plan)*)this->c2r_vorticity = FFTW(mpi_plan_many_dft_c2r)( \
+            3, sizes, 3, \
+            FFTW_MPI_DEFAULT_BLOCK, \
+            FFTW_MPI_DEFAULT_BLOCK, \
+            this->cvorticity, this->rvorticity, \
+            MPI_COMM_WORLD, \
+            FFTW_ESTIMATE); \
+ \
+    *(FFTW(plan)*)this->r2c_vorticity = FFTW(mpi_plan_many_dft_r2c)( \
+            3, sizes, 3, \
+            FFTW_MPI_DEFAULT_BLOCK, \
+            FFTW_MPI_DEFAULT_BLOCK, \
+            this->rvorticity, this->cvorticity, \
+            MPI_COMM_WORLD, \
+            FFTW_ESTIMATE); \
+ \
+    *(FFTW(plan)*)this->c2r_velocity = FFTW(mpi_plan_many_dft_c2r)( \
+            3, sizes, 3, \
+            FFTW_MPI_DEFAULT_BLOCK, \
+            FFTW_MPI_DEFAULT_BLOCK, \
+            this->cvelocity, this->rvelocity, \
+            MPI_COMM_WORLD, \
+            FFTW_ESTIMATE); \
+ \
+    *(FFTW(plan)*)this->r2c_velocity = FFTW(mpi_plan_many_dft_r2c)( \
+            3, sizes, 3, \
+            FFTW_MPI_DEFAULT_BLOCK, \
+            FFTW_MPI_DEFAULT_BLOCK, \
+            this->rvelocity, this->cvelocity, \
+            MPI_COMM_WORLD, \
+            FFTW_ESTIMATE); \
 } \
  \
 template<> \
 fluid_solver<R>::~fluid_solver() \
 { \
-    delete (X(plan)*)this->c2r_vorticity;\
-    delete (X(plan)*)this->r2c_vorticity;\
-    delete (X(plan)*)this->c2r_velocity ;\
-    delete (X(plan)*)this->r2c_velocity ;\
+    FFTW(destroy_plan)(*(FFTW(plan)*)this->c2r_vorticity);\
+    FFTW(destroy_plan)(*(FFTW(plan)*)this->r2c_vorticity);\
+    FFTW(destroy_plan)(*(FFTW(plan)*)this->c2r_velocity );\
+    FFTW(destroy_plan)(*(FFTW(plan)*)this->r2c_velocity );\
+ \
+    delete (FFTW(plan)*)this->c2r_vorticity;\
+    delete (FFTW(plan)*)this->r2c_vorticity;\
+    delete (FFTW(plan)*)this->c2r_velocity ;\
+    delete (FFTW(plan)*)this->r2c_velocity ;\
+ \
+    FFTW(free)(this->cvorticity);\
+    FFTW(free)(this->rvorticity);\
+    FFTW(free)(this->cvelocity);\
+    FFTW(free)(this->rvelocity);\
 } \
  \
 template<> \
@@ -70,10 +112,10 @@ void fluid_solver<R>::step() \
 
 /*****************************************************************************/
 /* now actually use the macro defined above                                  */
-FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_DOUBLE, double, fftw_complex)
 FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_FLOAT, float, fftwf_complex)
-FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_LONG_DOUBLE, long double, fftwl_complex)
-FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_QUAD, __float128, fftwq_complex)
+FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_DOUBLE, double, fftw_complex)
+//FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_LONG_DOUBLE, long double, fftwl_complex)
+//FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_QUAD, __float128, fftwq_complex)
 /*****************************************************************************/
 
 

src/fluid_solver.hpp

@@ -56,7 +56,13 @@ class fluid_solver
         void *c2r_velocity;
         void *r2c_velocity;
 
-        fluid_solver();
+        /* parameters */
+        rnumber nu;
+
+        /* methods */
+
+        fluid_solver(
+                int nx, int ny, int nz);
         ~fluid_solver();
 
         void step();