P3DFFT stands for Parallel Three-Dimensional Fast Fourier Transforms. It is a library for large-scale computer simulations on parallel platforms. It implements 3D FFT and related algorithms such as Chebyshev transform (an important class of algorithm for simulations in a wide range of fields). P3DFFT uses 2D, or pencil, decomposition. For more information:
Recently, I wrote a 3D spectral solver for solving the gravitational Poisson equation thanks to P3DFFT. The solver is written in C. My P3DFFT version is compiled with STRIDE1. The sample tests work fine.
For solving the Poisson equation in Fourier space I get inspiration from the driver_spec.c test file. You can find a minimal version of my code at the end of the post. I did a simple test that consists on transforming a real density distribution, divide it by the square of the wave vector and then do the backward transform. As you can see in the attached image (left initial density, right potential ), I do not recover the initial density profile (up to a factor). The result seems rather sheared. Do you have an idea from where this might come from ? I suspect it comes from the definition of the wave vectors but I cannot see where ...
Thank you for your help and valuable time,
Steven
static int real_start[3], sdims_real[3], real_end[3];
static int cmplx_end[3];
static unsigned char op_f[]="fft", op_b[]="tff";
/* Real and complex arrays */
static double *rho_real, *rho_cmplx;
static double *phi_real, *phi_cmplx;
static int N=0;
static int Nx, Ny, Nz, Nzh;
proc_dims[0]=c.bny;
proc_dims[1]=c.bnz;
//int real_start[3], sdims_real[3], real_end[3];
//int glob_start_cmplx[3], ldims_cmplx[3], cmplx_end[3];
// Initialise P3DFFT
Cp3dfft_setup(proc_dims, Nx, Ny, Nz, MPI_Comm_c2f(MPI_COMM_WORLD), Nx, Ny, Nz, 0, memsize);
/* Get dimensions for input array - real numbers, X-pencil shape.
Note that we are following the Fortran ordering, i.e.
the dimension with stride-1 is X. */
conf=1;
Cp3dfft_get_dims(real_start, real_end, sdims_real, conf);
/* Get dimensions for output array - complex numbers, Z-pencil shape.
Stride-1 dimension could be X or Z, depending on how the library
was compiled (stride1 option) */
conf=2;
Cp3dfft_get_dims(glob_start_cmplx, cmplx_end, ldims_cmplx, conf);
size_real = sdims_real[0] *sdims_real[1] *sdims_real[2];
size_cmplx = ldims_cmplx[0]*ldims_cmplx[1]*ldims_cmplx[2];
// Now allocate initial and final arrays in physical space as real-valued 1D
// storage containing a contiguous 3D local array
rho_real = (double *) malloc(sizeof(double) * size_real);
phi_real = (double *) malloc(sizeof(double) * size_real);
rho_cmplx = (double *) malloc(sizeof(double) * size_cmplx*2 );
phi_cmplx = (double *) malloc(sizeof(double) * size_cmplx*2 );
// -> Affect values to density
/* --- FORWARD TRANSFORM FOR DENSITY ------------------------------------ */
Cp3dfft_ftran_r2c(rho_real, rho_cmplx, op_f);
for (i=0; i < 2*ldims_cmplx[2]; i++)
{
II = i/2 + glob_start_cmplx[2] - 1;
kx = 2*M_PI*II/LX;
for (j=0; j < ldims_cmplx[1]; j++)
{
JJ = glob_start_cmplx[1] - 1 + j;
if (2*JJ > Ny){
JJ = JJ-Ny;
}
ky = 2*M_PI*JJ/LY;
for (k=0; k < ldims_cmplx[0]; k++)
{
KK = glob_start_cmplx[0] - 1 + k;
if (2*KK > Nz){
KK = KK-Nz;
}
kz = 2*M_PI*KK/LZ;
index_cmplx = (i*ldims_cmplx[1] + j)*ldims_cmplx[0] + k;
if (kx*kx + ky*ky + kz*kz==0.0)
{
/* The mean value of the potential is equal to 0 */
phi_cmplx[index_cmplx] = 0.0;
} else {
phi_cmplx[index_cmplx] = rho_cmplx[index_cmplx]
* (-pi4G) / (kx*kx + ky*ky + kz*kz);
}
}
}
}
/* --- BACKWARD TRANSFORMS TO REAL SPACE ------------------------------- */
Cp3dfft_btran_c2r(phi_cmplx, phi_real, op_b);
HI,
Recently, I wrote a 3D spectral solver for solving the gravitational Poisson equation thanks to P3DFFT. The solver is written in C. My P3DFFT version is compiled with STRIDE1. The sample tests work fine.
For solving the Poisson equation in Fourier space I get inspiration from the
driver_spec.ctest file. You can find a minimal version of my code at the end of the post. I did a simple test that consists on transforming a real density distribution, divide it by the square of the wave vector and then do the backward transform. As you can see in the attached image (left initial density, right potential ), I do not recover the initial density profile (up to a factor). The result seems rather sheared. Do you have an idea from where this might come from ? I suspect it comes from the definition of the wave vectors but I cannot see where ...Thank you for your help and valuable time, Steven