12. On the GPU

A large grid asks for millions of distances per regrid. The last chapter runs the batch of chapter 5 on the device with OpenACC and checks it against the host batch.

use fossil, only : surface_stl_object, SIGN_PSEUDO_NORMAL, SIGN_FAST_WINDING_NUMBER, STATUS_OK
use penf,   only : I4P, R8P
use vecfor, only : vector_R8P, ex_R8P, ey_R8P, ez_R8P

implicit none
integer(I4P), parameter       :: n = 24                    ! points per direction
type(surface_stl_object)      :: surface
type(vector_R8P), allocatable :: points(:)
real(R8P),        allocatable :: d_host(:), d_device(:)
character(14),   parameter    :: names(2) = ['pseudo-normal ', 'winding number']
integer(I4P)                  :: i, j, k, p, s, algorithm(2), status

call surface%load_from_file(file_name='dragon-placed.stl', guess_format=.true.)  ! SAH BVH, the device tree
allocate(points(n**3), d_host(n**3), d_device(n**3))
p = 0
do k = 1, n ; do j = 1, n ; do i = 1, n
   p = p + 1
   points(p) = (-0.6_R8P + 1.2_R8P * (i - 1) / (n - 1)) * ex_R8P + (-0.6_R8P + 1.2_R8P * (j - 1) / (n - 1)) * ey_R8P + &
               (-0.6_R8P + 1.2_R8P * (k - 1) / (n - 1)) * ez_R8P
enddo ; enddo ; enddo

call surface%enter_device                                   ! copy the BVH arrays to the device, once
algorithm = [SIGN_PSEUDO_NORMAL, SIGN_FAST_WINDING_NUMBER]
do s = 1, 2
   call surface%distance_many_device(points=points, distances=d_device, is_signed=.true., is_square_root=.true., &
                                     sign_algorithm=algorithm(s), status=status)
   if (status /= STATUS_OK) error stop 'distance_many_device failed'
   call surface%distance_many(points=points, distances=d_host, is_signed=.true., is_square_root=.true., &
                              sign_algorithm=algorithm(s))
   print '(A,A,I0,A,I0)', names(s), ': device vs host mismatches ', &
                             count(d_device /= d_host), ' of ', size(points)
enddo
call surface%exit_device                                    ! release the device copies

• enter_device copies the flat arrays the device kernels read (the BVH nodes, the triangle payload, the pseudo-normals, the winding-number cache) to the device, once. The surface must be analyzed with the SAH BVH, the default: on an octree enter_device refuses with a message on stderr and the surface stays on the host. A second call is a no-op.
• distance_many_device has the arguments and the defaults of distance_many (unsigned and squared unless asked). Each call moves only the points in and the distances out. Two signs run on the device, SIGN_PSEUDO_NORMAL and SIGN_FAST_WINDING_NUMBER; any other sign algorithm, or a surface not on the device, is refused with a message on stderr.
• exit_device releases the device copies. destroy does not: call it before the surface is destroyed or goes out of scope, and before analyze or anything else that rebuilds the tree. move (chapter 11) refreshes the device copy itself.
• query_many_device, box_classify_many_device, wall_distance_many_device and the curve batches follow the same pattern.

Running it

$ tut12_gpu
pseudo-normal : device vs host mismatches 0 of 13824
winding number: device vs host mismatches 0 of 13824

This output comes from a gfortran build. Without an OpenACC compiler the !$acc directives are comments, so the "device" batch runs on the host, through the same kernels, and matches the host batch bit for bit: 0 mismatches of 13824 for both signs. The program is still useful there, as a check that the device code path is correct.

The figure is drawn from that run: the device code path, executed on the host, not a GPU.

On a GPU, build FOSSIL with nvfortran 26.9 or later (NVIDIA HPC SDK). The static-nvidia mode builds the library, tests-nvidia the distance benchmark; the local_nvf varset sets the compute capability (cc89 by default, edit $NVF_CC for another GPU, e.g. cc80 for an A100):

fobis build --mode static-nvidia --varset local_nvf
fobis build --mode tests-nvidia --varset local_nvf
nvfortran -acc -gpu=cc89,nofma -Mnofma -I static/mod docs/examples/src/tut12_gpu.f90 static/libfossil.a -o tut12_gpu

The flags -gpu=nofma -Mnofma turn off fused multiply-add contraction on the device and on the host. FMA contraction differs between the two, and with it the device results are no longer bit-identical to the host ones.

What you learned

enter_device once, distance_many_device per batch, exit_device at the end; built with gfortran the same calls run on the host with identical results. Reference: GPU distance queries.

This is the end of the tutorial. The guide and the reference cover every method and argument.