diff --git a/bfps/cpp/fluid_solver.cpp b/bfps/cpp/fluid_solver.cpp
index ce7230e5b91e04bce51fb3c2a58d6651941562ce..4f282ce232e31f19c4124a64affb99d405834611 100644
--- a/bfps/cpp/fluid_solver.cpp
+++ b/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); \
} \
diff --git a/bfps/cpp/fluid_solver.hpp b/bfps/cpp/fluid_solver.hpp
index 9eead599b7421dbbe8171139d1533c4e9af2b096..883a3a9664730c970e8383b5a58aeb5709702507 100644
--- a/bfps/cpp/fluid_solver.hpp
+++ b/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();
diff --git a/bfps/cpp/fluid_solver_base.cpp b/bfps/cpp/fluid_solver_base.cpp
index bbdccad8f4b386007e42907bbe8ed5317c9728bd..548d5506d704c5f4383c9921d9688c8c04d9666f 100644
--- a/bfps/cpp/fluid_solver_base.cpp
+++ b/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; \
diff --git a/bfps/cpp/fluid_solver_base.hpp b/bfps/cpp/fluid_solver_base.hpp
index 64659a0488ffc576b3959970e002521dde28d8eb..1ea75eaff7ef527368194af22b6d0b8f818a6ae1 100644
--- a/bfps/cpp/fluid_solver_base.hpp
+++ b/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);
diff --git a/meta/Velocity gradient.ipynb b/meta/Velocity gradient.ipynb
index deaa7e596a699113671997b70b4d5e3e3517b722..deea6cb898882ec954e5ecf9b6ca42ba7214172a 100644
--- a/meta/Velocity gradient.ipynb
+++ b/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()"
]
},
{
diff --git a/setup.py b/setup.py
index 76de4ad3f5e3122ad46e3bd3c6702aeb5a47f83a..ae97c93421499a46f260781ed8a1b0722ec2b42a 100644
--- a/setup.py
+++ b/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',