From 8d35451e49a255d69f2c24fde1bcfd38088ebfb0 Mon Sep 17 00:00:00 2001
From: Chichi Lalescu <clalesc1@jhu.edu>
Date: Fri, 19 Jun 2015 21:03:03 +0200
Subject: [PATCH] split fluid_solver_base from fluid_solver
now it will be easier to handle functions like correlations etc.
---
makefile | 1 +
src/fluid_solver.cpp | 140 ++++-------------------------------
src/fluid_solver.hpp | 30 ++------
src/fluid_solver_base.cpp | 150 ++++++++++++++++++++++++++++++++++++++
src/fluid_solver_base.hpp | 115 +++++++++++++++++++++++++++++
5 files changed, 285 insertions(+), 151 deletions(-)
create mode 100644 src/fluid_solver_base.cpp
create mode 100644 src/fluid_solver_base.hpp
diff --git a/makefile b/makefile
index 3c163b9f..89b34da5 100644
--- a/makefile
+++ b/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})
diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp
index c6fde9c5..88c522af 100644
--- a/src/fluid_solver.cpp
+++ b/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,16 +257,18 @@ void fluid_solver<R>::step(double dt) \
/*****************************************************************************/
/* now actually use the macro defined above */
-FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_FLOAT,
- float,
- fftwf_complex,
- MPI_REAL4,
- MPI_COMPLEX8)
-FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_DOUBLE,
- double,
- fftw_complex,
- MPI_REAL8,
- MPI_COMPLEX16)
+FLUID_SOLVER_DEFINITIONS(
+ FFTW_MANGLE_FLOAT,
+ float,
+ fftwf_complex,
+ MPI_REAL4,
+ MPI_COMPLEX8)
+FLUID_SOLVER_DEFINITIONS(
+ FFTW_MANGLE_DOUBLE,
+ double,
+ fftw_complex,
+ MPI_REAL8,
+ MPI_COMPLEX16)
/*****************************************************************************/
diff --git a/src/fluid_solver.hpp b/src/fluid_solver.hpp
index d189ac4a..12f3a824 100644
--- a/src/fluid_solver.hpp
+++ b/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);
diff --git a/src/fluid_solver_base.cpp b/src/fluid_solver_base.cpp
new file mode 100644
index 00000000..fadf5d79
--- /dev/null
+++ b/src/fluid_solver_base.cpp
@@ -0,0 +1,150 @@
+/***********************************************************************
+*
+* 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>;
+/*****************************************************************************/
+
diff --git a/src/fluid_solver_base.hpp b/src/fluid_solver_base.hpp
new file mode 100644
index 00000000..dd284139
--- /dev/null
+++ b/src/fluid_solver_base.hpp
@@ -0,0 +1,115 @@
+/***********************************************************************
+*
+* 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
+
--
GitLab