Skip to content
Snippets Groups Projects
Commit 68ed2280 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

still debugging

parent 1b794b9e
Branches
No related tags found
1 merge request!23WIP: Feature/use cmake
Pipeline #
Hide whitespace changes
Inline Side-by-side
bfps/cpp/field.cpp
+ 117
− 45
View file @ 68ed2280
......@@ -73,11 +73,10 @@ field<rnumber, be, fc>::field(
nfftw[0] = nz;
nfftw[1] = ny;
nfftw[2] = nx;
//ptrdiff_t tmp_local_size;
ptrdiff_t tmp_local_size;
ptrdiff_t local_n0, local_0_start;
ptrdiff_t local_n1, local_1_start;
//tmp_local_size = fftw_mpi_local_size_many_transposed(
fftw_interface<rnumber>::mpi_local_size_many_transposed(
tmp_local_size = fftw_interface<rnumber>::mpi_local_size_many_transposed(
3, nfftw, ncomp(fc),
FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK, this->comm,
&local_n0, &local_0_start,
......@@ -86,6 +85,8 @@ field<rnumber, be, fc>::field(
sizes[0] = nz; sizes[1] = ny; sizes[2] = nx;
subsizes[0] = local_n0; subsizes[1] = ny; subsizes[2] = nx;
starts[0] = local_0_start; starts[1] = 0; starts[2] = 0;
DEBUG_MSG("local_0_start = %ld, local_1_start = %ld\n",
local_0_start, local_1_start);
this->rlayout = new field_layout<fc>(
sizes, subsizes, starts, this->comm);
this->npoints = this->rlayout->full_size / ncomp(fc);
......@@ -99,6 +100,11 @@ field<rnumber, be, fc>::field(
starts[0] = local_1_start; starts[1] = 0; starts[2] = 0;
this->clayout = new field_layout<fc>(
sizes, subsizes, starts, this->comm);
DEBUG_MSG("local_size = %ld, rlayout->local_size = %ld, rmemlayout->local_size = %ld, clayout->local_size = %ld\n",
tmp_local_size,
this->rlayout->local_size,
this->rmemlayout->local_size,
this->clayout->local_size);
this->data = fftw_interface<rnumber>::alloc_real(
this->rmemlayout->local_size);
memset(this->data, 0, sizeof(rnumber)*this->rmemlayout->local_size);
......@@ -1017,77 +1023,143 @@ void field<rnumber, be, fc>::symmetrize()
{
TIMEZONE("field::symmetrize");
assert(!this->real_space_representation);
ptrdiff_t ii, cc;
typename fftw_interface<rnumber>::complex *data = this->get_cdata();
//this->ift();
//this->dft();
//this->normalize();
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;
if (this->myrank == this->clayout->rank[0][0])
{
for (cc = 0; cc < ncomp(fc); cc++)
data[cc][1] = 0.0;
for (ii = 1; ii < ptrdiff_t(this->clayout->sizes[1]/2); ii++)
for (cc = 0; cc < ncomp(fc); cc++) {
( *(data + cc + ncomp(fc)*(this->clayout->sizes[1] - ii)*this->clayout->sizes[2]))[0] =
(*(data + cc + ncomp(fc)*( ii)*this->clayout->sizes[2]))[0];
( *(data + cc + ncomp(fc)*(this->clayout->sizes[1] - ii)*this->clayout->sizes[2]))[1] =
-(*(data + cc + ncomp(fc)*( ii)*this->clayout->sizes[2]))[1];
}
}
typename fftw_interface<rnumber>::complex *buffer;
buffer = fftw_interface<rnumber>::alloc_complex(ncomp(fc)*this->clayout->sizes[1]);
ptrdiff_t yy;
/*ptrdiff_t tindex;*/
int ranksrc, rankdst;
for (yy = 1; yy < ptrdiff_t(this->clayout->sizes[0]/2); yy++) {
ranksrc = this->clayout->rank[0][yy];
rankdst = this->clayout->rank[0][this->clayout->sizes[0] - yy];
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)
for (ii = 0; ii < ptrdiff_t(this->clayout->sizes[1]); ii++)
for (cc = 0; cc < ncomp(fc); cc++)
{
ptrdiff_t iyy = 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++)
(*(buffer + ncomp(fc)*ii+cc))[imag_comp] =
(*(data + ncomp(fc)*((yy - this->clayout->starts[0])*this->clayout->sizes[1] + ii)*this->clayout->sizes[2] + cc))[imag_comp];
(*(buffer + ncomp(fc)*iz+cc))[imag_comp] =
(*(cdata + ncomp(fc)*cindex + cc))[imag_comp];
}
}
if (ranksrc != rankdst)
{
if (this->clayout->myrank == ranksrc)
MPI_Send((void*)buffer,
ncomp(fc)*this->clayout->sizes[1], mpi_real_type<rnumber>::complex(), rankdst, yy,
this->clayout->comm);
ncomp(fc)*this->clayout->sizes[1],
mpi_real_type<rnumber>::complex(),
rankdst, iy,
this->clayout->comm);
if (this->clayout->myrank == rankdst)
MPI_Recv((void*)buffer,
ncomp(fc)*this->clayout->sizes[1], mpi_real_type<rnumber>::complex(), ranksrc, yy,
this->clayout->comm, mpistatus);
ncomp(fc)*this->clayout->sizes[1],
mpi_real_type<rnumber>::complex(),
ranksrc, iy,
this->clayout->comm,
mpistatus);
}
if (this->clayout->myrank == rankdst)
{
for (ii = 1; ii < ptrdiff_t(this->clayout->sizes[1]); ii++)
for (cc = 0; cc < ncomp(fc); cc++)
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++)
{
(*(data + ncomp(fc)*((this->clayout->sizes[0] - yy - this->clayout->starts[0])*this->clayout->sizes[1] + ii)*this->clayout->sizes[2] + cc))[0] =
(*(buffer + ncomp(fc)*(this->clayout->sizes[1]-ii)+cc))[0];
(*(data + ncomp(fc)*((this->clayout->sizes[0] - yy - this->clayout->starts[0])*this->clayout->sizes[1] + ii)*this->clayout->sizes[2] + cc))[1] =
-(*(buffer + ncomp(fc)*(this->clayout->sizes[1]-ii)+cc))[1];
(*(cdata + ncomp(fc)*cindex + cc))[0] = (*(buffer + ncomp(fc)*iz+cc))[0];
(*(cdata + ncomp(fc)*cindex + cc))[1] = -(*(buffer + ncomp(fc)*iz+cc))[1];
}
for (cc = 0; cc < ncomp(fc); cc++)
}
ptrdiff_t cindex = this->get_cindex(0, iyy, 0);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
{
(*((data + cc + ncomp(fc)*(this->clayout->sizes[0] - yy - this->clayout->starts[0])*this->clayout->sizes[1]*this->clayout->sizes[2])))[0] = (*(buffer + cc))[0];
(*((data + cc + ncomp(fc)*(this->clayout->sizes[0] - yy - this->clayout->starts[0])*this->clayout->sizes[1]*this->clayout->sizes[2])))[1] = -(*(buffer + cc))[1];
(*(cdata + cc + ncomp(fc)*cindex))[0] = (*(buffer + cc))[0];
(*(cdata + cc + ncomp(fc)*cindex))[1] = -(*(buffer + cc))[1];
}
}
}
fftw_interface<rnumber>::free(buffer);
delete mpistatus;
/* put asymmetric data to 0 */
// symmetrize kx = 0, ky = 0 line
if (this->clayout->myrank == this->clayout->rank[0][0])
{
for (ptrdiff_t iz = 1; iz < ptrdiff_t(this->clayout->sizes[1]/2); iz++)
{
ptrdiff_t cindex0 = this->get_cindex(0, 0, iz);
ptrdiff_t cindex1 = this->get_cindex(0, 0, this->clayout->sizes[1] - iz);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
{
rnumber rtmp = ((*(cdata + cc + ncomp(fc)*cindex0))[0] +
(*(cdata + cc + ncomp(fc)*cindex1))[0])/2;
rnumber ctmp = ((*(cdata + cc + ncomp(fc)*cindex0))[1] -
(*(cdata + cc + ncomp(fc)*cindex1))[1])/2;
(*(cdata + cc + ncomp(fc)*cindex0))[0] = rtmp;
(*(cdata + cc + ncomp(fc)*cindex1))[0] = rtmp;
(*(cdata + cc + ncomp(fc)*cindex0))[1] = ctmp;
(*(cdata + cc + ncomp(fc)*cindex1))[1] = -ctmp;
}
}
}
// make 0 mode real
if (this->myrank == this->clayout->rank[0][0])
{
for (ptrdiff_t cc = 0; cc < ncomp(fc); cc++)
cdata[cc][1] = 0.0;
}
// put kx = nx/2 modes to 0
for (ptrdiff_t iy = 0; iy < ptrdiff_t(this->clayout->subsizes[0]); iy++)
for (ptrdiff_t iz = 0; iz < ptrdiff_t(this->clayout->subsizes[1]); iz++)
{
ptrdiff_t cindex = this->get_cindex(this->clayout->sizes[2]-1, iy, iz);
for (int cc = 0; cc < int(ncomp(fc)); cc++) {
(*(cdata + cc + ncomp(fc)*cindex))[0] = 0.0;
(*(cdata + cc + ncomp(fc)*cindex))[1] = 0.0;
}
}
// put ky = ny/2 modes to 0
if (this->clayout->myrank == this->clayout->rank[0][this->clayout->sizes[0]/2])
{
ptrdiff_t tindex = ncomp(fc)*(this->clayout->sizes[0]/2 - this->clayout->starts[0])*this->clayout->sizes[1]*this->clayout->sizes[2];
for (ii = 0; ii < ptrdiff_t(this->clayout->sizes[1]); ii++)
for (ptrdiff_t iz = 0; iz < ptrdiff_t(this->clayout->subsizes[1]); iz++)
for (ptrdiff_t ix = 0; ix < ptrdiff_t(this->clayout->subsizes[2]); ix++)
{
std::fill_n((rnumber*)(data + tindex), ncomp(fc)*2*this->clayout->sizes[2], 0.0);
tindex += ncomp(fc)*this->clayout->sizes[2];
ptrdiff_t cindex = this->get_cindex(ix, this->clayout->sizes[0]/2-this->clayout->starts[0], iz);
for (int cc = 0; cc < int(ncomp(fc)); cc++) {
(*(cdata + cc + ncomp(fc)*cindex))[0] = 0.0;
(*(cdata + cc + ncomp(fc)*cindex))[1] = 0.0;
}
}
}
// put kz = nz/2 modes to 0
for (ptrdiff_t iy = 0; iy < ptrdiff_t(this->clayout->subsizes[0]); iy++)
for (ptrdiff_t ix = 0; ix < ptrdiff_t(this->clayout->subsizes[2]); ix++)
{
ptrdiff_t cindex = this->get_cindex(ix, iy, this->clayout->sizes[1]/2);
for (int cc = 0; cc < int(ncomp(fc)); cc++) {
(*(cdata + cc + ncomp(fc)*cindex))[0] = 0.0;
(*(cdata + cc + ncomp(fc)*cindex))[1] = 0.0;
}
}
//tindex = ncomp(fc)*();
//std::fill_n((rnumber*)(data + tindex), ncomp(fc)*2, 0.0);
///* put asymmetric data to 0 */
//if (this->clayout->myrank == this->clayout->rank[0][this->clayout->sizes[0]/2])
//{
// ptrdiff_t tindex = ncomp(fc)*(this->clayout->sizes[0]/2 - this->clayout->starts[0])*this->clayout->sizes[1]*this->clayout->sizes[2];
// for (ii = 0; ii < ptrdiff_t(this->clayout->sizes[1]); ii++)
// {
// std::fill_n((rnumber*)(data + tindex), ncomp(fc)*2*this->clayout->sizes[2], 0.0);
// tindex += ncomp(fc)*this->clayout->sizes[2];
// }
//}
////tindex = ncomp(fc)*();
////std::fill_n((rnumber*)(data + tindex), ncomp(fc)*2, 0.0);
}
template <typename rnumber,
......
bfps/cpp/full_code/field_test.cpp
+ 1
− 1
View file @ 68ed2280
......@@ -50,7 +50,7 @@ int field_test<rnumber>::do_work(void)
DEFAULT_FFTW_FLAG);
std::default_random_engine rgen;
std::normal_distribution<rnumber> rdist;
rgen.seed(1);
rgen.seed(2);
//auto gaussian = std::bind(rgen, rdist);
kspace<FFTW,SMOOTH> *kk = new kspace<FFTW, SMOOTH>(
scal_field->clayout, this->dkx, this->dky, this->dkz);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment