Merge branch 'develop'

I guess some people would call this a prerelease.
Most of the mechanisms are there, and I'm already using the code for
production, so we're very close to 1.0.
I still have a couple of easy todo items to touch before that.

.gitignore

 .ipynb_checkpoints/*
+meta/.ipynb_checkpoints/*
 .*.swp
 *.pyc
 .tox

README.rst

+================================
 Big Fluid and Particle Simulator
 ================================
 
@@ -12,26 +13,98 @@ in the repository), and it is working as expected. Parameters and
 statistics are stored in HDF5 format, together with code information,
 so simulation data should be "future proof".
 
+Users of this code are expected to either use `NavierStokes` objects
+directly, or construct their own class that inherits this class.
+The way I use it is to I inherit and add custom statistics as necessary; I
+also have private C++ code that can get added and used when needed.
+I plan on adding documentation on the procedure as when other people
+show interest in using the code, or when time permits.
+
+
+Installation
+------------
+
+**Postprocessing only**
+
+.. code:: bash
+
+    python setup.py install
+
+(add `--user` or `sudo` as appropriate).
+`setup.py` should tell you about the various packages you need.
+
+**Full installation**
+
+If you want to run simulations on the machine where you're installing,
+you will need to call `build` before installing.
+
+.. code:: bash
+
+    python setup.py build
+    python setup.py install
+
+The `build` command will most likely fail unless you modify
+`machine_settings.py` appropriately for your machine.
+Also, in order to run the C++ code you need to have an MPI compiler
+installed, the HDF5 C library as well as FFTW 3 (at least 3.3 I think).
+
+
 Comments
 --------
 
 * particles: initialization of multistep solvers is done with lower
   order methods, so don't be surprised if direct convergence tests fail.
 
-TODO
-----
+* I am using this code mainly with Python 3.4, but Python 2.7
+  compatibility should be kept since mayavi (well, vtk actually) only
+  works on Python 2.
+
 
-* add different numerical method, for testing purposes
+===========
+Development
+===========
 
-* test involving hydrodynamic similarity
+`bfps` is using `git` for version tracking (see https://git-scm.com/
+for description).
+The branching model described at
+http://nvie.com/posts/a-successful-git-branching-model/ should be
+adhered to as strictly as possible.
 
-* test anisotropic grids
+The usable `VERSION` number will be constructed according to semantic
+versioning rules (see http://semver.org/), `MAJOR.MINOR.PATCH`.
+In principle, if you use `bfps` and you've created children of the
+`NavierStokes` class, you should not need to rewrite your code unless
+the `MAJOR` version changes.
 
-* test non-cubic domains
+Versioning guidelines
+---------------------
 
-* use HDF5 io for fields
+There are 2 main branches, `master` and `develop`.
+`setup.py` will call `git` to read in the `VERSION`: it will get the
+latest available tag.
+If the active branch name contains either of the strings `develop`,
+`feature` or `bugfix`, then the full output of `git describe --tags`
+will be used;
+otherwise, only the string before the dash (if a dash exists) will be
+used.
 
-* try to make code more memory efficient
+At the moment the following rules seem adequate.
+Once I get feedback from someone who actually knows how to handle bigger
+projects, these may change.
 
-* make code py3 compatible
+1. New features are worked on in branches forked from `develop`, with
+   the branch name of the form `feature/bla`.
+   Feature branches are merged back into `develop` only after all tests
+   pass.
+2. Whenever the `develop` branch is merged into `master`, the last
+   commit on the `develop` branch is tagged with `MAJOR.MINOR`.
+3. Whenever a bug is discovered in version X.Y, a new branch called `vX.Y`
+   is forked from the corresponding `master` merge point.
+   A new bugfix branch is forked from `vX.Y`, the bug is fixed, and then
+   this bugfix branch is merged into all affected branches (this includes
+   `vX.Y`).
+   After merging, the respective merge points are tagged adequately (`vX.Y`
+   gets the tag X.Y.1).
+4. Whenever a bug is discovered in version X.Y.Z, a bugfix branch is
+   forked from `vX.Y`, and then rule 3 is adapted accordingly.
 

bfps/__init__.py

@@ -28,28 +28,20 @@ import os
 import subprocess
 import pickle
 
-from pkg_resources import get_distribution, DistributionNotFound
+import pkg_resources
 
-try:
-    _dist = get_distribution('bfps')
-    # Normalize case for Windows systems
-    dist_loc = os.path.normcase(os.path.realpath(_dist.location))
-    here = os.path.normcase(__file__)
-    if not here.startswith(os.path.join(dist_loc, 'bfps')):
-        # not installed, but there is another version that *is*
-        header_dir = os.path.join(os.path.dirname(here), 'cpp')
-        lib_dir = os.path.join(os.path.dirname(here), os.pardir)
-        raise DistributionNotFound
-    header_dir = os.path.join(os.path.join(dist_loc, 'bfps'), 'cpp')
-    lib_dir = _dist.location
-    __version__ = _dist.version
-except DistributionNotFound:
-    __version__ = ''
+__version__ = pkg_resources.require('bfps')[0].version
+
+_dist = pkg_resources.get_distribution('bfps')
+dist_loc = os.path.realpath(_dist.location)
+here = os.path.normcase(__file__)
+header_dir = os.path.join(os.path.join(dist_loc, 'bfps'), 'cpp')
+lib_dir = os.path.join(dist_loc, 'bfps')
 
 install_info = pickle.load(
         open(os.path.join(os.path.dirname(here),
                           'install_info.pickle'),
-             'r'))
+             'rb'))
 
 from .code import code
 from .fluid_converter import fluid_converter
@@ -85,7 +77,7 @@ def get_parser(base_class = NavierStokes,
     parser.add_argument('--njobs',
             type = int, dest = 'njobs',
             default = njobs)
-    c = base_class()
+    c = base_class(simname = simname)
     for k in sorted(c.parameters.keys()):
         parser.add_argument(
                 '--{0}'.format(k),

bfps/base.py

@@ -25,6 +25,7 @@
 
 
 import os
+import sys
 import numpy as np
 import h5py
 import bfps
@@ -48,8 +49,7 @@ class base(object):
         self.simname = simname
         return None
     def cdef_pars(self):
-        key = self.parameters.keys()
-        key.sort()
+        key = sorted(list(self.parameters.keys()))
         src_txt = ''
         for i in range(len(key)):
             if type(self.parameters[key[i]]) == int:
@@ -60,21 +60,18 @@ class base(object):
                 src_txt += 'double ' + key[i] + ';\n'
         return src_txt
     def cread_pars(self):
-        key = self.parameters.keys()
-        key.sort()
+        key = sorted(list(self.parameters.keys()))
         src_txt = ('int read_parameters(hid_t data_file_id)\n{\n'
                  + 'hid_t dset, memtype, space;\n'
                  + 'hsize_t dims[1];\n'
-                 + 'char *string_data;\n'
-                 + 'std::string tempstr;\n')
+                 + 'char *string_data;\n')
         for i in range(len(key)):
             src_txt += 'dset = H5Dopen(data_file_id, "parameters/{0}", H5P_DEFAULT);\n'.format(key[i])
             if type(self.parameters[key[i]]) == int:
                 src_txt += 'H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &{0});\n'.format(key[i])
             elif type(self.parameters[key[i]]) == str:
                 src_txt += ('space = H5Dget_space(dset);\n' +
-                            'memtype = H5Tcopy(H5T_C_S1);\n' +
-                            'H5Tset_size(memtype, H5T_VARIABLE);\n' +
+                            'memtype = H5Dget_type(dset);\n' +
                             'H5Sget_simple_extent_dims(space, dims, NULL);\n' +
                             'string_data = (char*)malloc(dims[0]*sizeof(char));\n' +
                             'H5Dread(dset, memtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &string_data);\n' +
@@ -88,8 +85,7 @@ class base(object):
         src_txt += 'return 0;\n}\n' # finishing read_parameters
         return src_txt
     def cprint_pars(self):
-        key = self.parameters.keys()
-        key.sort()
+        key = sorted(list(self.parameters.keys()))
         src_txt = ''
         for i in range(len(key)):
             if type(self.parameters[key[i]]) == int:
@@ -104,7 +100,14 @@ class base(object):
             os.makedirs(self.work_dir)
         ofile = h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'w-')
         for k in self.parameters.keys():
-            ofile['parameters/' + k] = self.parameters[k]
+            if (type(self.parameters[k]) == str) and (sys.version[0] == 3):
+                ofile.create_dataset('parameters/' + k,
+                                     (1,),
+                                     dtype = 'S10')
+                #ofile['parameters/' + k] = self.parameters[k].encode('ascii', 'ignore')
+                #print (ofile['parameters/' + k])
+            else:
+                ofile['parameters/' + k] = self.parameters[k]
         ofile['iteration'] = int(iter0)
         for k in bfps.install_info.keys():
             ofile['install_info/' + k] = str(bfps.install_info[k])

bfps/code.py

@@ -24,10 +24,11 @@
 
 
 
-import h5py
-import subprocess
 import os
+import sys
 import shutil
+import subprocess
+import h5py
 from datetime import datetime
 import math
 import bfps
@@ -139,7 +140,8 @@ class code(base):
             raise IOError('header not there:\n' +
                           '{0}\n'.format(os.path.join(bfps.header_dir, 'base.hpp')) +
                           '{0}\n'.format(bfps.dist_loc))
-        libraries = ['bfps'] + bfps.install_info['libraries']
+        libraries = ['bfps']
+        libraries += bfps.install_info['libraries']
 
         command_strings = ['g++']
         command_strings += [self.name + '.cpp', '-o', self.name]

bfps/cpp/Lagrange_polys.hpp

@@ -28,14 +28,14 @@
 
 #define LAGRANGE_POLYS
 
-void beta_Lagrange_n1(int deriv, double x, double *poly_val);
-void beta_Lagrange_n2(int deriv, double x, double *poly_val);
-void beta_Lagrange_n3(int deriv, double x, double *poly_val);
-void beta_Lagrange_n4(int deriv, double x, double *poly_val);
-void beta_Lagrange_n5(int deriv, double x, double *poly_val);
-void beta_Lagrange_n6(int deriv, double x, double *poly_val);
-void beta_Lagrange_n7(int deriv, double x, double *poly_val);
-void beta_Lagrange_n8(int deriv, double x, double *poly_val);
+void beta_Lagrange_n1(const int deriv, const double x, double *__restrict__ poly_val);
+void beta_Lagrange_n2(const int deriv, const double x, double *__restrict__ poly_val);
+void beta_Lagrange_n3(const int deriv, const double x, double *__restrict__ poly_val);
+void beta_Lagrange_n4(const int deriv, const double x, double *__restrict__ poly_val);
+void beta_Lagrange_n5(const int deriv, const double x, double *__restrict__ poly_val);
+void beta_Lagrange_n6(const int deriv, const double x, double *__restrict__ poly_val);
+void beta_Lagrange_n7(const int deriv, const double x, double *__restrict__ poly_val);
+void beta_Lagrange_n8(const int deriv, const double x, double *__restrict__ poly_val);
 
 #endif//LAGRANGE_POLYS
 

bfps/cpp/fftw_tools.cpp

@@ -28,9 +28,7 @@
 #include "base.hpp"
 #include "fftw_tools.hpp"
 
-//#ifdef NDEBUG
-//#undef NDEBUG
-//#endif//NDEBUG
+#define NDEBUG
 
 template <class rnumber>
 int clip_zero_padding(

bfps/cpp/field_descriptor.cpp

@@ -22,10 +22,9 @@
 *                                                                     *
 **********************************************************************/
 
-// code is generally compiled via setuptools, therefore NDEBUG is present
-//#ifdef NDEBUG
-//#undef NDEBUG
-//#endif//NDEBUG
+
+
+#define NDEBUG
 
 #include <stdlib.h>
 #include <algorithm>

bfps/cpp/fluid_solver.cpp

@@ -22,10 +22,9 @@
 *                                                                     *
 **********************************************************************/
 
-// code is generally compiled via setuptools, therefore NDEBUG is present
-//#ifdef NDEBUG
-//#undef NDEBUG
-//#endif//NDEBUG
+
+
+#define NDEBUG
 
 #include <cassert>
 #include <cmath>
@@ -71,8 +70,8 @@ fluid_solver<R>::fluid_solver( \
     this->cvorticity = FFTW(alloc_complex)(this->cd->local_size);\
     this->cvelocity  = FFTW(alloc_complex)(this->cd->local_size);\
     this->rvorticity = FFTW(alloc_real)(this->cd->local_size*2);\
-    this->rvelocity  = (R*)(this->cvelocity);\
-    /*this->rvelocity  = FFTW(alloc_real)(this->cd->local_size*2);*/\
+    /*this->rvelocity  = (R*)(this->cvelocity);*/\
+    this->rvelocity  = FFTW(alloc_real)(this->cd->local_size*2);\
  \
     this->ru = this->rvelocity;\
     this->cu = this->cvelocity;\
@@ -156,6 +155,7 @@ fluid_solver<R>::fluid_solver( \
     /* initialization of fields must be done AFTER planning */ \
     std::fill_n((R*)this->cvorticity, this->cd->local_size*2, 0.0); \
     std::fill_n((R*)this->cvelocity, this->cd->local_size*2, 0.0); \
+    std::fill_n(this->rvelocity, this->cd->local_size*2, 0.0); \
     std::fill_n(this->rvorticity, this->cd->local_size*2, 0.0); \
     std::fill_n((R*)this->cv[1], this->cd->local_size*2, 0.0); \
     std::fill_n(this->rv[1], this->cd->local_size*2, 0.0); \
@@ -191,6 +191,7 @@ fluid_solver<R>::~fluid_solver() \
     FFTW(free)(this->cvorticity);\
     FFTW(free)(this->rvorticity);\
     FFTW(free)(this->cvelocity);\
+    FFTW(free)(this->rvelocity);\
     DEBUG_MSG("freed arrays\n"); \
     DEBUG_MSG("exiting ~fluid_solver\n"); \
 } \
@@ -200,6 +201,7 @@ void fluid_solver<R>::compute_vorticity() \
 { \
     ptrdiff_t tindex; \
     CLOOP_K2( \
+            this, \
             tindex = 3*cindex; \
             if (k2 <= this->kM2) \
             { \
@@ -221,6 +223,7 @@ void fluid_solver<R>::compute_velocity(FFTW(complex) *vorticity) \
 { \
     ptrdiff_t tindex; \
     CLOOP_K2( \
+            this, \
             tindex = 3*cindex; \
             if (k2 <= this->kM2 && k2 > 0) \
             { \
@@ -288,6 +291,7 @@ void fluid_solver<R>::add_forcing(\
     { \
         double knorm; \
         CLOOP( \
+                this, \
                 knorm = sqrt(this->kx[xindex]*this->kx[xindex] + \
                              this->ky[yindex]*this->ky[yindex] + \
                              this->kz[zindex]*this->kz[zindex]); \
@@ -328,6 +332,7 @@ void fluid_solver<R>::omega_nonlin( \
     this->dealias(this->cu, 3); \
     /* $\imath k \times Fourier(u \times \omega)$ */ \
     CLOOP( \
+            this, \
             tindex = 3*cindex; \
             { \
                 tmp[0][0] = -(this->ky[yindex]*this->cu[tindex+2][1] - this->kz[zindex]*this->cu[tindex+1][1]); \
@@ -351,6 +356,7 @@ void fluid_solver<R>::step(double dt) \
     std::fill_n((R*)this->cv[1], this->cd->local_size*2, 0.0); \
     this->omega_nonlin(0); \
     CLOOP_K2( \
+            this, \
             if (k2 <= this->kM2) \
             { \
                 factor0 = exp(-this->nu * k2 * dt); \
@@ -362,6 +368,7 @@ void fluid_solver<R>::step(double dt) \
  \
     this->omega_nonlin(1); \
     CLOOP_K2( \
+            this, \
             if (k2 <= this->kM2) \
             { \
                 factor0 = exp(-this->nu * k2 * dt/2); \
@@ -375,6 +382,7 @@ void fluid_solver<R>::step(double dt) \
  \
     this->omega_nonlin(2); \
     CLOOP_K2( \
+            this, \
             if (k2 <= this->kM2) \
             { \
                 factor0 = exp(-this->nu * k2 * dt * 0.5); \
@@ -459,11 +467,13 @@ void fluid_solver<R>::compute_Eulerian_acceleration(R *acceleration) \
 { \
     this->omega_nonlin(0); \
     CLOOP_K2( \
+            this, \
             for (int cc=0; cc<3; cc++) \
                 for (int i=0; i<2; i++) \
                     this->cv[1][3*cindex+cc][i] = this->cu[3*cindex+cc][i] - this->nu*k2*this->cv[0][3*cindex+cc][i]; \
             ); \
     CLOOP_K2( \
+            this, \
             if (k2 <= this->kM2 && k2 > 0) \
             { \
                 this->cv[2][3*cindex+0][0] = -(this->ky[yindex]*this->cv[1][3*cindex+2][1] - this->kz[zindex]*this->cv[1][3*cindex+1][1]) / k2; \
@@ -500,6 +510,7 @@ void fluid_solver<R>::compute_pressure(FFTW(complex) *pressure) \
     FFTW(execute)(*((FFTW(plan)*)this->vr2c[1])); \
     this->dealias(this->cv[1], 3); \
     CLOOP_K2( \
+            this, \
             if (k2 <= this->kM2 && k2 > 0) \
             { \
                 tindex = 3*cindex; \
@@ -524,6 +535,7 @@ void fluid_solver<R>::compute_pressure(FFTW(complex) *pressure) \
     FFTW(execute)(*((FFTW(plan)*)this->vr2c[1])); \
     this->dealias(this->cv[1], 3); \
     CLOOP_K2( \
+            this, \
             if (k2 <= this->kM2 && k2 > 0) \
             { \
                 tindex = 3*cindex; \
@@ -539,41 +551,23 @@ void fluid_solver<R>::compute_pressure(FFTW(complex) *pressure) \
 } \
  \
 template<> \
-void fluid_solver<R>::compute_vel_gradient(FFTW(complex) *A) \
+void fluid_solver<R>::compute_gradient_statistics( \
+        FFTW(complex) *vec, \
+        double *gradu_moments, \
+        double *trS2QR_moments, \
+        ptrdiff_t *gradu_hist, \
+        ptrdiff_t *trS2QR_hist, \
+        ptrdiff_t *QR2D_hist, \
+        double trS2QR_max_estimates[], \
+        double gradu_max_estimates[], \
+        int nbins, \
+        int QR2D_nbins) \
 { \
-    ptrdiff_t tindex; \
-    this->compute_velocity(this->cvorticity); \
-    std::fill_n((R*)A, 3*2*this->cd->local_size, 0.0); \
-    FFTW(complex) *dx_u, *dy_u, *dz_u; \
-    dx_u = A; \
-    dy_u = A + this->cd->local_size; \
-    dz_u = A + 2*this->cd->local_size; \
-    CLOOP_K2( \
-            if (k2 <= this->kM2) \
-            { \
-                tindex = 3*cindex; \
-                for (int cc=0; cc<3; cc++) \
-                { \
-                    dx_u[tindex + cc][0] = -this->kx[xindex]*this->cu[tindex+cc][1]; \
-                    dx_u[tindex + cc][1] =  this->kx[xindex]*this->cu[tindex+cc][0]; \
-                    dy_u[tindex + cc][0] = -this->ky[yindex]*this->cu[tindex+cc][1]; \
-                    dy_u[tindex + cc][1] =  this->ky[yindex]*this->cu[tindex+cc][0]; \
-                    dz_u[tindex + cc][0] = -this->kz[zindex]*this->cu[tindex+cc][1]; \
-                    dz_u[tindex + cc][1] =  this->kz[zindex]*this->cu[tindex+cc][0]; \
-                } \
-            } \
-            ); \
-} \
- \
-template<> \
-void fluid_solver<R>::compute_trS2(R *trS2) \
-{ \
-    ptrdiff_t tindex; \
     FFTW(complex) *ca; \
     R *ra; \
     ca = FFTW(alloc_complex)(this->cd->local_size*3); \
     ra = (R*)(ca); \
-    this->compute_vel_gradient(ca); \
+    this->compute_vector_gradient(ca, vec); \
     for (int cc=0; cc<3; cc++) \
     { \
         std::copy( \
@@ -587,20 +581,93 @@ void fluid_solver<R>::compute_trS2(R *trS2) \
                 ra + cc*this->cd->local_size*2); \
     } \
     /* velocity gradient is now stored, in real space, in ra */ \
+    std::fill_n(this->rv[1], 2*this->cd->local_size, 0.0); \
     R *dx_u, *dy_u, *dz_u; \
     dx_u = ra; \
     dy_u = ra + 2*this->cd->local_size; \
     dz_u = ra + 4*this->cd->local_size; \
+    double binsize[2]; \
+    double tmp_max_estimate[3]; \
+    tmp_max_estimate[0] = trS2QR_max_estimates[0]; \
+    tmp_max_estimate[1] = trS2QR_max_estimates[1]; \
+    tmp_max_estimate[2] = trS2QR_max_estimates[2]; \
+    binsize[0] = 2*tmp_max_estimate[2] / QR2D_nbins; \
+    binsize[1] = 2*tmp_max_estimate[1] / QR2D_nbins; \
+    ptrdiff_t *local_hist = new ptrdiff_t[QR2D_nbins*QR2D_nbins]; \
+    std::fill_n(local_hist, QR2D_nbins*QR2D_nbins, 0); \
     RLOOP( \
             this, \
-            tindex = 3*rindex; \
-            trS2[rindex] = (dx_u[tindex+0]*dx_u[tindex+0] + \
-                            dy_u[tindex+1]*dy_u[tindex+1] + \
-                            dz_u[tindex+2]*dz_u[tindex+2] + \
-                          ((dx_u[tindex+1]+dy_u[tindex+0])*(dx_u[tindex+1]+dy_u[tindex+0]) + \
-                           (dy_u[tindex+2]+dz_u[tindex+1])*(dy_u[tindex+2]+dz_u[tindex+1]) + \
-                           (dz_u[tindex+0]+dx_u[tindex+2])*(dz_u[tindex+0]+dx_u[tindex+2]))/2) / this->normalization_factor; \
+            R AxxAxx; \
+            R AyyAyy; \
+            R AzzAzz; \
+            R AxyAyx; \
+            R AyzAzy; \
+            R AzxAxz; \
+            R Sxy; \
+            R Syz; \
+            R Szx; \
+            ptrdiff_t tindex = 3*rindex; \
+            AxxAxx = dx_u[tindex+0]*dx_u[tindex+0]; \
+            AyyAyy = dy_u[tindex+1]*dy_u[tindex+1]; \
+            AzzAzz = dz_u[tindex+2]*dz_u[tindex+2]; \
+            AxyAyx = dx_u[tindex+1]*dy_u[tindex+0]; \
+            AyzAzy = dy_u[tindex+2]*dz_u[tindex+1]; \
+            AzxAxz = dz_u[tindex+0]*dx_u[tindex+2]; \
+            this->rv[1][tindex+1] = - (AxxAxx + AyyAyy + AzzAzz)/2 - AxyAyx - AyzAzy - AzxAxz; \
+            this->rv[1][tindex+2] = - (dx_u[tindex+0]*(AxxAxx/3 + AxyAyx + AzxAxz) + \
+                                       dy_u[tindex+1]*(AyyAyy/3 + AxyAyx + AyzAzy) + \
+                                       dz_u[tindex+2]*(AzzAzz/3 + AzxAxz + AyzAzy) + \
+                                       dx_u[tindex+1]*dy_u[tindex+2]*dz_u[tindex+0] + \
+                                       dx_u[tindex+2]*dy_u[tindex+0]*dz_u[tindex+1]); \
+            int bin0 = int(floor((this->rv[1][tindex+2] + tmp_max_estimate[2]) / binsize[0])); \
+            int bin1 = int(floor((this->rv[1][tindex+1] + tmp_max_estimate[1]) / binsize[1])); \
+            if ((bin0 >= 0 && bin0 < QR2D_nbins) && \
+                (bin1 >= 0 && bin1 < QR2D_nbins)) \
+                local_hist[bin1*QR2D_nbins + bin0]++; \
+            Sxy = dx_u[tindex+1]+dy_u[tindex+0]; \
+            Syz = dy_u[tindex+2]+dz_u[tindex+1]; \
+            Szx = dz_u[tindex+0]+dx_u[tindex+2]; \
+            this->rv[1][tindex] = (AxxAxx + AyyAyy + AzzAzz + \
+                                   (Sxy*Sxy + Syz*Syz + Szx*Szx)/2); \
             ); \
+    MPI_Allreduce( \
+            local_hist, \
+            QR2D_hist, \
+            QR2D_nbins * QR2D_nbins, \
+            MPI_INT64_T, MPI_SUM, this->cd->comm); \
+    delete[] local_hist; \
+    this->compute_rspace_stats3( \
+            this->rv[1], \
+            trS2QR_moments, \
+            trS2QR_hist, \
+            tmp_max_estimate, \
+            nbins); \
+    double *tmp_moments = new double[10*3]; \
+    ptrdiff_t *tmp_hist = new ptrdiff_t[nbins*3]; \
+    for (int cc=0; cc<3; cc++) \
+    { \
+        tmp_max_estimate[0] = gradu_max_estimates[cc*3 + 0]; \
+        tmp_max_estimate[1] = gradu_max_estimates[cc*3 + 1]; \
+        tmp_max_estimate[2] = gradu_max_estimates[cc*3 + 2]; \
+        this->compute_rspace_stats3( \
+                dx_u, \
+                tmp_moments, \
+                tmp_hist, \
+                tmp_max_estimate, \
+                nbins); \
+        for (int n = 0; n < 10; n++) \
+        for (int i = 0; i < 3 ; i++) \
+        { \
+            gradu_moments[(n*3 + cc)*3 + i] = tmp_moments[n*3 + i]; \
+        } \
+        for (int n = 0; n < nbins; n++) \
+        for (int i = 0; i < 3; i++) \
+        { \
+            gradu_hist[(n*3 + cc)*3 + i] = tmp_hist[n*3 + i]; \
+        } \
+    } \
+    delete[] tmp_moments; \
+    delete[] tmp_hist; \
     FFTW(free)(ca); \
 } \
  \
@@ -612,34 +679,11 @@ void fluid_solver<R>::compute_Lagrangian_acceleration(R *acceleration) \
     pressure = FFTW(alloc_complex)(this->cd->local_size/3); \
     this->compute_velocity(this->cvorticity); \
     this->ift_velocity(); \
-    /*clip_zero_padding<R>(this->rd, this->rvelocity, 3); \
-    std::copy( \
-            this->rvelocity, \
-            this->rvelocity + this->rd->local_size, \
-            acceleration); \
-    return; */\
     this->compute_pressure(pressure); \
-    /*this->compute_trS2(this->ru); \
-    RLOOP( \
-            this, \
-            tindex = 3*rindex; \
-            this->rv[1][tindex+0] = this->ru[rindex]; \
-            this->rv[1][tindex+1] = (this->rv[0][tindex+0]*this->rv[0][tindex+0] + \
-                                     this->rv[0][tindex+1]*this->rv[0][tindex+1] + \
-                                     this->rv[0][tindex+2]*this->rv[0][tindex+2])/(2*this->normalization_factor); \
-             ); \
-    FFTW(execute)(*((FFTW(plan)*)this->vr2c[1])); \
-    CLOOP_K2( \
-            if (k2 < this->kM2 && k2 > 0) \
-            { \
-            tindex = 3*cindex; \
-            pressure[cindex][0] = (this->cv[1][tindex+0][0] - this->cv[1][tindex+1][0]) / k2; \
-            pressure[cindex][1] = (this->cv[1][tindex+0][1] - this->cv[1][tindex+1][1]) / k2; \
-            } \
-            );*/ \
     this->compute_velocity(this->cvorticity); \
     std::fill_n((R*)this->cv[1], 2*this->cd->local_size, 0.0); \
     CLOOP_K2( \
+            this, \
             if (k2 <= this->kM2) \
             { \
                 tindex = 3*cindex; \
@@ -668,51 +712,6 @@ void fluid_solver<R>::compute_Lagrangian_acceleration(R *acceleration) \
             this->rv[1], \
             this->rv[1] + 2*this->cd->local_size, \
             acceleration); \
-    /* ********* */ \
-    /* debugging */ \
-    /* check k2*p = -(trS2 - vorticity^2/2) */ \
-    /*this->compute_trS2(this->ru); \
-    RLOOP( \
-            this, \
-            tindex = 3*rindex; \
-            this->rv[1][tindex+0] = this->ru[rindex]; \
-            this->rv[1][tindex+1] = (this->rv[0][tindex+0]*this->rv[0][tindex+0] + \
-                                     this->rv[0][tindex+1]*this->rv[0][tindex+1] + \
-                                     this->rv[0][tindex+2]*this->rv[0][tindex+2])/(2*this->normalization_factor); \
-             ); \
-    this->dealias(this->cu, 3); \
-    FFTW(execute)(*((FFTW(plan)*)this->vr2c[1])); \