diff --git a/bfps/cpp/fluid_solver.cpp b/bfps/cpp/fluid_solver.cpp
index fb57f747f4262fb1df1ea27271561ab0cacba3b7..04b4aa5a1c88ee40bc70b0ad62b88f6d35935280 100644
--- a/bfps/cpp/fluid_solver.cpp
+++ b/bfps/cpp/fluid_solver.cpp
@@ -202,6 +202,7 @@ void fluid_solver<R>::compute_vorticity() \
{ \
ptrdiff_t tindex; \
CLOOP_K2( \
+ this, \
tindex = 3*cindex; \
if (k2 <= this->kM2) \
{ \
@@ -223,6 +224,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) \
{ \
@@ -290,6 +292,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]); \
@@ -330,6 +333,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]); \
@@ -353,6 +357,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); \
@@ -364,6 +369,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); \
@@ -377,6 +383,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); \
@@ -461,11 +468,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; \
@@ -502,6 +511,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; \
@@ -526,6 +536,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; \
@@ -655,6 +666,7 @@ void fluid_solver<R>::compute_Lagrangian_acceleration(R *acceleration) \
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; \
diff --git a/bfps/cpp/fluid_solver_base.cpp b/bfps/cpp/fluid_solver_base.cpp
index 9b5a23c19d745d4b01d3e597eb15460e0a456895..d570f80f98d3be1a990c56592b192fc97478e21a 100644
--- a/bfps/cpp/fluid_solver_base.cpp
+++ b/bfps/cpp/fluid_solver_base.cpp
@@ -70,6 +70,7 @@ void fluid_solver_base<rnumber>::cospectrum(cnumber *a, cnumber *b, double *spec
std::fill_n(cospec_local, this->nshells*9, 0);
int tmp_int;
CLOOP_K2_NXMODES(
+ this,
if (k2 <= this->kM2)
{
tmp_int = int(sqrt(k2)/this->dk)*9;
@@ -98,6 +99,7 @@ void fluid_solver_base<rnumber>::cospectrum(cnumber *a, cnumber *b, double *spec
double factor = 1;
int tmp_int;
CLOOP_K2_NXMODES(
+ this,
if (k2 <= this->kM2)
{
factor = nxmodes*pow(k2, k2exponent);
@@ -316,6 +318,7 @@ fluid_solver_base<R>::fluid_solver_base( \
std::fill_n(nshell_local, this->nshells, 0.0); \
double knorm; \
CLOOP_K2_NXMODES( \
+ this, \
if (k2 < this->kM2) \
{ \
knorm = sqrt(k2); \
@@ -366,6 +369,7 @@ void fluid_solver_base<R>::low_pass_Fourier(FFTW(complex) *a, const int howmany,
const int howmany2 = 2*howmany; \
/*DEBUG_MSG("entered low_pass_Fourier, kmax=%lg km2=%lg howmany2=%d\n", kmax, km2, howmany2);*/ \
CLOOP_K2( \
+ this, \
/*DEBUG_MSG("kx=%lg ky=%lg kz=%lg k2=%lg\n", \
this->kx[xindex], \
this->ky[yindex], \
@@ -386,6 +390,7 @@ void fluid_solver_base<R>::dealias(FFTW(complex) *a, const int howmany) \
} \
double tval; \
CLOOP_K2( \
+ this, \
tval = this->Fourier_filter[int(round(k2/this->dk2))]; \
for (int tcounter = 0; tcounter < howmany; tcounter++) \
for (int i=0; i<2; i++) \
@@ -398,6 +403,7 @@ void fluid_solver_base<R>::force_divfree(FFTW(complex) *a) \
{ \
FFTW(complex) tval; \
CLOOP_K2( \
+ this, \
if (k2 > 0) \
{ \
tval[0] = (this->kx[xindex]*((*(a + cindex*3 ))[0]) + \
@@ -428,6 +434,7 @@ void fluid_solver_base<R>::compute_vector_gradient(FFTW(complex) *A, FFTW(comple
dy_u = A + this->cd->local_size; \
dz_u = A + 2*this->cd->local_size; \
CLOOP_K2( \
+ this, \
if (k2 <= this->kM2) \
{ \
tindex = 3*cindex; \
diff --git a/bfps/cpp/fluid_solver_base.hpp b/bfps/cpp/fluid_solver_base.hpp
index 4d015b0bec2853c534cf0546c2d635cf18568e1a..4cf7c2725ecf887932973ec97322a1a160ac4ee2 100644
--- a/bfps/cpp/fluid_solver_base.hpp
+++ b/bfps/cpp/fluid_solver_base.hpp
@@ -110,32 +110,32 @@ class fluid_solver_base
/* macros for loops */
/* Fourier space loop */
-#define CLOOP(expression) \
+#define CLOOP(obj, expression) \
\
{ \
ptrdiff_t cindex = 0; \
- for (ptrdiff_t yindex = 0; yindex < this->cd->subsizes[0]; yindex++) \
- for (ptrdiff_t zindex = 0; zindex < this->cd->subsizes[1]; zindex++) \
- for (ptrdiff_t xindex = 0; xindex < this->cd->subsizes[2]; xindex++) \
+ for (ptrdiff_t yindex = 0; yindex < obj->cd->subsizes[0]; yindex++) \
+ for (ptrdiff_t zindex = 0; zindex < obj->cd->subsizes[1]; zindex++) \
+ for (ptrdiff_t xindex = 0; xindex < obj->cd->subsizes[2]; xindex++) \
{ \
expression; \
cindex++; \
} \
}
-#define CLOOP_NXMODES(expression) \
+#define CLOOP_NXMODES(obj, expression) \
\
{ \
ptrdiff_t cindex = 0; \
- for (ptrdiff_t yindex = 0; yindex < this->cd->subsizes[0]; yindex++) \
- for (ptrdiff_t zindex = 0; zindex < this->cd->subsizes[1]; zindex++) \
+ for (ptrdiff_t yindex = 0; yindex < obj->cd->subsizes[0]; yindex++) \
+ for (ptrdiff_t zindex = 0; zindex < obj->cd->subsizes[1]; zindex++) \
{ \
int nxmodes = 1; \
ptrdiff_t xindex = 0; \
expression; \
cindex++; \
nxmodes = 2; \
- for (xindex = 1; xindex < this->cd->subsizes[2]; xindex++) \
+ for (xindex = 1; xindex < obj->cd->subsizes[2]; xindex++) \
{ \
expression; \
cindex++; \
@@ -143,44 +143,44 @@ class fluid_solver_base
} \
}
-#define CLOOP_K2(expression) \
+#define CLOOP_K2(obj, expression) \
\
{ \
double k2; \
ptrdiff_t cindex = 0; \
- for (ptrdiff_t yindex = 0; yindex < this->cd->subsizes[0]; yindex++) \
- for (ptrdiff_t zindex = 0; zindex < this->cd->subsizes[1]; zindex++) \
- for (ptrdiff_t xindex = 0; xindex < this->cd->subsizes[2]; xindex++) \
+ for (ptrdiff_t yindex = 0; yindex < obj->cd->subsizes[0]; yindex++) \
+ for (ptrdiff_t zindex = 0; zindex < obj->cd->subsizes[1]; zindex++) \
+ for (ptrdiff_t xindex = 0; xindex < obj->cd->subsizes[2]; xindex++) \
{ \
- k2 = (this->kx[xindex]*this->kx[xindex] + \
- this->ky[yindex]*this->ky[yindex] + \
- this->kz[zindex]*this->kz[zindex]); \
+ k2 = (obj->kx[xindex]*obj->kx[xindex] + \
+ obj->ky[yindex]*obj->ky[yindex] + \
+ obj->kz[zindex]*obj->kz[zindex]); \
expression; \
cindex++; \
} \
}
-#define CLOOP_K2_NXMODES(expression) \
+#define CLOOP_K2_NXMODES(obj, expression) \
\
{ \
double k2; \
ptrdiff_t cindex = 0; \
- for (ptrdiff_t yindex = 0; yindex < this->cd->subsizes[0]; yindex++) \
- for (ptrdiff_t zindex = 0; zindex < this->cd->subsizes[1]; zindex++) \
+ for (ptrdiff_t yindex = 0; yindex < obj->cd->subsizes[0]; yindex++) \
+ for (ptrdiff_t zindex = 0; zindex < obj->cd->subsizes[1]; zindex++) \
{ \
int nxmodes = 1; \
ptrdiff_t xindex = 0; \
- k2 = (this->kx[xindex]*this->kx[xindex] + \
- this->ky[yindex]*this->ky[yindex] + \
- this->kz[zindex]*this->kz[zindex]); \
+ k2 = (obj->kx[xindex]*obj->kx[xindex] + \
+ obj->ky[yindex]*obj->ky[yindex] + \
+ obj->kz[zindex]*obj->kz[zindex]); \
expression; \
cindex++; \
nxmodes = 2; \
- for (xindex = 1; xindex < this->cd->subsizes[2]; xindex++) \
+ for (xindex = 1; xindex < obj->cd->subsizes[2]; xindex++) \
{ \
- k2 = (this->kx[xindex]*this->kx[xindex] + \
- this->ky[yindex]*this->ky[yindex] + \
- this->kz[zindex]*this->kz[zindex]); \
+ k2 = (obj->kx[xindex]*obj->kx[xindex] + \
+ obj->ky[yindex]*obj->ky[yindex] + \
+ obj->kz[zindex]*obj->kz[zindex]); \
expression; \
cindex++; \
} \