# Help modeling advection-diffusion equation with ModelingToolkit.jl

**URL:** <https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, modelingtoolkit\
**Created:** [October 16, 2023, 5:07pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028 "2023-10-16T17:07:30Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 16, 2023, 5:07pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/1 "2023-10-16T17:07:30Z")

</div>

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:

```julia
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:

```julia
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`.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 16, 2023, 5:56pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/2 "2023-10-16T17:56:20Z")

</div>

You forgot the discretization steps? See:

> **[Getting Started · MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/tutorials/brusselator/)**
>
> Documentation for MethodOfLines.jl.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 16, 2023, 6:09pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/3 "2023-10-16T18:09:11Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 16, 2023, 6:40pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/4 "2023-10-16T18:40:57Z")

</div>

But you never called discretize in your code.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 16, 2023, 6:51pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/5 "2023-10-16T18:51:29Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 16, 2023, 6:58pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/6 "2023-10-16T18:58:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 16, 2023, 7:38pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/7 "2023-10-16T19:38:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 16, 2023, 8:10pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/8 "2023-10-16T20:10:14Z")

</div>

@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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 16, 2023, 8:23pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/9 "2023-10-16T20:23:30Z")

</div>

> [@juliohm](#):
>
> 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.

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?

> **[Solving the Heat Equation · MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/tutorials/heat/)**
>
> Documentation for MethodOfLines.jl.

---

<div class="post-metadata">

**Author:** ![p\_f](https://avatars.discourse-cdn.com/v4/letter/p/45deac/32.png) [@p\_f](https://discourse.julialang.org/u/p_f)\
**Post date:** [October 16, 2023, 8:41pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/10 "2023-10-16T20:41:16Z")

</div>

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)?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 16, 2023, 11:28pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/11 "2023-10-16T23:28:38Z")

</div>

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

```julia
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))

```

```julia
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

```

---

<div class="post-metadata">

**Author:** ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@Sevi](https://discourse.julialang.org/u/Sevi)\
**Post date:** [October 17, 2023, 4:48am UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/12 "2023-10-17T04:48:09Z")

</div>

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

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

```

As described here

> **[Adding parameters · MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/tutorials/params/)**
>
> Documentation for MethodOfLines.jl.

The output of this should now list the parameters:

```julia
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

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

```

Not sure why the defaults are not recognized though…

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 17, 2023, 7:20am UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/13 "2023-10-17T07:20:37Z")

</div>

> [@Sevi](#):
>
> Not sure why the defaults are not recognized though…

@xtalax is that a known issue?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 17, 2023, 1:13pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/14 "2023-10-17T13:13:57Z")

</div>

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

```julia
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:

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

```

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 17, 2023, 5:08pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/15 "2023-10-17T17:08:48Z")

</div>

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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 17, 2023, 5:27pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/16 "2023-10-17T17:27:35Z")

</div>

Here is the working version:

```julia
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 17, 2023, 6:43pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/17 "2023-10-17T18:43:23Z")

</div>

> [@juliohm](#):
>
> 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.

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?

> [@juliohm](#):
>
> 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`.

That should get an error throw.

> [@juliohm](#):
>
> 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.

> [@juliohm](#):
>
> 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.

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

> [@juliohm](#):
>
> 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.

Can you point to which ones you ran into?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 17, 2023, 7:08pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/18 "2023-10-17T19:08:47Z")

</div>

> [@ChrisRackauckas](#):
>
> That should get an error throw.

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.

> [@ChrisRackauckas](#):
>
> Can you point to which ones you ran into?

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

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 17, 2023, 7:14pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/19 "2023-10-17T19:14:11Z")

</div>

> [@ChrisRackauckas](#):
>
> it’s not a good idea to pick a solver while debugging (especially the one that was chosen).

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

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [October 17, 2023, 8:05pm UTC](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028/20 "2023-10-17T20:05:11Z")

</div>

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.

[Next page](https://discourse.julialang.org/t/help-modeling-advection-diffusion-equation-with-modelingtoolkit-jl/105028.md?page=2)