For time-stepping large PDE problems, it’s common to read from one array (the previous time step), and write to a second array (the result at the new time step). Afterward the arrays are swapped. This is also known as the “flip-flop” approach, or ping-ponging. The separate arrays also make the update step easier to parallelise (no race conditions).
What is your favourite implementation approach? (Scroll down for examples)
- Separate arrays with parity flag
- Array with extra dimension
- Pointer swap
- Overwrite old array
- Other (describe in your reply)
Separate arrays with parity flag
real, allocatable :: f1(:), f2(:)
integer :: flipflop
! ... initialize f
flipflop = 1
do step = 1, nsteps
select case(flipflop)
case(1)
call sweep(fin=f1,fout=f2)
case(2)
call sweep(fin=f2,fout=f1)
end select
flipflop = 3 - flipflop
end do
If zero-based indexing is used, the parity flag can be inverted using
! 0-based indexing
flipflop = 1 - flipflop
Array with extra dimension
This one relies on array slicing or dummy argument association. Integer flags are used to denote the old and new arrays.
real, allocatable :: f(:,:)
integer :: tmp, old, new
allocate(f(n,2))
! ... initialize f
old = 1
new = 2
do step = 1, nsteps
call sweep(fin=f(:,kold),fout=f(:,knew))
tmp = old
old = new
new = tmp
end do
Alternatively, the integer flags can be calculated before the update like this:
old = mod(step+1,2) + 1
new = mod(step+2,2) + 1
or after the update, like this:
old = 3 - old
new = 3 - new
Pointer swap
real, allocatable, target :: f1(:), f2(:)
real, pointer :: fold(:), fnew(:), ftmp(:)
fold => f1
fnew => f2
do step = 1, nsteps
call sweep(fin=fold,fout=fnew)
ftmp => fold
fold => fnew
fnew => ftmp
end do
Overwrite old array
real, allocatable :: fold(:), fnew(:)
do step = 1, nsteps
call sweep(fin=fold,fout=fnew)
fold = fnew
end do