concept finite difference in category fortran

This is an excerpt from Manning's book Modern Fortran: Building efficient parallel applications MEAP V13. Login to get full access to this book.

If you want to delve deeper into the math behind this problem, head over to appendix B. There, I explain the gradient, which is the key concept behind advection, and how to express it in computer code using finite differences. This step is important, as it forms the foundation to express all other terms in the shallow water equations. Otherwise, if you want to skip the math and jump straight to programming, carry on!

Listing 8.9. One-dimensional version of the finite difference function

language="fortran" linenumbering="unnumbered">pure function diff(x) result(dx)
  real(real32), intent(in) :: x(:) #1
  real(real32) :: dx(size(x))      #2
  integer(int32) :: i, im
  im = size(x)
  dx = 0
  do concurrent(i = 2:im-1)         #3
    dx(i) = 0.5 * (x(i+1) - x(i-1)) #4
  end do
end function diff

To extend it to two dimensions, we need to adapt the declaration of the input array x and the result dx to both be two-dimensional arrays. Furthermore, for a finite difference in the x direction, we’ll calculate the difference between i+1 and i-1 elements for each i, applied to the whole array slice in the y direction, as shown in the following listing.

Listing 8.10. Two-dimensional version of the diff function, for differencing along the x direction

language="fortran" linenumbering="unnumbered">pure function diffx(x) result(dx)
  real(real32), intent(in) :: x(:,:) #1
  real(real32) :: dx(size(x, dim=1), size(x, dim=2)) #2
  integer(int32) :: i, im
  im = size(x, dim=1) #3
  dx = 0
  dx(2:im-1,:) = 0.5 * (x(3:im,:) - x(1:im-2,:)) #4
end function diffx

In listing 8.10, diffx expects a two-dimensional real array x as the only input argument. This is an assumed-shape array (x(:,:)), so the function will accept a two-dimensional real array of any size. We then use the size built-in function to declare the result dx to have exactly the same size as the input array x. To calculate the difference in x, we loop over all inner elements in the first dimension and apply a whole-array slice in the second dimension. Like before, the result of each element doesn’t depend on any other element, so we can evaluate dx in any order by applying a do concurrent loop. Note that now in the two-dimensional case, diffx only returns the finite difference along the array rows. We’ll also need a diffy function to compute the finite difference along the array columns. I leave this to you as an exercise (“Exercise 3” sidebar).

Exercise 3: Computing finite difference in the y direction

We just implemented a two-dimensional variant of the original diff function that we used in the one-dimensional solver, but only for calculating the difference in x direction (along rows). Using diffx from listing 8.10 as a template, can you implement the function diffy that will return the finite difference in y direction (along columns)? Is either one of these two functions likely to be more efficient than the other, and why?

Figure  B.3. Approximating partial derivatives in the advection equation with finite differences. The top equation approximates the change of u in time (tendency), and the bottom equation approximates the change of u in space (gradient). All terms on the far right side can be represented with variables in a computer program. Discrete time and space indices are shown as n and i, respectively.