Trixi.jl unable to convert elixir_advection_diffusion.jl 2d example into 1d example

Specific Domains
1

I am trying to create working 1d advection diffusion example with boundaries. I have first modified slightly code from Trixi.jl/examples/dgmulti_2d/elixir_advection_diffusion.jl at main · trixi-framework/Trixi.jl and confirmed that it is working correctly (code below). Next, I wanted to convert this 2d example into 1d example, but I get cryptic error message Cannot `convert` an object of type SVector{1, Float64} to an object of type Float64 that I cannot resolve.

1D failing example:

using OrdinaryDiffEq
using Trixi
using Trixi.SummationByPartsOperators

sbp_fd = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = 2, xmin = 0.0, xmax = 1.0, N = 11)

dg = DGMulti(polydeg = 3, element_type = Line(), approximation_type = sbp_fd,
             surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
             volume_integral = VolumeIntegralWeakForm())

equations = LinearScalarAdvectionEquation1D(1.5)
equations_parabolic = LaplaceDiffusion1D(5.0e-2, equations)

initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation1D) = SVector(0.0)
initial_condition = initial_condition_zero

# tag different boundary segments
left(x, tol = 50 * eps()) = abs(x[1] + 1) < tol
right(x, tol = 50 * eps()) = abs(x[1] - 1) < tol
is_on_boundary = Dict(:left => left, :right => right)

cells_per_dimension = (16,)
mesh = DGMultiMesh(dg, cells_per_dimension; is_on_boundary)

# BC types
boundary_condition_left = BoundaryConditionDirichlet((x, t, equations_parabolic) -> SVector(1 + 0.1 * x[1]))
boundary_condition_zero = BoundaryConditionDirichlet((x, t, equations_parabolic) -> SVector(0.0))
boundary_condition_neumann_zero = BoundaryConditionNeumann((x, t, equations_parabolic) -> SVector(0.0))

# define inviscid boundary conditions
boundary_conditions = (; :left => boundary_condition_left,
                       :right => boundary_condition_do_nothing)

# define viscous boundary conditions
boundary_conditions_parabolic = (; :left => boundary_condition_left,
                                 :right => boundary_condition_neumann_zero)

semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
                                             initial_condition, dg;
                                             boundary_conditions = (boundary_conditions,
                                                                    boundary_conditions_parabolic))

tspan = (0.0, 1.5)
ode = semidiscretize(semi, tspan)

time_int_tol = 1e-6
sol = solve(ode, Tsit5(); abstol = time_int_tol, reltol = time_int_tol,
            ode_default_options()..., callback = CallbackSet())

2D working example:

using OrdinaryDiffEq
using Trixi
using Trixi.SummationByPartsOperators

sbp_fd = derivative_operator(MattssonNordström2004(), derivative_order = 1, accuracy_order = 2, xmin = 0.0, xmax = 1.0, N = 11)

dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = sbp_fd,
             surface_integral = SurfaceIntegralWeakForm(flux_lax_friedrichs),
             volume_integral = VolumeIntegralWeakForm())

equations = LinearScalarAdvectionEquation2D(1.5, 1.0)
equations_parabolic = LaplaceDiffusion2D(5.0e-2, equations)

initial_condition_zero(x, t, equations::LinearScalarAdvectionEquation2D) = SVector(0.0)
initial_condition = initial_condition_zero

# tag different boundary segments
left(x, tol = 50 * eps()) = abs(x[1] + 1) < tol
right(x, tol = 50 * eps()) = abs(x[1] - 1) < tol
bottom(x, tol = 50 * eps()) = abs(x[2] + 1) < tol
top(x, tol = 50 * eps()) = abs(x[2] - 1) < tol
is_on_boundary = Dict(:left => left, :right => right, :top => top, :bottom => bottom)

cells_per_dimension = (16, 16)
mesh = DGMultiMesh(dg, cells_per_dimension; is_on_boundary)

# BC types
boundary_condition_left = BoundaryConditionDirichlet((x, t, equations_parabolic) -> SVector(1 + 0.1 * x[2]))
boundary_condition_zero = BoundaryConditionDirichlet((x, t, equations_parabolic) -> SVector(0.0))
boundary_condition_neumann_zero = BoundaryConditionNeumann((x, t, equations_parabolic) -> SVector(0.0))

# define inviscid boundary conditions
boundary_conditions = (; :left => boundary_condition_left,
                       :bottom => boundary_condition_zero,
                       :top => boundary_condition_do_nothing,
                       :right => boundary_condition_do_nothing)

# define viscous boundary conditions
boundary_conditions_parabolic = (; :left => boundary_condition_left,
                                 :bottom => boundary_condition_zero,
                                 :top => boundary_condition_zero,
                                 :right => boundary_condition_neumann_zero)

semi = SemidiscretizationHyperbolicParabolic(mesh, (equations, equations_parabolic),
                                             initial_condition, dg;
                                             boundary_conditions = (boundary_conditions,
                                                                    boundary_conditions_parabolic))

tspan = (0.0, 1.5)
ode = semidiscretize(semi, tspan)

time_int_tol = 1e-6
sol = solve(ode, Tsit5(); abstol = time_int_tol, reltol = time_int_tol,
            ode_default_options()..., callback = CallbackSet())

Error message for 1D example:

ERROR: MethodError: Cannot `convert` an object of type SVector{1, Float64} to an object of type Float64
The function `convert` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  convert(::Type{T}, ::VectorizationBase.AbstractSIMD) where T<:Union{Bool, Float16, Float32, Float64, Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8, SIMDTypes.Bit}
   @ VectorizationBase C:\Users\user\.julia\packages\VectorizationBase\7mwzi\src\base_defs.jl:201
  convert(::Type{T}, ::Union{Static.StaticBool{N}, Static.StaticFloat64{N}, Static.StaticInt{N}} where N) where T<:Number 
   @ Static C:\Users\user\.julia\packages\Static\YQhfT\src\Static.jl:434
  convert(::Type{T}, ::VectorizationBase.LazyMulAdd{M, O, I}) where {M, O, I, T<:Number}
   @ VectorizationBase C:\Users\user\.julia\packages\VectorizationBase\7mwzi\src\lazymul.jl:25
  ...

Stacktrace:
  [1] macro expansion
    @ C:\Users\user\.julia\packages\StaticArraysCore\7xxEJ\src\StaticArraysCore.jl:88 [inlined]
  [2] convert_ntuple
    @ C:\Users\user\.julia\packages\StaticArraysCore\7xxEJ\src\StaticArraysCore.jl:84 [inlined]
  [3] SVector{1, Float64}(x::Tuple{SVector{1, Float64}})
    @ StaticArraysCore C:\Users\user\.julia\packages\StaticArraysCore\7xxEJ\src\StaticArraysCore.jl:120
  [4] convert
    @ C:\Users\user\.julia\packages\StaticArrays\DsPgf\src\convert.jl:189 [inlined]
  [5] maybe_convert_elt(::Type{SVector{1, Float64}}, vals::SVector{1, SVector{1, Float64}})
    @ StructArrays C:\Users\user\.julia\packages\StructArrays\bgm9N\src\utils.jl:192
  [6] setindex!
    @ C:\Users\user\.julia\packages\StructArrays\bgm9N\src\structarray.jl:378 [inlined]
  [7] _setindex!
    @ .\abstractarray.jl:1466 [inlined]
  [8] setindex!
    @ .\abstractarray.jl:1443 [inlined]
  [9] macro expansion
    @ C:\Users\user\.julia\packages\Trixi\ZwFHs\src\solvers\dgmulti\dg_parabolic.jl:325 [inlined]
 [10] #calc_viscous_fluxes!##6
    @ C:\Users\user\.julia\packages\Polyester\almvr\src\closure.jl:309 [inlined]
 [11] macro expansion
    @ C:\Users\user\.julia\packages\Polyester\almvr\src\batch.jl:246 [inlined]
 [12] _batch_no_reserve
    @ C:\Users\user\.julia\packages\Polyester\almvr\src\batch.jl:168 [inlined]
 [13] batch
    @ C:\Users\user\.julia\packages\Polyester\almvr\src\batch.jl:334 [inlined]
 [14] macro expansion
    @ C:\Users\user\.julia\packages\Polyester\almvr\src\closure.jl:456 [inlined]
 [15] macro expansion
    @ C:\Users\user\.julia\packages\Trixi\ZwFHs\src\auxiliary\auxiliary.jl:213 [inlined]
 [16] calc_viscous_fluxes!(flux_viscous::SVector{…}, u::StructArray{…}, gradients::SVector{…}, mesh::DGMultiMesh{…}, equations::LaplaceDiffusion1D{…}, dg::DGMulti{…}, cache::@NamedTuple{…}, cache_parabolic::@NamedTuple{…})
    @ Trixi C:\Users\user\.julia\packages\Trixi\ZwFHs\src\solvers\dgmulti\dg_parabolic.jl:311
 [17] macro expansion
    @ C:\Users\user\.julia\packages\Trixi\ZwFHs\src\solvers\dgmulti\dg_parabolic.jl:491 [inlined]
 [18] macro expansion
    @ C:\Users\user\.julia\packages\TrixiBase\38UFo\src\trixi_timeit.jl:64 [inlined]
 [19] rhs_parabolic!(du::StructArray{…}, u::StructArray{…}, t::Float64, mesh::DGMultiMesh{…}, equations_parabolic::LaplaceDiffusion1D{…}, boundary_conditions::@NamedTuple{…}, source_terms::Nothing, dg::DGMulti{…}, parabolic_scheme::ViscousFormulationBassiRebay1, cache::@NamedTuple{…}, cache_parabolic::@NamedTuple{…})
    @ Trixi C:\Users\user\.julia\packages\Trixi\ZwFHs\src\solvers\dgmulti\dg_parabolic.jl:490
 [20] macro expansion
    @ C:\Users\user\.julia\packages\TrixiBase\38UFo\src\trixi_timeit.jl:64 [inlined]
 [21] rhs_parabolic!(du_ode::RecursiveArrayTools.VectorOfArray{…}, u_ode::RecursiveArrayTools.VectorOfArray{…}, semi::SemidiscretizationHyperbolicParabolic{…}, t::Float64)
    @ Trixi C:\Users\user\.julia\packages\Trixi\ZwFHs\src\semidiscretization\semidiscretization_hyperbolic_parabolic.jl:352
 [22] ODEFunction
    @ C:\Users\user\.julia\packages\SciMLBase\b4Q81\src\scimlfunctions.jl:2597 [inlined]
 [23] (::SplitFunction{…})(du::RecursiveArrayTools.VectorOfArray{…}, u::RecursiveArrayTools.VectorOfArray{…}, p::SemidiscretizationHyperbolicParabolic{…}, t::Float64)
    @ SciMLBase C:\Users\user\.julia\packages\SciMLBase\b4Q81\src\scimlfunctions.jl:2627
 [24] initialize!
    @ C:\Users\user\.julia\packages\OrdinaryDiffEqTsit5\o07IY\src\tsit_perform_step.jl:175 [inlined]
 [25] __init(prob::ODEProblem{…}, alg::Tsit5{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_discretes::Bool, save_start::Bool, save_end::Nothing, callback::CallbackSet{…}, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Float64, reltol::Float64, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias::ODEAliasSpecifier, initializealg::DiffEqBase.DefaultInit, kwargs::@Kwargs{})    
    @ OrdinaryDiffEqCore C:\Users\user\.julia\packages\OrdinaryDiffEqCore\GMkz9\src\solve.jl:578
 [26] __init (repeats 2 times)
    @ C:\Users\user\.julia\packages\OrdinaryDiffEqCore\GMkz9\src\solve.jl:11 [inlined]
 [27] __solve(::ODEProblem{…}, ::Tsit5{…}; kwargs::@Kwargs{…})
    @ OrdinaryDiffEqCore C:\Users\user\.julia\packages\OrdinaryDiffEqCore\GMkz9\src\solve.jl:6
 [28] __solve
    @ C:\Users\user\.julia\packages\OrdinaryDiffEqCore\GMkz9\src\solve.jl:1 [inlined]
 [29] #solve_call#23
    @ C:\Users\user\.julia\packages\DiffEqBase\p82Yh\src\solve.jl:127 [inlined]
 [30] solve_call
    @ C:\Users\user\.julia\packages\DiffEqBase\p82Yh\src\solve.jl:84 [inlined]
 [31] #solve_up#30
    @ C:\Users\user\.julia\packages\DiffEqBase\p82Yh\src\solve.jl:563 [inlined]
 [32] solve_up
    @ C:\Users\user\.julia\packages\DiffEqBase\p82Yh\src\solve.jl:540 [inlined]
 [33] #solve#29
    @ C:\Users\user\.julia\packages\DiffEqBase\p82Yh\src\solve.jl:530 [inlined]
 [34] top-level scope
    @ c:\Code\MWE\Trixi_SBP_2Dto1D\example1d.jl:49
Some type information was truncated. Use `show(err)` to see complete types.
2

Yeah that is a bug in the code.

I outline the fix in the corresponding Trixi PR, hopefully it will get merged soon.

For you

function Trixi.flux(u, gradients, orientation::Integer, equations_parabolic::LaplaceDiffusion1D)
    dudx = gradients[1] # Extract first (and only) component from gradients
    # orientation == 1
    return equations_parabolic.diffusivity * dudx
end

Should solve the issue.

1 Like
3

Btw: Maybe open for stuff like this directly an issue at Trixi.jl, odds are much higher you get a quick response over there.

4

The PR created by @DanDoe has been merged and a new version of Trixi.jl including it should be released soon (see New version: Trixi v0.13.11 by JuliaRegistrator · Pull Request #140903 · JuliaRegistries/General · GitHub).