Discretization

ODEFunctionExpr(pdesys, discretization)

Generate an expression for the ODE function produced by method-of-lines discretization.

generate_code(pdesys, discretization[, filename])

Write generated discretized ODE function code for pdesys to filename.

get_discrete(pdesys, discretization)

Return a map from symbolic variables to the grids and discrete variables generated by discretization.

discretize(pdesys, discretization;
           analytic = nothing, checks = true, fallback = true, kwargs...)

Discretize pdesys and return a problem ready to solve.

For a time-dependent system this builds a DAEProblem. MethodOfLines emits residuals of the form D(u) - f ~ 0, which are already implicit-DAE form, so no mtkcompile is needed and the array (slice-form) equations survive at the System level. Residual code generation expands one entry per unknown, so discretize time and the first-call compile of prob.f grow with resolution (see MethodOfLines.jl#691; ModelingToolkit.jl#5139). Calling solve(prob) selects the default DAE algorithm.

A few systems cannot be posed as a first-order DAE — those second order in time, and those whose initialization equations BrownFullBasicInit would not honour. Those fall back to mtkcompile plus an ODEProblem, which scalarizes the array equations at the System level. Pass fallback = false to make that an error instead.

Supplying analytic selects the compiled ODEProblem path because analytic solutions are attached through the compiled ODEFunction.

Time-independent systems have no derivative to keep implicit and discretize to a NonlinearProblem as before.

Explicit Runge–Kutta methods such as Tsit5() solve ODEProblems, not the DAEProblem returned by this method. To use one, start from symbolic_discretize and compile the discretized system:

sys, tspan = symbolic_discretize(pdesys, discretization)
prob = ODEProblem(mtkcompile(sys), nothing, tspan)
sol = solve(prob, Tsit5())

chebyspace(N, dom)

Construct N Chebyshev-spaced grid points over the domain dom.

WENOScheme(; epsilon = 1.0e-6)

Jiang-Shu WENO-5 advection scheme for uniform and non-uniform grids.

Keyword Arguments

• epsilon: A quantity used to prevent vanishing denominators in the scheme, defaults to 1e-6. More sensitive problems will benefit from a smaller value. It is defined as a functional scheme.

MOLFiniteDifference(dxs, time=nothing;
                    approx_order = 2, advection_scheme = UpwindScheme(),
                    grid_align = CenterAlignedGrid(), kwargs...)

A discretization algorithm.

Arguments

• dxs: A vector of pairs of parameters to the grid step in this dimension, i.e. [x => 0.2, y => 0.1]. For a non-uniform rectilinear grid, replace any or all of the step sizes with the grid to use with that variable. It must be an AbstractVector, but not a StepRangeLen.
• time: The continuous variable, usually time. If time = nothing, discretization yields a NonlinearProblem. Defaults to nothing.

Keywords

• approx_order: The order of the derivative approximation.
• advection_scheme: The scheme used to discretize first-order spatial derivatives and associated coefficients. Defaults to UpwindScheme(). WENOScheme() is more stable and accurate at the cost of complexity.
• grid_align: The grid alignment value. Use center_align, edge_align, or StaggeredGrid() as appropriate for the discretization.
• edge_aligned_var: Required with StaggeredGrid(). The dependent variable, or a collection of dependent variables, that is edge aligned. Every other dependent variable is center aligned. An entry u => x or u => (x, y) edge aligns u in those directions only, as a Yee grid needs (Ex(t, x, y) => x).
• kwargs: Additional keyword arguments passed to the generated problem.

Fields

• dxs: A dictionary mapping each discretized independent variable to an integer grid size, a spacing, or an explicit grid.
• time: The independent variable left undiscretized, or nothing for a fully discretized system.
• approx_order: The requested finite-difference approximation order.
• advection_scheme: The scheme used for first-order spatial derivatives.
• grid_align: The grid alignment marker.
• should_transform: Whether supported symbolic transformations are applied before discretization.
• useIR: Whether ModelingToolkit's intermediate representation is used.
• callbacks: Symbolic discretization callbacks.
• kwargs: Additional keyword arguments forwarded to the generated problem.

Example

using ModelingToolkit
using MethodOfLines

@parameters t x
discretization = MOLFiniteDifference([x => 0.1], t)

ChebyshevCollocation(n; fft = true)

Chebyshev–Lobatto pseudospectral collocation with n grid points on the domain [a, b], including both endpoints. Spatial derivatives of order d are discretized with the differentiation matrix of the degree n - 1 polynomial interpolant through the grid values (equivalently, the maximal-stencil Fornberg weights on the Lobatto nodes). When an AbstractFFTs backend such as FFTW is loaded and fft = true, whole-direction derivatives in the array form are computed through the Chebyshev transform (a DCT-I via the FFT of the even extension) and the coefficient recurrence instead of the dense matrix, at O(n log n) per application; the two agree to roundoff. Pass fft = false to always use the dense matrix.

Use this scheme for directions with non-periodic (truncating) boundary conditions such as Dirichlet or Neumann conditions.

FourierCollocation(n; fft = true)

Fourier pseudospectral collocation with n distinct equispaced grid points on the periodic domain [a, b). The grid stores n + 1 points, with the upper endpoint identified with the lower endpoint by the periodic boundary condition u(t, a) ~ u(t, b).

Spatial derivatives of order d are discretized with the differentiation matrix of the trigonometric interpolant, computed as the dth power of the first-derivative matrix. When an AbstractFFTs backend such as FFTW is loaded and fft = true, whole-direction derivatives in the array form are applied with real FFTs instead of the dense matrix, O(n log n) per application rather than O(n^2); the two agree to roundoff. Pass fft = false to always use the dense matrix.

This scheme requires a periodic boundary condition in the associated direction; conversely, periodic directions must use FourierCollocation.

PseudospectralDiscretization(dxs, time = nothing; kwargs...)

A pseudospectral (collocation) discretization algorithm.

Each spatial independent variable is discretized on a collocation grid, and every spatial derivative Differential(x)^d is replaced by the dense differentiation matrix of the spectral interpolant in that direction, applied in physical (grid) space; no transform to spectral coefficients is taken. Boundary conditions replace the equation at the boundary nodes, exactly as in MOLFiniteDifference, with derivatives in a condition taken from the boundary row of the differentiation matrix. Nested derivative terms such as Differential(x)(a(u) * Differential(x)(u)) are evaluated directly by collocation, so no PDE-system transformation is needed.

The interior of each PDE is emitted as one symbolic array equation over slices of the discretized variables, with each derivative an opaque operator holding its differentiation matrix, so the number of symbolic equations and the size of the symbolic expressions are independent of the resolution. Residual code generation expands one entry per unknown (see MethodOfLines.jl#691; ModelingToolkit.jl#5139). Derivatives are applied with FFTs in both Fourier and Chebyshev directions when an AbstractFFTs backend such as FFTW is loaded. See the pseudospectral documentation page for the boundary conditions each grid type accepts and the scaling.

Arguments

• dxs: A vector of pairs mapping each independent variable to a collocation specification:
  • x => n::Integer or x => ChebyshevCollocation(n) discretizes x on n Chebyshev–Lobatto points. Use for non-periodic directions.
  • x => FourierCollocation(n) discretizes x on n distinct equispaced points and requires a periodic boundary condition u(t, a) ~ u(t, b) in that direction. Higher-order periodic matching conditions such as Differential(x)(u(t, a)) ~ Differential(x)(u(t, b)) are redundant - they are satisfied identically by the trigonometric interpolant - and are skipped.
  • x => grid::AbstractVector uses a custom set of n collocation nodes (polynomial-interpolation derivatives via Fornberg weights).
• time: The continuous variable, usually time. If time = nothing, discretization yields a NonlinearProblem. Defaults to nothing.

Keywords

• kwargs: Additional keyword arguments passed to the generated problem.

Example

using ModelingToolkit, MethodOfLines

@parameters t x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2

# Chebyshev collocation on 32 points
discretization = PseudospectralDiscretization([x => 32], t)

# Fourier collocation on 64 distinct points, for a periodic domain
discretization = PseudospectralDiscretization([x => FourierCollocation(64)], t)

Limitations

• Multi-domain (interface) boundary conditions are not supported; only same-variable periodic conditions may be used with FourierCollocation.
• Integral terms are not supported.
• Unlike MOLFiniteDifference, there is no upwinding or special scheme selection: all derivative orders in a direction share the same spectral differentiation matrix.
• Derivatives are dense matrix products, O(n^2) per application in one dimension, and the Jacobian is dense.

MOLDiscCallback(f, disc_ps)

Wrap a discretization callback function and its parameters for symbolic callback handling.

center_align

Grid alignment value for center-aligned finite difference grids.

edge_align

Grid alignment value for edge-aligned finite difference grids.

FunctionalScheme

F = FunctionalScheme{interior_points, boundary_points}(interior, lower, upper, is_nonuniform, parameters; name)

A user definable scheme that takes a set of functions as input. The functions define the derivative at the interior, lower boundary, and upper boundary.

lower and upper should be vectors of functions. In general, upper and lower must be at least floor(interior_points/2) long. Where you have no good approximation for a derivative at the boundary, you can use nothing as a placeholder. MethodOfLines will then attempt to use an extrapolation here where necessary. Be warned that this can lead to instability.

The boundary functions define the derivative at their index in the function vector, numbering from the boundary. For example, if boundary_points = 3, the first function in the vector will define the derivative at the boundary, the second at the boundary plus one step, and the third at the boundary plus two steps.

The functions making up the scheme take the following inputs:

Functions must be of the form f(u, p, t, deriv_iv, d_iv).

Where the function would branch on a value of u, p, t, or d_iv, use ifelse instead of standard conditionals.

For the interior, u takes a vector of dependent variable values in the direction of the derivative of length interior_points. interior_points must be odd, as this function defines the derivative at the center of the input points.

For the lower and upper boundaries, u takes a vector of dependent variable values of length boundary_points. This will be the boundary_points number of points closest to the lower and upper boundary respectively. p will take all parameter values in the order specified in the PDESystem, with the scheme's parameters prepended to the list.

deriv_iv takes a vector of independent variable values of the same support as for u, for the independent variable in the direction of the derivative.

If is_nonuniform is false, d_iv will take a scalar value of the stepsize between the points used to call the function in u and deriv_iv.

If is_nonuniform is true, the scheme must be able to accept d_iv as a vector of stepsizes between the points used to call the function in u and deriv_iv, therefore of length length(u)-1. A method should also be defined for the case where d_iv is a scalar, in which case the stepsizes are assumed to be uniform.

UpwindScheme([approx_order])

Upwind finite difference scheme for advection terms.