diff --git a/makefile b/makefile
index 9eb470ce27518da157a908ec4e0e5455977e8bf2..3c163b9fc93aa176c525e323a442be7673e08cb3 100644
--- a/makefile
+++ b/makefile
@@ -61,6 +61,7 @@ base_files := \
fftw_tools \
Morton_shuffler \
p3DFFT_to_iR \
+ vector_field \
fluid_solver
#headers := $(patsubst %, ./src/%.hpp, ${base_files})
diff --git a/py_template/test_FFT.cpp b/py_template/test_FFT.cpp
index 34d6f95453898b065654cea34a5f67f5a8dfbd0e..23d7e104d5b2456c5187c811f923aefd1eceafcb 100644
--- a/py_template/test_FFT.cpp
+++ b/py_template/test_FFT.cpp
@@ -2,15 +2,22 @@ fluid_solver<float> *fs;
fs = new fluid_solver<float>(32, 32, 32);
DEBUG_MSG("fluid_solver object created\n");
-fs->fc->read(
+vector_field<float> cv(fs->cd, fs->cvorticity);
+vector_field<float> rv(fs->cd, fs->rvorticity);
+
+fs->cd->read(
"Kdata0",
(void*)fs->cvorticity);
fftwf_execute(*(fftwf_plan*)fs->c2r_vorticity);
+//rv*(1. / (fs->rd->sizes[0]*fs->rd->sizes[1]*fs->rd->sizes[2]));
fftwf_execute(*(fftwf_plan*)fs->r2c_vorticity);
-fs->fc->write(
+cv = cv*(1. / (fs->rd->sizes[0]*fs->rd->sizes[1]*fs->rd->sizes[2]));
+fs->cd->write(
"Kdata1",
(void*)fs->cvorticity);
+DEBUG_MSG("full size is %ld\n", fs->rd->full_size);
+
delete fs;
DEBUG_MSG("fluid_solver object deleted\n");
diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp
index e24827858a9a4dbb878967886afe32f5b56686d4..d10e6960bb72b301ee9b71bc056b6d87c6894e1c 100644
--- a/src/fluid_solver.cpp
+++ b/src/fluid_solver.cpp
@@ -28,7 +28,7 @@
/*****************************************************************************/
/* macro for specializations to numeric types compatible with FFTW */
-#define FLUID_SOLVER_DEFINITIONS(FFTW, R, C) \
+#define FLUID_SOLVER_DEFINITIONS(FFTW, R, C, MPI_RNUM, MPI_CNUM) \
\
template<> \
fluid_solver<R>::fluid_solver( \
@@ -36,20 +36,29 @@ fluid_solver<R>::fluid_solver( \
int ny, \
int nz) \
{ \
- get_descriptors_3D<R>(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);\
+ 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[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);\
+ this->rvelocity = FFTW(alloc_real)(this->cd->local_size*2);\
\
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]};\
+ ptrdiff_t sizes[] = {nz, \
+ ny, \
+ nx};\
\
*(FFTW(plan)*)this->c2r_vorticity = FFTW(mpi_plan_many_dft_c2r)( \
3, sizes, 3, \
@@ -101,6 +110,9 @@ fluid_solver<R>::~fluid_solver() \
FFTW(free)(this->rvorticity);\
FFTW(free)(this->cvelocity);\
FFTW(free)(this->rvelocity);\
+ \
+ delete this->cd; \
+ delete this->rd; \
} \
\
template<> \
@@ -112,10 +124,16 @@ void fluid_solver<R>::step() \
/*****************************************************************************/
/* now actually use the macro defined above */
-FLUID_SOLVER_DEFINITIONS(FFTW_MANGLE_FLOAT, float, fftwf_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)
+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 06e34240d06419cb510da0b884cf88004cce6ba0..0104e6bf9917ba099826b6c9ff2de7be8de01e77 100644
--- a/src/fluid_solver.hpp
+++ b/src/fluid_solver.hpp
@@ -22,6 +22,7 @@
#include <stdlib.h>
#include <iostream>
#include "field_descriptor.hpp"
+#include "vector_field.hpp"
#ifndef FLUID_SOLVER
@@ -41,7 +42,7 @@ class fluid_solver
private:
typedef rnumber cnumber[2];
public:
- field_descriptor<rnumber> *fc, *fr;
+ field_descriptor<rnumber> *cd, *rd;
/* fields */
rnumber *rvorticity;
diff --git a/src/vector_field.cpp b/src/vector_field.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3a9412522fef80ec024d0cb7ec1c01df97cb4dde
--- /dev/null
+++ b/src/vector_field.cpp
@@ -0,0 +1,42 @@
+#include "vector_field.hpp"
+
+/* destructor doesn't actually do anything */
+template <class rnumber>
+vector_field<rnumber>::~vector_field()
+{}
+
+
+template <class rnumber>
+vector_field<rnumber>::vector_field(
+ field_descriptor<rnumber> *d,
+ rnumber *data)
+{
+ this->is_real = true;
+ this->cdata = (rnumber (*)[2])(data);
+ this->rdata = data;
+ this->descriptor = d;
+}
+
+template <class rnumber>
+vector_field<rnumber>::vector_field(
+ field_descriptor<rnumber> *d,
+ rnumber (*data)[2])
+{
+ this->is_real = false;
+ this->rdata = (rnumber*)(&data[0][0]);
+ this->cdata = data;
+ this->descriptor = d;
+}
+
+template <class rnumber>
+vector_field<rnumber>& vector_field<rnumber>::operator*(rnumber factor)
+{
+ ptrdiff_t i;
+ for (i = 0;
+ i < this->descriptor->local_size * 2;
+ i++)
+ *(this->rdata + i) *= factor;
+ return *this;
+}
+
+template class vector_field<float>;
diff --git a/src/vector_field.hpp b/src/vector_field.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0d8aa70c13cd063fbe8bb01adf2f32bfee52c5ac
--- /dev/null
+++ b/src/vector_field.hpp
@@ -0,0 +1,19 @@
+#include "field_descriptor.hpp"
+
+template <class rnumber>
+class vector_field
+{
+ private:
+ field_descriptor<rnumber> *descriptor;
+ rnumber *rdata;
+ rnumber (*cdata)[2];
+ bool is_real;
+ public:
+ vector_field(field_descriptor<rnumber> *d, rnumber *data);
+ vector_field(field_descriptor<rnumber> *d, rnumber (*data)[2]);
+ ~vector_field();
+
+ /* various operators */
+ vector_field &operator*(rnumber factor);
+};
+
diff --git a/test.py b/test.py
old mode 100644
new mode 100755
index bc3d31ff2443e747768397a246e5183f4419b00a..94cfa417882591616e8b10afccec856a980c94ee
--- a/test.py
+++ b/test.py
@@ -1,6 +1,9 @@
+#! /usr/bin/env python2
+
import numpy as np
import subprocess
import pyfftw
+import matplotlib.pyplot as plt
def run_test(
test_name = 'test_FFT',
@@ -46,13 +49,29 @@ def generate_data_3D(
a[ii] = 0
return a
-Kdata0 = generate_data_3D(32, p = 2).astype(np.complex64)
+n = 32
+Kdata00 = generate_data_3D(n, p = 2).astype(np.complex64)
+Kdata01 = generate_data_3D(n, p = 2).astype(np.complex64)
+Kdata02 = generate_data_3D(n, p = 2).astype(np.complex64)
+Kdata0 = np.zeros(
+ Kdata00.shape + (3,),
+ Kdata00.dtype)
+Kdata0[..., 0] = Kdata00
+Kdata0[..., 1] = Kdata01
+Kdata0[..., 2] = Kdata02
Kdata0.tofile("Kdata0")
run_test('test_FFT')
Kdata1 = np.fromfile('Kdata1', dtype = np.complex64).reshape(Kdata0.shape)
print np.max(np.abs(Kdata0 - Kdata1))
+print np.max(np.abs(Kdata0))
+fig = plt.figure(figsize=(12, 6))
+a = fig.add_subplot(121)
+a.imshow(abs(Kdata0[4, :, :, 2]), interpolation = 'none')
+a = fig.add_subplot(122)
+a.imshow(abs(Kdata1[4, :, :, 2]), interpolation = 'none')
+fig.savefig('tmp.pdf', format = 'pdf')