You didn’t supply initial conditions.
Are you sure this is the problem? As it seems, I don’t need to give the initial conditions explicity when defining the ODEProblem probDAE_MTK1 either. And this is in line with the following tutorial:
https://mtk.sciml.ai/stable/mtkitize_tutorials/modelingtoolkitize/
Nevertheless, I tried the following:
# Now try to structural_simplify the system before solving it
de2 = structural_simplify(de)
u0reduced = [u0[1:200]; u0[202:end]]
probDAE_MTK2 = ODEProblem(de2,u0reduced,(0.0,1500),p_noDualCache,jac=true,sparse=true)
sol_MTK2 = solve(probDAE_MTK2,QBDF())
This still gives the same error, which I find relieving, because defining u0reduced forced me to make assumptions about state / equation ordering in system de2 after structural_simplify, which can very well be wrong.
So, what is my problem?
Can you share what the error is?
I get the following error: UndefVarError: x_{201} not defined
Complete code to reproduce:
using DifferentialEquations
using DifferentialEquations.DiffEqBase: dualcache
using LinearAlgebra
using ModelingToolkit
using SparseArrays
using BenchmarkTools
# Generate the constants
const N = 100
dx = 0.5 / N
q = dualcache(zeros(N,1))
rR = dualcache(zeros(N,1))
ϱu = 0
dp = 100e2
Tin = 200
# generate index lists to find the solution variables easily in solution vector
const indices_Tgas = 1:N
const indices_Tsol = indices_Tgas .+ N
const indices_mflux = indices_Tsol[end] + 1
const indices_RConcentrationSolid = indices_mflux[end] .+ (1:N)
# indices = (indices_Tgas, indices_Tsol, indices_mflux, indices_RConcentrationSolid)
p_noDualCache = (dx, # control volume length
dp, # boundary condition: pressure differential across domain
Tin, # boundary condition: Gas inlet temperature
q, # caching variable for q
rR, # caching variable for rR
ϱu)
function DiscretizedSystemMassMatrixForm_noDualcache!(du,u,p,t)
# N,dx, dp, Tin, q, rR, ϱu, indices_Tgas, indices_Tsol, indices_mflux, indices_RConcentrationSolid = p
dx, dp, Tin, q, rR, ϱu = p
Tgas = @view u[indices_Tgas]
Tsol = @view u[indices_Tsol]
dTgas = @view du[indices_Tgas]
dTsol = @view du[indices_Tsol]
mflux = @view u[indices_mflux]
dmflux = @view du[indices_mflux]
RConcsol = @view u[indices_RConcentrationSolid]
dRConcsol = @view du[indices_RConcentrationSolid]
q = 1e4 * (Tgas - Tsol)
ϱu = mflux[1]
rR = 5.5e6 .* exp.(.-5680 ./(273 .+ Tsol)) .* (1 .- exp.(.-RConcsol./50))
# gas energy balance
@inbounds begin
dTgas[1] = ϱu * 1000 * (Tin - Tgas[1]) - q[1] * dx
end
@inbounds for k in 2:N
dTgas[k] = ϱu * 1000 * (Tgas[k-1] - Tgas[k]) - q[k] * dx
end
# solids energy balance
@inbounds for k in 1:N
dTsol[k] = q[k] / (1000 * 2000)
end
# gas momentum balance
@inbounds begin
dmflux[1] = - mflux[1] + 5e-3 * sqrt(dp)
end
# reactant concentration equations
@inbounds for k in 1:N
dRConcsol[k] = -rR[k]
end
end
# build initial conditions
u0 = zeros(N*3 + 1)
u0[indices_Tgas] .= Tin
u0[indices_Tsol] .= 20
u0[indices_mflux] = 0.5
u0[indices_RConcentrationSolid] .= 3000
# build mass matrix
Ii = [indices_Tsol; indices_RConcentrationSolid];
Jj = [indices_Tsol; indices_RConcentrationSolid];
V = ones(length(Ii));
Msparse = sparse(Ii,Jj,V, 3*N+1, 3*N+1)
f_DAE = ODEFunction(DiscretizedSystemMassMatrixForm_noDualcache!, mass_matrix=Msparse)
probDAE = ODEProblem(f_DAE,u0,(0.0,1500),p_noDualCache)
de = modelingtoolkitize(probDAE)
# MTKized system back to numerical problem:
probDAE_MTK1 = ODEProblem(de,[],(0.0,1500),p_noDualCache,jac=true,sparse=true)
# solve the numerical problem - works
sol_MTK1 = solve(probDAE_MTK1,QBDF())
# pretty much as fast as solving non-MTKized system with Jacobian sparsity, as in above posts.
@benchmark solve(probDAE_MTK1,QBDF())
# Now try to structural_simplify the system before solving it
de2 = structural_simplify(de)
# system de2 only has 200 states after structural simplification
probDAE_MTK2 = ODEProblem(de2,[],(0.0,1500),p_noDualCache,jac=true,sparse=true)
# Error thrown: UndefVarError: x_{201} not defined
sol_MTK2 = solve(probDAE_MTK2,QBDF())
If simplicity of syntax and ease of use are what you are after, check out MethodOfLines.jl, with which you can specify, discretize and solve PDE systems in close to purely mathematical syntax.
Compile times are quite high right now, but the solves themselves are comparable to manual semidiscretization, so still useful for parameter sweeps, ensembles and optimization problems.
Changes are coming which will improve compile times by a lot.
Thank you, I am more or less aware of the beauty the usage of MethodOfLines and ModelingToolkit offers, and I will definitely make use of them, when (if) I make the decision to use Julia in professional and productive settings. I am currently evaluating this.
For that reason, I wanted to get a feeling how much the performance of above codes can be pushed with some effort in coding.