Make ArrayPartitions of GPU arrays work in DifferentialEquations solvers - #416

Checklist

• Appropriate tests were added
• [Tests Pass] Any code changes were done in a way that does not break public API
• [NA] All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• [NA] Any new documentation only uses public API

Additional context

The only function that was preventing ArrayPartitions of GPU arrays to work in DifferentialEquations solvers was recursivefill!, which was deferring to Base.fill! with scalar indexing. This PR makes fixes two things:

• recursivefill! calls fill! on each array within the ArrayPartitions.x tuple - which now defer to the right GPU kernels.
• Filling does not require index recalculations for each setindex! - which was unnecessary work.

Example Code (similar to what was added to the tests):

using RecursiveArrayTools
using Metal

a = (Metal.rand(5), Metal.rand(5))
pA = ArrayPartition(a)

# This now works
RecursiveArrayTools.recursivefill!(pA, false)

Possible Issues

• I added an unrolled_foreach(f, x::Tuple) implementation for type stability when calling fill! on each subarray in the ArrayPartition.x tuple. This now uses a closure for f which updates each subarray inline; it is type-stable when checking with Cthulhu, and it does not update any variable outside the closure; I see the SciML guidelines do not recommend closures - should I rewrite that as a struct and all that?
• I added a specialised implementation of recursivefill! for ArrayPartitions; all tests seem to work, but I don't have the "bigger picture" of how this is used across the ecosystem - is this okay / non-ambiguous? Should the other recursivefill! implementations in utils.jl be changed?