Help modeling advection-diffusion equation with ModelingToolkit.jl

Specific Domains Modelling & Simulations
,
1

First time trying to use ModelingToolkit.jl in a project. Appreciate if someone can help with a simple 1D example of advection-diffusion:

R \frac{\partial c}{\partial t} = D \frac{\partial^2 c}{\partial x^2} - v \frac{\partial c}{\partial x}

with c(x,t) the concentration satisfying the initial condition c(x,0) = c_o.

I tried the following based on the documentation:

using ModelingToolkit
using DifferentialEquations

@variables x t c(x,t)

@parameters R=1.0 D=1.0 v=1.0 cₒ=1.0

∂t = Differential(t)
∂x = Differential(x)

# advection-diffusion equation
eqn = R * ∂t(c) ~ D * ∂x(c)^2 - v * ∂x(c)

# initial and boundary conditions
bcs = [c(x,0) ~ cₒ]

# space and time domains
domains = [x ∈ (0, 1), t ∈ (0, 1)]

# define the PDE system
@named sys = PDESystem(eqn, bcs, domains, [x, t], [c(x,t)])

# simplify the PDE system
simpsys = structural_simplify(sys)

but the code errors in the initial condition saying:

ERROR: Sym c(x, t) is not callable. Use @syms c(x, t)(var1, var2,...) to create it as a callable.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] (::SymbolicUtils.BasicSymbolic{Real})(::SymbolicUtils.BasicSymbolic{Real}, ::Vararg{Any})
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/ssQsQ/src/types.jl:867
 [3] (::Num)(::Num, ::Vararg{Any})
   @ Symbolics ~/.julia/packages/Symbolics/gBKZv/src/num.jl:19
 [4] top-level scope
   @ ~/Desktop/transport/main.jl:16

I tried using the macro @syms without success.

Can you please confirm that the macros @variables and @parameters are used correctly in the above example? I’ve encountered examples in the documentation that use @parameters for t and x and examples that use @variables for t and x.

2

You forgot the discretization steps? See:

3

The error occurs before the discretization steps, in the creation of the bcs variable. Probably related to the @parameters versus @variables usage?

4

But you never called discretize in your code.

5

Can you help with the missing steps? I don’t understand why we need to discretize. The code above doesn’t even define the boundary conditions successfully.

6

@variables c(..). This is all covered in the first tutorial that is linked.

7

BTW I’m curious, what tutorial or part of the docs were you starting from? You did something odd here, so I presume something in the docs misled you. It would be good to know what your reading source was that led you here.

8

@ChrisRackauckas unfortunately I cannot make the code work after couple attempts and variations. Can you please clarify the difference between @parameters and @variables?

Regarding the documentation, I can try to find the links again, but many examples switch between @parameters and @variables as I mentioned even when dealing with x and t which are space and time variables.

9

That shouldn’t matter, though @xtalax we should clean that up.

I’m curious though, can you show me a code where you do the discretization call? The code above has no discretize call at all.

Did you start from the heat equation codes or the getting started?

10

I had this error yesterday, replacing the variable declaration c(x,t) with c(..) should fix it

also in your eqn do you need to say c(x,t)?

11

Here is the updated MWE. Now the error occurs in the discretization step:

using ModelingToolkit
using DifferentialEquations
using MethodOfLines

@parameters x t
@variables c(..)

∂t = Differential(t)
∂x = Differential(x)

# advection-diffusion constants
@parameters R=1.0 D=1.0 v=1.0 cₒ=1.0

# advection-diffusion equation
eqn = R * ∂t(c(x, t)) ~ D * ∂x(c(x, t))^2 - v * ∂x(c(x, t))

# initial and boundary conditions
bcs = [c(x, 0) ~ cₒ]

# space and time domains
dom = [x ∈ (0, 1), t ∈ (0, 1)]

# define PDE system
@named sys = PDESystem(eqn, bcs, dom, [x, t], [c(x, t)])

# convert the PDE system into an ODE problem
prob = discretize(sys, MOLFiniteDifference([x => 0.1], t))

ERROR: MethodError: no method matching AbstractFloat(::Type{SymbolicUtils.BasicSymbolic{Real}})

Closest candidates are:
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number}
   @ Base char.jl:50
  (::Type{T})(::SymbolicUtils.Symbolic) where T<:Union{AbstractFloat, Integer, Complex{<:AbstractFloat}, Complex{<:Integer}}
   @ Symbolics ~/.julia/packages/Symbolics/gBKZv/src/Symbolics.jl:148
  (::Type{T})(::Base.TwicePrecision) where T<:Number
   @ Base twiceprecision.jl:266
  ...

Stacktrace:
  [1] float(x::Type)
    @ Base ./float.jl:294
  [2] promote_to_concrete(vs::Vector{SymbolicUtils.BasicSymbolic{Real}}; tofloat::Bool, use_union::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oIgbi/src/utils.jl:659
  [3] varmap_to_vars(varmap::Vector{Pair}, varlist::Vector{Any}; defaults::Dict{Any, Any}, check::Bool, toterm::Function, promotetoconcrete::Nothing, tofloat::Bool, use_union::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oIgbi/src/variables.jl:93
  [4] get_u0_p(sys::ODESystem, u0map::Vector{Pair}, parammap::SciMLBase.NullParameters; use_union::Bool, tofloat::Bool, symbolic_u0::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oIgbi/src/systems/diffeqs/abstractodesystem.jl:786
  [5] get_u0_p
    @ ~/.julia/packages/ModelingToolkit/oIgbi/src/systems/diffeqs/abstractodesystem.jl:768 [inlined]
  [6] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Pair}, parammap::SciMLBase.NullParameters; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, use_union::Bool, tofloat::Bool, symbolic_u0::Bool, u0_constructor::typeof(identity), kwargs::Base.Pairs{Symbol, Integer, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:t, :has_difference, :check_length), Tuple{Int64, Bool, Bool}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oIgbi/src/systems/diffeqs/abstractodesystem.jl:811
  [7] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Pair}, tspan::Tuple{Int64, Int64}, parammap::SciMLBase.NullParameters; callback::Nothing, check_length::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oIgbi/src/systems/diffeqs/abstractodesystem.jl:936
  [8] ODEProblem
    @ ~/.julia/packages/ModelingToolkit/oIgbi/src/systems/diffeqs/abstractodesystem.jl:929 [inlined]
  [9] (ODEProblem{true, SciMLBase.AutoSpecialize})(sys::ODESystem, u0map::Vector{Pair}, tspan::Tuple{Int64, Int64})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oIgbi/src/systems/diffeqs/abstractodesystem.jl:929
 [10] (ODEProblem{true})(::ODESystem, ::Vector{Pair}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oIgbi/src/systems/diffeqs/abstractodesystem.jl:916
 [11] (ODEProblem{true})(::ODESystem, ::Vector{Pair}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oIgbi/src/systems/diffeqs/abstractodesystem.jl:915
 [12] ODEProblem(::ODESystem, ::Vector{Pair}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oIgbi/src/systems/diffeqs/abstractodesystem.jl:912
 [13] ODEProblem(::ODESystem, ::Vector{Pair}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oIgbi/src/systems/diffeqs/abstractodesystem.jl:911
 [14] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ PDEBase ~/.julia/packages/PDEBase/aRCcl/src/discretization_state.jl:74
 [15] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ PDEBase ~/.julia/packages/PDEBase/aRCcl/src/discretization_state.jl:55
 [16] top-level scope
    @ ~/Desktop/transport/main.jl:27
12

I think PDESystem doesn’t automatically extract the default parameters from @parameters. You can put them into the PDESystem directly instead.

PDESystem(eqn, bcs, dom, [x, t], [c(x, t)], [R=>1.0, D=>1.0, v=>1.0, cₒ=>1.0])

As described here

The output of this should now list the parameters:

PDESystem
Equations: Equation[R*Differential(t)(c(x, t)) ~ -v*Differential(x)(c(x, t)) + D*(Differential(x)(c(x, t))^2)]
Boundary Conditions: Equation[c(x, 0) ~ cₒ]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(x, 0 .. 1), Symbolics.VarDomainPairing(t, 0.0 .. 1.0)]
Dependent Variables: Num[c(x, t)]
Independent Variables: Num[x, t]
Parameters: Pair{Num, Float64}[R => 1.0, D => 1.0, v => 1.0, cₒ => 1.0]
Default Parameter ValuesDict{Any, Any}()

whereas before it was

...
Parameters: SciMLBase.NullParameters()
Default Parameter ValuesDict{Any, Any}()

Not sure why the defaults are not recognized though…

13

@xtalax is that a known issue?

14

Updated the MWE. The example now fails to compile the problem:

using ModelingToolkit
using DifferentialEquations
using MethodOfLines

@variables x t
@variables c(..)

∂t = Differential(t)
∂x = Differential(x)

# advection-diffusion constants
R = 1.0
D = 1.0
v = 1.0
cₒ = 1.0

# advection-diffusion equation
eqn = R * ∂t(c(x, t)) ~ D * ∂x(c(x, t))^2 - v * ∂x(c(x, t))

# initial and boundary conditions
bcs = [c(x, 0) ~ cₒ]

# space and time domains
dom = [x ∈ (0, 1), t ∈ (0, 1)]

# define PDE system
@named sys = PDESystem(eqn, bcs, dom, [x, t], [c(x, t)], [R=>R, D=>D, v=>v, cₒ=>cₒ])

# convert the PDE into an ODE problem
prob = discretize(sys, MOLFiniteDifference([x => 0.1], t))

# solve the problem
solve(prob)

I also get a couple of warnings:

┌ Warning: : no method matching get_unit for arguments (Pair{Float64, Float64},).
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oIgbi/src/systems/validation.jl:154
15

It all boils down to this @parameters versus @variables macros and how they interact with the rest of the system. I really appreciate if you can clarify the meaning of each of these macros and when/how to use them.

If I completely remove the @parameters and treat R, D, v, co as normal Julia variables, the discretization works without warnings, but the solve call returns an error:

ERROR: MethodError: no method matching SciMLBase.PDETimeSeriesSolution(::SciMLBase.PDETimeSeriesSolution{Float64, 1, Dict{Num, Matrix{Float64}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ODESystem, Dict{Any, Any}, Vector{Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}}, CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, TRBDF2{1, false, LinearSolve.DefaultLinearSolver, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, TRBDF2{0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, OrdinaryDiffEq.CompositeInterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ODESystem, Dict{Any, Any}, Vector{Any}}, Nothing, ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.CompositeCache{Tuple{OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.TRBDF2Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.TRBDF2Tableau{Float64, Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ODESystem, Dict{Any, Any}, Vector{Any}}, Nothing, ODESystem}, Float64, Vector{Float64}}, FiniteDiff.JacobianCache{Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, UnitRange{Int64}, Nothing, Val{:forward}(), Float64}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.DefaultLinearSolver, LinearSolve.DefaultLinearSolverInit{LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, Tuple{LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, Vector{Int64}}, Nothing, Nothing, Nothing, LinearAlgebra.SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}, LinearAlgebra.Cholesky{Float64, Matrix{Float64}}, LinearAlgebra.Cholesky{Float64, Matrix{Float64}}, Tuple{LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int32}}, Base.RefValue{Int32}}, Tuple{LinearAlgebra.LU{Float64, Matrix{Float64}, Vector{Int64}}, Base.RefValue{Int64}}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, Bool}}}}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, TRBDF2{0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}}, DiffEqBase.Stats, Vector{Int64}}, Nothing, Vector{Float64}, Tuple{Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Vector{SymbolicUtils.BasicSymbolic{Real}}, Vector{Num}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#637#generated_observed#555"{Bool, ODESystem, Dict{Any, Any}, Vector{Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}}, CompositeAlgorithm{Tuple{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, TRBDF2{1, false, LinearSolve.DefaultLinearSolver, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}}, OrdinaryDiffEq.AutoSwitchCache{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, TRBDF2{0, false, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, Rational{Int64}, Int64}}, Dict{Num, Interpolations.GriddedInterpolation{Float64, 2, Matrix{Float64}, Interpolations.Gridded{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Vector{Float64}, Vector{Float64}}}}, DiffEqBase.Stats}, ::MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization})

Closest candidates are:
  SciMLBase.PDETimeSeriesSolution(::SciMLBase.AbstractODESolution{T}, ::MethodOfLines.MOLMetadata) where T
   @ MethodOfLines ~/.julia/packages/MethodOfLines/8J3AD/src/interface/solution/timedep.jl:2
16

Here is the working version:

using ModelingToolkit
using DifferentialEquations
using MethodOfLines

@variables x t
@variables c(..)

∂t = Differential(t)
∂x = Differential(x)

# advection-diffusion constants
R = 1.0
D = 1.0
v = 1.0
cₒ = 1.0

# advection-diffusion equation
eqn = R * ∂t(c(x, t)) ~ D * ∂x(c(x, t))^2 - v * ∂x(c(x, t))

# initial and boundary conditions
bcs = [c(x, 0) ~ cₒ]

# space and time domains
dom = [x ∈ (0, 1), t ∈ (0, 1)]

# define PDE system
@named sys = PDESystem(eqn, bcs, dom, [x, t], [c(x, t)])

# convert the PDE into an ODE problem
prob = discretize(sys, MOLFiniteDifference([x => 100], t))

# solve the problem
sol = solve(prob, Tsit5(), saveat=0.2)

Feedback:

  • The docs assume that the reader is familiar with @variables @parameters and related symbolic manipulations. I would add comments in the quickstart explaining the role of these macros, and a special page dedicated to advanced use cases with parameter initialization, etc.
  • The syntax Interval(0, 1) can be replaced with (0, 1) in the domains. That should be preferred in simple examples to avoid importing the ModelingToolkit.Interval.
  • It would be nice to expand on the post-procesing steps in the quickstart, explaining what one can do with the solution object. Plotting the solution is just one possible route, and the supported syntax with Plots.jl is not clear either.
  • Some of the error messages along the way were completely unrelated to the root of the issue. I know this issue is hard to tackle, but given the wide use of the framework, it is worth the effort.
  • The solve function didn’t work without an explicit solver. Perhaps add an observation somewhere in the docs that this automatic solver selection doesn’t work always and that it is good practice to pick a solver while debugging.

Now I will start to modify the parameters to see if the behavior matches the expected behavior. Thanks so far for the help.

2 Likes
17

That’s not quite true though. This seems to be an oversight that it’s not tested, and I’ll get that fixed up today, but generally it’s not a good idea to pick a solver while debugging (especially the one that was chosen).

@xtalax can we fix up the Getting Started to just use the default?

That should get an error throw.

@xtalax I’m going to do a docs pass on that. The docs are much further behind than they should be.

Can you point to which ones you ran into?

18

You mean that using t \in (0, 1) is unsupported? The code worked fine for me without the Interval constructor. It is a nice feature for those interested in “rectangular” domains made of Cartesian products of intervals.

I think some of them are documented here along the corresponding MWE.

19

Can you share some recommendation on a good solver for this family of equations?

20

Still on documentation…

Did you consider writing an intro book on the diffeq ecosystem to non-diffeq people? Which parts of the software already converged in terms of design?

From personal experience writing Geospatial Data Science with Julia, the book format is the most effective for communicating vision and software organization with a broader scientific community.