add: PP option to write real pressure to binary

fix: pressure representation status in vorticity_equation

TurTLE/PP.py

@@ -449,13 +449,17 @@ class PP(_code):
         parser_resize = subparsers.add_parser(
                 'resize',
                 help = 'get joint acceleration and velocity statistics')
+        parser_write_rpressure = subparsers.add_parser(
+                'write_rpressure',
+                help = 'write real pressure field to binary')
         for pp_type in [
                 'resize',
                 'joint_acc_vel_stats',
                 'bandpass_stats',
                 'get_rfields',
                 'field_single_to_double',
-                'native_binary_to_hdf5']:
+                'native_binary_to_hdf5',
+                'write_rpressure']:
             eval('self.simulation_parser_arguments(parser_' + pp_type + ')')
             eval('self.job_parser_arguments(parser_' + pp_type + ')')
             eval('self.parameters_to_parser_arguments(parser_' + pp_type + ')')

bfps/cpp/full_code/write_rpressure.cpp

+#include <string>
+#include <cmath>
+#include "write_rpressure.hpp"
+#include "scope_timer.hpp"
+
+
+template <typename rnumber>
+int write_rpressure<rnumber>::initialize(void)
+{
+    
+    // initialize
+    this->NSVE_field_stats<rnumber>::initialize();   
+    
+    //iteration list needed by postprocess.cpp for the main loop
+    hid_t parameter_file = H5Fopen(
+            (this->simname + std::string("_post.h5")).c_str(),
+            H5F_ACC_RDONLY,
+            H5P_DEFAULT);
+    this->iteration_list = hdf5_tools::read_vector<int>(
+            parameter_file,
+            "/write_rpressure/parameters/iteration_list");
+    H5Fclose(parameter_file);
+    if (this->myrank==0) DEBUG_MSG("Iteration list[0] = %d \n", this->iteration_list[0]);
+    
+    //get the pressure from here
+    this->ve = new vorticity_equation<rnumber, FFTW>(
+            this->simname.c_str(),
+            this->nx,
+            this->ny,
+            this->nz,
+            this->dkx,
+            this->dky,
+            this->dkz,
+            this->vorticity->fftw_plan_rigor);
+    
+    // output field
+    this->pressure = new field<rnumber, FFTW, ONE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            this->vorticity->fftw_plan_rigor);
+    
+    // write to binary 
+    this->bin_IO_scalar = new field_binary_IO<rnumber, REAL, ONE>(
+                this->pressure->rlayout->sizes,
+                this->pressure->rlayout->subsizes,
+                this->pressure->rlayout->starts,
+                this->pressure->rlayout->comm);
+    
+    if (this->myrank==0) DEBUG_MSG("Initialize end\n");
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int write_rpressure<rnumber>::clip_zero_padding(
+    field <rnumber, FFTW, ONE>* f)
+{
+    rnumber *a = f->get_rdata();
+    rnumber *b = f->get_rdata();
+    ptrdiff_t copy_size = f->rlayout->sizes[2] * ncomp(ONE);
+    ptrdiff_t skip_size = copy_size + 2*ncomp(ONE);
+    for (int i0 = 0; i0 < f->rlayout->subsizes[0]; i0++)
+        for (int i1 = 0; i1 < f->rlayout->sizes[1]; i1++)
+        {
+            std::copy(a, a + copy_size, b);
+            a += skip_size;
+            b += copy_size;
+        }
+    return EXIT_SUCCESS;
+}
+
+
+template <typename rnumber>
+int write_rpressure<rnumber>::work_on_current_iteration(void)
+{
+    if (this->myrank==0) DEBUG_MSG("\nComputing the pressure\n");
+    //compute_pressure    
+    this->read_current_cvorticity();
+    *this->ve->cvorticity = this->vorticity->get_cdata();
+    
+    this->ve->compute_velocity(this->ve->cvorticity);
+    this->ve->cvelocity->ift();
+    this->ve->compute_pressure(this->pressure);
+    
+    
+    if (this->myrank==0) DEBUG_MSG("Computed the pressure\n");
+    
+    if (this->myrank==0) DEBUG_MSG("vorticity in ve fourier = [%lG, %lG; ...; %lG, %lG;  %lG, %lG;.....; %lG, %lG;] \n", this->vorticity->get_cdata()[0][0], this->vorticity->get_cdata()[0][1], this->vorticity->get_cdata()[1][0], this->vorticity->get_cdata()[1][1], this->vorticity->get_cdata()[2][0], this->vorticity->get_cdata()[2][1],this->vorticity->get_cdata()[65][0], this->vorticity->get_cdata()[65][1]); 
+   
+    if (this->myrank==0) DEBUG_MSG("vorticity in ve fourier = [%lG, %lG; ...; %lG, %lG;  %lG, %lG;.....; %lG, %lG;] \n", this->ve->cvorticity->get_cdata()[0][0], this->ve->cvorticity->get_cdata()[0][1], this->ve->cvorticity->get_cdata()[1][0], this->ve->cvorticity->get_cdata()[1][1], this->ve->cvorticity->get_cdata()[2][0], this->ve->cvorticity->get_cdata()[2][1],this->ve->cvorticity->get_cdata()[65][0], this->ve->cvorticity->get_cdata()[65][1]);
+    this->ve->compute_velocity(this->vorticity);
+    
+    if (this->myrank==0) DEBUG_MSG("velocity in ve fourier = [%lG, %lG; ...; %lG, %lG;  %lG, %lG;.....; %lG, %lG;] \n", this->ve->u->get_cdata()[0][0], this->ve->u->get_cdata()[0][1], this->ve->u->get_cdata()[1][0], this->ve->u->get_cdata()[1][1], this->ve->u->get_cdata()[2][0], this->ve->u->get_cdata()[2][1],this->ve->u->get_cdata()[65][0], this->ve->u->get_cdata()[65][1]);
+    
+    if (this->myrank==0) DEBUG_MSG("pressure fourier = [%lG, %lG; ...; %lG, %lG;  %lG, %lG;.....; %lG, %lG;] \n", this->pressure->get_cdata()[0][0], this->pressure->get_cdata()[0][1], this->pressure->get_cdata()[1][0], this->pressure->get_cdata()[1][1], this->pressure->get_cdata()[2][0], this->pressure->get_cdata()[2][1],this->pressure->get_cdata()[65][0], this->pressure->get_cdata()[65][1]);
+     
+    //transform to real space
+    if (this->pressure->real_space_representation == false)  
+        this->pressure->ift();
+    
+    if (this->myrank==0) DEBUG_MSG("Real space\n");
+    
+    if (this->myrank==0) DEBUG_MSG("pressure before clipping = [%lG; ...; %lG;.....; %lG;] \n", this->pressure->get_rdata()[0], this->pressure->get_rdata()[2], this->pressure->get_rdata()[13549]);
+    
+    //write out to file
+    this->clip_zero_padding(pressure);
+    if (this->myrank==0) DEBUG_MSG("pressure after clipping = [%lG; ...; %lG;.....; %lG;] \n", this->pressure->get_rdata()[0], this->pressure->get_rdata()[2], this->pressure->get_rdata()[13549]);
+    
+    //write to file -- THIS IS OK
+    char itername[16];
+    sprintf(itername, "i%.5x", this->iteration);
+    std::string native_binary_fname;
+
+    native_binary_fname = (
+            this->simname +
+            std::string("_rpressure_") +
+            std::string(itername));
+    if (this->myrank ==0) DEBUG_MSG("%s\n", native_binary_fname.c_str());
+    
+    this->bin_IO_scalar->write(
+            native_binary_fname,
+            this->pressure->get_rdata());
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int write_rpressure<rnumber>::finalize(void)
+{
+     
+    delete this->bin_IO_scalar;
+    delete this->pressure;
+    delete this->ve;   
+    this->NSVE_field_stats<rnumber>::finalize();
+    return EXIT_SUCCESS;
+}
+
+template class write_rpressure<float>;
+template class write_rpressure<double>;
+

bfps/cpp/full_code/write_rpressure.hpp

+/**********************************************************************
+*                                                                     *
+*  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 WRITE_RPRESSURE
+#define WRITE_RPRESSURE
+
+#include <cstdlib>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <vector>
+#include "base.hpp"
+#include "field.hpp"
+#include "field_binary_IO.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/NSVE_field_stats.hpp"
+
+/** generate and write out pressure field in real space **/
+
+template <typename rnumber>
+class write_rpressure: public NSVE_field_stats<rnumber>
+{
+    public:
+        int niter_out;
+        
+        vorticity_equation<rnumber, FFTW> *ve;        
+        field<rnumber, FFTW, ONE> *pressure;
+        field_binary_IO <rnumber, REAL, ONE> *bin_IO_scalar;        
+        int clip_zero_padding(field <rnumber, FFTW, ONE>* );
+        
+        write_rpressure(
+                const MPI_Comm COMMUNICATOR,
+                const std::string &simulation_name):
+            NSVE_field_stats<rnumber>(
+                    COMMUNICATOR,
+                    simulation_name){}
+        virtual ~write_rpressure(){}
+
+        int initialize(void);
+        int work_on_current_iteration(void);
+        int finalize(void);
+};
+
+
+#endif//WRITE_PRESSURE
+

cpp/vorticity_equation.cpp

@@ -595,7 +595,9 @@ void vorticity_equation<rnumber, be>::compute_pressure(field<rnumber, be, ONE> *
     TIMEZONE("vorticity_equation::compute_pressure");
     pressure->real_space_representation = false;
     /* assume velocity is already in real space representation */
-
+    /* set the pressure representation to Fourier space */
+    pressure->real_space_representation = false;
+    
     this->v[1]->real_space_representation = true;
     /* diagonal terms 11 22 33 */
     this->v[1]->RLOOP (