sdsc/p3dfft

Problem with definition of wave vectors

Opened this issue · 0 comments

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.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);
![sheared potential](https://github.com/sdsc/p3dfft/assets/63067010/2546ead5-52fc-4f14-b0a5-81587c91b409)