Fortran 2D grid w/ Periodic Boundary Conditions - efficient implementation

Viewed 141

I am using Fortran to solve some coupled PDEs using finite-differences and spectral methods. At a certain point those grids are fed into FFTW. I assume the domain is periodic and I am struggling to find the most efficient way to deal with the problem of indexing. Options:

1 - Ghost Cells:

I define my grids as

real :: y(0:nx+1,0:ny+1)

(do finite differences)

call update_ghost(y)

Problem : when I use FFTW I want the input to be the grid without the ghost cells, so I do:

call fftw(phi(1:nx,1:ny),phi_k)

and this leads to temporary array creation since the data is not contiguous.

2 - Indirect Adressing:

real :: y(nx,ny)

I then calculate vectors with the coordinates of the neighbours taking into account PBC then do

do j=1,ny
   do i=1,nx
      laplacian(i,j) = y(xnext(i),j) + y(xprev(i),j) + ...
   enddo
enddo

Problem: Indirect Adressing is expensive

3 - Calculate indices inside the loop

do j=1,ny
   do i=1,nx
      xnext = modulo(i+2,nx) + 1
      ynext = modulo(i,nx) + 1
      xprev = modulo(j+2,ny) + 1
      yprev = modulo(j,ny) + 1
      laplacian(i,j) = y(xnext,j) + y(xprev(i),j) + ...
   enddo
enddo

Problem: extra calculations

Does anyone have any other options? Right now I use option 1 but it kills me to think of all the temporary array creation that is happening behind the scenes.

0 Answers
Related