"ExtraVariablesSystemException: The system is unbalanced" Error while solving a PDE!

[aelmokadem]

Hi,

I am trying to solve this PDE system and getting the error below. Any ideas what I am missing here?!

Thank you.

Modeling code:

using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets

@parameters t r 
Dt = Differential(t)
Dr = Differential(r)
Drr = Differential(r)^2

Rcap = 8.0
Rkrogh = 72.0
Rec=1e5 
Av=6.023e23 
D=1.3e-7 
P=2.8e-7 
Γ=1.0 
kon=6e-3 
KD=0.1 
ϵ=0.24 
Rs=1.0 
kinT=0.01 
kinB=0.01 
keT=0.03 
keB=0.03 
feT=0.5 
feB=0.5 
Rcell=10.0 
A=60.0 
B=40.0 
ka=0.6*60 
kb=0.17*24*60 
Cplasma0=1*0.02/2

# derived quantities
Vcell = (4.0/3.0)*pi*(Rcell^3)
koff = KD*kon
Cplasma = Cplasma0 * ((A*exp(-ka*t)) + B*exp(-kb*t))

# variables
@variables u(..) T(t)=Rec/(Av * Vcell) DT(t)=0.0 Te(t)=0.0 DTe(t)=0.0 

Rxn = kon*(u(t, r)/ϵ)*((T/ϵ)*Γ) - koff*DT
Cplasma = Cplasma0 * ((A*exp(-ka*t)) + B*exp(-kb*t))

eqs  = [
    Dt(u(t, r)) ~ D*(1/r*Dr(u(t, r)) + Drr(u(t, r))) - Rxn,
    Dt(T) ~ Rs - Rxn - kinT*T + keT*feT*Te,
    Dt(DT) ~ Rxn - kinB*DT + keB*feB*DTe,
    Dt(Te) ~ kinT*T - keT*Te,
    Dt(DTe) ~ kinB*DT - keB*DTe
]

# initial and boundary conditions
# Space and time domains
domains = [t ∈ Interval(0.0, 50.0*24.0*60.0),
           r ∈ Interval(Rcap, Rkrogh)]

bcs = [u(0.0, r) ~ 0.0,
       D*Dr(u(t, Rcap)) ~ P*(Cplasma - (u(t, Rcap)/ϵ)),
       Dr(u(t, Rkrogh)) ~ 0.0]

# define PDE system
@named pdesys = PDESystem(eqs, bcs, domains, [t, r], [u(t, r), T, DT, Te, DTe])

# Method of lines discretization
dr = 8.0
discretization = MOLFiniteDifference([r => dr], t)

# Convert the PDE problem into an ODE problem and solve
prob = discretize(pdesys, discretization)
sol = solve(prob, Tsit5())

Error message:

The system of equations is:
Equation[-1.3e-7(0.015625(u(t))[1] + 0.0625(-0.125(u(t))[1] + 0.125(u(t))[2]) - 0.03125(u(t))[2] + 0.0156
25(u(t))[3]) - 0.0006000000000000001DT(t) + Differential(t)((u(t))[2]) + 0.10416666666666669(u(t))[2]*T(t) ~ 0, -1.3e-7(0.041666666666666664(-0.125(u(t))[2] + 0.125(u(t))[3]) + 0.015625(u(t))[2] - 0.03125(u(t))[3] + 0.015625(u(t))[4]) - 0.0006000000000000001DT(t) + Differential(t)((u(t))[3]) + 0.10416666666666669(u(t))[3]*T(t) ~ 0, -1.3e-7(0.015625(u(t))[3] + 0.03125(-0.125(u(t))[3] + 0.125(u(t))[4]) - 0.03125(u(t))[4] + 0.015625(u(t))[5]) - 0.0006000000000000001DT(t) + Differential(t)((u(t))[4]) + 0.10416666666666669(u(t))[4]*T(t) ~ 0, -1.3e-7(0.025(-0.125(u(t))[4] + 0.125(u(t))[5]) + 0.015625(u(t))[4] - 0.03125(u(t))[5] + 0.015625(u(t))[6]) + Differential(t)((u(t))[5]) - 0.0006000000000000001DT(t) + 0.10416666666666669(u(t))[5]*T(t) ~ 0, -1.3e-7(0.015625(u(t))[5] + 0.020833333333333332(-0.125(u(t))[5] + 0.125(u(t))[6]) - 0.03125(u(t))[6] + 0.015625(u(t))[7]) - 0.0006000000000000001DT(t) + Differential(t)((u(t))[6]) + 0.10416666666666669(u(t))[6]*T(t) ~ 0, -1.3e-7(0.017857142857142856(-0.125(u(t))[6] + 0.125(u(t))[7]) + 0.015625(u(t))[6] - 0.03125(u(t))[7] + 0.015625(u(t))[8]) - 0.0006000000000000001DT(t) + Differential(t)((u(t))[7]) + 0.10416666666666669(u(t))[7]*T(t) ~ 0, -1.3e-7(0.015625(-0.125(u(t))[7] + 0.125(u(t))[8]) + 0.015625(u(t))[7] - 0.03125(u(t))[8] + 0.015625(u(t))[9]) - 0.0006000000000000001DT(t) + Differential(t)((u(t))[8]) + 0.10416666666666669(u(t))[8]*T(t) ~ 0, -1.0 - 0.0006000000000000001DT(t) + 0.01T(t) - 0.015Te(t) + Differential(t)(T(t)) + 0.10416666666666669T(t)*u(t, r) ~ 0, 0.0106DT(t) + Differential(t)(DT(t)) - 0.015DTe(t) - 0.10416666666666669T(t)*u(t, r) ~ 0, -0.01T(t) + 0.03Te(t) + Differential(t)(Te(t)) ~ 0, Differential(t)(DTe(t)) - 0.01DT(t) + 0.03DTe(t) ~ 0, 1.3e-7(-0.1875(u(t))[1] + 0.25(u(t))[2] - 0.0625(u(t))[3]) ~ 2.8e-7(-4.166666666666667(u(t))[1] + 0.01(40.0exp(-244.8t) + 60.0exp(-36.0t))), 0.0625(u(t))[7] - 0.25(u(t))[8] + 0.1875(u(t))[9] ~ 0.0]                                                                             
There are 13 variables and 13 equations.

There are 11 time derivatives.

The variables without time derivatives are:
SymbolicUtils.BasicSymbolic{Real}[(u(t))[1], (u(t))[9]]

The equations without time derivatives are:
Equation[1.3e-7(-0.1875(u(t))[1] + 0.25(u(t))[2] - 0.0625(u(t))[3]) ~ 2.8e-7(-4.166666666666667(u(t))[1] 
+ 0.01(40.0exp(-244.8t) + 60.0exp(-36.0t))), 0.0625(u(t))[7] - 0.25(u(t))[8] + 0.1875(u(t))[9] ~ 0.0]    ERROR: ExtraVariablesSystemException: The system is unbalanced. There are 14 highest order derivative var
iables and 13 equations.                                                                                 More variables than equations, here are the potential extra variable(s):
 (u(t))[2]
 (u(t))[3]
 (u(t))[4]
 (u(t))[5]
 (u(t))[6]
 (u(t))[7]
 (u(t))[8]
 T(t)
 DT(t)
 Te(t)
 DTe(t)
 u(t, r)
Stacktrace:
  [1] error_reporting(state::TearingState{…}, bad_idxs::Vector{…}, n_highest_vars::Int64, iseqs::Bool, or
ig_inputs::Set{…})                                                                                           @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/dCa81/src/structural_tr
ansformation/utils.jl:50                                                                                   [2] check_consistency(state::TearingState{ODESystem}, orig_inputs::Set{Any})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/dCa81/src/structural_tr
ansformation/utils.jl:85                                                                                   [3] _structural_simplify!(state::TearingState{…}, io::Nothing; simplify::Bool, check_consistency::Bool,
 fully_determined::Bool, kwargs::@Kwargs{})                                                                  @ ModelingToolkit ~/.julia/packages/ModelingToolkit/dCa81/src/systems/systemstructure.jl:609
  [4] _structural_simplify!
    @ ~/.julia/packages/ModelingToolkit/dCa81/src/systems/systemstructure.jl:597 [inlined]
  [5] structural_simplify!(state::TearingState{…}, io::Nothing; simplify::Bool, check_consistency::Bool, 
fully_determined::Bool, kwargs::@Kwargs{})                                                                   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/dCa81/src/systems/systemstructure.jl:563
  [6] __structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/dCa81/src/systems/systems.jl:56
  [7] __structural_simplify
    @ ~/.julia/packages/ModelingToolkit/dCa81/src/systems/systems.jl:36 [inlined]
  [8] structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, kwargs::@Kwargs{})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/dCa81/src/systems/systems.jl:21
  [9] structural_simplify
    @ ~/.julia/packages/ModelingToolkit/dCa81/src/systems/systems.jl:19 [inlined]
 [10] structural_simplify(sys::ODESystem)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/dCa81/src/systems/systems.jl:19
 [11] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, kwargs::@K
wargs{})                                                                                                     @ PDEBase ~/.julia/packages/PDEBase/Aqj4G/src/discretization_state.jl:60
 [12] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
    @ PDEBase ~/.julia/packages/PDEBase/Aqj4G/src/discretization_state.jl:55
 [13] top-level scope
    @ REPL[93]:1
Some type information was truncated. Use `show(err)` to see complete types.

[ChrisRackauckas]

I think this is a bug in the parsing because u(t, r) shouldn't end up in the final states. It's not entirely clear why that's happening though.

[aelmokadem]

Thanks @ChrisRackauckas for the reply. So, is there something that I can do to help sort it out?

[ChrisRackauckas]

If you're willing to dive into the parser

[aelmokadem]

@ChrisRackauckas I think that would be beyond me. I am just curious about what in this particular problem that is causing this and if there's something I can change in the problem to avoid the issue! I have other examples that work fine!