Vectorized Bruss

 kernel_u! = let N=N, xyd=xyd_brusselator, dx=step(xyd_brusselator)
    @inline function (du, u, A, B, α, II, I, t)
      i, j = Tuple(I)
      x = xyd[I[1]]
      y = xyd[I[2]]
      ip1 = limit(i+1, N); im1 = limit(i-1, N)
      jp1 = limit(j+1, N); jm1 = limit(j-1, N)
      du[II[i,j,1]] = α[II[i,j,1]]*(u[II[im1,j,1]] + u[II[ip1,j,1]] + u[II[i,jp1,1]] + u[II[i,jm1,1]] - 4u[II[i,j,1]])/dx^2 +
      B[II[i,j,1]] + u[II[i,j,1]]^2*u[II[i,j,2]] - (A[II[i,j,1]] + 1)*u[II[i,j,1]] + brusselator_f(x, y, t)
    end
  end
  kernel_v! = let N=N, xyd=xyd_brusselator, dx=step(xyd_brusselator)
    @inline function (du, u, A, B, α, II, I, t)
      i, j = Tuple(I)
      ip1 = limit(i+1, N)
      im1 = limit(i-1, N)
      jp1 = limit(j+1, N)
      jm1 = limit(j-1, N)
      du[II[i,j,2]] = α[II[i,j,2]]*(u[II[im1,j,2]] + u[II[ip1,j,2]] + u[II[i,jp1,2]] + u[II[i,jm1,2]] - 4u[II[i,j,2]])/dx^2 +
      A[II[i,j,1]]*u[II[i,j,1]] - u[II[i,j,1]]^2*u[II[i,j,2]]
    end
  end
  brusselator_2d = let N=N, xyd=xyd_brusselator, dx=step(xyd_brusselator)
    function (du, u, p, t)
      @inbounds begin
        ii1 = N^2
        ii2 = ii1+N^2
        ii3 = ii2+2(N^2)
        A = @view p[1:ii1]
        B = @view p[ii1+1:ii2]
        α = @view p[ii2+1:ii3]
        II = LinearIndices((N, N, 2))
        kernel_u!.(Ref(du), Ref(u), Ref(A), Ref(B), Ref(α), Ref(II), CartesianIndices((N, N)), t)
        kernel_v!.(Ref(du), Ref(u), Ref(A), Ref(B), Ref(α), Ref(II), CartesianIndices((N, N)), t)
        return nothing
      end
    end
  end