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.
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
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.