split fluid_solver_base from fluid_solver

now it will be easier to handle functions like correlations etc.

makefile

@@ -62,6 +62,7 @@ base_files := \
 	Morton_shuffler \
 	p3DFFT_to_iR \
 	vector_field \
+	fluid_solver_base \
 	fluid_solver
 
 #headers := $(patsubst %, ./src/%.hpp, ${base_files})

src/fluid_solver.cpp

@@ -25,47 +25,6 @@
 #include "fftw_tools.hpp"
 
 
-
-/*****************************************************************************/
-/* macros for loops                                                          */
-
-/* Fourier space loop */
-#define CLOOP(expression) \
- \
-{ \
-    int cindex; \
-    for (int yindex = 0; yindex < this->cd->subsizes[0]; yindex++) \
-    for (int zindex = 0; zindex < this->cd->subsizes[1]; zindex++) \
-    { \
-        cindex = (yindex * this->cd->sizes[1] + zindex)*this->cd->sizes[2]; \
-    for (int xindex = 0; xindex < this->cd->subsizes[2]; xindex++) \
-        { \
-            expression; \
-            cindex++; \
-        } \
-    } \
-}
-
-/* real space loop */
-#define RLOOP(expression) \
- \
-{ \
-    int rindex; \
-    for (int zindex = 0; zindex < this->rd->subsizes[0]; zindex++) \
-    for (int yindex = 0; yindex < this->rd->subsizes[1]; yindex++) \
-    { \
-        rindex = (zindex * this->rd->sizes[1] + yindex)*this->rd->sizes[2]; \
-    for (int xindex = 0; xindex < this->rd->subsizes[2]; xindex++) \
-        { \
-            expression; \
-            rindex++; \
-        } \
-    } \
-}
-
-/*****************************************************************************/
-
-
 /*****************************************************************************/
 /* macro for specializations to numeric types compatible with FFTW           */
 
@@ -78,20 +37,9 @@ fluid_solver<R>::fluid_solver( \
         int nz, \
         double DKX, \
         double DKY, \
-        double DKZ) \
+        double DKZ) : fluid_solver_base<R>(nx , ny , nz, \
+                                           DKX, DKY, DKZ) \
 { \
-    int ntmp[4]; \
-    ntmp[0] = nz; \
-    ntmp[1] = ny; \
-    ntmp[2] = nx; \
-    ntmp[3] = 3; \
-    this->rd = new field_descriptor<R>( \
-            4, ntmp, MPI_RNUM, MPI_COMM_WORLD);\
-    ntmp[0] = ny; \
-    ntmp[1] = nz; \
-    ntmp[2] = nx/2 + 1; \
-    this->cd = new field_descriptor<R>( \
-            4, ntmp, MPI_CNUM, MPI_COMM_WORLD);\
     this->cvorticity = FFTW(alloc_complex)(this->cd->local_size);\
     this->cvelocity  = FFTW(alloc_complex)(this->cd->local_size);\
     this->rvorticity = FFTW(alloc_real)(this->cd->local_size*2);\
@@ -172,70 +120,11 @@ fluid_solver<R>::fluid_solver( \
     /* ``physical'' parameters etc */ \
  \
     this->nu = 0.1; \
-    this->dkx = DKX; \
-    this->dky = DKY; \
-    this->dkz = DKZ; \
-    this->kx = new double[this->cd->sizes[2]]; \
-    this->ky = new double[this->cd->subsizes[0]]; \
-    this->kz = new double[this->cd->sizes[1]]; \
-    this->knullx = new bool[this->cd->sizes[2]]; \
-    this->knully = new bool[this->cd->subsizes[0]]; \
-    this->knullz = new bool[this->cd->sizes[1]]; \
-    this->nonzerokx = int(this->rd->sizes[2] / 3); \
-    this->kMx = this->dkx*(this->nonzerokx-1); \
-    this->nonzeroky = int(this->rd->sizes[1] / 3); \
-    this->kMy = this->dky*(this->nonzeroky-1); \
-    this->nonzeroky = 2*this->nonzeroky - 1; \
-    this->nonzerokz = int(this->rd->sizes[0] / 3); \
-    this->kMz = this->dkz*(this->nonzerokz-1); \
-    this->nonzerokz = 2*this->nonzerokz - 1; \
-    int i, ii; \
-    for (i = 0; i<this->cd->sizes[2]; i++) \
-    { \
-        this->kx[i] = i*this->dkx; \
-        if (i < this->nonzerokx) \
-            this->knullx[i] = false; \
-        else \
-            this->knullx[i] = true; \
-    } \
-    for (i = 0; i<this->cd->subsizes[0]; i++) \
-    { \
-        int tval = (this->nonzeroky+1)/2; \
-        ii = i + this->cd->starts[0]; \
-        if (ii <= this->rd->sizes[1]/2) \
-            this->ky[i] = this->dky*ii; \
-        else \
-            this->ky[i] = this->dky*(ii - this->rd->sizes[1]); \
-        if (ii < tval || (this->rd->sizes[1] - ii) < tval) \
-            this->knully[i] = false; \
-        else \
-            this->knully[i] = true; \
-    } \
-    for (i = 0; i<this->cd->sizes[1]; i++) \
-    { \
-        int tval = (this->nonzerokz+1)/2; \
-        if (i <= this->rd->sizes[0]/2) \
-            this->kz[i] = this->dkz*i; \
-        else \
-            this->kz[i] = this->dkz*(i - this->rd->sizes[0]); \
-        if (i < tval || (this->rd->sizes[0] - i) < tval) \
-            this->knullz[i] = false; \
-        else \
-            this->knullz[i] = true; \
-    } \
 } \
  \
 template<> \
 fluid_solver<R>::~fluid_solver() \
 { \
- \
-    delete this->kx;\
-    delete this->ky;\
-    delete this->kz;\
-    delete this->knullx;\
-    delete this->knully;\
-    delete this->knullz;\
- \
     FFTW(destroy_plan)(*(FFTW(plan)*)this->c2r_vorticity);\
     FFTW(destroy_plan)(*(FFTW(plan)*)this->r2c_vorticity);\
     FFTW(destroy_plan)(*(FFTW(plan)*)this->c2r_velocity );\
@@ -260,9 +149,6 @@ fluid_solver<R>::~fluid_solver() \
     FFTW(free)(this->rvorticity);\
     FFTW(free)(this->cvelocity);\
     FFTW(free)(this->rvelocity);\
- \
-    delete this->cd; \
-    delete this->rd; \
 } \
  \
 template<> \
@@ -371,12 +257,14 @@ void fluid_solver<R>::step(double dt) \
 
 /*****************************************************************************/
 /* now actually use the macro defined above                                  */
-FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_FLOAT,
+FLUID_SOLVER_DEFINITIONS(
+        FFTW_MANGLE_FLOAT,
         float,
         fftwf_complex,
         MPI_REAL4,
         MPI_COMPLEX8)
-FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_DOUBLE,
+FLUID_SOLVER_DEFINITIONS(
+        FFTW_MANGLE_DOUBLE,
         double,
         fftw_complex,
         MPI_REAL8,

src/fluid_solver.hpp

@@ -23,6 +23,7 @@
 #include <iostream>
 #include "field_descriptor.hpp"
 #include "vector_field.hpp"
+#include "fluid_solver_base.hpp"
 
 #ifndef FLUID_SOLVER
 
@@ -37,22 +38,18 @@ extern int myrank, nprocs;
  * */
 
 template <class rnumber>
-class fluid_solver
+class fluid_solver:public fluid_solver_base<rnumber>
 {
-    private:
-        typedef rnumber cnumber[2];
     public:
-        field_descriptor<rnumber> *cd, *rd;
-
         /* fields */
         rnumber *rvorticity;
         rnumber *rvelocity ;
-        cnumber *cvorticity;
-        cnumber *cvelocity ;
+        typename fluid_solver_base<rnumber>::cnumber *cvorticity;
+        typename fluid_solver_base<rnumber>::cnumber *cvelocity ;
 
         /* short names for velocity, and 4 vorticity fields */
         rnumber *ru, *rv[4];
-        cnumber *cu, *cv[4];
+        typename fluid_solver_base<rnumber>::cnumber *cu, *cv[4];
 
         /* plans */
         void *c2r_vorticity;
@@ -62,22 +59,8 @@ class fluid_solver
         void *uc2r, *ur2c;
         void *vr2c[3], *vc2r[3];
 
-        /* simulation parameters */
-        int iteration;
-
         /* physical parameters */
         rnumber nu;
-        rnumber dkx, dky, dkz, dk;
-
-        /* mode and dealiasing information */
-        double kMx, kMy, kMz, kM, kM2;
-        double *kx, *ky, *kz;
-        bool *knullx, *knully, *knullz;
-        int nonzerokx, nonzeroky, nonzerokz;
-        double *kshell;
-        int64_t *nshell;
-        int nshells;
-
 
         /* methods */
         fluid_solver(
@@ -87,9 +70,6 @@ class fluid_solver
                 double DKX = 1.0,
                 double DKY = 1.0,
                 double DKZ = 1.0);
-
-
-
         ~fluid_solver();
 
         void omega_nonlin(int src);

src/fluid_solver_base.cpp

+/***********************************************************************
+*
+*  Copyright 2015 Max Planck Institute for Dynamics and SelfOrganization
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*
+* Contact: Cristian.Lalescu@ds.mpg.de
+*
+************************************************************************/
+
+
+#include <cassert>
+#include <cmath>
+#include "fluid_solver_base.hpp"
+#include "fftw_tools.hpp"
+
+
+
+/*****************************************************************************/
+/* macro for specializations to numeric types compatible with FFTW           */
+
+#define FLUID_SOLVER_BASE_DEFINITIONS(FFTW, R, C, MPI_RNUM, MPI_CNUM) \
+ \
+template<> \
+fluid_solver_base<R>::fluid_solver_base( \
+        int nx, \
+        int ny, \
+        int nz, \
+        double DKX, \
+        double DKY, \
+        double DKZ) \
+{ \
+    int ntmp[4]; \
+    ntmp[0] = nz; \
+    ntmp[1] = ny; \
+    ntmp[2] = nx; \
+    ntmp[3] = 3; \
+    this->rd = new field_descriptor<R>( \
+            4, ntmp, MPI_RNUM, MPI_COMM_WORLD);\
+    ntmp[0] = ny; \
+    ntmp[1] = nz; \
+    ntmp[2] = nx/2 + 1; \
+    this->cd = new field_descriptor<R>( \
+            4, ntmp, MPI_CNUM, MPI_COMM_WORLD);\
+ \
+    this->dkx = DKX; \
+    this->dky = DKY; \
+    this->dkz = DKZ; \
+    this->kx = new double[this->cd->sizes[2]]; \
+    this->ky = new double[this->cd->subsizes[0]]; \
+    this->kz = new double[this->cd->sizes[1]]; \
+    this->knullx = new bool[this->cd->sizes[2]]; \
+    this->knully = new bool[this->cd->subsizes[0]]; \
+    this->knullz = new bool[this->cd->sizes[1]]; \
+    this->nonzerokx = int(this->rd->sizes[2] / 3); \
+    this->kMx = this->dkx*(this->nonzerokx-1); \
+    this->nonzeroky = int(this->rd->sizes[1] / 3); \
+    this->kMy = this->dky*(this->nonzeroky-1); \
+    this->nonzeroky = 2*this->nonzeroky - 1; \
+    this->nonzerokz = int(this->rd->sizes[0] / 3); \
+    this->kMz = this->dkz*(this->nonzerokz-1); \
+    this->nonzerokz = 2*this->nonzerokz - 1; \
+    int i, ii; \
+    for (i = 0; i<this->cd->sizes[2]; i++) \
+    { \
+        this->kx[i] = i*this->dkx; \
+        if (i < this->nonzerokx) \
+            this->knullx[i] = false; \
+        else \
+            this->knullx[i] = true; \
+    } \
+    for (i = 0; i<this->cd->subsizes[0]; i++) \
+    { \
+        int tval = (this->nonzeroky+1)/2; \
+        ii = i + this->cd->starts[0]; \
+        if (ii <= this->rd->sizes[1]/2) \
+            this->ky[i] = this->dky*ii; \
+        else \
+            this->ky[i] = this->dky*(ii - this->rd->sizes[1]); \
+        if (ii < tval || (this->rd->sizes[1] - ii) < tval) \
+            this->knully[i] = false; \
+        else \
+            this->knully[i] = true; \
+    } \
+    for (i = 0; i<this->cd->sizes[1]; i++) \
+    { \
+        int tval = (this->nonzerokz+1)/2; \
+        if (i <= this->rd->sizes[0]/2) \
+            this->kz[i] = this->dkz*i; \
+        else \
+            this->kz[i] = this->dkz*(i - this->rd->sizes[0]); \
+        if (i < tval || (this->rd->sizes[0] - i) < tval) \
+            this->knullz[i] = false; \
+        else \
+            this->knullz[i] = true; \
+    } \
+} \
+ \
+template<> \
+fluid_solver_base<R>::~fluid_solver_base() \
+{ \
+ \
+    delete this->kx;\
+    delete this->ky;\
+    delete this->kz;\
+    delete this->knullx;\
+    delete this->knully;\
+    delete this->knullz;\
+ \
+    delete this->cd; \
+    delete this->rd; \
+}
+
+/*****************************************************************************/
+
+
+
+/*****************************************************************************/
+/* now actually use the macro defined above                                  */
+FLUID_SOLVER_BASE_DEFINITIONS(
+        FFTW_MANGLE_FLOAT,
+        float,
+        fftwf_complex,
+        MPI_REAL4,
+        MPI_COMPLEX8)
+FLUID_SOLVER_BASE_DEFINITIONS(
+        FFTW_MANGLE_DOUBLE,
+        double,
+        fftw_complex,
+        MPI_REAL8,
+        MPI_COMPLEX16)
+/*****************************************************************************/
+
+
+
+/*****************************************************************************/
+/* finally, force generation of code for single precision                    */
+template class fluid_solver_base<float>;
+/*****************************************************************************/
+

src/fluid_solver_base.hpp

+/***********************************************************************
+*
+*  Copyright 2015 Max Planck Institute for Dynamics and SelfOrganization
+*
+* Licensed under the Apache License, Version 2.0 (the "License");
+* you may not use this file except in compliance with the License.
+* You may obtain a copy of the License at
+*
+*     http://www.apache.org/licenses/LICENSE-2.0
+*
+* Unless required by applicable law or agreed to in writing, software
+* distributed under the License is distributed on an "AS IS" BASIS,
+* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+* See the License for the specific language governing permissions and
+* limitations under the License.
+*
+* Contact: Cristian.Lalescu@ds.mpg.de
+*
+************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include "field_descriptor.hpp"
+
+#ifndef FLUID_SOLVER_BASE
+
+#define FLUID_SOLVER_BASE
+
+extern int myrank, nprocs;
+
+
+/* container for field descriptor, fields themselves, parameters, etc
+ * using the same big macro idea that they're using in fftw3.h
+ * I feel like I should quote:  Ugh.
+ * */
+
+template <class rnumber>
+class fluid_solver_base
+{
+    protected:
+        typedef rnumber cnumber[2];
+    public:
+        field_descriptor<rnumber> *cd, *rd;
+
+        /* simulation parameters */
+        int iteration;
+
+        /* physical parameters */
+        rnumber dkx, dky, dkz, dk;
+
+        /* mode and dealiasing information */
+        double kMx, kMy, kMz, kM, kM2;
+        double *kx, *ky, *kz;
+        bool *knullx, *knully, *knullz;
+        int nonzerokx, nonzeroky, nonzerokz;
+        double *kshell;
+        int64_t *nshell;
+        int nshells;
+
+
+        /* methods */
+        fluid_solver_base(
+                int nx,
+                int ny,
+                int nz,
+                double DKX = 1.0,
+                double DKY = 1.0,
+                double DKZ = 1.0);
+        ~fluid_solver_base();
+};
+
+
+
+/*****************************************************************************/
+/* macros for loops                                                          */
+
+/* Fourier space loop */
+#define CLOOP(expression) \
+ \
+{ \
+    int cindex; \
+    for (int yindex = 0; yindex < this->cd->subsizes[0]; yindex++) \
+    for (int zindex = 0; zindex < this->cd->subsizes[1]; zindex++) \
+    { \
+        cindex = (yindex * this->cd->sizes[1] + zindex)*this->cd->sizes[2]; \
+    for (int xindex = 0; xindex < this->cd->subsizes[2]; xindex++) \
+        { \
+            expression; \
+            cindex++; \
+        } \
+    } \
+}
+
+/* real space loop */
+#define RLOOP(expression) \
+ \
+{ \
+    int rindex; \
+    for (int zindex = 0; zindex < this->rd->subsizes[0]; zindex++) \
+    for (int yindex = 0; yindex < this->rd->subsizes[1]; yindex++) \
+    { \
+        rindex = (zindex * this->rd->sizes[1] + yindex)*this->rd->sizes[2]; \
+    for (int xindex = 0; xindex < this->rd->subsizes[2]; xindex++) \
+        { \
+            expression; \
+            rindex++; \
+        } \
+    } \
+}
+
+/*****************************************************************************/
+
+#endif//FLUID_SOLVER_BASE
+