Tips for bookkeeping cell- versus node-centered data in finite difference codes

Help
1

Howdy, this question is not Fortran-specific, but it may have good Fortran-specific solutions. In finite difference schemes it is reasonably common to have quantities tabulated on offset meshes. For instance in 1D, it is common to consider both a grid at “cell centers” and also one with points on “cell edges”. When implementing difference schemes, it is important to keep straight which quantities are defined on which grid, especially when it comes to handling boundary conditions. Doing this bookkeeping mentally gets tedious and error-prone very quickly (especially as one goes to higher dimensions). I’m wondering if there are tried-and-true tricks or patterns for abstracting some of these bookkeeping details away.

2 Likes
2

One trick I have used in the past is to use custom lower and upper bounds. The cell-centered quantities use the default indexing from 1, while quantities at the edges are numbered from 0:

integer, parameter :: n = 100
real :: conc(n)    ! cell centers
real :: flux(0:n)  ! edges

Another trick more typical of node-based methods is to extend your arrays with ghost layers. In this case you could use a numbering like:

real :: conc(0:n+1)  ! node-centered

where nodes 0 and n+1 are the ghost nodes. You can use these efficiently to get periodic or symmetry boundaries, while preserving a “central” stencil, i.e.

integer, parameter :: n = 21 ! number of nodes
real :: conc(0:n+1) ! n+1 remains unused in this example
real, parameter :: tolerance = 1.e-6, th  = 1.0
real :: dx, factor, tmp, sumerr
integer :: i
dx = 1.0/real(n-1)
factor = 1.0/((th*dx)**2 + 2)
conc(n) = 1.0       ! fixed concentration at right-most node

outer_loop: do
  conc(0) = conc(2) ! leads to symmetry at conc(1)
  sumerr = 0
  do i = 1, n-1
     tmp = factor*(conc(i+1) + conc(i-1))
     sumerr = sumerr + abs(conc(i) - tmp)
     conc(i) = tmp
  end do
  if (sumerr < tolerance) exit outer_loop
end do outer_loop

If your are writing a parallel code, the ghost layers will be used to perform halo exchanges. If your stencils require more than one level of neighbor access, you will need more ghost nodes.

This is not really related to finite difference methods, but for good performance you might also want to use blocking (meaning you declare your arrays to some multiple of the register size) in order to improve cache utilization.

6 Likes