Help modeling advection-diffusion equation with ModelingToolkit.jl

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