diff --git a/DifferentiationInterface/ext/DifferentiationInterfaceForwardDiffExt/DifferentiationInterfaceForwardDiffExt.jl b/DifferentiationInterface/ext/DifferentiationInterfaceForwardDiffExt/DifferentiationInterfaceForwardDiffExt.jl
index 9a7428ae5..36e34aac3 100644
--- a/DifferentiationInterface/ext/DifferentiationInterfaceForwardDiffExt/DifferentiationInterfaceForwardDiffExt.jl
+++ b/DifferentiationInterface/ext/DifferentiationInterfaceForwardDiffExt/DifferentiationInterfaceForwardDiffExt.jl
@@ -21,6 +21,7 @@ using DifferentiationInterface:
     outer,
     unwrap,
     with_contexts
+import ForwardDiff.DiffResults as DR
 using ForwardDiff.DiffResults:
     DiffResults, DiffResult, GradientResult, HessianResult, MutableDiffResult
 using ForwardDiff:
diff --git a/DifferentiationInterface/ext/DifferentiationInterfaceForwardDiffExt/onearg.jl b/DifferentiationInterface/ext/DifferentiationInterfaceForwardDiffExt/onearg.jl
index ee4328da4..9e3a7c6b9 100644
--- a/DifferentiationInterface/ext/DifferentiationInterfaceForwardDiffExt/onearg.jl
+++ b/DifferentiationInterface/ext/DifferentiationInterfaceForwardDiffExt/onearg.jl
@@ -165,9 +165,11 @@ function DI.value_and_gradient!(
     f::F, grad, ::AutoForwardDiff, x, contexts::Vararg{Context,C}
 ) where {F,C}
     fc = with_contexts(f, contexts...)
-    result = MutableDiffResult(zero(eltype(x)), (grad,))
+    result = DiffResult(zero(eltype(x)), (grad,))
     result = gradient!(result, fc, x)
-    return DiffResults.value(result), DiffResults.gradient(result)
+    y = DR.value(result)
+    grad === DR.gradient(result) || copyto!(grad, DR.gradient(result))
+    return y, grad
 end
 
 function DI.value_and_gradient(
@@ -176,7 +178,7 @@ function DI.value_and_gradient(
     fc = with_contexts(f, contexts...)
     result = GradientResult(x)
     result = gradient!(result, fc, x)
-    return DiffResults.value(result), DiffResults.gradient(result)
+    return DR.value(result), DR.gradient(result)
 end
 
 function DI.gradient!(
@@ -213,9 +215,11 @@ function DI.value_and_gradient!(
     contexts::Vararg{Context,C},
 ) where {F,C}
     fc = with_contexts(f, contexts...)
-    result = MutableDiffResult(zero(eltype(x)), (grad,))
+    result = DiffResult(zero(eltype(x)), (grad,))
     result = gradient!(result, fc, x, prep.config)
-    return DiffResults.value(result), DiffResults.gradient(result)
+    y = DR.value(result)
+    grad === DR.gradient(result) || copyto!(grad, DR.gradient(result))
+    return y, grad
 end
 
 function DI.value_and_gradient(
@@ -224,7 +228,7 @@ function DI.value_and_gradient(
     fc = with_contexts(f, contexts...)
     result = GradientResult(x)
     result = gradient!(result, fc, x, prep.config)
-    return DiffResults.value(result), DiffResults.gradient(result)
+    return DR.value(result), DR.gradient(result)
 end
 
 function DI.gradient!(
@@ -255,9 +259,11 @@ function DI.value_and_jacobian!(
 ) where {F,C}
     fc = with_contexts(f, contexts...)
     y = fc(x)
-    result = MutableDiffResult(y, (jac,))
+    result = DiffResult(y, (jac,))
     result = jacobian!(result, fc, x)
-    return DiffResults.value(result), DiffResults.jacobian(result)
+    y = DR.value(result)
+    jac === DR.jacobian(result) || copyto!(jac, DR.jacobian(result))
+    return y, jac
 end
 
 function DI.value_and_jacobian(
@@ -302,9 +308,11 @@ function DI.value_and_jacobian!(
 ) where {F,C}
     fc = with_contexts(f, contexts...)
     y = fc(x)
-    result = MutableDiffResult(y, (jac,))
+    result = DiffResult(y, (jac,))
     result = jacobian!(result, fc, x, prep.config)
-    return DiffResults.value(result), DiffResults.jacobian(result)
+    y = DR.value(result)
+    jac === DR.jacobian(result) || copyto!(jac, DR.jacobian(result))
+    return y, jac
 end
 
 function DI.value_and_jacobian(
@@ -457,11 +465,12 @@ function DI.value_gradient_and_hessian!(
     f::F, grad, hess, ::AutoForwardDiff, x, contexts::Vararg{Context,C}
 ) where {F,C}
     fc = with_contexts(f, contexts...)
-    result = MutableDiffResult(one(eltype(x)), (grad, hess))
+    result = DiffResult(one(eltype(x)), (grad, hess))
     result = hessian!(result, fc, x)
-    return (
-        DiffResults.value(result), DiffResults.gradient(result), DiffResults.hessian(result)
-    )
+    y = DR.value(result)
+    grad === DR.gradient(result) || copyto!(grad, DR.gradient(result))
+    hess === DR.hessian(result) || copyto!(hess, DR.hessian(result))
+    return (y, grad, hess)
 end
 
 function DI.value_gradient_and_hessian(
@@ -470,9 +479,7 @@ function DI.value_gradient_and_hessian(
     fc = with_contexts(f, contexts...)
     result = HessianResult(x)
     result = hessian!(result, fc, x)
-    return (
-        DiffResults.value(result), DiffResults.gradient(result), DiffResults.hessian(result)
-    )
+    return (DR.value(result), DR.gradient(result), DR.hessian(result))
 end
 
 ### Prepared
@@ -527,11 +534,12 @@ function DI.value_gradient_and_hessian!(
     contexts::Vararg{Context,C},
 ) where {F,C}
     fc = with_contexts(f, contexts...)
-    result = MutableDiffResult(one(eltype(x)), (grad, hess))
+    result = DiffResult(one(eltype(x)), (grad, hess))
     result = hessian!(result, fc, x, prep.manual_result_config)
-    return (
-        DiffResults.value(result), DiffResults.gradient(result), DiffResults.hessian(result)
-    )
+    y = DR.value(result)
+    grad === DR.gradient(result) || copyto!(grad, DR.gradient(result))
+    hess === DR.hessian(result) || copyto!(hess, DR.hessian(result))
+    return (y, grad, hess)
 end
 
 function DI.value_gradient_and_hessian(
@@ -540,7 +548,5 @@ function DI.value_gradient_and_hessian(
     fc = with_contexts(f, contexts...)
     result = HessianResult(x)
     result = hessian!(result, fc, x, prep.auto_result_config)
-    return (
-        DiffResults.value(result), DiffResults.gradient(result), DiffResults.hessian(result)
-    )
+    return (DR.value(result), DR.gradient(result), DR.hessian(result))
 end
diff --git a/DifferentiationInterface/ext/DifferentiationInterfaceReverseDiffExt/DifferentiationInterfaceReverseDiffExt.jl b/DifferentiationInterface/ext/DifferentiationInterfaceReverseDiffExt/DifferentiationInterfaceReverseDiffExt.jl
index bc691fd67..a6f2a1969 100644
--- a/DifferentiationInterface/ext/DifferentiationInterfaceReverseDiffExt/DifferentiationInterfaceReverseDiffExt.jl
+++ b/DifferentiationInterface/ext/DifferentiationInterfaceReverseDiffExt/DifferentiationInterfaceReverseDiffExt.jl
@@ -15,7 +15,9 @@ using DifferentiationInterface:
     NoPullbackPrep,
     unwrap,
     with_contexts
-using ReverseDiff.DiffResults: DiffResults, DiffResult, GradientResult, MutableDiffResult
+import ReverseDiff.DiffResults as DR
+using ReverseDiff.DiffResults:
+    DiffResults, DiffResult, GradientResult, HessianResult, MutableDiffResult
 using LinearAlgebra: dot, mul!
 using ReverseDiff:
     CompiledGradient,
diff --git a/DifferentiationInterface/ext/DifferentiationInterfaceReverseDiffExt/onearg.jl b/DifferentiationInterface/ext/DifferentiationInterfaceReverseDiffExt/onearg.jl
index 2cd21d07c..a316d1ca6 100644
--- a/DifferentiationInterface/ext/DifferentiationInterfaceReverseDiffExt/onearg.jl
+++ b/DifferentiationInterface/ext/DifferentiationInterfaceReverseDiffExt/onearg.jl
@@ -84,10 +84,12 @@ end
 function DI.value_and_gradient!(
     f, grad, prep::ReverseDiffGradientPrep, ::AutoReverseDiff, x
 )
-    y = f(x)  # TODO: remove once ReverseDiff#251 is fixed
-    result = MutableDiffResult(y, (grad,))
+    y = f(x)  # TODO: ReverseDiff#251
+    result = DiffResult(y, (grad,))
     result = gradient!(result, prep.tape, x)
-    return DiffResults.value(result), DiffResults.derivative(result)
+    y = DR.value(result)
+    grad === DR.gradient(result) || copyto!(grad, DR.gradient(result))
+    return y, grad
 end
 
 function DI.value_and_gradient(
@@ -115,10 +117,12 @@ function DI.value_and_gradient!(
     f, grad, ::NoGradientPrep, ::AutoReverseDiff, x, contexts::Vararg{Context,C}
 ) where {C}
     fc = with_contexts(f, contexts...)
-    y = fc(x)  # TODO: remove once ReverseDiff#251 is fixed
-    result = MutableDiffResult(y, (grad,))
+    y = fc(x)  # TODO: ReverseDiff#251
+    result = DiffResult(y, (grad,))
     result = gradient!(result, fc, x)
-    return DiffResults.value(result), DiffResults.derivative(result)
+    y = DR.value(result)
+    grad === DR.gradient(result) || copyto!(grad, DR.gradient(result))
+    return y, grad
 end
 
 function DI.value_and_gradient(
@@ -162,9 +166,11 @@ function DI.value_and_jacobian!(
     f, jac, prep::ReverseDiffOneArgJacobianPrep, ::AutoReverseDiff, x
 )
     y = f(x)
-    result = MutableDiffResult(y, (jac,))
+    result = DiffResult(y, (jac,))
     result = jacobian!(result, prep.tape, x)
-    return DiffResults.value(result), DiffResults.derivative(result)
+    y = DR.value(result)
+    jac === DR.jacobian(result) || copyto!(jac, DR.jacobian(result))
+    return y, jac
 end
 
 function DI.value_and_jacobian(f, prep::ReverseDiffOneArgJacobianPrep, ::AutoReverseDiff, x)
@@ -190,9 +196,11 @@ function DI.value_and_jacobian!(
 ) where {C}
     fc = with_contexts(f, contexts...)
     y = fc(x)
-    result = MutableDiffResult(y, (jac,))
+    result = DiffResult(y, (jac,))
     result = jacobian!(result, fc, x)
-    return DiffResults.value(result), DiffResults.derivative(result)
+    y = DR.value(result)
+    jac === DR.jacobian(result) || copyto!(jac, DR.jacobian(result))
+    return y, jac
 end
 
 function DI.value_and_jacobian(
@@ -220,46 +228,49 @@ end
 
 ### Without contexts
 
-struct ReverseDiffHessianPrep{T} <: HessianPrep
-    tape::T
+struct ReverseDiffHessianGradientPrep{GT,HT} <: HessianPrep
+    gradient_tape::GT
+    hessian_tape::HT
 end
 
 function DI.prepare_hessian(f, ::AutoReverseDiff{Compile}, x) where {Compile}
-    tape = HessianTape(f, x)
+    gradient_tape = GradientTape(f, x)
+    hessian_tape = HessianTape(f, x)
     if Compile
-        tape = compile(tape)
+        gradient_tape = compile(gradient_tape)
+        hessian_tape = compile(hessian_tape)
     end
-    return ReverseDiffHessianPrep(tape)
+    return ReverseDiffHessianGradientPrep(gradient_tape, hessian_tape)
 end
 
-function DI.hessian!(_f, hess, prep::ReverseDiffHessianPrep, ::AutoReverseDiff, x)
-    return hessian!(hess, prep.tape, x)
+function DI.hessian!(_f, hess, prep::ReverseDiffHessianGradientPrep, ::AutoReverseDiff, x)
+    return hessian!(hess, prep.hessian_tape, x)
 end
 
-function DI.hessian(_f, prep::ReverseDiffHessianPrep, ::AutoReverseDiff, x)
-    return hessian!(prep.tape, x)
+function DI.hessian(_f, prep::ReverseDiffHessianGradientPrep, ::AutoReverseDiff, x)
+    return hessian!(prep.hessian_tape, x)
 end
 
 function DI.value_gradient_and_hessian!(
-    f, grad, hess, prep::ReverseDiffHessianPrep, ::AutoReverseDiff, x
+    f, grad, hess, prep::ReverseDiffHessianGradientPrep, ::AutoReverseDiff, x
 )
-    y = f(x)  # TODO: remove once ReverseDiff#251 is fixed
-    result = MutableDiffResult(y, (grad, hess))
-    result = hessian!(result, prep.tape, x)
-    return (
-        DiffResults.value(result), DiffResults.gradient(result), DiffResults.hessian(result)
-    )
+    y = f(x)  # TODO: ReverseDiff#251
+    result = DiffResult(y, (grad, hess))
+    result = hessian!(result, prep.hessian_tape, x)
+    y = DR.value(result)
+    # grad === DR.gradient(result) || copyto!(grad, DR.gradient(result))
+    grad = gradient!(grad, prep.gradient_tape, x)  # TODO: ReverseDiff#251
+    hess === DR.hessian(result) || copyto!(hess, DR.hessian(result))
+    return y, grad, hess
 end
 
 function DI.value_gradient_and_hessian(
-    f, prep::ReverseDiffHessianPrep, ::AutoReverseDiff, x
+    f, prep::ReverseDiffHessianGradientPrep, ::AutoReverseDiff, x
 )
     y = f(x)  # TODO: remove once ReverseDiff#251 is fixed
-    result = MutableDiffResult(y, (similar(x), similar(x, length(x), length(x))))
-    result = hessian!(result, prep.tape, x)
-    return (
-        DiffResults.value(result), DiffResults.gradient(result), DiffResults.hessian(result)
-    )
+    result = DiffResult(y, (similar(x), similar(x, length(x), length(x))))
+    result = hessian!(result, prep.hessian_tape, x)
+    return (DR.value(result), DR.gradient(result), DR.hessian(result))
 end
 
 ### With contexts
@@ -286,22 +297,22 @@ function DI.value_gradient_and_hessian!(
     f, grad, hess, ::NoHessianPrep, ::AutoReverseDiff, x, contexts::Vararg{Context,C}
 ) where {C}
     fc = with_contexts(f, contexts...)
-    y = fc(x)  # TODO: remove once ReverseDiff#251 is fixed
-    result = MutableDiffResult(y, (grad, hess))
+    y = fc(x)  # TODO: ReverseDiff#251
+    result = DiffResult(y, (grad, hess))
     result = hessian!(result, fc, x)
-    return (
-        DiffResults.value(result), DiffResults.gradient(result), DiffResults.hessian(result)
-    )
+    y = DR.value(result)
+    # grad === DR.gradient(result) || copyto!(grad, DR.gradient(result))
+    grad = gradient!(grad, fc, x)  # TODO: ReverseDiff#251
+    hess === DR.hessian(result) || copyto!(hess, DR.hessian(result))
+    return y, grad, hess
 end
 
 function DI.value_gradient_and_hessian(
     f, ::NoHessianPrep, ::AutoReverseDiff, x, contexts::Vararg{Context,C}
 ) where {C}
     fc = with_contexts(f, contexts...)
-    y = fc(x)  # TODO: remove once ReverseDiff#251 is fixed
-    result = MutableDiffResult(y, (similar(x), similar(x, length(x), length(x))))
+    y = fc(x)  # TODO: ReverseDiff#251
+    result = HessianResult(x)
     result = hessian!(result, fc, x)
-    return (
-        DiffResults.value(result), DiffResults.gradient(result), DiffResults.hessian(result)
-    )
+    return (DR.value(result), DR.gradient(result), DR.hessian(result))
 end
diff --git a/DifferentiationInterface/test/Back/Zygote/test.jl b/DifferentiationInterface/test/Back/Zygote/test.jl
index f75cd5ba7..fe868626e 100644
--- a/DifferentiationInterface/test/Back/Zygote/test.jl
+++ b/DifferentiationInterface/test/Back/Zygote/test.jl
@@ -40,7 +40,7 @@ test_differentiation(
 if VERSION >= v"1.10"
     test_differentiation(
         AutoZygote(),
-        vcat(component_scenarios(), gpu_scenarios(), static_scenarios());
+        vcat(component_scenarios(), gpu_scenarios());
         second_order=false,
         logging=LOGGING,
     )
diff --git a/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestComponentArraysExt/DifferentiationInterfaceTestComponentArraysExt.jl b/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestComponentArraysExt/DifferentiationInterfaceTestComponentArraysExt.jl
index d15960d62..ff3f7bf15 100644
--- a/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestComponentArraysExt/DifferentiationInterfaceTestComponentArraysExt.jl
+++ b/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestComponentArraysExt/DifferentiationInterfaceTestComponentArraysExt.jl
@@ -25,9 +25,7 @@ function comp_to_num_pullback(x, dy)
 end
 
 function comp_to_num_scenarios_onearg(x::ComponentVector; dx::AbstractVector, dy::Number)
-    nb_args = 1
     f = comp_to_num
-    y = f(x)
     dy_from_dx = comp_to_num_pushforward(x, dx)
     dx_from_dy = comp_to_num_pullback(x, dy)
     grad = comp_to_num_gradient(x)
diff --git a/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestJLArraysExt/DifferentiationInterfaceTestJLArraysExt.jl b/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestJLArraysExt/DifferentiationInterfaceTestJLArraysExt.jl
index 6b99e4d6f..278d899fe 100644
--- a/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestJLArraysExt/DifferentiationInterfaceTestJLArraysExt.jl
+++ b/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestJLArraysExt/DifferentiationInterfaceTestJLArraysExt.jl
@@ -1,54 +1,42 @@
 module DifferentiationInterfaceTestJLArraysExt
 
+using Compat
+import DifferentiationInterface as DI
 using DifferentiationInterfaceTest
 import DifferentiationInterfaceTest as DIT
 using JLArrays: JLArray, jl
 using Random: AbstractRNG, default_rng
 
-num_to_arr_jlvector(x) = DIT.num_to_arr(x, JLArray{Float64,1})
-num_to_arr_jlmatrix(x) = DIT.num_to_arr(x, JLArray{Float64,2})
-
-DIT.pick_num_to_arr(::Type{<:JLArray{<:Real,1}}) = num_to_arr_jlvector
-DIT.pick_num_to_arr(::Type{<:JLArray{<:Real,2}}) = num_to_arr_jlmatrix
-
-function DIT.gpu_scenarios(rng::AbstractRNG=default_rng(); linalg=true)
-    x_ = rand(rng)
-    dx_ = rand(rng)
-    dy_ = rand(rng)
-
-    x_6 = jl(rand(rng, 6))
-    dx_6 = jl(rand(rng, 6))
-
-    x_2_3 = jl(rand(rng, 2, 3))
-    dx_2_3 = jl(rand(rng, 2, 3))
-
-    dy_12 = jl(rand(rng, 12))
-    dy_6_2 = jl(rand(rng, 6, 2))
-    dy_6 = jl(rand(rng, 6))
-    dy_2_3 = jl(rand(rng, 2, 3))
-
-    V = typeof(dy_6)
-    M = typeof(dy_2_3)
+myjl(f::Function) = f
+function myjl(::DIT.NumToArr{A}) where {T,N,A<:AbstractArray{T,N}}
+    return DIT.NumToArr(JLArray{T,N})
+end
 
-    scens = vcat(
-        # one argument
-        DIT.num_to_arr_scenarios_onearg(x_, V; dx=dx_, dy=dy_6),
-        DIT.num_to_arr_scenarios_onearg(x_, M; dx=dx_, dy=dy_2_3),
-        DIT.arr_to_num_scenarios_onearg(x_6; dx=dx_6, dy=dy_, linalg),
-        DIT.arr_to_num_scenarios_onearg(x_2_3; dx=dx_2_3, dy=dy_, linalg),
-        DIT.vec_to_vec_scenarios_onearg(x_6; dx=dx_6, dy=dy_12),
-        DIT.vec_to_mat_scenarios_onearg(x_6; dx=dx_6, dy=dy_6_2),
-        DIT.mat_to_vec_scenarios_onearg(x_2_3; dx=dx_2_3, dy=dy_12),
-        DIT.mat_to_mat_scenarios_onearg(x_2_3; dx=dx_2_3, dy=dy_6_2),
-        # two arguments
-        DIT.num_to_arr_scenarios_twoarg(x_, V; dx=dx_, dy=dy_6),
-        DIT.num_to_arr_scenarios_twoarg(x_, M; dx=dx_, dy=dy_2_3),
-        DIT.vec_to_vec_scenarios_twoarg(x_6; dx=dx_6, dy=dy_12),
-        DIT.vec_to_mat_scenarios_twoarg(x_6; dx=dx_6, dy=dy_6_2),
-        DIT.mat_to_vec_scenarios_twoarg(x_2_3; dx=dx_2_3, dy=dy_12),
-        DIT.mat_to_mat_scenarios_twoarg(x_2_3; dx=dx_2_3, dy=dy_6_2),
+myjl(f::DIT.MultiplyByConstant) = f
+myjl(f::DIT.WritableClosure) = f
+
+myjl(x::Number) = x
+myjl(x::AbstractArray) = jl(x)
+myjl(x::Tuple) = map(myjl, x)
+myjl(x::DI.Constant) = DI.Constant(myjl(DI.unwrap(x)))
+myjl(::Nothing) = nothing
+
+function myjl(scen::Scenario{op,pl_op,pl_fun}) where {op,pl_op,pl_fun}
+    @compat (; f, x, y, tang, contexts, res1, res2) = scen
+    return Scenario{op,pl_op,pl_fun}(
+        myjl(f);
+        x=myjl(x),
+        y=myjl(y),
+        tang=myjl(tang),
+        contexts=myjl(contexts),
+        res1=myjl(res1),
+        res2=myjl(res2),
     )
-    return scens
+end
+
+function DIT.gpu_scenarios(args...; kwargs...)
+    scens = DIT.default_scenarios(args...; kwargs...)
+    return myjl.(scens)
 end
 
 end
diff --git a/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestStaticArraysExt/DifferentiationInterfaceTestStaticArraysExt.jl b/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestStaticArraysExt/DifferentiationInterfaceTestStaticArraysExt.jl
index d2bbd6e6a..34e58d68d 100644
--- a/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestStaticArraysExt/DifferentiationInterfaceTestStaticArraysExt.jl
+++ b/DifferentiationInterfaceTest/ext/DifferentiationInterfaceTestStaticArraysExt/DifferentiationInterfaceTestStaticArraysExt.jl
@@ -1,69 +1,79 @@
 module DifferentiationInterfaceTestStaticArraysExt
 
+using Compat
+import DifferentiationInterface as DI
 using DifferentiationInterfaceTest
 import DifferentiationInterfaceTest as DIT
 using Random: AbstractRNG, default_rng
 using SparseArrays: SparseArrays, SparseMatrixCSC, nnz, spdiagm
 using StaticArrays: MArray, MMatrix, MVector, SArray, SMatrix, SVector
 
-num_to_arr_svector(x) = DIT.num_to_arr(x, SVector{6,Float64})
-num_to_arr_smatrix(x) = DIT.num_to_arr(x, SMatrix{2,3,Float64,6})
-
-DIT.pick_num_to_arr(::Type{<:SVector}) = num_to_arr_svector
-DIT.pick_num_to_arr(::Type{<:SMatrix}) = num_to_arr_smatrix
-
-function DIT.static_scenarios(rng::AbstractRNG=default_rng(); linalg=true)
-    x_ = rand(rng)
-    dx_ = rand(rng)
-    dy_ = rand(rng)
-
-    x_6 = rand(rng, 6)
-    dx_6 = rand(rng, 6)
-
-    x_2_3 = rand(rng, 2, 3)
-    dx_2_3 = rand(rng, 2, 3)
-
-    dy_6 = rand(rng, 6)
-    dy_12 = rand(rng, 12)
-    dy_2_3 = rand(rng, 2, 3)
-    dy_6_2 = rand(rng, 6, 2)
-
-    SV_6 = SVector{6}
-    MV_6 = MVector{6}
-    SV_12 = SVector{12}
-    MV_12 = MVector{12}
-
-    SM_2_3 = SMatrix{2,3}
-    MM_2_3 = MMatrix{2,3}
-    SM_6_2 = SMatrix{6,2}
-    MM_6_2 = MMatrix{6,2}
-
-    scens = vcat(
-        # one argument
-        DIT.num_to_arr_scenarios_onearg(x_, SV_6; dx=dx_, dy=SV_6(dy_6)),
-        DIT.num_to_arr_scenarios_onearg(x_, SM_2_3; dx=dx_, dy=SM_2_3(dy_2_3)),
-        DIT.arr_to_num_scenarios_onearg(SV_6(x_6); dx=SV_6(dx_6), dy=dy_, linalg),
-        DIT.arr_to_num_scenarios_onearg(SM_2_3(x_2_3); dx=SM_2_3(dx_2_3), dy=dy_, linalg),
-        DIT.vec_to_vec_scenarios_onearg(SV_6(x_6); dx=SV_6(dx_6), dy=SV_12(dy_12)),
-        DIT.vec_to_mat_scenarios_onearg(SV_6(x_6); dx=SV_6(dx_6), dy=SM_6_2(dy_6_2)),
-        DIT.mat_to_vec_scenarios_onearg(SM_2_3(x_2_3); dx=SM_2_3(dx_2_3), dy=SV_12(dy_12)),
-        DIT.mat_to_mat_scenarios_onearg(
-            SM_2_3(x_2_3); dx=SM_2_3(dx_2_3), dy=SM_6_2(dy_6_2)
-        ),
-        # two arguments
-        DIT.num_to_arr_scenarios_twoarg(x_, MV_6; dx=dx_, dy=MV_6(dy_6)),
-        DIT.num_to_arr_scenarios_twoarg(x_, MM_2_3; dx=dx_, dy=MM_2_3(dy_2_3)),
-        DIT.vec_to_vec_scenarios_twoarg(MV_6(x_6); dx=MV_6(dx_6), dy=MV_12(dy_12)),
-        DIT.vec_to_mat_scenarios_twoarg(MV_6(x_6); dx=MV_6(dx_6), dy=MM_6_2(dy_6_2)),
-        DIT.mat_to_vec_scenarios_twoarg(MM_2_3(x_2_3); dx=MM_2_3(dx_2_3), dy=MV_12(dy_12)),
-        DIT.mat_to_mat_scenarios_twoarg(
-            MM_2_3(x_2_3); dx=MM_2_3(dx_2_3), dy=MM_6_2(dy_6_2)
-        ),
+mySArray(f::Function) = f
+myMArray(f::Function) = f
+
+mySArray(::DIT.NumToArr{A}) where {T,A<:AbstractVector{T}} = DIT.NumToArr(SVector{6,T})
+myMArray(::DIT.NumToArr{A}) where {T,A<:AbstractVector{T}} = DIT.NumToArr(MVector{6,T})
+
+mySArray(::DIT.NumToArr{A}) where {T,A<:AbstractMatrix{T}} = DIT.NumToArr(SMatrix{2,3,T,6})
+myMArray(::DIT.NumToArr{A}) where {T,A<:AbstractMatrix{T}} = DIT.NumToArr(MMatrix{2,3,T,6})
+
+mySArray(f::DIT.MultiplyByConstant) = f
+myMArray(f::DIT.MultiplyByConstant) = f
+
+mySArray(f::DIT.WritableClosure) = f
+myMArray(f::DIT.WritableClosure) = f
+
+mySArray(x::Number) = x
+myMArray(x::Number) = x
+
+mySArray(x::AbstractVector{T}) where {T} = convert(SVector{length(x),T}, x)
+myMArray(x::AbstractVector{T}) where {T} = convert(MVector{length(x),T}, x)
+
+function mySArray(x::AbstractMatrix{T}) where {T}
+    return convert(SMatrix{size(x, 1),size(x, 2),T,length(x)}, x)
+end
+function myMArray(x::AbstractMatrix{T}) where {T}
+    return convert(MMatrix{size(x, 1),size(x, 2),T,length(x)}, x)
+end
+
+mySArray(x::Tuple) = map(mySArray, x)
+myMArray(x::Tuple) = map(myMArray, x)
+
+mySArray(x::DI.Constant) = DI.Constant(mySArray(DI.unwrap(x)))
+myMArray(x::DI.Constant) = DI.Constant(myMArray(DI.unwrap(x)))
+
+mySArray(::Nothing) = nothing
+myMArray(::Nothing) = nothing
+
+function mySArray(scen::Scenario{op,pl_op,pl_fun}) where {op,pl_op,pl_fun}
+    @compat (; f, x, y, tang, contexts, res1, res2) = scen
+    return Scenario{op,pl_op,pl_fun}(
+        mySArray(f);
+        x=mySArray(x),
+        y=pl_fun == :in ? myMArray(y) : mySArray(y),
+        tang=mySArray(tang),
+        contexts=mySArray(contexts),
+        res1=mySArray(res1),
+        res2=mySArray(res2),
     )
-    scens = filter(scens) do s
-        DIT.operator_place(s) == :out || s.x isa Union{Number,MArray}
-    end
-    return scens
+end
+
+function myMArray(scen::Scenario{op,pl_op,pl_fun}) where {op,pl_op,pl_fun}
+    @compat (; f, x, y, tang, contexts, res1, res2) = scen
+    return Scenario{op,pl_op,pl_fun}(
+        myMArray(f);
+        x=myMArray(x),
+        y=pl_fun == :in ? myMArray(y) : myMArray(y),
+        tang=myMArray(tang),
+        contexts=myMArray(contexts),
+        res1=myMArray(res1),
+        res2=myMArray(res2),
+    )
+end
+
+function DIT.static_scenarios(args...; kwargs...)
+    scens = DIT.default_scenarios(args...; kwargs...)
+    return vcat(mySArray.(scens), myMArray.(scens))
 end
 
 end
diff --git a/DifferentiationInterfaceTest/src/scenarios/allocfree.jl b/DifferentiationInterfaceTest/src/scenarios/allocfree.jl
index c2b095491..4707031f8 100644
--- a/DifferentiationInterfaceTest/src/scenarios/allocfree.jl
+++ b/DifferentiationInterfaceTest/src/scenarios/allocfree.jl
@@ -1,6 +1,5 @@
 function identity_scenarios(x::Number; dx::Number, dy::Number)
     f = identity
-    y = f(x)
     dy_from_dx = dx
     dx_from_dy = dy
     der = one(x)
@@ -14,7 +13,6 @@ end
 
 function sum_scenarios(x::AbstractArray; dx::AbstractArray, dy::Number)
     f = sum
-    y = f(x)
     dy_from_dx = sum(dx)
     dx_from_dy = (similar(x) .= dy)
     grad = similar(x)
@@ -48,7 +46,7 @@ end
 Create a vector of [`Scenario`](@ref)s with functions that do not allocate.
 
 !!! warning
-    At the moment, some second-order scenarios are excluded.
+    At the moment, second-order scenarios are excluded.
 """
 function allocfree_scenarios(rng::AbstractRNG=default_rng())
     x_ = rand(rng)
diff --git a/DifferentiationInterfaceTest/src/scenarios/default.jl b/DifferentiationInterfaceTest/src/scenarios/default.jl
index 99d466d0a..25b85eda6 100644
--- a/DifferentiationInterfaceTest/src/scenarios/default.jl
+++ b/DifferentiationInterfaceTest/src/scenarios/default.jl
@@ -71,22 +71,20 @@ end
 
 ## Number to array
 
-multiplicator(::Type{V}) where {V<:AbstractVector} = convert(V, float.(1:6))
-multiplicator(::Type{M}) where {M<:AbstractMatrix} = convert(M, float.(reshape(1:6, 2, 3)))
+multiplicator(::Type{A}) where {A<:AbstractVector} = convert(A, float.(1:6))
+multiplicator(::Type{A}) where {A<:AbstractMatrix} = convert(A, reshape(float.(1:6), 2, 3))
 
-function num_to_arr(x::Number, ::Type{A}) where {A<:AbstractArray}
+struct NumToArr{A} end
+NumToArr(::Type{A}) where {A} = NumToArr{A}()
+Base.eltype(::NumToArr{A}) where {A} = eltype(A)
+
+function (f::NumToArr{A})(x::Number) where {A}
     a = multiplicator(A)
     return sin.(x .* a)
 end
 
-num_to_arr_vector(x) = num_to_arr(x, Vector{Float64})
-num_to_arr_matrix(x) = num_to_arr(x, Matrix{Float64})
-
-pick_num_to_arr(::Type{<:Vector}) = num_to_arr_vector
-pick_num_to_arr(::Type{<:Matrix}) = num_to_arr_matrix
-
-function num_to_arr!(y::AbstractArray, x::Number)::Nothing
-    a = multiplicator(typeof(y))
+function (f!::NumToArr{A})(y::AbstractArray, x::Number)::Nothing where {A}
+    a = multiplicator(A)
     y .= sin.(x .* a)
     return nothing
 end
@@ -114,8 +112,7 @@ end
 function num_to_arr_scenarios_onearg(
     x::Number, ::Type{A}; dx::Number, dy::AbstractArray
 ) where {A<:AbstractArray}
-    f = pick_num_to_arr(A)
-    y = f(x)
+    f = NumToArr(A)
     dy_from_dx = num_to_arr_pushforward(x, dx, A)
     dx_from_dy = num_to_arr_pullback(x, dy, A)
     der = num_to_arr_derivative(x, A)
@@ -142,8 +139,8 @@ end
 function num_to_arr_scenarios_twoarg(
     x::Number, ::Type{A}; dx::Number, dy::AbstractArray
 ) where {A<:AbstractArray}
-    f! = num_to_arr!
-    y = similar(num_to_arr(x, A))
+    f! = NumToArr(A)
+    y = similar(NumToArr(A)(x))
     f!(y, x)
     dy_from_dx = num_to_arr_pushforward(x, dx, A)
     dx_from_dy = num_to_arr_pullback(x, dy, A)
@@ -170,9 +167,12 @@ end
 
 ## Array to number
 
-arr_to_num_aux_linalg(x; α, β) = sum(vec(x .^ α) .* transpose(vec(x .^ β)))
+const α = 4
+const β = 6
+
+arr_to_num_linalg(x::AbstractArray) = sum(vec(x .^ α) .* transpose(vec(x .^ β)))
 
-function arr_to_num_aux_no_linalg(x; α, β)
+function arr_to_num_no_linalg(x::AbstractArray)
     n = length(x)
     s = zero(eltype(x))
     for i in 1:n, j in 1:n
@@ -181,7 +181,7 @@ function arr_to_num_aux_no_linalg(x; α, β)
     return s
 end
 
-function arr_to_num_aux_gradient(x0; α, β)
+function arr_to_num_gradient(x0)
     x = Array(x0)  # GPU arrays don't like indexing
     g = similar(x)
     for k in eachindex(g, x)
@@ -194,7 +194,7 @@ function arr_to_num_aux_gradient(x0; α, β)
     return convert(typeof(x0), g)
 end
 
-function arr_to_num_aux_hessian(x0; α, β)
+function arr_to_num_hessian(x0)
     x = Array(x0)  # GPU arrays don't like indexing
     H = similar(x, length(x), length(x))
     for k in axes(H, 1), l in axes(H, 2)
@@ -211,16 +211,6 @@ function arr_to_num_aux_hessian(x0; α, β)
     return convert(typeof(similar(x0, length(x0), length(x0))), H)
 end
 
-const DEFAULT_α = 4
-const DEFAULT_β = 6
-
-arr_to_num_linalg(x::AbstractArray)::Number =
-    arr_to_num_aux_linalg(x; α=DEFAULT_α, β=DEFAULT_β)
-arr_to_num_no_linalg(x::AbstractArray)::Number =
-    arr_to_num_aux_no_linalg(x; α=DEFAULT_α, β=DEFAULT_β)
-
-arr_to_num_gradient(x) = arr_to_num_aux_gradient(x; α=DEFAULT_α, β=DEFAULT_β)
-arr_to_num_hessian(x) = arr_to_num_aux_hessian(x; α=DEFAULT_α, β=DEFAULT_β)
 arr_to_num_pushforward(x, dx) = dot(arr_to_num_gradient(x), dx)
 arr_to_num_pullback(x, dy) = arr_to_num_gradient(x) .* dy
 arr_to_num_hvp(x, dx) = reshape(arr_to_num_hessian(x) * vec(dx), size(x))
@@ -229,7 +219,6 @@ function arr_to_num_scenarios_onearg(
     x::AbstractArray; dx::AbstractArray, dy::Number, linalg=true
 )
     f = linalg ? arr_to_num_linalg : arr_to_num_no_linalg
-    y = f(x)
     dy_from_dx = arr_to_num_pushforward(x, dx)
     dx_from_dy = arr_to_num_pullback(x, dy)
     grad = arr_to_num_gradient(x)
@@ -305,7 +294,6 @@ function vec_to_vec_scenarios_onearg(
     x::AbstractVector; dx::AbstractVector, dy::AbstractVector
 )
     f = vec_to_vec
-    y = f(x)
     dy_from_dx = vec_to_vec_pushforward(x, dx)
     dx_from_dy = vec_to_vec_pullback(x, dy)
     jac = vec_to_vec_jacobian(x)
@@ -344,7 +332,6 @@ function vec_to_mat_scenarios_onearg(
     x::AbstractVector; dx::AbstractVector, dy::AbstractMatrix
 )
     f = vec_to_mat
-    y = f(x)
     dy_from_dx = vec_to_mat_pushforward(x, dx)
     dx_from_dy = vec_to_mat_pullback(x, dy)
     jac = vec_to_mat_jacobian(x)
@@ -391,7 +378,6 @@ function mat_to_vec_scenarios_onearg(
     x::AbstractMatrix; dx::AbstractMatrix, dy::AbstractVector
 )
     f = mat_to_vec
-    y = f(x)
     dy_from_dx = mat_to_vec_pushforward(x, dx)
     dx_from_dy = mat_to_vec_pullback(x, dy)
     jac = mat_to_vec_jacobian(x)
@@ -437,7 +423,6 @@ function mat_to_mat_scenarios_onearg(
     x::AbstractMatrix; dx::AbstractMatrix, dy::AbstractMatrix
 )
     f = mat_to_mat
-    y = f(x)
     dy_from_dx = mat_to_mat_pushforward(x, dx)
     dx_from_dy = mat_to_mat_pullback(x, dy)
     jac = mat_to_mat_jacobian(x)
diff --git a/DifferentiationInterfaceTest/src/scenarios/scenario.jl b/DifferentiationInterfaceTest/src/scenarios/scenario.jl
index 1692057a2..30484c469 100644
--- a/DifferentiationInterfaceTest/src/scenarios/scenario.jl
+++ b/DifferentiationInterfaceTest/src/scenarios/scenario.jl
@@ -68,7 +68,6 @@ function Scenario{op,pl_op}(
 ) where {op,pl_op}
     @assert op in ALL_OPS
     @assert pl_op in (:in, :out)
-    f!(y, x, map(unwrap, contexts)...)
     return Scenario{op,pl_op,:in}(f!; x, y, tang, contexts, res1, res2)
 end
 
diff --git a/DifferentiationInterfaceTest/test/weird.jl b/DifferentiationInterfaceTest/test/weird.jl
index 2b19baeb2..dc3f17567 100644
--- a/DifferentiationInterfaceTest/test/weird.jl
+++ b/DifferentiationInterfaceTest/test/weird.jl
@@ -21,11 +21,20 @@ using Zygote: Zygote
 
 LOGGING = get(ENV, "CI", "false") == "false"
 
-## Weird arrays
+## Generate all scenarios
 
-test_differentiation(
-    AutoForwardDiff(), vcat(component_scenarios(), static_scenarios()); logging=LOGGING
+gpu_scenarios(;
+    include_constantified=true, include_closurified=true, include_batchified=true
 )
+static_scenarios(;
+    include_constantified=true, include_closurified=true, include_batchified=true
+)
+
+## Weird arrays
+
+test_differentiation(AutoForwardDiff(), static_scenarios(); logging=LOGGING)
+
+test_differentiation(AutoForwardDiff(), component_scenarios(); logging=LOGGING)
 
 test_differentiation(AutoZygote(), gpu_scenarios(); second_order=false, logging=LOGGING)