In addition to an extra equation and a finite difference term for each of the x and y axes, we need specific implementations of diff for each of them. diffx will return a difference along array rows, and diffy along array columns.
The function in the following listing served us well in the one-dimensional solver.
Listing 8.9 One-dimensional version of the finite difference function
pure function diff(x) result(dx)
real(real32), intent(in) :: x(:) ❶
real(real32) :: dx(size(x)) ❷
integer(int32) :: i, im
im = size(x)
dx = 0
do concurrent(i = 2:im-1) ❸
dx(i) = 0.5 * (x(i+1) - x(i-1)) ❹
end do
end function diff❸ Can be evaluated in any order
❹ Difference between neighboring cells
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 diff function for differencing in the x direction
pure function diffx(x) result(dx)
real(real32), intent(in) :: x(:,:) ❶
real(real32) :: dx(size(x, dim=1), size(x, dim=2)) ❷
integer(int32) :: i, im
im = size(x, dim=1) ❸
dx = 0
dx(2:im-1,:) = 0.5 * (x(3:im,:) - x(1:im-2,:)) ❹
end function diffx❶ The input array is now two-dimensional.
❸ Length of the first dimension
❹ Finite difference with whole-array slices along the y direction
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 apply a whole-array slice in both dimensions. 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 whole-array slice. Note that 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?
You can find the solution in the “Answer key” section near the end of this chapter.
If you follow along with the code from GitHub, you can find these functions in tsunami/src/ch08/mod_diff.f 90.
However, our job here is not done yet. We still need to be able to pass the Field instance to diffx and diffy, which expect real arrays as input. In the next subsection, I describe how you can make a wrapper function that will accept a Field instance as an input argument and perform a finite difference calculation on its data component.