From 1590ae680cc7a6da8f20d3ad1cb2523258a1680b Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Fri, 17 Jul 2015 15:55:00 +0200
Subject: [PATCH] fix continuation inconsistency

---
 src/field_descriptor.cpp  | 31 ++++++++++++++++++++++---------
 src/fluid_solver.cpp      | 17 +++++++++++++----
 src/fluid_solver_base.cpp | 10 ++++++++--
 src/fluid_solver_base.hpp |  1 +
 4 files changed, 44 insertions(+), 15 deletions(-)

diff --git a/src/field_descriptor.cpp b/src/field_descriptor.cpp
index 23e681a3..65e65b55 100644
--- a/src/field_descriptor.cpp
+++ b/src/field_descriptor.cpp
@@ -91,7 +91,7 @@ field_descriptor<rnumber>::field_descriptor(
     tsizes[ndims-1] *= sizeof(rnumber);
     tsubsizes[ndims-1] *= sizeof(rnumber);
     tstarts[ndims-1] *= sizeof(rnumber);
-    if (this->mpi_dtype == MPI_COMPLEX8 ||
+    if (this->mpi_dtype == MPI_COMPLEX ||
         this->mpi_dtype == MPI_COMPLEX16)
     {
         tsizes[ndims-1] *= 2;
@@ -224,20 +224,28 @@ field_descriptor<rnumber>::~field_descriptor()
     delete[] this->rank;
 }
 
-template<>
-int field_descriptor<float>::read(
+template<class rnumber>
+int field_descriptor<rnumber>::read(
         const char *fname,
         void *buffer)
 {
+    MPI_Datatype ttype;
+    if (sizeof(rnumber)==4)
+        ttype = MPI_COMPLEX;
+    else if (sizeof(rnumber)==8)
+        ttype = MPI_DOUBLE_COMPLEX;
+    DEBUG_MSG("aloha 00\n");
     if (this->subsizes[0] > 0)
     {
         MPI_Info info;
         MPI_Info_create(&info);
         MPI_File f;
-        int read_size = this->local_size*sizeof(float);
+        ptrdiff_t read_size = this->local_size*sizeof(rnumber);
+        DEBUG_MSG("aloha %ld\n", read_size);
         char ffname[200];
-        if (this->mpi_dtype == MPI_COMPLEX8)
+        if (this->mpi_dtype == ttype)
             read_size *= 2;
+        DEBUG_MSG("aloha %ld\n", read_size);
         sprintf(ffname, "%s", fname);
 
         MPI_File_open(
@@ -264,19 +272,24 @@ int field_descriptor<float>::read(
     return EXIT_SUCCESS;
 }
 
-template<>
-int field_descriptor<float>::write(
+template<class rnumber>
+int field_descriptor<rnumber>::write(
         const char *fname,
         void *buffer)
 {
+    MPI_Datatype ttype;
+    if (sizeof(rnumber)==4)
+        ttype = MPI_COMPLEX;
+    else if (sizeof(rnumber)==8)
+        ttype = MPI_DOUBLE_COMPLEX;
     if (this->subsizes[0] > 0)
     {
         MPI_Info info;
         MPI_Info_create(&info);
         MPI_File f;
-        int read_size = this->local_size*sizeof(float);
+        ptrdiff_t read_size = this->local_size*sizeof(rnumber);
         char ffname[200];
-        if (this->mpi_dtype == MPI_COMPLEX8)
+        if (this->mpi_dtype == ttype)
             read_size *= 2;
         sprintf(ffname, "%s", fname);
 
diff --git a/src/fluid_solver.cpp b/src/fluid_solver.cpp
index 9d8618a0..075ba616 100644
--- a/src/fluid_solver.cpp
+++ b/src/fluid_solver.cpp
@@ -389,11 +389,11 @@ void fluid_solver<R>::step(double dt) \
     this->symmetrize(this->cu, 3); \
     this->force_divfree(this->cv[0]); \
     this->low_pass_Fourier(this->cv[0], 3, this->kM); \
-    this->write_base("cvorticity", this->cvorticity); \
+    /*this->write_base("cvorticity", this->cvorticity); \
     this->read_base("cvorticity", this->cvorticity); \
         this->low_pass_Fourier(this->cvorticity, 3, this->kM); \
         this->force_divfree(this->cvorticity); \
-        this->symmetrize(this->cvorticity, 3); \
+        this->symmetrize(this->cvorticity, 3);*/ \
     /*this->ift_vorticity(); \
     this->dft_vorticity(); \
     CLOOP( \
@@ -408,14 +408,19 @@ void fluid_solver<R>::step(double dt) \
 template<> \
 int fluid_solver<R>::read(char field, char representation) \
 { \
+    char fname[512]; \
     int read_result; \
     if (field == 'v') \
     { \
         if (representation == 'c') \
         { \
-            read_result = this->read_base("cvorticity", this->cvorticity); \
+            this->fill_up_filename("cvorticity", fname); \
+            read_result = this->cd->read(fname, (void*)this->cvorticity); \
             if (read_result != EXIT_SUCCESS) \
+            { \
+                DEBUG_MSG("read error"); \
                 return read_result; \
+            } \
         } \
         if (representation == 'r') \
         { \
@@ -440,8 +445,12 @@ int fluid_solver<R>::read(char field, char representation) \
 template<> \
 int fluid_solver<R>::write(char field, char representation) \
 { \
+    char fname[512]; \
     if ((field == 'v') && (representation == 'c')) \
-        return this->write_base("cvorticity", this->cvorticity); \
+    { \
+        this->fill_up_filename("cvorticity", fname); \
+        return this->cd->write(fname, (void*)this->cvorticity); \
+    } \
     if ((field == 'v') && (representation == 'r')) \
     { \
         FFTW(execute)(*((FFTW(plan)*)this->c2r_vorticity )); \
diff --git a/src/fluid_solver_base.cpp b/src/fluid_solver_base.cpp
index adea136a..033dcce4 100644
--- a/src/fluid_solver_base.cpp
+++ b/src/fluid_solver_base.cpp
@@ -30,11 +30,17 @@
 #include "fftw_tools.hpp"
 
 
+template <class rnumber>
+void fluid_solver_base<rnumber>::fill_up_filename(const char *base_name, char *destination)
+{
+    sprintf(destination, "%s_%s_i%.5x", this->name, base_name, this->iteration); \
+}
+
 template <class rnumber>
 void fluid_solver_base<rnumber>::clean_up_real_space(rnumber *a, int howmany)
 {
-//    for (ptrdiff_t rindex = 0; rindex < this->cd->local_size*2; rindex += howmany*(this->rd->subsizes[2]+2))
-//        std::fill_n(a+rindex+this->rd->subsizes[2]*howmany, 2*howmany, 0.0);
+    for (ptrdiff_t rindex = 0; rindex < this->cd->local_size*2; rindex += howmany*(this->rd->subsizes[2]+2))
+        std::fill_n(a+rindex+this->rd->subsizes[2]*howmany, 2*howmany, 0.0);
 }
 
 template <class rnumber>
diff --git a/src/fluid_solver_base.hpp b/src/fluid_solver_base.hpp
index b948ab5d..a85c4446 100644
--- a/src/fluid_solver_base.hpp
+++ b/src/fluid_solver_base.hpp
@@ -80,6 +80,7 @@ class fluid_solver_base
         rnumber correl_vec(cnumber *a, cnumber *b);
         void cospectrum(cnumber *a, cnumber *b, double *spec, const double k2exponent = 0.0);
         void write_spectrum(const char *fname, cnumber *a, const double k2exponent = 0.0);
+        void fill_up_filename(const char *base_name, char *full_name);
         int read_base(const char *fname, rnumber *data);
         int read_base(const char *fname, cnumber *data);
         int write_base(const char *fname, rnumber *data);
-- 
GitLab