broken --- add gradient statistics

at the moment it's broken because I needed to change the prototype of
rspace_stats.

bfps/cpp/fluid_solver.cpp

@@ -539,41 +539,20 @@ 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 *moments, \
+        ptrdiff_t *hist, \
+        ptrdiff_t *QR2D_hist, \
+        double max_estimate, \
+        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,21 +566,85 @@ 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 */ \
+    R *vals_trS2 = this->rv[1]; \
+    R *vals_Q = this->rv[1] + 2*this->cd->local_size/3; \
+    R *vals_R = this->rv[1] + 4*this->cd->local_size/3; \
+    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 norm3half = this->normalization_factor*sqrt(this->normalization_factor); \
+    double binsize[2]; \
+    double tmp_max_estimate[4]; \
+    tmp_max_estimate[0] = max_estimate; \
+    tmp_max_estimate[1] = max_estimate; \
+    tmp_max_estimate[2] = max_estimate*sqrt(max_estimate); \
+    tmp_max_estimate[3] = 1.0; /* 4th component is gonna be disregarded anyway... */ \
+    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]; \
+            vals_Q[rindex] = (- ((AxxAxx + AyyAyy + AzzAzz)/2 - AxyAyx - AyzAzy - AzxAxz) / \
+                                this->normalization_factor); \
+            vals_R[rindex] = - (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]) / norm3half; \
+            int bin0 = int((vals_R[rindex] + tmp_max_estimate[2]) / binsize[0]); \
+            int bin1 = int((vals_Q[rindex] + 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]; \
+            vals_trS2[rindex] = ((AxxAxx + AyyAyy + AzzAzz + \
+                                  (Sxy*Sxy + Syz*Syz + Szx*Szx)/2) / \
+                                 this->normalization_factor); \
             ); \
     FFTW(free)(ca); \
+    MPI_Allreduce( \
+            local_hist, \
+            QR2D_hist, \
+            QR2D_nbins * QR2D_nbins, \
+            MPI_INT64_T, MPI_SUM, this->cd->comm); \
+    delete[] local_hist; \
+    double *tmp_moments = new double[4*10]; \
+    ptrdiff_t *tmp_hist = new ptrdiff_t[4*nbins]; \
+    this->compute_rspace_stats( \
+            this->rv[1], \
+            tmp_moments, \
+            tmp_hist, \
+            tmp_max_estimate, \
+            nbins); \
+    for (int i=0; i<10; i++) \
+        for (int j=0; j<3; j++) \
+            moments[i*3+j] = tmp_moments[i*4+j]; \
+    delete[] tmp_moments; \
+    for (int i=0; i<nbins; i++) \
+        for (int j=0; j<3; j++) \
+            hist[i*3+j] = tmp_hist[i*4+j]; \
+    delete[] tmp_hist; \
 } \
  \
 template<> \
@@ -612,31 +655,7 @@ 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( \
@@ -668,51 +687,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])); \
-    CLOOP_K2( \
-            if (k2 < this->kM2 && k2 > 0) \
-            { \
-                tindex = 3*cindex; \
-                this->cv[1][tindex+2][0] = k2*pressure[cindex][0]; \
-                this->cv[1][tindex+2][1] = k2*pressure[cindex][1]; \
-            } \
-            );*/ \
-    /*char fname[256]; \
-    this->fill_up_filename("pressure_trS2_enstrophy", fname); \
-    this->cd->write(fname, this->cv[1]);*/ \
-    /*double difference = 0.0; \
-    double itval, rtval; \
-    CLOOP_K2_NXMODES( \
-            if (k2 < this->kM2 && k2 > 0) \
-            { \
-            tindex = 3*cindex; \
-            rtval = k2*pressure[cindex][0] - this->cv[1][tindex+0][0] + this->cv[1][tindex+1][0]; \
-            itval = k2*pressure[cindex][1] - this->cv[1][tindex+0][1] + this->cv[1][tindex+1][1]; \
-            difference += nxmodes*(itval*itval + rtval*rtval); \
-            } \
-            ); \
-    itval = difference; \
-    MPI_Allreduce( \
-            &itval, \
-            &difference, \
-            1, \
-            MPI_DOUBLE, MPI_SUM, this->cd->comm); \
-    if (myrank == 0) DEBUG_MSG("difference is %g\n", difference);*/ \
-    /* debugging */ \
-    /* ********* */ \
     FFTW(free)(pressure); \
 } \
 

bfps/cpp/fluid_solver.hpp

@@ -82,11 +82,18 @@ class fluid_solver:public fluid_solver_base<rnumber>
                 unsigned FFTW_PLAN_RIGOR = FFTW_MEASURE);
         ~fluid_solver(void);
 
+        void compute_gradient_statistics(
+                rnumber (*vec)[2],
+                double *moments,
+                ptrdiff_t *histograms_1D,
+                ptrdiff_t *histogram_QR2D,
+                double trS2_max_estimate = 1.0,
+                int nbins_1D = 256,
+                int nbins_2D = 64);
+
         void compute_vorticity(void);
         void compute_velocity(rnumber (*vorticity)[2]);
         void compute_pressure(rnumber (*pressure)[2]);
-        void compute_vel_gradient(rnumber (*A)[2]);
-        void compute_trS2(rnumber *trS2);
         void compute_Eulerian_acceleration(rnumber *dst);
         void compute_Lagrangian_acceleration(rnumber *dst);
         void ift_velocity();

bfps/cpp/fluid_solver_base.cpp

@@ -130,17 +130,18 @@ void fluid_solver_base<rnumber>::compute_rspace_stats(
         rnumber *a,
         double *moments,
         ptrdiff_t *hist,
-        double max_estimate,
+        double max_estimate[],
         int nbins)
 {
     double *local_moments = fftw_alloc_real(10*4);
-    double val_tmp[4], binsize, pow_tmp[4];
+    double val_tmp[4], binsize[4], pow_tmp[4];
     ptrdiff_t *local_hist = new ptrdiff_t[nbins*4];
     int bin;
-    binsize = 2*max_estimate / nbins;
+    for (int i=0; i<4; i++)
+        binsize[i] = 2*max_estimate[i] / nbins;
     std::fill_n(local_hist, nbins*4, 0);
     std::fill_n(local_moments, 10*4, 0);
-    local_moments[3] = max_estimate;
+    local_moments[3] = max_estimate[3];
     RLOOP(
         this,
         std::fill_n(pow_tmp, 4, 1.0);
@@ -155,7 +156,7 @@ void fluid_solver_base<rnumber>::compute_rspace_stats(
             local_moments[0*4+3] = val_tmp[3];
         if (val_tmp[3] > local_moments[9*4+3])
             local_moments[9*4+3] = val_tmp[3];
-        bin = int(val_tmp[3]*2/binsize);
+        bin = int(val_tmp[3]*2/binsize[3]);
         if (bin >= 0 && bin < nbins)
             local_hist[bin*4+3]++;
         for (int i=0; i<3; i++)
@@ -164,7 +165,7 @@ void fluid_solver_base<rnumber>::compute_rspace_stats(
                 local_moments[0*4+i] = val_tmp[i];
             if (val_tmp[i] > local_moments[9*4+i])
                 local_moments[9*4+i] = val_tmp[i];
-            bin = int((val_tmp[i] + max_estimate) / binsize);
+            bin = int((val_tmp[i] + max_estimate[i]) / binsize[i]);
             if (bin >= 0 && bin < nbins)
                 local_hist[bin*4+i]++;
         }
@@ -418,6 +419,32 @@ void fluid_solver_base<R>::force_divfree(FFTW(complex) *a) \
 } \
  \
 template<> \
+void fluid_solver_base<R>::compute_vector_gradient(FFTW(complex) *A, FFTW(complex) *cvec) \
+{ \
+    ptrdiff_t tindex; \
+    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]*cvec[tindex+cc][1]; \
+                    dx_u[tindex + cc][1] =  this->kx[xindex]*cvec[tindex+cc][0]; \
+                    dy_u[tindex + cc][0] = -this->ky[yindex]*cvec[tindex+cc][1]; \
+                    dy_u[tindex + cc][1] =  this->ky[yindex]*cvec[tindex+cc][0]; \
+                    dz_u[tindex + cc][0] = -this->kz[zindex]*cvec[tindex+cc][1]; \
+                    dz_u[tindex + cc][1] =  this->kz[zindex]*cvec[tindex+cc][0]; \
+                } \
+            } \
+            ); \
+} \
+ \
+template<> \
 void fluid_solver_base<R>::symmetrize(FFTW(complex) *data, const int howmany) \
 { \
     ptrdiff_t ii, cc; \

bfps/cpp/fluid_solver_base.hpp

@@ -93,8 +93,9 @@ class fluid_solver_base
         void compute_rspace_stats(rnumber *a,
                                   double *moments,
                                   ptrdiff_t *hist,
-                                  double max_estimate = 1.0,
+                                  double max_estimate[4],
                                   int nbins = 256);
+        void compute_vector_gradient(rnumber (*A)[2], rnumber(*source)[2]);
         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);

meta/Velocity gradient.ipynb

@@ -125,7 +125,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 34,
+   "execution_count": 37,
    "metadata": {
     "collapsed": false
    },
@@ -136,7 +136,9 @@
      "text": [
       "0\n",
       "0\n",
-      "0\n"
+      "0\n",
+      "-Axx*(Axx**2/3 + Axy*Ayx + Axz*Azx) - Axy*Ayz*Azx - Axz*Ayx*Azy - Ayy*(Axy*Ayx + Ayy**2/3 + Ayz*Azy) - Azz*(Axz*Azx + Ayz*Azy + Azz**2/3)\n",
+      "Axx*Ayy*Azz - Axx*Ayz*Azy - Axy*Ayx*Azz + Axy*Ayz*Azx + Axz*Ayx*Azy - Axz*Ayy*Azx\n"
      ]
     }
    ],
@@ -156,11 +158,15 @@
     "          A[0, 1]*A[1, 2]*A[2, 0] +\n",
     "          A[0, 2]*A[1, 0]*A[2, 1])\n",
     "print sp.simplify(Ralt - R)\n",
+    "#print sp.simplify(Ralt + A.det())\n",
     "trS2alt = (AxxAxx + AyyAyy + AzzAzz +\n",
     "           ((A[0, 1] + A[1, 0])**2 +\n",
     "            (A[1, 2] + A[2, 1])**2 +\n",
     "            (A[2, 0] + A[0, 2])**2)/2)\n",
-    "print sp.simplify(trS2alt - trS2)"
+    "print sp.simplify(trS2alt - trS2)\n",
+    "\n",
+    "print Ralt\n",
+    "print A.det()"
    ]
   },
   {

setup.py

@@ -46,11 +46,11 @@ except:
 VERSION = date_name
 
 src_file_list = ['field_descriptor',
+                 'fluid_solver_base',
+                 'fluid_solver',
                  'interpolator',
                  'particles',
                  'fftw_tools',
-                 'fluid_solver_base',
-                 'fluid_solver',
                  'spline_n1',
                  'spline_n2',
                  'spline_n3',