diff --git a/CMakeLists.txt b/CMakeLists.txt
index aa5e8c13fe89ef2610e5bcf1a1581089737485ad..3b062abff905e526a6d6e657e1d11f96c449c816 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -272,6 +272,7 @@ set(cpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/symmetrize_test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_output_test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/get_rfields.cpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/write_rpressure.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/bandpass_stats.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_single_to_double.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/resize.cpp
@@ -321,6 +322,7 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/symmetrize_test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_output_test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/get_rfields.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/write_rpressure.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/bandpass_stats.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_single_to_double.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/resize.hpp
@@ -454,6 +456,10 @@ if (BUILD_TESTING)
NAME test_pp_get_rfields
COMMAND turtle PP get_rfields --simname dns_nsveparticles --iter0 0 --iter1 64
WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_pp_write_rpressure
+ COMMAND turtle PP write_rpressure --simname dns_nsveparticles --iter0 0 --iter1 64
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
add_test(
NAME test_pp_joint_acc_vel_stats
COMMAND turtle PP joint_acc_vel_stats --simname dns_nsveparticles --iter0 0 --iter1 64
diff --git a/TurTLE/PP.py b/TurTLE/PP.py
index e1cebad2bc079fcca2843d2e6be63e2d9228b4bd..acf3b003c0881f89fbed97fffb67c897a6b00786 100644
--- a/TurTLE/PP.py
+++ b/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 + ')')
diff --git a/cpp/full_code/write_rpressure.cpp b/cpp/full_code/write_rpressure.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ceeca61c127f6d31424e805095cd73d5a0e1eb08
--- /dev/null
+++ b/cpp/full_code/write_rpressure.cpp
@@ -0,0 +1,138 @@
+#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>;
+
diff --git a/cpp/full_code/write_rpressure.hpp b/cpp/full_code/write_rpressure.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e7e2d6d9c86c16bf3297fb7c5856ecfd86c23824
--- /dev/null
+++ b/cpp/full_code/write_rpressure.hpp
@@ -0,0 +1,68 @@
+/**********************************************************************
+* *
+* 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
+
diff --git a/cpp/vorticity_equation.cpp b/cpp/vorticity_equation.cpp
index f062667104303cf50b15c903eca319eddc880595..31a46cddeb18b77f4d4049f14350e07bb66472bb 100644
--- a/cpp/vorticity_equation.cpp
+++ b/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 (