From 8d88006bbcbcb9d05e1b25f6a7f2b5b268bc0346 Mon Sep 17 00:00:00 2001
From: Chichi Lalescu <chichilalescu@gmail.com>
Date: Thu, 25 May 2017 21:45:53 +0200
Subject: [PATCH] real space conversion works

---
 bfps/PP.py                                   | 10 ++-
 bfps/cpp/full_code/NSVE_field_stats.cpp      | 94 ++++++++++++++++++++
 bfps/cpp/full_code/NSVE_field_stats.hpp      | 63 +++++++++++++
 bfps/cpp/full_code/get_rfields.cpp           | 93 +++++++++++++++++++
 bfps/cpp/full_code/get_rfields.hpp           | 61 +++++++++++++
 bfps/cpp/full_code/native_binary_to_hdf5.cpp |  1 +
 bfps/cpp/full_code/postprocess.cpp           |  3 +-
 setup.py                                     |  2 +
 8 files changed, 325 insertions(+), 2 deletions(-)
 create mode 100644 bfps/cpp/full_code/NSVE_field_stats.cpp
 create mode 100644 bfps/cpp/full_code/NSVE_field_stats.hpp
 create mode 100644 bfps/cpp/full_code/get_rfields.cpp
 create mode 100644 bfps/cpp/full_code/get_rfields.hpp

diff --git a/bfps/PP.py b/bfps/PP.py
index ffa0b2da..71c67286 100644
--- a/bfps/PP.py
+++ b/bfps/PP.py
@@ -451,6 +451,12 @@ class PP(_code):
         self.simulation_parser_arguments(parser_native_binary_to_hdf5)
         self.job_parser_arguments(parser_native_binary_to_hdf5)
         self.parameters_to_parser_arguments(parser_native_binary_to_hdf5)
+        parser_get_rfields = subparsers.add_parser(
+                'get_rfields',
+                help = 'get real space velocity field')
+        self.simulation_parser_arguments(parser_get_rfields)
+        self.job_parser_arguments(parser_get_rfields)
+        self.parameters_to_parser_arguments(parser_get_rfields)
         return None
     def prepare_launch(
             self,
@@ -671,6 +677,8 @@ class PP(_code):
                 njobs = opt.njobs,
                 hours = opt.minutes // 60,
                 minutes = opt.minutes % 60,
-                no_submit = opt.no_submit)
+                no_submit = opt.no_submit,
+                err_file = 'err_' + self.dns_type,
+                out_file = 'out_' + self.dns_type)
         return None
 
diff --git a/bfps/cpp/full_code/NSVE_field_stats.cpp b/bfps/cpp/full_code/NSVE_field_stats.cpp
new file mode 100644
index 00000000..344cf2ce
--- /dev/null
+++ b/bfps/cpp/full_code/NSVE_field_stats.cpp
@@ -0,0 +1,94 @@
+#include <string>
+#include <cmath>
+#include "NSVE_field_stats.hpp"
+#include "scope_timer.hpp"
+
+
+template <typename rnumber>
+int NSVE_field_stats<rnumber>::initialize(void)
+{
+    this->postprocess::read_parameters();
+    this->vorticity = new field<rnumber, FFTW, THREE>(
+            nx, ny, nz,
+            this->comm,
+            DEFAULT_FFTW_FLAG);
+    this->vorticity->real_space_representation = false;
+    hid_t parameter_file = H5Fopen(
+            (this->simname + std::string(".h5")).c_str(),
+            H5F_ACC_RDONLY,
+            H5P_DEFAULT);
+    if (!H5Lexists(parameter_file, "field_dtype", H5P_DEFAULT))
+        this->bin_IO = NULL;
+    else
+    {
+        hid_t dset = H5Dopen(parameter_file, "field_dtype", H5P_DEFAULT);
+        hid_t space = H5Dget_space(dset);
+        hid_t memtype = H5Dget_type(dset);
+        char *string_data = (char*)malloc(256);
+        H5Dread(dset, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &string_data);
+        // check that we're using the correct data type
+        // field_dtype SHOULD be something like "<f4", "<f8", ">f4", ">f8"
+        // first character is ordering, which is machine specific
+        // for the other two I am checking that they have the correct values
+        assert(string_data[1] == 'f');
+        assert(string_data[2] == '0' + sizeof(rnumber));
+        free(string_data);
+        H5Sclose(space);
+        H5Tclose(memtype);
+        H5Dclose(dset);
+        this->bin_IO = new field_binary_IO<rnumber, COMPLEX, THREE>(
+                this->vorticity->clayout->sizes,
+                this->vorticity->clayout->subsizes,
+                this->vorticity->clayout->starts,
+                this->vorticity->clayout->comm);
+    }
+    H5Fclose(parameter_file);
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_field_stats<rnumber>::read_current_cvorticity(void)
+{
+    this->vorticity->real_space_representation = false;
+    if (this->bin_IO != NULL)
+    {
+        char itername[16];
+        sprintf(itername, "i%.5x", this->iteration);
+        std::string native_binary_fname = (
+                this->simname +
+                std::string("_cvorticity_") +
+                std::string(itername));
+        this->bin_IO->read(
+                native_binary_fname,
+                this->vorticity->get_cdata());
+    }
+    else
+    {
+        this->vorticity->io(
+                this->simname + std::string("_fields.h5"),
+                "vorticity",
+                this->iteration,
+                true);
+    }
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_field_stats<rnumber>::finalize(void)
+{
+    if (this->bin_IO != NULL)
+        delete this->bin_IO;
+    delete this->vorticity;
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_field_stats<rnumber>::work_on_current_iteration(void)
+{
+    DEBUG_MSG("entered NSVE_field_stats::work_on_current_iteration\n");
+    return EXIT_SUCCESS;
+}
+
+template class NSVE_field_stats<float>;
+template class NSVE_field_stats<double>;
+
diff --git a/bfps/cpp/full_code/NSVE_field_stats.hpp b/bfps/cpp/full_code/NSVE_field_stats.hpp
new file mode 100644
index 00000000..d544c0c7
--- /dev/null
+++ b/bfps/cpp/full_code/NSVE_field_stats.hpp
@@ -0,0 +1,63 @@
+/**********************************************************************
+*                                                                     *
+*  Copyright 2017 Max Planck Institute                                *
+*                 for Dynamics and Self-Organization                  *
+*                                                                     *
+*  This file is part of bfps.                                         *
+*                                                                     *
+*  bfps is free software: you can redistribute it and/or modify       *
+*  it under the terms of the GNU General Public License as published  *
+*  by the Free Software Foundation, either version 3 of the License,  *
+*  or (at your option) any later version.                             *
+*                                                                     *
+*  bfps is distributed in the hope that it will be useful,            *
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of     *
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      *
+*  GNU General Public License for more details.                       *
+*                                                                     *
+*  You should have received a copy of the GNU General Public License  *
+*  along with bfps.  If not, see <http://www.gnu.org/licenses/>       *
+*                                                                     *
+* Contact: Cristian.Lalescu@ds.mpg.de                                 *
+*                                                                     *
+**********************************************************************/
+
+
+
+#ifndef NSVE_FIELD_STATS_HPP
+#define NSVE_FIELD_STATS_HPP
+
+#include <cstdlib>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <vector>
+#include "base.hpp"
+#include "field.hpp"
+#include "field_binary_IO.hpp"
+#include "full_code/postprocess.hpp"
+
+template <typename rnumber>
+class NSVE_field_stats: public postprocess
+{
+    private:
+        field_binary_IO<rnumber, COMPLEX, THREE> *bin_IO;
+    public:
+        field<rnumber, FFTW, THREE> *vorticity;
+
+        NSVE_field_stats(
+                const MPI_Comm COMMUNICATOR,
+                const std::string &simulation_name):
+            postprocess(
+                    COMMUNICATOR,
+                    simulation_name){}
+        virtual ~NSVE_field_stats(){}
+
+        virtual int initialize(void);
+        virtual int work_on_current_iteration(void);
+        virtual int finalize(void);
+
+        int read_current_cvorticity(void);
+};
+
+#endif//NSVE_FIELD_STATS_HPP
+
diff --git a/bfps/cpp/full_code/get_rfields.cpp b/bfps/cpp/full_code/get_rfields.cpp
new file mode 100644
index 00000000..0df8b564
--- /dev/null
+++ b/bfps/cpp/full_code/get_rfields.cpp
@@ -0,0 +1,93 @@
+#include <string>
+#include <cmath>
+#include "get_rfields.hpp"
+#include "scope_timer.hpp"
+
+
+template <typename rnumber>
+int get_rfields<rnumber>::initialize(void)
+{
+    this->NSVE_field_stats<rnumber>::initialize();
+    this->kk = new kspace<FFTW, SMOOTH>(
+            this->vorticity->clayout, this->dkx, this->dky, this->dkz);
+    hid_t parameter_file = H5Fopen(
+            (this->simname + std::string(".h5")).c_str(),
+            H5F_ACC_RDONLY,
+            H5P_DEFAULT);
+    hid_t dset = H5Dopen(parameter_file, "/parameters/niter_out", H5P_DEFAULT);
+    H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->niter_out);
+    H5Dclose(dset);
+    if (H5Lexists(parameter_file, "/parameters/checkpoints_per_file", H5P_DEFAULT))
+    {
+        dset = H5Dopen(parameter_file, "/parameters/checkpoints_per_file", H5P_DEFAULT);
+        H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->checkpoints_per_file);
+        H5Dclose(dset);
+    }
+    else
+        this->checkpoints_per_file = 1;
+    this->iteration_list = hdf5_tools::read_vector<int>(
+            parameter_file,
+            "/get_rfields/iteration_list");
+    H5Fclose(parameter_file);
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int get_rfields<rnumber>::work_on_current_iteration(void)
+{
+    DEBUG_MSG("entered get_rfields::work_on_current_iteration\n");
+    this->read_current_cvorticity();
+    field<rnumber, FFTW, THREE> *vel = new field<rnumber, FFTW, THREE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            this->vorticity->fftw_plan_rigor);
+
+    vel->real_space_representation = false;
+    this->kk->CLOOP_K2(
+                [&](ptrdiff_t cindex,
+                    ptrdiff_t xindex,
+                    ptrdiff_t yindex,
+                    ptrdiff_t zindex,
+                    double k2){
+        if (k2 <= this->kk->kM2 && k2 > 0)
+        {
+            vel->cval(cindex,0,0) = -(this->kk->ky[yindex]*this->vorticity->cval(cindex,2,1) - this->kk->kz[zindex]*this->vorticity->cval(cindex,1,1)) / k2;
+            vel->cval(cindex,0,1) =  (this->kk->ky[yindex]*this->vorticity->cval(cindex,2,0) - this->kk->kz[zindex]*this->vorticity->cval(cindex,1,0)) / k2;
+            vel->cval(cindex,1,0) = -(this->kk->kz[zindex]*this->vorticity->cval(cindex,0,1) - this->kk->kx[xindex]*this->vorticity->cval(cindex,2,1)) / k2;
+            vel->cval(cindex,1,1) =  (this->kk->kz[zindex]*this->vorticity->cval(cindex,0,0) - this->kk->kx[xindex]*this->vorticity->cval(cindex,2,0)) / k2;
+            vel->cval(cindex,2,0) = -(this->kk->kx[xindex]*this->vorticity->cval(cindex,1,1) - this->kk->ky[yindex]*this->vorticity->cval(cindex,0,1)) / k2;
+            vel->cval(cindex,2,1) =  (this->kk->kx[xindex]*this->vorticity->cval(cindex,1,0) - this->kk->ky[yindex]*this->vorticity->cval(cindex,0,0)) / k2;
+        }
+        else
+            std::fill_n((rnumber*)(vel->get_cdata()+3*cindex), 6, 0.0);
+    }
+    );
+    vel->symmetrize();
+    vel->ift();
+
+    std::string fname = (
+            this->simname +
+            std::string("_checkpoint_") +
+            std::to_string(this->iteration / (this->niter_out*this->checkpoints_per_file)) +
+            std::string(".h5"));
+    vel->io(
+            fname,
+            "velocity",
+            this->iteration,
+            false);
+
+    delete vel;
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int get_rfields<rnumber>::finalize(void)
+{
+    delete this->kk;
+    this->NSVE_field_stats<rnumber>::finalize();
+    return EXIT_SUCCESS;
+}
+
+template class get_rfields<float>;
+template class get_rfields<double>;
+
diff --git a/bfps/cpp/full_code/get_rfields.hpp b/bfps/cpp/full_code/get_rfields.hpp
new file mode 100644
index 00000000..ae669aa3
--- /dev/null
+++ b/bfps/cpp/full_code/get_rfields.hpp
@@ -0,0 +1,61 @@
+/**********************************************************************
+*                                                                     *
+*  Copyright 2017 Max Planck Institute                                *
+*                 for Dynamics and Self-Organization                  *
+*                                                                     *
+*  This file is part of bfps.                                         *
+*                                                                     *
+*  bfps is free software: you can redistribute it and/or modify       *
+*  it under the terms of the GNU General Public License as published  *
+*  by the Free Software Foundation, either version 3 of the License,  *
+*  or (at your option) any later version.                             *
+*                                                                     *
+*  bfps is distributed in the hope that it will be useful,            *
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of     *
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the      *
+*  GNU General Public License for more details.                       *
+*                                                                     *
+*  You should have received a copy of the GNU General Public License  *
+*  along with bfps.  If not, see <http://www.gnu.org/licenses/>       *
+*                                                                     *
+* Contact: Cristian.Lalescu@ds.mpg.de                                 *
+*                                                                     *
+**********************************************************************/
+
+
+
+#ifndef GET_RFIELDS_HPP
+#define GET_RFIELDS_HPP
+
+#include <cstdlib>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <vector>
+#include "base.hpp"
+#include "field.hpp"
+#include "field_binary_IO.hpp"
+#include "full_code/NSVE_field_stats.hpp"
+
+template <typename rnumber>
+class get_rfields: public NSVE_field_stats<rnumber>
+{
+    public:
+        int checkpoints_per_file;
+        int niter_out;
+        kspace<FFTW, SMOOTH> *kk;
+
+        get_rfields(
+                const MPI_Comm COMMUNICATOR,
+                const std::string &simulation_name):
+            NSVE_field_stats<rnumber>(
+                    COMMUNICATOR,
+                    simulation_name){}
+        virtual ~get_rfields(){}
+
+        int initialize(void);
+        int work_on_current_iteration(void);
+        int finalize(void);
+};
+
+#endif//GET_RFIELDS_HPP
+
diff --git a/bfps/cpp/full_code/native_binary_to_hdf5.cpp b/bfps/cpp/full_code/native_binary_to_hdf5.cpp
index 78a29add..7774e2de 100644
--- a/bfps/cpp/full_code/native_binary_to_hdf5.cpp
+++ b/bfps/cpp/full_code/native_binary_to_hdf5.cpp
@@ -61,6 +61,7 @@ int native_binary_to_hdf5<rnumber>::read_parameters(void)
     this->iteration_list = hdf5_tools::read_vector<int>(
             parameter_file,
             "/native_binary_to_hdf5/iteration_list");
+    H5Fclose(parameter_file);
     return EXIT_SUCCESS;
 }
 
diff --git a/bfps/cpp/full_code/postprocess.cpp b/bfps/cpp/full_code/postprocess.cpp
index 193dafeb..05242f39 100644
--- a/bfps/cpp/full_code/postprocess.cpp
+++ b/bfps/cpp/full_code/postprocess.cpp
@@ -8,6 +8,8 @@
 
 int postprocess::main_loop(void)
 {
+    DEBUG_MSG("entered postprocess::main_loop\n");
+    DEBUG_MSG("interation_list.size = %d\n", iteration_list.size());
     this->start_simple_timer();
     for (unsigned int iteration_counter = 0;
          iteration_counter < iteration_list.size();
@@ -32,7 +34,6 @@ int postprocess::read_parameters()
     hid_t parameter_file;
     hid_t dset, memtype, space;
     char fname[256];
-    hsize_t dims[1];
     char *string_data;
     sprintf(fname, "%s.h5", this->simname.c_str());
     parameter_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
diff --git a/setup.py b/setup.py
index cd2dd02a..1178f547 100644
--- a/setup.py
+++ b/setup.py
@@ -89,6 +89,8 @@ print('This is bfps version ' + VERSION)
 
 ### lists of files and MANIFEST.in
 src_file_list = ['hdf5_tools',
+                 'full_code/get_rfields',
+                 'full_code/NSVE_field_stats',
                  'full_code/native_binary_to_hdf5',
                  'full_code/postprocess',
                  'full_code/code_base',
-- 
GitLab