Hacky way of obtaining `n_vars` is dangerous when using different data types

During my recent ventures into floating point data types beyond `Float64` I encountered crashing simulations/callbacks when not being careful.

In particular,

```julia
# Reinterpret the solution array as an array of conservative variables,
# compute the solution variables via broadcasting, and reinterpret the
# result as a plain array of floating point numbers
data = Array(reinterpret(eltype(u),
solution_variables.(reinterpret(SVector{nvariables(equations),
eltype(u)}, u),
Ref(equations))))

# Find out variable count by looking at output from `solution_variables` function
n_vars = size(data, 1)
```

is dangerous if datatypes are not consistent across solver, equations, mesh, ode, callbacks etc.

We use this hack at different places: https://github.com/search?q=repo%3Atrixi-framework%2FTrixi.jl+%22n_vars+%3D+size%28%22&type=code

The issue points to multiple repository locations using `n_vars = size(...)` and shows a `reinterpret`/`solution_variables` path; begin with the linked code search and trace those call sites through solver, equations, mesh, ODE, and callback setup. Compare how variable counts and element types are obtained at each site. Done means the dangerous inference is removed consistently and mixed floating-point types no longer cause simulation or callback crashes.

Written by the indexing model from the issue text.