From efb8a0abfa9ecaf794a72121bd03e3c1b57f1707 Mon Sep 17 00:00:00 2001
From: Chichi Lalescu <clalesc1@jhu.edu>
Date: Fri, 19 Jun 2015 21:28:30 +0200
Subject: [PATCH] add preliminary correlation function

---
 py_template/test_step.cpp |  6 ++----
 src/fluid_solver.cpp      |  1 +
 src/fluid_solver_base.cpp | 23 +++++++++++++++++++++++
 src/fluid_solver_base.hpp |  2 ++
 test.py                   | 16 +---------------
 5 files changed, 29 insertions(+), 19 deletions(-)

diff --git a/py_template/test_step.cpp b/py_template/test_step.cpp
index f76cda73..deca5c8c 100644
--- a/py_template/test_step.cpp
+++ b/py_template/test_step.cpp
@@ -6,11 +6,9 @@ fs->cd->read(
         "Kdata0",
         (void*)fs->cvorticity);
 DEBUG_MSG("field read\n");
+DEBUG_MSG("######### %d %g\n", fs->iteration, fs->correl_vec(fs->cvorticity, fs->cvorticity));
 fs->step(0.0001);
-DEBUG_MSG("after time step\n");
-fs->cd->write(
-        "Kdata1",
-        (void*)fs->cvorticity);
+DEBUG_MSG("######### %d %g\n", fs->iteration, fs->correl_vec(fs->cvorticity, fs->cvorticity));
 
 delete fs;
 DEBUG_MSG("fluid_solver object deleted\n");
diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp
index 88c522af..fbfe6ce5 100644
--- a/src/fluid_solver.cpp
+++ b/src/fluid_solver.cpp
@@ -250,6 +250,7 @@ void fluid_solver<R>::step(double dt) \
             this->cv[3][cindex*3+2][1] = (this->cv[0][cindex*3+2][1]*factor0 + 2*(this->cv[1][cindex*3+2][1] + dt*this->cu[cindex*3+2][1])*factor1)/3; \
             ); \
  \
+    this->iteration++; \
 }
 /*****************************************************************************/
 
diff --git a/src/fluid_solver_base.cpp b/src/fluid_solver_base.cpp
index fadf5d79..645aa675 100644
--- a/src/fluid_solver_base.cpp
+++ b/src/fluid_solver_base.cpp
@@ -40,6 +40,8 @@ fluid_solver_base<R>::fluid_solver_base( \
         double DKY, \
         double DKZ) \
 { \
+    this->iteration = 0; \
+ \
     int ntmp[4]; \
     ntmp[0] = nz; \
     ntmp[1] = ny; \
@@ -119,6 +121,27 @@ fluid_solver_base<R>::~fluid_solver_base() \
  \
     delete this->cd; \
     delete this->rd; \
+} \
+ \
+template<> \
+R fluid_solver_base<R>::correl_vec(C *a, C *b) \
+{ \
+    double val_process = 0.0, val; \
+    double k2; \
+    int factor = 1;\
+    CLOOP( \
+            k2 = (this->kx[xindex]*this->kx[xindex] + \
+                  this->ky[yindex]*this->ky[yindex] + \
+                  this->kz[zindex]*this->kz[zindex]); \
+            factor = (xindex == 0) ? 1 : 2; \
+            val_process += factor * ((*(a + cindex))[0] * (*(b + cindex))[0] + \
+                                     (*(a + cindex))[1] * (*(b + cindex))[1]); \
+            );\
+    MPI_Allreduce( \
+            (void*)(&val_process), \
+            (void*)(&val), \
+            1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD); \
+    return R(val); \
 }
 
 /*****************************************************************************/
diff --git a/src/fluid_solver_base.hpp b/src/fluid_solver_base.hpp
index dd284139..f0001e53 100644
--- a/src/fluid_solver_base.hpp
+++ b/src/fluid_solver_base.hpp
@@ -68,6 +68,8 @@ class fluid_solver_base
                 double DKY = 1.0,
                 double DKZ = 1.0);
         ~fluid_solver_base();
+
+        rnumber correl_vec(cnumber *a, cnumber *b);
 };
 
 
diff --git a/test.py b/test.py
index b2bdd40f..48432b61 100755
--- a/test.py
+++ b/test.py
@@ -60,19 +60,5 @@ Kdata0[..., 0] = Kdata00
 Kdata0[..., 1] = Kdata01
 Kdata0[..., 2] = Kdata02
 Kdata0.tofile("Kdata0")
-run_test('test_step')
-Kdata1 = np.fromfile('Kdata1', dtype = np.complex64).reshape(Kdata0.shape)
-
-#print np.max(np.abs(Kdata0 - Kdata1))
-print np.max(np.abs(Kdata0))
-print np.max(np.abs(Kdata1))
-#
-#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')
-
-
+run_test('test_step', ncpu = 1)
 
-- 
GitLab