diff --git a/TurTLE/test/test_turtle_NSVEparticles.py b/TurTLE/test/test_turtle_NSVEparticles.py
index 8329a0b5ceebdb0f8477fbb45dc05645891a774e..4ed1d9968ae0288fdcb05f26dd072575435770c4 100644
--- a/TurTLE/test/test_turtle_NSVEparticles.py
+++ b/TurTLE/test/test_turtle_NSVEparticles.py
@@ -89,6 +89,7 @@ def main():
y0 = f0['tracers0/rhs/{0}'.format(iteration)][...]
y1 = f1['tracers0/rhs/{0}'.format(iteration)][...]
rhs_error = np.max(np.abs(y0 - y1))
+ print(field_error, traj_error, rhs_error)
assert(field_error < 1e-5)
assert(traj_error < 1e-5)
assert(rhs_error < 1e-5)
diff --git a/cpp/field.cpp b/cpp/field.cpp
index 2fd79ce4fe7d25f2f7f8bdae0a500ee40e0a8330..a5779752ee4ca232f4883042c3171ec4c45ab18a 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -1278,16 +1278,20 @@ void field<rnumber, be, fc>::symmetrize()
typename fftw_interface<rnumber>::complex *cdata = this->get_cdata();
// symmetrize kx = 0 plane, line by line, for ky != 0
MPI_Status *mpistatus = new MPI_Status;
- typename fftw_interface<rnumber>::complex *buffer = new typename fftw_interface<rnumber>::complex[ncomp(fc)*this->clayout->sizes[1]];
- //typename fftw_interface<rnumber>::complex *buffer;
- //buffer = fftw_interface<rnumber>::alloc_complex(ncomp(fc)*this->clayout->sizes[1]);
+ // bufferp will hold data from "plus", i.e. iy
+ // bufferm will hold data from "minus", i.e. ny - iy
+ typename fftw_interface<rnumber>::complex *bufferp = new typename fftw_interface<rnumber>::complex[ncomp(fc)*this->clayout->sizes[1]];
+ typename fftw_interface<rnumber>::complex *bufferm = new typename fftw_interface<rnumber>::complex[ncomp(fc)*this->clayout->sizes[1]];
/*ptrdiff_t tindex;*/
- int ranksrc, rankdst;
+ int rankp, rankm;
+ // for each ky slice except ky=0
for (ptrdiff_t iy = 1; iy < ptrdiff_t(this->clayout->sizes[0]/2); iy++)
{
- ranksrc = this->clayout->rank[0][iy];
- rankdst = this->clayout->rank[0][this->clayout->sizes[0] - iy];
- if (this->clayout->myrank == ranksrc)
+ // read rank plus and rank minus
+ rankp = this->clayout->rank[0][iy];
+ rankm = this->clayout->rank[0][this->clayout->sizes[0] - iy];
+ // if my rank is rank plus, I should fill out the plus buffer
+ if (this->clayout->myrank == rankp)
{
ptrdiff_t iyy = iy - this->clayout->starts[0];
for (ptrdiff_t iz = 0; iz < ptrdiff_t(this->clayout->sizes[1]); iz++)
@@ -1295,50 +1299,103 @@ void field<rnumber, be, fc>::symmetrize()
ptrdiff_t cindex = this->get_cindex(0, iyy, iz);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
for (int imag_comp=0; imag_comp<2; imag_comp++)
- (*(buffer + ncomp(fc)*iz+cc))[imag_comp] =
+ (*(bufferp + ncomp(fc)*iz+cc))[imag_comp] =
(*(cdata + ncomp(fc)*cindex + cc))[imag_comp];
}
}
- if (ranksrc != rankdst)
+ // if my rank is rank minus, I should fill out the minus buffer
+ if (this->clayout->myrank == rankm)
{
- if (this->clayout->myrank == ranksrc)
- MPI_Send((void*)buffer,
+ ptrdiff_t iyy = (this->clayout->sizes[0] - iy) - this->clayout->starts[0];
+ for (ptrdiff_t iz = 0; iz < ptrdiff_t(this->clayout->sizes[1]); iz++)
+ {
+ ptrdiff_t cindex = this->get_cindex(0, iyy, iz);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++)
+ for (int imag_comp=0; imag_comp<2; imag_comp++)
+ (*(bufferm + ncomp(fc)*iz+cc))[imag_comp] =
+ (*(cdata + ncomp(fc)*cindex + cc))[imag_comp];
+ }
+ }
+ // if ranks are different, send and receive
+ if (rankp != rankm)
+ {
+ // if my rank is rank plus, send buffer plus, receive buffer minus
+ if (this->clayout->myrank == rankp)
+ {
+ MPI_Send((void*)bufferp,
ncomp(fc)*this->clayout->sizes[1],
mpi_real_type<rnumber>::complex(),
- rankdst, iy,
+ rankm, 2*iy+0,
this->clayout->comm);
- if (this->clayout->myrank == rankdst)
- MPI_Recv((void*)buffer,
+ MPI_Recv((void*)bufferm,
ncomp(fc)*this->clayout->sizes[1],
mpi_real_type<rnumber>::complex(),
- ranksrc, iy,
+ rankm, 2*iy+1,
this->clayout->comm,
mpistatus);
+ }
+ // if my rank is rank minus, receive buffer plus, send buffer minus
+ if (this->clayout->myrank == rankm)
+ {
+ MPI_Recv((void*)bufferp,
+ ncomp(fc)*this->clayout->sizes[1],
+ mpi_real_type<rnumber>::complex(),
+ rankp, 2*iy+0,
+ this->clayout->comm,
+ mpistatus);
+ MPI_Send((void*)bufferm,
+ ncomp(fc)*this->clayout->sizes[1],
+ mpi_real_type<rnumber>::complex(),
+ rankp, 2*iy+1,
+ this->clayout->comm);
+ }
+ }
+ // if I my rank is either plus or minus, I should update my slice
+ if (this->clayout->myrank == rankp)
+ {
+ ptrdiff_t iyy = iy - this->clayout->starts[0];
+ for (ptrdiff_t iz = 1; iz < ptrdiff_t(this->clayout->sizes[1]); iz++)
+ {
+ ptrdiff_t izz = (this->clayout->sizes[1] - iz);
+ ptrdiff_t cindex = this->get_cindex(0, iyy, izz);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++)
+ {
+ (*(cdata + ncomp(fc)*cindex + cc))[0] = ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*iz+cc))[0])/2;
+ (*(cdata + ncomp(fc)*cindex + cc))[1] = ((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferp + ncomp(fc)*iz+cc))[1])/2;
+ }
+ }
+ ptrdiff_t cindex = this->get_cindex(0, iyy, 0);
+ for (int cc = 0; cc < int(ncomp(fc)); cc++)
+ {
+ (*(cdata + cc + ncomp(fc)*cindex))[0] = ((*(bufferp + cc))[0] + (*(bufferm + cc))[0])/2;
+ (*(cdata + cc + ncomp(fc)*cindex))[1] = ((*(bufferp + cc))[1] - (*(bufferm + cc))[1])/2;
+ }
}
- if (this->clayout->myrank == rankdst)
+ // if I my rank is either plus or minus, I should update my slice
+ if (this->clayout->myrank == rankm)
{
ptrdiff_t iyy = (this->clayout->sizes[0] - iy) - this->clayout->starts[0];
for (ptrdiff_t iz = 1; iz < ptrdiff_t(this->clayout->sizes[1]); iz++)
{
ptrdiff_t izz = (this->clayout->sizes[1] - iz);
ptrdiff_t cindex = this->get_cindex(0, iyy, izz);
- //DEBUG_MSG("iy = %ld, iz = %ld\n", iy, iz);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
{
- (*(cdata + ncomp(fc)*cindex + cc))[0] = (*(buffer + ncomp(fc)*iz+cc))[0];
- (*(cdata + ncomp(fc)*cindex + cc))[1] = -(*(buffer + ncomp(fc)*iz+cc))[1];
+ (*(cdata + ncomp(fc)*cindex + cc))[0] = ((*(bufferp + ncomp(fc)*iz+cc))[0] + (*(bufferm + ncomp(fc)*iz+cc))[0])/2;
+ (*(cdata + ncomp(fc)*cindex + cc))[1] = -((*(bufferp + ncomp(fc)*iz+cc))[1] - (*(bufferp + ncomp(fc)*iz+cc))[1])/2;
}
}
ptrdiff_t cindex = this->get_cindex(0, iyy, 0);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
{
- (*(cdata + cc + ncomp(fc)*cindex))[0] = (*(buffer + cc))[0];
- (*(cdata + cc + ncomp(fc)*cindex))[1] = -(*(buffer + cc))[1];
+ (*(cdata + cc + ncomp(fc)*cindex))[0] = ((*(bufferp + cc))[0] + (*(bufferm + cc))[0])/2;
+ (*(cdata + cc + ncomp(fc)*cindex))[1] = -((*(bufferp + cc))[1] - (*(bufferm + cc))[1])/2;
}
}
}
//fftw_interface<rnumber>::free(buffer);
- delete[] buffer;
+ delete[] bufferp;
+ delete[] bufferm;
delete mpistatus;
// symmetrize kx = 0, ky = 0 line
if (this->clayout->myrank == this->clayout->rank[0][0])
@@ -1349,8 +1406,12 @@ void field<rnumber, be, fc>::symmetrize()
ptrdiff_t cindex1 = this->get_cindex(0, 0, this->clayout->sizes[1] - iz);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
{
- (*(cdata + cc + ncomp(fc)*cindex1))[0] = (*(cdata + cc + ncomp(fc)*cindex0))[0];
- (*(cdata + cc + ncomp(fc)*cindex1))[1] = -(*(cdata + cc + ncomp(fc)*cindex0))[1];
+ double re = ((*(cdata + cc + ncomp(fc)*cindex0))[0] + (*(cdata + cc + ncomp(fc)*cindex1))[0])/2;
+ double im = ((*(cdata + cc + ncomp(fc)*cindex0))[1] - (*(cdata + cc + ncomp(fc)*cindex1))[1])/2;
+ (*(cdata + cc + ncomp(fc)*cindex0))[0] = re;
+ (*(cdata + cc + ncomp(fc)*cindex0))[1] = im;
+ (*(cdata + cc + ncomp(fc)*cindex1))[0] = re;
+ (*(cdata + cc + ncomp(fc)*cindex1))[1] = -im;
}
}
}