Catalyst with units

I originally asked this on Slack, but I was looking for examples with Catalyst and units. It wasn't immediately clear how to do this, so improved documentation would be appreciated. Looking at some [tests](https://github.com/SciML/Catalyst.jl/blob/a058ce8f315af10b19c966fd4b84b11b20716c8e/test/miscellaneous_tests/units.jl) I arrived at: ``` using Catalyst using DynamicQuantities using OrdinaryDiffEqDefault ``` ```julia @parameters t [unit=u"s"] @species X(t) [unit=u"mol/m^3"] @parameters d [unit=u"1/s"] rxs = [ Reaction(d, [X], []) ] example = ReactionSystem(rxs, t, name=:example) |> complete ``` But now the question arises: How do I run this model. My goal here is that `u0` and `p` are also provided in units and that the solution also remembers the units (and the plotting and so forth) ``` tspan_example = (0.0, 14*24.0) .* u"s" p_example = [ :d => -1.0 * u"s^-1" ] u0_example = [ :X => 10 * u"mol/m^3" ] ode_example = ODEProblem(example, u0_example, tspan_example, p_example) ``` fails with: ``` Cannot create an additive identity from `Type{<:DynamicQuantities.Quantity}`, as the dimensions are unknown. Please use `zero(::DynamicQuantities.Quantity)` instead. ``` Without units ``` tspan_example = (0.0, 14*60) p_example = [ :d => -1.0 ] u0_example = [ :X => 10 ] ode_example = ODEProblem(example, u0_example, tspan_example, p_example) sol_example = solve(ode_example, saveat=0.1) ``` does work. ``` function ustrip(varmap) map(varmap) do (name, value) name => DynamicQuantities.ustrip( ModelingToolkit.get_unit(name), value) end end function from_units(rn, u, tspan, p) u = symmap_to_varmap(rn, u) p = symmap_to_varmap(rn, p) tspan = map(tspan) do t DynamicQuantities.ustrip(ModelingToolkit.get_unit(rn.t), t) end ustrip(u), tspan, ustrip(p) end ``` ``` ode_example = ODEProblem(example, from_units(example, u0_example, tspan_example, p_example)...) ``` Allows me to do 50% of what I want, but I still need to add the units manually back on to the solution, and of course ideally Catalyst would handle that for me. ``` julia> sol_example[:X] |> typeof Vector{Float64} (alias for Array{Float64, 1}) ```