  Numerical methods such as the Finite Element Method (FEM) have been
successfully adapted to utilize the computational power of GPU accelerators.
However, much of the effort around applying FEM to GPU's has been focused on
high-order FEM due to higher arithmetic intensity and order of accuracy. For
applications such as the simulation of subsurface processes, high levels of
heterogeneity results in high-resolution grids characterized by highly
discontinuous (cell-wise) material property fields. Moreover, due to the
significant uncertainties in the characterization of the domain of interest,
e.g. geologic reservoirs, the benefits of high order accuracy are reduced, and
low-order methods are typically employed. In this study, we present a strategy
for implementing highly performant low-order matrix-free FEM operator kernels
in the context of the conjugate gradient (CG) method. Performance results of
matrix-free Laplace and isotropic elasticity operator kernels are presented and
are shown to compare favorably to matrix-based SpMV operators on V100, A100,
and MI250X GPUs.
