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.