16. On the GPU

A large grid asks for millions of distances per regrid. This chapter runs the queries of the earlier chapters on the device with OpenACC, one batch per call, and checks every result against the host routine it replaces: the distance batch of chapter 5, the immersed-boundary queries of chapter 8, the box classification of chapter 9 and the 2-D section of chapter 7.

The pages are built with gfortran. Without an OpenACC compiler the !$acc directives are comments, so every "device" batch below runs on the host, through the same kernels the GPU runs. The runs therefore show the API and the bit-exactness check, never a speed-up; the same programs built with nvfortran run on the GPU (last section).

One batch of distances

The surface goes to the device once, then each batch moves only its points in and its distances out.

use fossil, only : surface_stl_object, AABB_TREE_OCTREE, SIGN_PSEUDO_NORMAL, SIGN_FAST_WINDING_NUMBER, &
                   SIGN_RAY_INTERSECTIONS, 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, octree
type(vector_R8P), allocatable :: points(:)
real(R8P),        allocatable :: d_pn(:), d_fwn(:), d_host(:)
integer(I4P)                  :: i, j, k, p, s, status

call surface%load_from_file(file_name='dragon-placed.stl', guess_format=.true.)  ! SAH BVH: the device tree
call surface%enter_device                                   ! copy the flat BVH arrays to the device, once
call surface%enter_device                                   ! a second call is a no-op
print '(A,L1)', 'on the device: ', surface%is_on_device

• enter_device copies the flat arrays the device kernels read: the BVH nodes, the triangle payload, the pseudo-normals, the winding-number cache, the bodies and the exact facet vertices (plus the feature-edge and thin trees of chapter 12 if they were built). A second call is a no-op; is_on_device tells whether the copy exists.
• The surface must be analyzed with the SAH BVH, the default tree of load_from_file.

allocate(points(n**3), d_pn(n**3), d_fwn(n**3), d_host(n**3))
p = 0
do k = 1, n ; do j = 1, n ; do i = 1, n                     ! x fastest, z slowest: slice k is contiguous
   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%distance_many_device(points=points, distances=d_pn, is_signed=.true., is_square_root=.true., &
                                  sign_algorithm=SIGN_PSEUDO_NORMAL, status=status)
call surface%distance_many(points=points, distances=d_host, is_signed=.true., is_square_root=.true., &
                           sign_algorithm=SIGN_PSEUDO_NORMAL)
print '(A,I0,A,I0)', 'pseudo-normal : device vs host mismatches ', count(d_pn /= d_host), ' of ', size(points)
call surface%distance_many_device(points=points, distances=d_fwn, is_signed=.true., is_square_root=.true., &
                                  sign_algorithm=SIGN_FAST_WINDING_NUMBER, status=status)
call surface%distance_many(points=points, distances=d_host, is_signed=.true., is_square_root=.true., &
                           sign_algorithm=SIGN_FAST_WINDING_NUMBER)
print '(A,I0,A,I0)', 'winding number: device vs host mismatches ', count(d_fwn /= d_host), ' of ', size(points)
print '(A,I0,A,I0)', 'the two signs disagree at ', count(sign(1._R8P, d_pn) /= sign(1._R8P, d_fwn)), &
                     ' points; |d| differs at ', count(abs(d_pn) /= abs(d_fwn))

• distance_many_device has the arguments and the defaults of distance_many: unsigned and squared unless is_signed=.true. and is_square_root=.true. are passed. The 24³ points are a lattice over [-0.6, 0.6]³, stored x fastest, so the plane k = const is a contiguous slice of the array.
• Each sign algorithm runs once on the device and once on the host, and the program counts the distances that differ in any bit.

The figure is one plane of that batch: the distance field around the dragon's section, the device code path executed on the host.

One slice at a time, or one batch

A solver that fills its grid plane by plane can call the host routine once per plane; on the device the whole grid goes in one call, to keep the GPU busy and pay the transfer latency once. The results must not depend on the split.

do k = 1, n                                                 ! the same points, one host call per slice z = const
   s = (k - 1) * n**2
   call surface%distance_many(points=points(s + 1:s + n**2), distances=d_host(s + 1:s + n**2), is_signed=.true., &
                              is_square_root=.true., sign_algorithm=SIGN_FAST_WINDING_NUMBER)
enddo
print '(A,I0,A,I0)', '24 host calls of one slice vs one device batch: mismatches ', count(d_host /= d_fwn), &
                     ' of ', size(points)

• Each point is an independent query: its traversal of the tree and its sign test do not read the other points of the batch, so splitting the batch into slices cannot change a bit. The run checks it: 24 host calls of one slice each against one device call of 13824 points, 0 mismatches.

Two signs, and what the device refuses

The device runs two of the four sign algorithms of chapter 4, and only on the BVH. Anything else is refused, loudly, instead of silently falling back to the host.

call surface%distance_many_device(points=points, distances=d_host, is_signed=.true., &
                                  sign_algorithm=SIGN_RAY_INTERSECTIONS, status=status)  ! not ported: refused
print '(A,L1)', 'ray-intersection sign on the device refused: ', status /= STATUS_OK
call octree%load_from_file(file_name='dragon-placed.stl', guess_format=.true., aabb_tree_kind=AABB_TREE_OCTREE)
call octree%enter_device                                    ! the device traverses the SAH BVH only: refused
print '(A,L1)', 'octree on the device: ', octree%is_on_device
call surface%exit_device                                    ! release the device copies, before destroy or analyze
print '(A,L1)', 'after exit_device: ', surface%is_on_device

$ tut16_a_batch
on the device: T
pseudo-normal : device vs host mismatches 0 of 13824
winding number: device vs host mismatches 0 of 13824
the two signs disagree at 180 points; |d| differs at 0
24 host calls of one slice vs one device batch: mismatches 0 of 13824
fossil_surface_stl_object%distance_many_device: refused -- sign_algorithm must be SIGN_PSEUDO_NORMAL or SIGN_FAST_WINDING_NUMBER on the device.
ray-intersection sign on the device refused: T
fossil_surface_stl_object%enter_device: refused -- device path requires AABB_TREE_SAH_BVH.
octree on the device: F
after exit_device: F

• SIGN_PSEUDO_NORMAL (the default) and SIGN_FAST_WINDING_NUMBER run on the device, both bit-identical to the host: 0 mismatches of 13824 each. The FWN branch calls the same bvh_winding_number_oac routine with the same constants as the host path, which is what keeps it bit-exact.
• SIGN_RAY_INTERSECTIONS and SIGN_SOLID_ANGLE need a ray traversal or an O(N) sum that is not ported: the call returns status /= STATUS_OK and writes the reason on stderr (the first refused line of the run).
• An octree (aabb_tree_kind=AABB_TREE_OCTREE) stays on the host: enter_device writes a message and is_on_device stays false. The device scope is the SAH BVH only.
• 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.
• The two signs agree on |d| everywhere and disagree on the sign at 180 of the 13824 points. The repaired dragon is not closed (chapter 2). The pseudo-normal decides from the closest feature alone and assumes a closed, outward oriented surface; the winding number integrates the whole surface and tolerates holes (chapter 4). The figure shows where they differ.

For a surface that is not closed, pass sign_algorithm=SIGN_FAST_WINDING_NUMBER on the device as on the host.

Rich queries on the device

The immersed-boundary queries of chapter 8 have device twins with the same arguments.

use fossil,     only : surface_stl_object, rigid_pose_t, pose_from_axis_angle, FEATURE_FACE, FEATURE_EDGE, &
                       FEATURE_VERTEX, SIGN_FAST_WINDING_NUMBER
use penf,       only : I4P, R8P
use vecfor,     only : vector_R8P, ex_R8P, ey_R8P, ez_R8P

implicit none
integer(I4P), parameter       :: n = 12                    ! probes per direction
real(R8P),    parameter       :: pi = 4._R8P * atan(1._R8P)
type(surface_stl_object)      :: surface
type(rigid_pose_t)            :: pose
type(vector_R8P), allocatable :: points(:)
real(R8P),        allocatable :: d(:, :), cp(:, :, :), nm(:, :, :)
integer(I4P),     allocatable :: facet(:, :), feature(:, :), body(:, :)
integer(I4P)                  :: i, j, k, p

call surface%load_from_file(file_name='dragon-placed.stl', guess_format=.true.)
call surface%enter_device
allocate(points(n**3), d(n**3, 2), cp(3, n**3, 2), nm(3, n**3, 2), facet(n**3, 2), feature(n**3, 2), body(n**3, 2))
p = 0
do k = 1, n ; do j = 1, n ; do i = 1, n
   p = p + 1
   points(p) = (-0.55_R8P + 1.1_R8P * (i - 1) / (n - 1)) * ex_R8P + &
               (-0.55_R8P + 1.1_R8P * (j - 1) / (n - 1)) * ey_R8P + (-0.55_R8P + 1.1_R8P * (k - 1) / (n - 1)) * ez_R8P
enddo ; enddo ; enddo
call surface%query_many_device(points=points, distances=d(:, 1), closest_points=cp(:, :, 1), facets=facet(:, 1), &
                               features=feature(:, 1), normals=nm(:, :, 1), bodies=body(:, 1), &
                               sign_algorithm=SIGN_FAST_WINDING_NUMBER)                  ! column 1: device
call surface%query_many(points=points, distances=d(:, 2), closest_points=cp(:, :, 2), facets=facet(:, 2), &
                        features=feature(:, 2), normals=nm(:, :, 2), bodies=body(:, 2), &
                        sign_algorithm=SIGN_FAST_WINDING_NUMBER)                         ! column 2: host
print '(A,I0,A)', 'query_many_device vs query_many, ', size(points), ' probes, mismatches:'
print '(A,I0,A,I0,A,I0)', '  distances ', count(d(:, 1) /= d(:, 2)), ', closest points ', &
      count(any(cp(:, :, 1) /= cp(:, :, 2), dim=1)), ', facets ', count(facet(:, 1) /= facet(:, 2))
print '(A,I0,A,I0,A,I0)', '  features ', count(feature(:, 1) /= feature(:, 2)), ', normals ', &
      count(any(nm(:, :, 1) /= nm(:, :, 2), dim=1)), ', bodies ', count(body(:, 1) /= body(:, 2))
print '(A,I0,A,I0,A,I0)', 'feet on a face ', count(feature(:, 1) == FEATURE_FACE), ', on an edge ', &
      count(feature(:, 1) == FEATURE_EDGE), ', on a vertex ', count(feature(:, 1) == FEATURE_VERTEX)

• query_many_device returns what query_many returns: distance, closest point, closest facet, the feature the foot lies on (FEATURE_FACE, FEATURE_EDGE, FEATURE_VERTEX), the pseudo-normal and the body. The kernel computes every output and copies out the requested ones. The CPU-only extras (ambiguous, second_foot_points, curvatures) have no device twin.

call surface%wall_distance_many_device(points=points, distances=d(:, 1), bodies=body(:, 1), cutoff=0.05_R8P)
call surface%wall_distance_many(points=points, distances=d(:, 2), bodies=body(:, 2), cutoff=0.05_R8P)
print '(A,I0,A,I0,A,I0)', 'wall distance, cutoff 0.05: mismatches ', count(d(:, 1) /= d(:, 2)), &
      ', probes nearer than the cutoff ', count(body(:, 1) /= 0), ' of ', size(points)

• wall_distance_many_device(points, distances, bodies, cutoff=) is the unsigned wall distance of chapter 8 with an early cutoff: the traversal starts at cutoff², so a probe with no wall closer gets cutoff and body 0.

pose = pose_from_axis_angle(axis=[0._R8P, 1._R8P, 0._R8P], angle=pi / 18._R8P, translation=[0.02_R8P, 0._R8P, 0._R8P])
call surface%query_many_device(points=points, distances=d(:, 1), facets=facet(:, 1), pose=pose)  ! chapter 11's step
call surface%query_many(points=points, distances=d(:, 2), facets=facet(:, 2), pose=pose)
print '(A,I0,A,I0)', 'posed query: distance mismatches ', count(d(:, 1) /= d(:, 2)), ', facet mismatches ', &
      count(facet(:, 1) /= facet(:, 2))
call surface%exit_device

• pose= (chapter 11) maps the points into the body frame on the host and queries the tree as entered: a moving body needs no new enter_device per step.

$ tut16_b_queries
query_many_device vs query_many, 1728 probes, mismatches:
  distances 0, closest points 0, facets 0
  features 0, normals 0, bodies 0
feet on a face 128, on an edge 477, on a vertex 1123
wall distance, cutoff 0.05: mismatches 0, probes nearer than the cutoff 153 of 1728
posed query: distance mismatches 0, facet mismatches 0

Every output matches the host bit for bit on the 12³ probes. Most feet fall on a vertex (1123) or an edge (477), not inside a face (128): seen from probes spread over the box, the closest point of a coarse mesh is most often one of its corners or edges.

Boxes and curves

AMR tagging (chapter 9) classifies boxes, and the 2-D section (chapter 7) has its own segment tree: both have device batches.

use fossil,     only : surface_stl_object, curve_object, BOX_CUT, BOX_NEAR, BOX_INSIDE, BOX_OUTSIDE
use penf,       only : I4P, R8P
use vecfor,     only : vector_R8P, ex_R8P, ey_R8P, ez_R8P

implicit none
integer(I4P), parameter   :: nb = 16, nc = 32, np = 64     ! 3-D boxes, 2-D boxes, 2-D points per direction
type(surface_stl_object)  :: surface
type(curve_object)        :: section
real(R8P)                 :: lo(3, nb**3), hi(3, nb**3), lo2(2, nc**2), hi2(2, nc**2), pt(2, np**2), d(np**2, 2), h
integer(I4P)              :: codes(nb**3, 2), codes2(nc**2, 2), i, j, k, b

call surface%load_from_file(file_name='dragon-placed.stl', guess_format=.true.)
call surface%enter_device
h = 1.2_R8P / nb
b = 0
do k = 1, nb ; do j = 1, nb ; do i = 1, nb
   b = b + 1
   lo(:, b) = -0.6_R8P + [i - 1, j - 1, k - 1] * h
   hi(:, b) = lo(:, b) + h
enddo ; enddo ; enddo
call surface%box_classify_many_device(lo=lo, hi=hi, codes=codes(:, 1), radius=0.05_R8P)  ! winding-number sign
call surface%box_classify_many(lo=lo, hi=hi, codes=codes(:, 2), radius=0.05_R8P)
print '(A,I0,A,I0)', 'surface boxes: device vs host mismatches ', count(codes(:, 1) /= codes(:, 2)), ' of ', nb**3
print '(A,I0,A,I0,A,I0,A,I0)', '  cut ', count(codes(:, 1) == BOX_CUT), ', near ', count(codes(:, 1) == BOX_NEAR), &
      ', inside ', count(codes(:, 1) == BOX_INSIDE), ', outside ', count(codes(:, 1) == BOX_OUTSIDE)
call surface%exit_device

• box_classify_many_device(lo=, hi=, codes=, radius=) gives the codes of box_classify_many: BOX_CUT, BOX_NEAR, BOX_INSIDE, BOX_OUTSIDE. The contact test is the same bvh_box_contact_oac kernel on both sides; the centre sign of a box not cut uses the winding number by default.
• On 16³ boxes of side 0.075 with a near radius 0.05, every box that reaches into the dragon is cut or near, so none is inside.

call section%load_from_file(file_name='dragon-section.dat')         ! chapter 7: the section z = 0
call section%enter_device                                           ! the segment BVH and its caches
h = 1.2_R8P / nc
b = 0
do j = 1, nc ; do i = 1, nc
   b = b + 1
   lo2(:, b) = -0.6_R8P + [i - 1, j - 1] * h
   hi2(:, b) = lo2(:, b) + h
enddo ; enddo
call section%box_classify_many_device(lo=lo2, hi=hi2, codes=codes2(:, 1), radius=0.02_R8P)
call section%box_classify_many(lo=lo2, hi=hi2, codes=codes2(:, 2), radius=0.02_R8P)
print '(A,I0,A,I0)', 'curve boxes: device vs host mismatches ', count(codes2(:, 1) /= codes2(:, 2)), ' of ', nc**2
b = 0
do j = 1, np ; do i = 1, np
   b = b + 1
   pt(:, b) = -0.6_R8P + 1.2_R8P * [i - 1, j - 1] / (np - 1._R8P)
enddo ; enddo
call section%distance_many_device(points=pt, distances=d(:, 1))     ! signed, negative inside
call section%distance_many(points=pt, distances=d(:, 2))
print '(A,I0,A,I0,A,I0)', 'curve distances: device vs host mismatches ', count(d(:, 1) /= d(:, 2)), ' of ', np**2, &
      ', inside ', count(d(:, 1) < 0._R8P)
call section%exit_device

• A curve_object has its own enter_device, exit_device, distance_many_device (points as (2, n) arrays, signed by default, negative inside), query_many_device and box_classify_many_device, with the same rules.

$ tut16_c_boxes
surface boxes: device vs host mismatches 0 of 4096
  cut 374, near 426, inside 0, outside 3296
curve boxes: device vs host mismatches 0 of 1024
curve distances: device vs host mismatches 0 of 4096, inside 205

Building for 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/tut16_a_batch.f90 static/libfossil.a -o tut16_a_batch

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 (1-ulp distances, sign and facet flips at near-ties). The figure programs of this chapter also need docs/examples/lib/fossil_fig.f90 compiled in.

What you learned

enter_device once, a *_device batch per call with the arguments of its host twin, exit_device before the tree changes or the surface goes away. Two signs (pseudo-normal, fast winding number) and the SAH BVH only; anything else is refused with a status. A batch gives the same bits as any split of it, and built with gfortran the same calls run on the host with identical results. Reference: GPU queries.

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