Open
Description
I am constructing an ODEProblem with u0 defined as dictionary, see function build_problem in below code. When doing so, solving the constructed problem errors with no method matching oneunit(::Type{Any}). If instead I build u0 as a plain, boring vector of Floats, solution of system works.
using OrdinaryDiffEq, ModelingToolkit
using DifferentialEquations
function build_system(n)
@parameters t
vars = @variables (Tgas[1:n](t) = fill(100,n)),
(Tsol[1:n](t) = fill(20, n)),
(cR[1:n](t) = fill(3000, n)),
(mflux(t) = 0.5)
@parameters params[1:3]
D = Differential(t)
dp, Tin, h = params
dx = h / n
############################################################
# Build equations ##########################################
############################################################
eqns = Array{Equation}(undef, 3 * n + 1)
count = 1
# specification as function or intermediate variabel does not matter. I could also make this an
# actual @variable, then I could retrieve it later.
rR = 5.5e6 .* exp.(.-5680 ./(273 .+ Tsol)) .* (1 .- exp.(.-cR./50))
q = 1e4 * (Tgas - Tsol)
# gas energy balance
eqns[count] = 0 ~ mflux * 1000 * (Tin - Tgas[1]) - q[1] * dx
count += 1
for i in 2:n
#eqns[count] = 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q(Tsol[i], Tgas[i]) * dx
eqns[count] = 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q[i] * dx
count += 1
end
# solids energy balance
for i in 1:n
#eqns[count] = D(Tsol[i]) ~ q(Tsol[i], Tgas[i]) / (1000 * 2000)
eqns[count] = D(Tsol[i]) ~ q[i] / (1000 * 2000)
count += 1
end
# gas momentum balance
eqns[count] = 0 ~ - mflux[1] + 5e-3 * sqrt(dp)
count += 1
# reactant concentration equations
for i in 1:n
#eqns[count] = D(cR[i]) ~ -rR(Tsol[i], cR[i])
eqns[count] = D(cR[i]) ~ -rR[i]
count += 1
end
@named sys = ODESystem(eqns, t, vcat(vars...), params)
return sys
end
function build_problem(sys; tspan = (0.0, 1.5), param_values = [100e2, 200, 0.5])
u0map = [
(sys.Tgas .=> fill(100.0, N))...
(sys.Tsol .=> fill(20.0, N))...
(sys.cR .=> fill(3000.0, N))...
(sys.mflux => 0.5)
]
prob = ODEProblem(sys, u0map, tspan, param_values, jac=true, sparse=true)
return prob
end
dp = 100e2
Tin = 200
h = 0.5
parameters = [dp, Tin, h]
N = 100
tspan = (0.0, 1500.0)
sys = build_system(N)
prob = build_problem(sys, tspan=tspan, param_values = parameters)
# error: no method matching oneunit(::Type{Any})
# error does not happen if a plain vector is used as u0 in ODEProblem construction.
sol = solve(prob, QBDF(), reltol = 1e-3)