Cannot use implicit time integration with automatic differentiation

Automatic differentiation is mentioned as one of the key features in the README.md, so I would expect this to work.
```julia
julia> using Trixi, OrdinaryDiffEqSDIRK

julia> trixi_include(default_example(), sol = nothing, stepsize_callback = nothing);
[ Info: You just called `trixi_include`. Julia may now compile the code, please be patient.

julia> sol = solve(ode, ImplicitEuler(); ode_default_options()..., callback=callbacks)

████████╗██████╗ ██╗██╗ ██╗██╗
╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
██║ ██████╔╝██║ ╚███╔╝ ██║
██║ ██╔══██╗██║ ██╔██╗ ██║
██║ ██║ ██║██║██╔╝ ██╗██║
╚═╝ ╚═╝ ╚═╝╚═╝╚═╝ ╚═╝╚═╝

[...]

────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 3.87500000e-06 s
Δt: 0.00000000e+00 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: NaN s
PID: Inf s
#DOFs per field: 4096 alloc'd memory: 361.143 MiB
#elements: 256

Variable: scalar
L2 error: 3.94120825e-06
Linf error: 1.23670425e-05
∑∂S/∂U ⋅ Uₜ : -2.19124035e-16
────────────────────────────────────────────────────────────────────────────────────────────────────

ERROR: First call to automatic differentiation for the Jacobian
failed. This means that the user `f` function is not compatible
with automatic differentiation. Methods to fix this include:

1. Turn off automatic differentiation (e.g. Rosenbrock23() becomes
Rosenbrock23(autodiff = AutoFiniteDiff())). More details can befound at
https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/
2. Improving the compatibility of `f` with ForwardDiff.jl automatic
differentiation (using tools like PreallocationTools.jl). More details
can be found at https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers
3. Defining analytical Jacobians. More details can be
found at https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/#SciMLBase.ODEFunction

Note: turning off automatic differentiation tends to have a very minimal
performance impact (for this use case, because it's forward mode for a
square Jacobian. This is different from optimization gradient scenarios).
However, one should be careful as some methods are more sensitive to
accurate gradients than others. Specifically, Rodas methods like `Rodas4`
and `Rodas5P` require accurate Jacobians in order to have good convergence,
while many other methods like BDF (`QNDF`, `FBDF`), SDIRK (`KenCarp4`),
and Rosenbrock-W (`Rosenbrock23`) do not. Thus if using an algorithm which
is sensitive to autodiff and solving at a low tolerance, please change the
algorithm as well.

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:265
(::Type{T})(::T) where T<:Number
@ Core boot.jl:900
Float64(::Irrational{:mad_constant})
@ StatsBase irrationals.jl:251
...

Stacktrace:
[1] jacobian!(J::Matrix{…}, f::Function, x::Vector{…}, fx::Vector{…}, integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, jac_config::Tuple{…})
@ OrdinaryDiffEqDifferentiation ~/.julia/packages/OrdinaryDiffEqDifferentiation/M5647/src/derivative_wrappers.jl:223
[2] calc_J!
@ ~/.julia/packages/OrdinaryDiffEqDifferentiation/M5647/src/derivative_utils.jl:222 [inlined]
[3] calc_W!(W::Matrix{…}, integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, nlsolver::OrdinaryDiffEqNonlinearSolve.NLSolver{…}, cache::OrdinaryDiffEqSDIRK.ImplicitEulerCache{…}, dtgamma::Float64, repeat_step::Bool, newJW::Nothing)
@ OrdinaryDiffEqDifferentiation ~/.julia/packages/OrdinaryDiffEqDifferentiation/M5647/src/derivative_utils.jl:627
[4] update_W!
@ ~/.julia/packages/OrdinaryDiffEqDifferentiation/M5647/src/derivative_utils.jl:719 [inlined]
[5] update_W!
@ ~/.julia/packages/OrdinaryDiffEqDifferentiation/M5647/src/derivative_utils.jl:718 [inlined]
[6] nlsolve!(nlsolver::OrdinaryDiffEqNonlinearSolve.NLSolver{…}, integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqSDIRK.ImplicitEulerCache{…}, repeat_step::Bool)
@ OrdinaryDiffEqNonlinearSolve ~/.julia/packages/OrdinaryDiffEqNonlinearSolve/e4hoO/src/nlsolve.jl:27
[7] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqSDIRK.ImplicitEulerCache{…}, repeat_step::Bool)
@ OrdinaryDiffEqSDIRK ~/.julia/packages/OrdinaryDiffEqSDIRK/YAA9M/src/sdirk_perform_step.jl:94
[8] perform_step!
@ ~/.julia/packages/OrdinaryDiffEqSDIRK/YAA9M/src/sdirk_perform_step.jl:79 [inlined]
[9] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:620
[10] __solve(::ODEProblem{…}, ::ImplicitEuler{…}; kwargs::@Kwargs{…})
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:7
[11] __solve
@ ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:1 [inlined]
[12] #solve_call#35
@ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:635 [inlined]
[13] solve_call
@ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:592 [inlined]
[14] #solve_up#44
@ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1142 [inlined]
[15] solve_up
@ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1120 [inlined]
[16] #solve#42
@ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1057 [inlined]
[17] top-level scope
@ REPL[7]:1

caused by: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
(::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
@ Base rounding.jl:265
(::Type{T})(::T) where T<:Number
@ Core boot.jl:900
Float64(::Irrational{:mad_constant})
@ StatsBase irrationals.jl:251
...

Stacktrace:
[1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12})
@ Base ./number.jl:7
[2] setindex!
@ ./array.jl:994 [inlined]
[3] macro expansion
@ ~/git/Trixi.jl/src/solvers/dgsem_tree/dg_2d.jl:557 [inlined]
[4] macro expansion
@ ~/.julia/packages/Polyester/eqrC9/src/closure.jl:525 [inlined]
[5] macro expansion
@ ~/git/Trixi.jl/src/auxiliary/auxiliary.jl:213 [inlined]
[6] prolong2interfaces!(cache::@NamedTuple{…}, u::Array{…}, mesh::TreeMesh{…}, equations::LinearScalarAdvectionEquation2D{…}, surface_integral::SurfaceIntegralWeakForm{…}, dg::DGSEM{…})
@ Trixi ~/git/Trixi.jl/src/solvers/dgsem_tree/dg_2d.jl:550
[7] macro expansion
@ ~/git/Trixi.jl/src/solvers/dgsem_tree/dg_2d.jl:130 [inlined]
[8] macro expansion
@ ~/git/TrixiBase.jl/src/trixi_timeit.jl:64 [inlined]
[9] rhs!(du::Array{…}, u::Array{…}, t::Float64, mesh::TreeMesh{…}, equations::LinearScalarAdvectionEquation2D{…}, boundary_conditions::Trixi.BoundaryConditionPeriodic, source_terms::Nothing, dg::DGSEM{…}, cache::@NamedTuple{…})
@ Trixi ~/git/Trixi.jl/src/solvers/dgsem_tree/dg_2d.jl:129
[10] macro expansion
@ ~/git/TrixiBase.jl/src/trixi_timeit.jl:64 [inlined]
[11] rhs!(du_ode::Vector{…}, u_ode::Vector{…}, semi::SemidiscretizationHyperbolic{…}, t::Float64)
@ Trixi ~/git/Trixi.jl/src/semidiscretization/semidiscretization_hyperbolic.jl:421
[12] ODEFunction
@ ~/.julia/packages/SciMLBase/c6Noy/src/scimlfunctions.jl:2470 [inlined]
[13] UJacobianWrapper
@ ~/.julia/packages/SciMLBase/c6Noy/src/function_wrappers.jl:32 [inlined]
[14] chunk_mode_jacobian!(result::Matrix{…}, f!::SciMLBase.UJacobianWrapper{…}, y::Vector{…}, x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
@ ForwardDiff ~/.julia/packages/ForwardDiff/UBbGT/src/jacobian.jl:187
[15] jacobian!
@ ~/.julia/packages/ForwardDiff/UBbGT/src/jacobian.jl:84 [inlined]
[16] jacobian!(::SciMLBase.UJacobianWrapper{…}, ::Vector{…}, ::Matrix{…}, ::DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{…}, ::ADTypes.AutoForwardDiff{…}, ::Vector{…})
@ DifferentiationInterfaceForwardDiffExt ~/.julia/packages/DifferentiationInterface/Yk2Kt/ext/DifferentiationInterfaceForwardDiffExt/twoarg.jl:489
[17] jacobian!(J::Matrix{…}, f::Function, x::Vector{…}, fx::Vector{…}, integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, jac_config::Tuple{…})
@ OrdinaryDiffEqDifferentiation ~/.julia/packages/OrdinaryDiffEqDifferentiation/M5647/src/derivative_wrappers.jl:221
[18] calc_J!
@ ~/.julia/packages/OrdinaryDiffEqDifferentiation/M5647/src/derivative_utils.jl:222 [inlined]
[19] calc_W!(W::Matrix{…}, integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, nlsolver::OrdinaryDiffEqNonlinearSolve.NLSolver{…}, cache::OrdinaryDiffEqSDIRK.ImplicitEulerCache{…}, dtgamma::Float64, repeat_step::Bool, newJW::Nothing)
@ OrdinaryDiffEqDifferentiation ~/.julia/packages/OrdinaryDiffEqDifferentiation/M5647/src/derivative_utils.jl:627
[20] update_W!
@ ~/.julia/packages/OrdinaryDiffEqDifferentiation/M5647/src/derivative_utils.jl:719 [inlined]
[21] update_W!
@ ~/.julia/packages/OrdinaryDiffEqDifferentiation/M5647/src/derivative_utils.jl:718 [inlined]
[22] nlsolve!(nlsolver::OrdinaryDiffEqNonlinearSolve.NLSolver{…}, integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqSDIRK.ImplicitEulerCache{…}, repeat_step::Bool)
@ OrdinaryDiffEqNonlinearSolve ~/.julia/packages/OrdinaryDiffEqNonlinearSolve/e4hoO/src/nlsolve.jl:27
[23] perform_step!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…}, cache::OrdinaryDiffEqSDIRK.ImplicitEulerCache{…}, repeat_step::Bool)
@ OrdinaryDiffEqSDIRK ~/.julia/packages/OrdinaryDiffEqSDIRK/YAA9M/src/sdirk_perform_step.jl:94
[24] perform_step!
@ ~/.julia/packages/OrdinaryDiffEqSDIRK/YAA9M/src/sdirk_perform_step.jl:79 [inlined]
[25] solve!(integrator::OrdinaryDiffEqCore.ODEIntegrator{…})
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:620
[26] __solve(::ODEProblem{…}, ::ImplicitEuler{…}; kwargs::@Kwargs{…})
@ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:7
[27] __solve
@ ~/.julia/packages/OrdinaryDiffEqCore/bMOsj/src/solve.jl:1 [inlined]
[28] #solve_call#35
@ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:635 [inlined]
[29] solve_call
@ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:592 [inlined]
[30] #solve_up#44
@ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1142 [inlined]
[31] solve_up
@ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1120 [inlined]
[32] #solve#42
@ ~/.julia/packages/DiffEqBase/zYZst/src/solve.jl:1057 [inlined]
[33] top-level scope
@ REPL[7]:1
Some type information was truncated. Use `show(err)` to see complete types.
```

Reproduce the example with `ImplicitEuler()` and automatic differentiation, then inspect `src/solvers/dgsem_tree/dg_2d.jl` around line 550 and `src/semidiscretization/semidiscretization_hyperbolic.jl` around line 421, where the stack trace reaches the failing Float64 conversion. Done means the reported implicit solve completes with automatic differentiation enabled and the relevant regression coverage passes.

Written by the indexing model from the issue text.