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

[DanielDoehring]

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

In particular,

# 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

[jlchan]

I'm trying to see the issue - is n_vars computed incorrectly if eltype(u) is not Float64 or Float32? I would have guessed that eltype(u) deals with changes in the data type.

[DanielDoehring]

Yeah I think it is not obvious. The way I encountered this issue was when I changed solver and mesh datatype to Float32 but forgot to change equation datatype as well.

In that case I think the Ref(equations) part is causing the trouble, as we have for that scenario:

eltype(u) Float32
Ref(equations) Base.RefValue{IdealGlmMhdEquations2D{Float64}}(IdealGlmMhdEquations2D with 9 variables)

which then spawned an 18-dimensional array.

[ranocha]

The fix is to materialize

solution_variables.(reinterpret(SVector{nvariables(equations), eltype(u)}, u),
                    Ref(equations))))

first and use the eltype of its eltype for reinterpret. Right now, the code assumes that you do not change some types.

Either way, mixing difference precisions of the solver, the mesh, and the equations is not good.