Combine PETSc.jl with KernelAbstractions/CUDA.jl

[boriskaus]

PETSc has GPU support since a few years, which allows putting matrices and vectors to the GPU using command-line options (provided that PETSc is configured with CUDA support). Yet, to get the best performance, also the residual kernels will have to be rewritten for GPU, which requires writing native CUDA code (for NVIDIA systems), or use Kokkos ("dependency hell" as someone told me).

Julia, on the other hand, is well-known for its excellent GPU support, for example, through packages such as KernelAbstractions, which allows compiling the same vcode for NVIDIA, AMD, or Mac hardware in a straightforward manner.

Here, the PETSc ex19 example (typically used to test PETSc installations) was translated from C to Julia and combined with KernelAbstractions such that residual routines also run on the GPU. It uses coloring and finite differences to approximate the Jacobian and works with multigrid preconditioners.

The documentation has been updated to show scalability tests on GPU vs. CPU (on 1 or 32 cores).

As a quick summary, it works and already shows great potential but can perhaps be further improved:

KSPSolve:

Resolution	GPU time (s)	GPU (GFlop/s)	CPU-1 time (s)	CPU-1 (GFlop/s)	CPU-32 time (s)	CPU-32 (GFlop/s)
513²	0.116	144.3	4.337	4.1	0.297	61.0
1025²	0.299	249.5	19.10	3.9	1.196	61.7
2049²	1.118	295.5	89.76	3.6	5.757	56.6
4097²	4.540	312.4	422.2	3.3	28.30	49.4

SNESSolve:

Resolution	GPU time (s)	GPU (GFlop/s)	CPU-1 time (s)	CPU-1 (GFlop/s)	CPU-32 time (s)	CPU-32 (GFlop/s)
513²	3.645	6.7	8.041	3.2	2.118	12.3
1025²	5.127	20.5	32.75	3.2	3.698	28.2
2049²	11.37	39.7	144.1	3.1	10.57	42.3
4097²	36.88	47.7	658.3	2.9	43.58	43.2

---

Disclaimer: Claude Sonnet 4.6 was used in preparing this; Valentin helped steer it back in place...

[vchuravy]

There should also be tests for withlocalarray(Array, ...) and co