Document solving a stiff PDE on the GPU with a sparse Jacobian

[AJ0070]

Adds a page on solving a stiff PDE on the GPU with a sparse Jacobian, using the Brusselator from #1566 as a complete worked example.

Three separate things have to line up for this to work, and all three report the same Scalar indexing is disallowed, which makes them hard to tell apart. The page covers each and ends with a short guide to reading the stack trace to decide which one you have hit.

The one that is easy to miss is that the right-hand side has to execute on the device. In the #1566 script every array is Ref-wrapped, so the index set is the only broadcastable argument and it alone decides where the loop runs. With a host CartesianIndices the whole RHS runs on the CPU and reads the device arrays element by element, before autodiff or the linear solver are involved. The page also notes that mul! against a CuSparseMatrixCSR has no ForwardDiff.Dual kernel, which is a second way to write an RHS that cannot be differentiated on the device.

The example was run on a T4: Rosenbrock23, Rosenbrock23 with KrylovJL_GMRES, and FBDF all match a CPU reference to about 1e-10. It needs #4452 for the Rosenbrock23 results shown, so this is best merged after that one.

AI Disclosure: Used Opus 5