Allow different numbers of variables in `PlotData2D`

[ranocha]

With StructArrays.jl v0.6.8:

julia> using Trixi, OrdinaryDiffEq, Plots

julia> trixi_include("examples/structured_2d_dgsem/elixir_advection_extended.jl")

julia> PlotData2D(sol) # fine
PlotData2DTriangulated{Matrix{Float64}, StructArray{SVector{1, Float64}, 2, Tuple{Matrix{Float64}}, Int64}, Matrix{Float64}, StructArray{SVector{1, Float64}, 2, Tuple{Matrix{Float64}}, Int64}, SVector{1, String}}(<x>, <y>, <data>, <plot_triangulation>, <x_face>, <y_face>, <face_data>, <variable_names>)

julia> cons2strange(u, ::LinearScalarAdvectionEquation2D) = SVector(u[1], u[1])
cons2strange (generic function with 1 method)

julia> Trixi.varnames(::typeof(cons2strange), ::LinearScalarAdvectionEquation2D) = ("a", "b")

julia> PlotData2D(sol; solution_variables=cons2strange) # fails
ERROR: MethodError: Cannot `convert` an object of type Tuple{Float64, Float64} to an object of type Float64
Closest candidates are:
  convert(::Type{T}, ::Static.StaticFloat64) where T<:AbstractFloat at ~/.julia/packages/Static/NDBmr/src/float.jl:22
  convert(::Type{T}, ::LLVM.GenericValue, ::LLVM.LLVMType) where T<:AbstractFloat at ~/.julia/packages/LLVM/szqwr/src/execution.jl:39
  convert(::Type{T}, ::LLVM.ConstantFP) where T<:AbstractFloat at ~/.julia/packages/LLVM/szqwr/src/core/value/constant.jl:111
  ...
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/StaticArrays/58yy1/src/util.jl:17 [inlined]
  [2] convert_ntuple
    @ ~/.julia/packages/StaticArrays/58yy1/src/util.jl:13 [inlined]
  [3] SVector{1, Float64}(x::Tuple{Tuple{Float64, Float64}})
    @ StaticArrays ~/.julia/packages/StaticArrays/58yy1/src/SArray.jl:28
  [4] StaticArray
    @ ~/.julia/packages/StaticArrays/58yy1/src/convert.jl:4 [inlined]
  [5] convert
    @ ~/.julia/packages/StaticArrays/58yy1/src/convert.jl:10 [inlined]
  [6] maybe_convert_elt(#unused#::Type{SVector{1, Float64}}, vals::SVector{2, Float64})
    @ StructArrays ~/.julia/packages/StructArrays/rICDm/src/utils.jl:197
  [7] setindex!
    @ ~/.julia/packages/StructArrays/rICDm/src/structarray.jl:363 [inlined]
  [8] transform_to_solution_variables!(u::StructArray{SVector{1, Float64}, 2, Tuple{Matrix{Float64}}, Int64}, solution_variables::typeof(cons2strange), equations::LinearScalarAdvectionEquation2D{Float64})
    @ Trixi ~/.julia/dev/Trixi/src/visualization/utilities.jl:53
  [9] Trixi.PlotData2DTriangulated(u::Array{Float64, 4}, mesh::StructuredMesh{2, Float64}, equations::LinearScalarAdvectionEquation2D{Float64}, dg::DGSEM{LobattoLegendreBasis{Float64, 4, SVector{4, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, Trixi.LobattoLegendreMortarL2{Float64, 4, Matrix{Float64}, Matrix{Float64}}, SurfaceIntegralWeakForm{FluxLaxFriedrichs{typeof(max_abs_speed_naive)}}, VolumeIntegralWeakForm}, cache::NamedTuple{(:elements,), Tuple{Trixi.ElementContainer{2, Float64, Float64, 3, 4, 5}}}; solution_variables::Function, nvisnodes::Int64)
    @ Trixi ~/.julia/dev/Trixi/src/visualization/types.jl:398
 [10] #PlotData2D#1103
    @ ~/.julia/dev/Trixi/src/visualization/types.jl:242 [inlined]
 [11] PlotData2D(u_ode::Vector{Float64}, semi::SemidiscretizationHyperbolic{StructuredMesh{2, Float64}, LinearScalarAdvectionEquation2D{Float64}, typeof(initial_condition_convergence_test), Trixi.BoundaryConditionPeriodic, Nothing, DGSEM{LobattoLegendreBasis{Float64, 4, SVector{4, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, Trixi.LobattoLegendreMortarL2{Float64, 4, Matrix{Float64}, Matrix{Float64}}, SurfaceIntegralWeakForm{FluxLaxFriedrichs{typeof(max_abs_speed_naive)}}, VolumeIntegralWeakForm}, NamedTuple{(:elements,), Tuple{Trixi.ElementContainer{2, Float64, Float64, 3, 4, 5}}}}; kwargs::Base.Pairs{Symbol, typeof(cons2strange), Tuple{Symbol}, NamedTuple{(:solution_variables,), Tuple{typeof(cons2strange)}}})
    @ Trixi ~/.julia/dev/Trixi/src/visualization/types.jl:229
 [12] #PlotData2D#1105
    @ ~/.julia/dev/Trixi/src/visualization/types.jl:285 [inlined]
 [13] top-level scope
    @ REPL[7]:1

With StructArrays.jl v0.6.7:

julia> using Trixi, OrdinaryDiffEq, Plots

julia> trixi_include("examples/structured_2d_dgsem/elixir_advection_extended.jl")

julia> PlotData2D(sol) # fine
PlotData2DTriangulated{Matrix{Float64}, StructArray{SVector{1, Float64}, 2, Tuple{Matrix{Float64}}, Int64}, Matrix{Float64}, StructArray{SVector{1, Float64}, 2, Tuple{Matrix{Float64}}, Int64}, SVector{1, String}}(<x>, <y>, <data>, <plot_triangulation>, <x_face>, <y_face>, <face_data>, <variable_names>)

julia> cons2strange(u, ::LinearScalarAdvectionEquation2D) = SVector(u[1], u[1])
cons2strange (generic function with 1 method)

julia> Trixi.varnames(::typeof(cons2strange), ::LinearScalarAdvectionEquation2D) = ("a", "b")

julia> pd = PlotData2D(sol; solution_variables=cons2strange) # works now
PlotData2DTriangulated{Matrix{Float64}, StructArray{SVector{1, Float64}, 2, Tuple{Matrix{Float64}}, Int64}, Matrix{Float64}, StructArray{SVector{1, Float64}, 2, Tuple{Matrix{Float64}}, Int64}, SVector{2, String}}(<x>, <y>, <data>, <plot_triangulation>, <x_face>, <y_face>, <face_data>, <variable_names>)

julia> plot(pd["a"]) # fine, as in our tutorial

This let's our CI of tutorial 8 fail right now, e.g. https://github.com/trixi-framework/Trixi.jl/pull/1140

CC @jlchan

[jlchan]

Hm...this is probably due to some default "convert on setindex" behavior introduced by Tim Holy (see https://github.com/JuliaArrays/StructArrays.jl/pull/227 and https://github.com/JuliaArrays/StructArrays.jl/issues/216 for an earlier discussion).

My guess is that PlotData2D creates a plotting array using similar(u) and doesn't take into account the transformation. I'll try to get to it later today.

[ranocha]

The issue is that we assume that solution_variables will have the same number of variables as the conserved variables when creating a PlotData2D object.

• Is this reasonable? If so, we should document it.
• Can we generalize it to allow another number of variables after the transformation by solution_variables? If we do that, we should update the part changed in #1143.

[jlchan]

We should be able to allow any number of variables (just call solution_variables on a single element of u to get the number of outputs and initialize uplot accordingly). We'd probably also need to change the construction of ufp (the plotting solution on mesh faces).

[jlchan]

IMO it's easier to assume solution_variables has the same number of variables, and to use something like ScalarPlotData2D for other fields, but I'm open to both options.

[ranocha]

The problem I see with ScalarPlotData2D right now is that you need to know the memory layout to generate the scalar data, which is kind of hacky for non-DGMulti solvers... Of course, this could also be seen as an issue with the non-DGMulti solvers :sweat_smile:

[ranocha]

Let's turn this issue into a note on a possibly interesting new feature for our visualization pipeline: It would be nice to allow different numbers of variables when processing the conserved variables with solution_variables in PlotData2D and friends.

[jlchan]

Sounds good to me