Updating the finite difference calculation

As you probably noticed in the previous subsection, in the new version of the tsunami simulator we’ll get to apply the diff function not once, but three times. It looks like our work with procedures and modules will finally pay off. There’s just one more thing we need to do, and that is make our finite difference function a tad more general. Recall that our original finite difference function calculated the upwind difference between the current and previous element:

dx(2:im) = x(2:im) - x(1:im-1)      ❶

❶ Whole-array difference

It turns out that this difference pattern works well only if the motion in the system is moving from left to right–that is, from lower to higher array indices. This is also why it worked well for our simple advection simulator, where the constant background flow was configured from left to right.

However, in a more realistic configuration, where water velocity is coupled with water height and water is allowed to move in both directions, we need a centered finite difference that can account for changes coming from either direction:

dx(2:im-1) = 0.5 * (x(3:im) - x(1:im-2))

In contrast to the original finite difference, here we’re calculating the difference between the next and previous element and dividing by 2 to account for the fact that we’re calculating the difference over two increments. Figure 4.7 illustrates the difference in calculating these differences.

The following listing provides the complete finite difference function, centered in space.

Listing 4.9 Finite difference of a one-dimensional array, centered in space

pure function diff_centered(x) result(dx)
  real(real32), intent(in) :: x(:)
  real(real32) :: dx(size(x))
  integer(int32) :: im
  im = size(x)
  dx(1) = x(2) - x(im)                  ❶
  dx(im) = x(1) - x(im-1)               ❷
  dx(2:im-1) = x(3:im) - x(1:im-2)      ❸
  dx = 0.5 * dx                         ❹
end function diff_centered

❶ Calculates the boundary value on the left

❷ Calculates the boundary value on the right

❸ Calculates the difference in the interior

❹ Divides all elements by 2

Figure 4.7 Illustration of upwind and centered finite difference calculations

We now have two different versions of the finite difference function: one for the centered difference in space, and another for the upstream difference in space. They both receive the same input, a one-dimensional array x, and they each return a result that has the same meaning–it’s just calculated differently. However, it’s important to use the centered difference function now to allow the water to move in both directions. Notice that in listing 4.9, I multiply dx by 0.5 rather than dividing it by 2. In theory, this is preferred, because for floating-point (real) numbers, multiplication is a much faster operation than division. In practice, however, a good compiler is typically able to optimize floating-point division for you, so this becomes a matter of style.