How it works

MethodOfLines.jl makes heavy use of Symbolics.jl and SymbolicUtils.jl, namely it's rule matching features to recognize terms which require particular discretizations.

See here for the highest level overview of the algorithm.

Given your discretization and PDESystem, we take each independent variable defined on the space to be discretized and create a corresponding range. We then take each dependant variable and create an array of symbolic variables to represent it in its discretized form. This is stored in a DiscreteSpace object, a useful abstraction.

We recognize boundary conditions, i.e whether they are on the upper or lower ends of the domain, or periodic here, and use this information to construct the interior of the domain for each equation here. Each PDE is matched to each dependant variable in this step by which variable is highest order in each PDE, with precedance given to time derivatives. This dictates which boundary conditions reduce the size of the interior for which PDE. This is done to ensure that there will be the same number of equations as discrete variable states, so that the system of equations is balanced.

Next, the boundary conditions are discretized, creating an equation for each point on the boundary in terms of the discretized variables, replacing any space derivatives in the direction of the boundary with their upwind finite difference expressions. This is the place to look to see how this happens.

After that, the system of PDEs is discretized creating a finite difference equation for each point in their interior. Specific terms are recognized, and the best implemented scheme for these terms dispatched. For example advection terms are discretized with the upwind scheme. There are also special schemes for the nonlinear laplacian and spherical laplacian. See here for how this term matching occurs, note that the order the generated rules are applied is important, with more specific rules applied first to avoid their terms being matched incorrectly by more general rules. The SymbolicUtils.jl docs are a useful companion here. See here for the practical implementation of the finite difference schemes.

Now we have a system of equations which are either ODEs, linear, or nonlinear equations and an equal number of unknowns. See here for the system that is generated for the Brusselator at low point count. The structure of the system is simplified with ModelingToolkit.structural_simplify, and then either an ODEProblem or NonlinearProblem is returned. Under the hood, the ODEProblem generates a fast semidiscretization, written in Julia with RuntimeGeneratedFunctions. See here for an example of the generated code for the Brusselator system at low point count.