Discretization

center_align

Grid alignment value for center-aligned finite difference grids.

edge_align

Grid alignment value for edge-aligned finite difference grids.

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.

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.

MOLDiscCallback(f, disc_ps)

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

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.
• 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)

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 generated code are independent of the resolution. 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.

UpwindScheme([approx_order])

Upwind finite difference scheme for advection terms.

ODEFunctionExpr(pdesys, discretization)

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

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.

chebyspace(N, dom)

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

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 reach the generated code intact. 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. 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())

BrownFullBasicInitUnsafeError(offenders)

Raised by DAEProblem(::PDESystem, ::MOLFiniteDifference) when the discretized system carries initialization equations that BrownFullBasicInit() would silently discard, so no default initialization algorithm can be chosen safely. offenders pairs each such equation with the reason it is not honored.

CenterAlignedGrid()

Grid alignment strategy that places the first and last grid points on the domain boundaries. Use the singleton center_align as the grid_align value in MOLFiniteDifference.

This is a stateless marker type and has no fields.

DiscreteSpace(domain, depvars, indepvars, discretization::MOLFiniteDifference)

A type that stores information about the discretized space. It takes each independent variable defined on the space to be discretized and create a corresponding range. It then takes each dependent variable and create an array of symbolic variables to represent it in its discretized form.

Arguments

• domain: The domain of the space.
• vars: A VariableMap object that contains the dependent and independent variables and other important values.
• discretization: The discretization algorithm.

Properties

• ū: The vector of dependent variables.
• args: The dictionary of the operations of dependent variables and the corresponding arguments, which include the time variable if given.
• discvars: The dictionary of dependent variables and the discrete symbolic representation of them. Note that this includes the boundaries. See the example below.
• time: The time variable. nothing for steady state problems.
• x̄: The vector of symbolic spatial variables.
• axies: The dictionary of symbolic spatial variables and their numerical discretizations.
• grid: Same as axies if CenterAlignedGrid is used. For EdgeAlignedGrid, interpolation will need to be defined ±dx/2 above and below the edges of the simulation domain, where dx is the step size in the direction of that edge.
• dxs: The discretization of symbolic spatial variables and their step sizes.
• Iaxies: The dictionary of the dependent variables and their CartesianIndices of the discretization.
• Igrid: Same as axies if CenterAlignedGrid is used. For EdgeAlignedGrid, one more index will be needed for extrapolation.
• x2i: The dictionary of symbolic spatial variables and their ordering.

Examples

julia> using MethodOfLines, DomainSets, ModelingToolkit
julia> using MethodOfLines:DiscreteSpace

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

julia> eq  = [Dt(u(t, x)) ~ Dxx(u(t, x))]
julia> bcs = [u(0, x) ~ cos(x),
              u(t, 0) ~ exp(-t),
              u(t, 1) ~ exp(-t) * cos(1)]

julia> domain = [t ∈ Interval(0.0, 1.0),
                 x ∈ Interval(0.0, 1.0)]

julia> dx = 0.1
julia> discretization = MOLFiniteDifference([x => dx], t)
julia> ds = DiscreteSpace(domain, [u(t,x).val], [x.val], discretization)

julia> ds.discvars[u(t,x)]
11-element Vector{Num}:
  u[1](t)
  u[2](t)
  u[3](t)
  u[4](t)
  u[5](t)
  u[6](t)
  u[7](t)
  u[8](t)
  u[9](t)
 u[10](t)
 u[11](t)

julia> ds.axies
Dict{Sym{Real, Base.ImmutableDict{DataType, Any}}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}} with 1 entry:
  x => 0.0:0.1:1.0

EdgeAlignedGrid()

Grid alignment strategy that places grid points half a spacing from the domain boundaries. Use the singleton edge_align as the grid_align value in MOLFiniteDifference, for example when edge-centered values improve Neumann boundary accuracy.

This is a stateless marker type and has no fields.

MOLMetadata

A type used to store data about a PDESystem, and how it was discretized by MethodOfLines.jl. Used to unpack the solution.

• discretespace: a DiscreteSpace object, used in the discretization.
• disc: a Discretization object, used in the discretization. Usually a MOLFiniteDifference object.
• pdesys: a PDESystem object, used in the discretization.

A matched nonlinear laplacian Dx(expr * Dx(u)): term is the matched additive term, pre any prefactor product and div any divisor (nothing when absent).

SpectralApply{J, S}

Symbolic array operator applying a differentiation operator (a dense block, or a whole-direction FFT operator) along axis J of its argument through apply_along, producing an array of size S. The operator is a field rather than a literal in the expression tree, so generated code references it as a constant instead of spelling out every entry.

SpectralDerivativeOperator

Dense differentiation matrix acting along one independent variable. mat[i, :] is the differentiation row for matrix-index i; rowmap[g] maps grid index g to its matrix row; taps lists the grid indices the columns act on. fast is an operator equivalent to mat over the whole direction that the array form prefers (an FFT operator, see fast_spectral_operator), or nothing.

The DifferentialDiscretizer analogue for PseudospectralDiscretization: maps Differential(x)^d to the corresponding SpectralDerivativeOperator.

A matched spherical laplacian Dx(x^2 * expr * Dx(u)) / x^2: term is the matched additive term, expr the inner coefficient besides x^2 and pre any prefactor product (nothing when absent).

StaggeredGrid()

Grid alignment strategy for variables whose grid locations are staggered relative to the primary grid. Pass an instance as grid_align to MOLFiniteDifference and specify the corresponding variable alignment when constructing the discretization.

This is a stateless marker type and has no fields.

DAEProblem(pdesys::PDESystem, discretization::MOLFiniteDifference; kwargs...)

Discretize pdesys and build a DAEProblem from the residuals MethodOfLines emits, without running mtkcompile. The array equations reach the generated code intact, which mtkcompile would undo.

initializealg defaults to BrownFullBasicInit(), which is the only algorithm that reproduces the discretize result here, and is chosen only when the discretized system's initialization equations are ones it honors; otherwise a BrownFullBasicInitUnsafeError is raised naming the offending equations. Passing initializealg explicitly overrides both the default and that check.

The solution is a PDETimeSeriesSolution, the same wrapper discretize produces, so it is indexed and interpolated by the PDESystem's own variables: sol[u(t, x)], sol(t, x).

A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the centered scheme.

A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the half offset centered scheme. See table 2 in https://web.njit.edu/~jiang/math712/fornberg.pdf

A helper function to compute the coefficients of a derivative operator including the boundary coefficients in the upwind scheme.

Idx(II::CartesianIndex, s::DiscreteSpace, u, indexmap)

Here indexmap maps the arguments of u in s to their ordering. Return a subindex of II that corresponds to only the spatial arguments of u.

axis_cumsum_pad(A, dim)

cumsum along dim after a leading zero slab. The result is one longer than A and starts at 0. Per-axis wrappers exist because @register_array_symbolic evaluates size with a symbolic dim.

axis_sum(A, dim)

dropdims(sum(A; dims = dim); dims = dim). Base sum(; dims =) on a MethodOfLines slice can leave a broken size.

_dae_problem(sys, tspan, discretization; kwargs...)

Build the DAEProblem from an already discretized system. Shared by DAEProblem(::PDESystem, ::MOLFiniteDifference) and discretize so that the system is discretized once.

apply_along(op, X, ::Val{J})

Apply the whole-direction derivative operator op along axis J of X. For a matrix op this is op * X contracted over axis J; the AbstractFFTs extension adds FFT operators.

Align a slice of u (in u's IV order) to the equation axes in ranges. Prefix axes are left as-is; any other subset is permuted and reshaped with singleton dimensions — the array form of Idx.

Decompose the interior into the boxes on which one array equation is valid, as the index ranges to take in each dimension (the boxes are their cartesian product). Returns those ranges and, per dimension, which of them holds points whose stencils do not wrap.

In a direction with no interface boundaries this is the single subbox of the interior on which every derivative resolves to the translation-invariant interior stencil, mirroring the branch conditions in central_difference_weights_and_stencil and _upwind_difference; interior points outside it — the frame — are discretized pointwise.

A wrapping end — self-periodic or two-domain — has no such boundary branch: the interior stencil applies through that end, but for points within a stencil of the seam some taps wrap and the wrapped tap is not contiguous with the rest. Splitting those points off as one range each keeps every tap of every box a single contiguous slice, at the cost of a handful of extra equations — as many as the stencil is wide, so the count still does not depend on the grid resolution. A direction may wrap at one end only; the other end still shrinks as a regular boundary.

nllap_orders maps each direction carrying a nonlinear laplacian to the coefficient's derivative orders above one; those directions additionally take the half-offset branch conditions and tap extents of array_nonlinlap_constraints. sph_orders does the same for spherical laplacians via array_spherical_constraints, and additionally keeps any r ≈ 0 points out of the core (the pointwise path treats them with a separate branch). mixedorders maps each direction a mixed derivative reaches along to the centered orders used there. Those use the centered stencil of that order, which for an odd mixed order is wider than the winding stencil the same entry of pdeorders would select — at first order and approx_order 4/6 this is the original mixed first-order case.

array_bc_eqs(s, boundary, interiormap, derivweights, bcmap)

Generate the equations for boundary as a single symbolic array equation over the face it occupies, rather than one scalar equation per point on that face.

This works because the index along the boundary's own direction is fixed across the face, so every point on it selects the same stencil weights and tap offsets — the same translation invariance the interior exploits, applied one dimension down. Without this, boundary equations stay pointwise and dominate the equation count in 2D and 3D, where they scale with the surface (O(n) and O(n^2)) while the interior collapses to one.

Throws ArrayFormFallback for boundaries with no slice representation, in which case the caller emits the pointwise form.

array_bc_eqs(s, boundary::HigherOrderInterfaceBoundary, interiormap, derivweights, bcmap)

Flux (or other higher-order) interface condition as one array equation over the face. boundary.u is discretized on this edge; boundary.u2 on the partner face at index 1. A self-periodic flux reuses one discrete variable on both ends; each call is bound to its own face. A 1D (single-point) face falls back.

array_bc_eqs(s, boundary::InterfaceBoundary, interiormap, derivweights, bcmap)

Equate the two faces an interface (periodic) boundary joins as a single array equation, the slice form of the disc1[II] ~ disc2[II + Ioffset] the pointwise path emits per point.

As in the pointwise path only the lower boundary of the pair carries the equations; the upper one repeats the same relation and contributes none.

Bind numeric spatial arguments of time-literal v into formula.

The IC is a continuous expression, so v(0, 1) substitutes x => 1, not a grid index. Free spatial arguments are left for the outer view to evaluate.

Interfaces for a boundary-value derivative in x: wrap if periodic, else [].

The slice form of Dop applied to u on the face at II0, for the boundary direction x_ (equation axis j of N): the weights and taps the pointwise path would use at a representative point on the face, with the taps as shifted slices.

array_boundary_derivative_expr for a spectral operator: the boundary row of the differentiation matrix applied to the full slice along x_.

Resolve a numeric boundary argument to a 1-based grid index at a domain edge, matching the arithmetic in newindex.

depvarderivbcmaps stencil for (Differential(x)^d)(u_) as a scalar or face slice over the core box.

Rules mapping each handleable boundary-value derivative in pde to its stencil.

Handleable (expr, x, d, u_) records in pde, e.g. Dx(u(t, 1)).

Index for u_: edge on pinned arguments, first core index on free ones.

Substitution rules mapping each boundary value in pde to its array element, face slice, or edge-aligned interpolant over the core box described by ranges. Call after array_validate_boundary_values.

Returns a Vector{<:Pair} even when empty, so vcat into ArrayifyContext.rules stays well-typed.

Numeric-arg depvars whose time slot is still s.time. Time-literals are collected by array_time_literal_terms so u(0, x) is not read as u(t, 0).

The array element (every spatial argument fixed) or face slice (some arguments free over the core box) corresponding to a boundary value like u(t, 1) or u(t, 0, y). On a center-aligned grid this is a direct index; on an edge-aligned grid it is the interpmap interpolant at the domain edge. Fixed dimensions use a singleton k:k range so the result broadcasts against the core slice; a fully-fixed reference is a scalar, which broadcasts as well.

Array form of central_difference on the core region for the even order derivative (Differential(x)^d)(u).

Per-slot coefficient arrays of a split functional scheme: getcoeffs(i) is evaluated once per point of rng and its nslots entries sliced into broadcastable arrays along dimension j of N, like array_weight_vals.

The array equation for one box of the interior, given as ranges (dimension => index range). depvars are the sliceable fields; allvars includes rank-dropping integrands. bcmap builds pinned-direction derivatives of interior boundary values.

array_corner_eqs(s, interiormap, u, N)

Equations for the points that lie outside the interior in two or more dimensions — the corners in 2D, and the edges as well as the corners in 3D — as one array equation per contiguous box instead of one scalar equation per point.

generate_corner_eqs! builds this region with setdiff, which yields a bag of indices and loses the structure. The region is in fact a union of boxes: along each dimension a point lies in the lower band, the interior band, or the upper band, and this region is exactly the combinations with at least two non-interior bands. Enumerating those gives 3^N - 2N - 1 boxes — none in 1D, the 4 corners in 2D, and the 12 edges plus 8 corners in 3D — a count that does not depend on the grid resolution.

Without this the 3D edges stay pointwise and cost 12n - 10 equations, which is what keeps 3D at O(n) once the faces are sliced.

Trapezoidal running integral of u along x as a slice over the core described by ranges. Reads the full axis of u in x.

The interpolated domain-edge value of u_ on an EdgeAlignedGrid.

A single pinned argument (u(t, 1), u(t, 0, y)) is the interpmap stencil along that axis — a scalar when every spatial argument is fixed (or the remaining ranges are length 1), a face of singleton-range taps otherwise. Several pinned arguments (u(t, 0, 0)) are the tensor product of those interpolators, which is the value at the domain corner the one-axis pointwise map does not form.

Resolve a numeric boundary argument to the interpolation index newindex uses with shift = true on an EdgeAlignedGrid: the lower domain edge is 1, the upper is n - 1 so half_offset_centered_difference sits on the midpoint between that node and the next.

array_boundary_edge_index is the unshifted (n) form and must stay that way: center-aligned views and any later derivative path index the last node directly.

Array form of function_scheme on the interior branch, for the first derivative of u in x: the shared per-direction trace with its taps replaced by shifted slices of u and its coefficient slots, if any, by the numeric arrays in crules.

Trace of a functional advection scheme for one direction x, shared by every variable advected in x: the traced expression, tap offsets, placeholder taps usyms, and on nonuniform grids the coefficient slots csyms with their split.

The interior is translation invariant, so the scheme traces once and its taps are later replaced by shifted slices; solution-dependent weights (WENO's smoothness indicators) are tap arithmetic and broadcast like any other term. Schemes that read the grid coordinate fall back: the pointwise path folds those numeric coordinates, which a trace would rebuild reassociated. Nonuniform grids trace the apply kernel of the scheme's array_scheme_split instead; schemes without a split fall back. Periodic nonuniform directions share the split kernel, with coefficient windows unwrapped across the seam by array_scheme_coeff_rules.

Whether the derivative of order d is discretized by a functional advection scheme (WENO and friends) rather than by the winding rules, mirroring the branch on the advection scheme in generate_finite_difference_rules: only the first derivative is handled by the scheme, the higher odd orders still wind.

The numeric grid values of x over the core region, shaped to broadcast along the dimension of x in an N-dimensional array expression.

Slice form of the half-offset operator D applied to v at offset o from each core point (interior branch of get_half_offset_weights_and_stencil). On nonuniform grids the weights vary per point, indexed as stencil_coefs[i + o - boundary_point_count].

True when a time-literal sits under any Differential. Dt of a frozen field is zero and would drop an initial-velocity residual; Dx(u(0, x)) stays fallback so this does not overlap #664.

Evaluate an IC formula at spatial_subs and the initial time.

Look up the order-0 IC formula for canonical u.

Canonical dependent variable => order-0 IC formula after symbolic_linear_solve. Time-literals evaluate this field; the discrete unknown is U(t), not U(t0).

Index ranges of u for an integral along x: the full axis in x, the equation core (or the full axis) in every other independent variable.

The interior branch of central_difference_weights_and_stencil for Dop along dimension j of N over the index range rng, as (weights, taps). On a nonuniform grid the interior weights vary from point to point, so each returned weight is the broadcastable numeric array of that tap's weight over rng rather than a scalar.

Numeric grid values of x interpolated to the half-offset point at offset o, shaped to broadcast along the dimension of x (slice form of map_ivs_to_interpolated, with _wrapperiodic-style wrapping).

Weights and tap indices of derivweights.interpmap[x] at the edge-aligned interpolation index II_j, from get_half_offset_weights_and_stencil with bs = [] — the same call boundary_value_maps makes for EdgeAlignedGrid.

True when expr is (Differential(x)^d)(u_) with x numeric on u_.

True when the argument of u_ for x is numeric.

Array form of mixed_central_difference on the core region for (Differential(x)^m * Differential(y)^n)(u).

The scalar scheme is the tensor product of the two centered stencils: a sum over the taps in x of a sum over the taps in y of wx*wy*u[II + kx + ky]. Every point of the core takes the interior branch of both, so the whole thing is one broadcasted sum of slices shifted along two axes at once — array_central_difference with a second shifted axis. The weights come from the same DerivativeOperators the pointwise path uses and their products are numeric, so the two agree term by term. The first-order case Dx(Dy(u)) is m = n = 1.

The (u, x, m, y, n) records for which the mixed derivative (Differential(x)^m * Differential(y)^n)(u) occurs in pde.

Two distinct spatial variables and a dependent variable is the only mixed family the pointwise path has a scheme for (generate_mixed_rules), and so the only one with a slice form here. Three-or-more spatial directions and mixed derivatives of anything other than a dependent variable reach arrayify with a spatial differential still in place and fall back.

Slice form of cartesian_nonlinear_laplacian for a matched Dx(expr * Dx(u)): at each half-offset point of the outer stencil, expr * Dx(u) is rebuilt over shifted slices (dependent variables interpolated, same-x derivatives via the half-offset operators, x interpolated numerically), then combined with the outer weights. Unhandled patterns in the coefficient surface as ArrayFormFallback from arrayify.

Band bounds and tap extents (lo, hi, mintap, maxtap) a nonlinear laplacian imposes in direction x: the interior branch conditions of get_half_offset_weights_and_stencil for the outer operator (applied at II - 1 on the clipped length n - 1) and, at each of its half-offset taps, for the interpolator and the inner half-offset derivatives.

Rules binding each matched term to its slice form: the laplacian from array_nonlinear_laplacian, grid-constant prefactors broadcast on, divisors discretized with the base rules (mirrors replacevals).

True when an occurrence of a higher-dimensional variable leaves an extra IV free. A face reference such as φ(t, x, 1.0) is not free.

Dependent-variable occurrences on both sides of pde.

permutedims as a symbolic array term. Promotion methods keep the axes.

Pointwise depvarderivbcmaps value for (Differential(x)^d)(u_) at this box.

Original and grid-valued keys for each handleable boundary-value derivative at this singleton box. Applied before expand_derivatives in the pointwise path.

reshape as a symbolic array term. Promotion methods keep the axes.

IC formula for u, with nested time-literals resolved.

u(0,x) ~ v(0,x) substitutes v's formula. A fixed spatial argument is bound in continuous coordinates (v(0, 1) becomes the formula at x = 1), not replaced by the whole field of v. Remaining depvars, including live fields left by parse_bcs, are replaced by that variable's IC.

Rules binding the coefficient slots of a split scheme trace to their numeric per-point arrays over the core box; empty on uniform grids. The windows fed to coeffs are the same grid windows the pointwise path sees; in a periodic direction taps beyond either end take the periodically shifted coordinate bcoord would produce.

array_scheme_split(F::FunctionalScheme)

Coefficient split of a functional scheme for nonuniform grids; nothing (the default) falls back to the pointwise path there. Not exported: a scheme opts in by defining a method on MethodOfLines.array_scheme_split, as WENO does.

A split (coeffs, apply, nslots) factors interior into grid geometry and solution arithmetic: coeffs(xwindow) maps the interior_points-long window of grid coordinates to nslots numbers, and apply(u, p, t, c) recombines them with the taps so that apply(u, p, t, coeffs(xwindow)) equals interior(u, p, t, xwindow, dxwindow) up to reassociation. apply is traced once on placeholder symbols and must be branch free; the slots are evaluated numerically per grid point. If each slot holds exactly the constant the scalar trace folds at that spot and enters apply linearly, the array form matches the pointwise path bitwise (WENO's split does, pinned by its property test); a split that reassociates the fold agrees to ~1e-15 relative instead.

Placeholder symbols standing in for one argument vector of a functional scheme while it is traced. They are substituted away before the expression leaves array_function_scheme, so they never reach the discretized system; the ## prefix keeps them clear of user names.

A slice of the array variable for u over the core region, shifted by offsets[j] in each dimension j it names and wrapped around the seam where that dimension wraps (see wrap_tap_range). A two-domain wrap reads the partner array; a self-periodic wrap stays on u. Dimensions absent from offsets are taken unshifted and unwrapped, which is what the pointwise path does with the dimensions a stencil does not reach along.

Naming a dimension with offset 0 is not the same as omitting it: _wrapperiodic maps the first index of a periodic dimension onto the last for every tap of a stencil that reaches along it, the centre tap included. A destination wrap is applied once: the partner array is not wrapped again, matching wrapinterface on a RefCartesianIndex.

A slice of the array variable for u over the core region, optionally shifted by offset in the dimension of shiftx, wrapped around the seam if that dimension is periodic (see wrap_periodic_range). A 0D variable is the matching scalar. The slice is aligned to the equation axes.

Dependent variables whose spatial IVs are a subset of args (including 0D). Rank-dropping integrands are excluded; array_integral_rules handles them.

Band bounds and tap extents (lo, hi, mintap, maxtap) a spherical laplacian imposes in direction x: those of the nonlinear laplacian it contains, plus the interior branch of the centered first derivative the scheme adds.

Slice form of spherical_diffusion for a matched Dx(x^2 * expr * Dx(u)) / x^2, away from x ≈ 0 (scheme 1 in appendix A of the paper referenced there): the coefficient at the grid points times the centered first derivative divided by the grid values of x plus the nonlinear laplacian of the inner coefficient.

Rules binding each matched spherical term to its slice form, grid-constant prefactors broadcast on (mirrors the splicing in generate_spherical_diffusion_rules).

Array form of the interior branch of the staggered generate_cartesian_rules: the same windmap weights, with the two taps fixed by each variable's alignment — (0, +1) for a center-aligned variable, (-1, 0) for an edge-aligned one — constant across the core.

Broadcasted weighted sum of shifted slices: the array-form analogue of sym_dot.

Substitute boundary values and time-literals into an already pointwise-discretized equation at the single point described by ranges (all singleton). Used for size-1 wrap boxes and frame points, where discretize_equation_at_point leaves interface-face and free-standing-corner boundary values symbolic. valmaps has already replaced free spatial arguments with their grid values inside those leftover terms (u(t, 0, y) appears as e.g. u(t, 0, 0.2), u(0, x) as u(0, x_i)), so each rule keys on that grid-valued form and maps to a scalar array element (center-aligned) or interpolant (edge-aligned).

Boundary-value derivatives are already replaced via extra_rules before expand_derivatives; this leftover sweep is for leftover values and time-literals. A time-literal's time slot is not substituted; the value is the IC formula.

The most negative and most positive tap offsets any derivative in the equation applies in direction x, mirroring the interior branches of central_difference_weights_and_stencil, _upwind_difference and get_f_taps_coords.

On a staggered grid the interior taps are (0, +1) or (-1, 0) depending on each variable's alignment; the union over both alignments is taken, so at worst one extra point per end lands in the pointwise frame.

A mixed derivative reaches along x with the centered stencil of the mixed order in that direction, rather than the winding one pdeorders would select for an odd order, so those orders take their centered taps too. For order 1 this is the first-order centered stencil, which at approximation orders 4/6 is wider than the winding stencil.

Spatial indices of a time-literal over ranges. Free arguments take the core range (or its first index when singleton); fixed arguments are domain-edge indices. The time slot is skipped.

Evaluate the IC formula of u_ on the spatial view. Arrays are one symbolic term so treesize does not grow with the grid. Fixed spatial arguments become length-1 axes to broadcast against the core, as in array_boundary_value_view.

Numeric literal in the time slot of u_, or nothing.

Aligns arguments(u_) with s.args[operation(u_)] so u(x, 0) is recognized when time is not first. A bare Number test would read u(0, x) as the spatial edge u(t, 0) whenever 0 is the left end of x.

The canonical fields named by occs, uniqued by isequal.

Array form of upwind_difference on the core region, for one winding direction.

Equation-level checks for interior boundary values and time-literals. Throws ArrayFormFallback when there is no slice form (staggered spatial edges, edge-aligned mixed time+spatial edges, edge-aligned boundary-value derivatives, non-IC times, off-edge samples).

Pure time-literals are valid on every grid alignment. Spatial boundary values are accepted on CenterAlignedGrid and EdgeAlignedGrid (the latter via array_edge_aligned_boundary_value). Mixed time+spatial literals require CenterAlignedGrid.

Equation-level check: every dependent-variable occurrence is a subset of the equation axes (including 0D), an integral rank drop (array_compatible_depvar), or a boundary value that fixes every extra IV.

The underlying (unscalarized) array variable of which s.discvars[u] holds the elements.

Per-point stencil weights as a broadcastable numeric array along dimension j of N, for nonuniform grids where the interior weights vary from point to point. getweights(i) returns the weight SVector at grid index i.

Whole-domain trapezoidal integral of u along x. Rank drops in x: a 1D integrand becomes a scalar; a higher-dimensional integrand keeps its other axes.

Array form of the winding selection for an odd derivative multiplied by expression expr, mirroring the pointwise path's ifelse(coef > 0, coef*pos, coef*neg).

When the coefficient does not vary over the grid — a literal, a parameter, or any expression of time alone — the wind direction is one scalar condition for the whole slice, so ifelse broadcasts and reproduces the pointwise path exactly.

A grid-varying coefficient needs a per-point condition, and ifelse cannot be broadcast over a symbolic array condition (the elementwise comparison carries symtype Any rather than Bool). Those use max(coef, 0)*pos + min(coef, 0)*neg, which agrees with ifelse on finite values but yields NaN rather than the finite branch when the unselected stencil is Inf/NaN.

Coordinate of raw tap index i in a wrapping direction. Taps at or below the first point take the lower-end destination chart, taps past the last point the upper-end destination. Self-periodic (dest === :self) is a single add/subtract of the period; a two-domain interface uses the partner grid and a contiguous (zero) shift when the physical edges coincide. Mirrors _wrapcoord bit for bit — keep them in lockstep.

array_periodic_coord is the self-periodic special case.

The wrap destination of interface boundary b on u along x: :self when the join is the same variable at the other end of the same independent variable, otherwise the partner's discrete variable, independent variable and grid length.

Throws when the join is the same end of both domains, or when the partner's array layout is not the same CartesianIndex the pointwise wrapinterface writes into.

The directions in which every dependent variable wraps — self-periodic or two-domain interface — mapped to a wrap specification.

A self-periodic direction is one whose interface boundaries join a variable to itself at the other end of the same independent variable, which is what u(t, 0) ~ u(t, 1) parses to. A two-domain interface joins different variables (typically at one end only); haslowerupper still reports that end as an interface, so the interior stencil applies there with taps wrapped onto the partner array by bwrap, which wrap_tap_range reproduces on slices.

A nonuniform wrapping direction is admitted: operators whose seam form the pointwise path cannot build (linear stencils, half-offset operators) throw at their own sites instead.

Wrap specification for one direction: source length n and, per end, where a tap that leaves the grid is read from.

lower / upper are nothing (that end is a regular boundary; shrink the core), :self (self-periodic: wrap onto the same array with the same length), or (; u, x, n) naming the partner variable a two-domain interface joins.

arrayify(expr, ctx)

Broadcast-aware substitution: rebuild expr bottom-up, replacing any subterm that matches a rule in ctx.rules (first match wins) and broadcasting any operation that receives an array-valued argument. Time differentials are applied directly to their (array-valued) arguments; any spatial differential that survives the rules means the expression contains a scheme this path does not support, so fall back.

Calls f(u, θ) map over the field slice and f([u, v], θ) over the stacked slices, scalars in the literal broadcast onto the slice; indexed calls select each point's output. Networks flagged by batched_callable are evaluated in one batched call.

A function that returns what to replace independent variables with in boundary equations

batched_callable(f)

Whether the symbolic callable parameter f evaluates all grid points in one call, one column per point, as Lux networks do. The default applies f point by point; the ModelingToolkitNeuralNets extension enables batching for its networks.

bcoord(I, bs, s, jx)

Physical coordinate of raw tap index I in the differentiated grid's chart. Wrapped taps use the exact interface chart transition: periodic shift = period length, contiguous shift = 0. Result is strictly increasing across the seam.

Mirrored bit for bit by array_periodic_coord for self-periodic directions and by array_wrap_coord for two-domain interfaces; keep them in lockstep.

brown_init_offenders(sys)

Return the initialization equations of sys that BrownFullBasicInit() would not honor, each paired with the reason, and an empty vector when the algorithm is safe for sys.

BrownFullBasicInit() holds the differential variables at their given values and solves for everything else, so an initialization equation survives it only when it fixes a single differential unknown to a value involving no other unknown. Everything else is reported, including equations this predicate cannot classify.

No stencil extents: spectral rows always span the entire direction, so no ghost-point padding is ever required.

cartesian_nonlinear_laplacian

Differential(x)(expr(x)*Differential(x)(u(x)))

Given an internal multiplying expression expr, return the correct finite difference equation for the nonlinear laplacian at the location in the grid given by II.

The inner derivative is discretized with the half offset centered scheme, giving the derivative at interpolated grid points offset by dx/2 from the regular grid.

The outer derivative is discretized with the centered scheme, giving the nonlinear laplacian at the grid point II. For first order returns something like this: d/dx( a du/dx ) ~ (a(x+1/2) * (u[i+1] - u[i]) - a(x-1/2) * (u[i] - u[i-1]) / dx^2

For 4th order, returns something like this:

first_finite_diffs = [a(x-3/2)*finitediff(u, i-3/2),
                      a(x-1/2)*finitediff(u, i-1/2),
                      a(x+1/2)*finitediff(u, i+1/2),
                      a(x+3/2)*finitediff(u, i+3/2)]

dot(central_finite_diff_weights, first_finite_diffs)

where finitediff(u, i) is the finite difference at the interpolated point i in the grid.

And so on.

central_difference_weights_and_stencil for SpectralDerivativeOperator: the stencil of an interior point is the full matrix row over all collocation points of the direction.

Performs a centered difference in x centered at index II of u ufunc is a function that returns the correct discretization indexed at Itap, it is designed this way to allow for central differences of arbitrary expressions which may be needed in some schemes

Check that term is a compatible derivative argument, and return the term if it is not and whether to expand.

Check that a derivative matching condition across a periodic seam holds identically under FourierCollocation, where both ends of the seam share a differentiation row. Evaluated at a single edge point: the row selection is the same across the face.

spectral_clip_interior!!: higher-order periodic matching conditions are satisfied identically by the trigonometric interpolant, so they must not truncate the interior - the corresponding equations are skipped in discretize_equation! as well.

Whole-domain trapezoidal integral as a scalar at the transverse location of II.

Turn term in to an auxiliary variable, and replace it in the equations and boundary conditions.

Modified copilot explanation: #= Here is the explanation for the code above:

1. First we create a new variable to represent our term, and replace the term with the new variable.
2. Then we generate the equation for the new variable, and add it to the equation list.
3. Then we generate the replacement rules for the boundary conditions, and substitute them into the new equation to infer auxiliary bcs.
4. Finally we add the new boundary conditions to the boundarymap and pmap. =#

Finds incompatible terms in the equations and returns them with the incompatible part and whether to expand the term.

map dependent variables

discretize_equation_array_form(pde, interior, s, depvars, derivweights, bcmap,
                               eqvar, indexmap, boundaryvalfuncs)

Discretize the interior of pde as symbolic array equations over slices of the discretized dependent variables — one per box of the decomposition built by array_bands, which is a single box unless the equation has periodic directions — plus pointwise scalar equations for any interior points whose stencils differ from the translation-invariant interior stencil. Throws ArrayFormFallback when pde contains a pattern that cannot be represented this way.

Discretize space

discretize_spectral_array_form(pde, interior, s, depvars, derivweights, eqvar, indexmap)

The interior of pde as a single symbolic array equation over the interior box, with each spatial derivative a dense differentiation block applied along its axis. Throws ArrayFormFallback for patterns without a slice form here: stationary systems and whatever arrayify declines.

fast_spectral_operator(spec, grid, d, mat)

An operator applying the order-d derivative over the whole direction with FFTs, equivalent to the dense matrix mat, or nothing when no FFT backend is available. Defined by the AbstractFFTs extension for FourierCollocation and ChebyshevCollocation; the default is nothing.

Returns the term if it is incompatible, and whether to expand the term.

Create a vector containing indices of the corners of the domain.

This is a catch all ruleset, as such it does not use @rule. Any even ordered derivative may be adequately approximated by these.

generate dxs

generate_extrap_eqs

Pads the boundaries with extrapolation equations, extrapolated with 6th order lagrangian polynomials. Reuses central_difference as this already dispatches the correct stencil, given a DerivativeOperator which contains the correct weights.

generate_finite_difference_rules

Generate a vector of finite difference rules to dictate what to replace variables in the pde with at the gridpoint II.

Care is taken to make sure that the rules only use points that are actually in the discretized grid by progressively up/downwinding the stencils when the gridpoint II is close to the boundary.

There is a general catch all ruleset that uses the cartesian centered difference scheme for derivatives, and simply the discretized variable at the given gridpoint for particular variables.

There are of course more specific schemes that are used to improve stability/speed/accuracy when particular forms are encountered in the PDE. These rules are applied first to override the general ruleset.

##Currently implemented special cases are as follows: - Spherical derivatives - Nonlinear laplacian uses a half offset centered scheme for the inner derivative to improve stability - Spherical nonlinear laplacian. - Upwind schemes to be used for odd ordered derivatives multiplied by a coefficient, downwinding when the coefficient is positive, and upwinding when the coefficient is negative.

Please submit an issue if you know of any special cases which impact stability or accuracy that are not implemented, with links to papers and/or code that demonstrates the special case.

generate_finite_difference_rules for SpectralDifferentialDiscretizer.

Emits a substitution rule for every Differential(x)^d)(u) appearing, plus a rule for each term headed by a spatial derivative whose argument is not a plain dependent variable call (nested derivatives such as Dx(u * Dx(u)), mixed derivatives Dx(Dy(u)), and boundary-value derivatives like Dx(u(t, 0))).

Get the weights and stencil for the inner half offset centered difference for the nonlinear laplacian for a given index and differentiating variable.

Does not discretize so that the weights can be used in a replacement rule TODO: consider refactoring this to harmonize with centered difference

Each index corresponds to the weights and index for the derivative at index i+1/2

Creates a ranking of the variables in the term, based on their derivative order. The heuristic that should work is, if there's a time derivative then use that variable, otherwise use the highest derivative for that variable. If there are two with the highest derivative, pick first from the list that hasn't been chosen for another equation

interface_layout_compatible(s, b, j)

Whether interface b may write a CartesianIndex of b.u into b.u2.

The pointwise wrap and face equations index the partner at the same slot j plus a shift along that axis. That is valid only when b.x2 occupies argument position j in b.u2, the two variables have the same number of spatial arguments, and every non-interface axis has the same discrete length. A different layout is not remapped.

ivs(u, s::DiscreteSpace)

Filter out the time variable and get the spatial variables of u in s.

Match the five @rule patterns of generate_nonlinlap_rules against the additive terms, keeping the last match per term (mirrors the pointwise path's Dict semantics). Grid-varying prefactors throw: the pointwise path leaves them undiscretized, so no slice form can reproduce them.

Match the three @rule patterns of generate_spherical_diffusion_rules against the additive terms, keeping the last match per term. Terms carrying a nonlinear laplacian match are skipped: the pointwise path keys both rulesets by the same term and the nonlinear laplacian entry, added later, wins its rules Dict. Grid-varying prefactors throw for the same reason as in match_nonlinlap_terms.

Performs a mixed centered difference in x centered at index II of u ufunc is a function that returns the correct discretization indexed at Itap, it is designed this way to allow for central differences of arbitrary expressions which may be needed in some schemes

The substitution key for the two-direction mixed derivative (Differential(x)^m * Differential(y)^n)(u).

This is the form both generate_mixed_rules and the array path match. The historical m = n = 1 spelling (Differential(x) * Differential(y))(u) is isequal to Differential(x)^1 * Differential(y)^1 in current Symbolics; mixed_derivative_keys still emits both spellings so a mismatch cannot drop a (1, 1) term.

Substitution keys for (Differential(x)^m * Differential(y)^n)(u), including the unpowered (Differential(x) * Differential(y))(u) spelling when m = n = 1.

Map each spatial direction to the centered mixed-derivative orders that reach along it.

Check if term is a compatible part of a nonlinear laplacian, including spherical laplacian, and return the argument to the innermost derivative if it is.

Derivative orders above one that the coefficient applies along m.x; order one is always contributed by the inner derivative. Accepts NonlinlapMatch and SphericalMatch.

Gets the parameter symbols of the system

Evaluate ex at the grid point II by substituting dependent variables with their discrete counterparts and independent variables with their grid values. Rules that replace a whole matched term must pass every captured factor through this, as the general variable/grid substitution rules never reach a rule's output.

spectral_apply(D, inner, II, s, derivweights, indexmap, x)

Differential(x)^d(inner) at II as a differentiation-matrix row dot. When inner is a dependent variable call the literal arguments pin the row (so Dx(u(t, a)) is the derivative at the boundary), otherwise inner is recursively evaluated at every tap of the row.

mat applied along axis j of the symbolic array X, as a SpectralApply term.

spectral_array_derivative(term, ranges, c::SpectralArrayContext)

(Differential(x)^d)(inner) over ranges as the differentiation block for the rows of ranges along x applied to inner evaluated on the full tap range of that direction. A dependent variable call with a literal argument pins that axis: a literal in the x slot selects the boundary row (Dx(u(t, a, y))), a literal elsewhere restricts the slice (Dx(u(t, x, b))). Terms that do not vary along x differentiate to zero.

spectral_arrayify(ex, ranges, c::SpectralArrayContext)

The slice form of ex over the box ranges (equation axis => index range). Every outermost spatial-derivative subterm becomes a SpectralApply term via spectral_array_derivative; everything else is handled by arrayify with the same slice, boundary-value, time-literal and grid rules the finite difference array form uses.

Whether ex varies along x: x itself appears, or a dependent variable is called with x (rather than a literal) in that slot.

spectral_diff_matrix(spec, grid, order)

Differentiation matrix of order order for the collocation scheme spec on grid. Fourier directions use the explicit trigonometric-interpolant first derivative matrix (Trefethen, Spectral Methods in MATLAB) raised to order. Chebyshev–Lobatto directions use the Weideman–Reddy chebdif recursion, which builds every order directly with the negative-sum trick and stays finite for thousands of nodes. Custom grids use the maximal-stencil Fornberg weights, i.e. the polynomial-interpolant differentiation matrix, which overflow beyond roughly a thousand nodes.

spectral_eval(ex, II, s, derivweights, indexmap)

Evaluate the expression ex at grid point II under spectral discretization. Spatial differentials become differentiation-matrix row dots; dependent variable calls become the matching discrete unknown (boundary literals resolve to the boundary index via newindex); independent variables become their grid value. This is the collocation analogue of the substitution rules used by the finite difference path, and transparently handles nested derivative terms such as Dx(a(u) * Dx(u)) by evaluating the argument at every tap of the outer derivative row.

spectral_grid(spec, a, b)

The collocation grid for spec on the domain [a, b].

spectral_taps_rowmap(spec, n)

The tapped grid indices and the grid-index-to-matrix-row map for spec on a grid of n points.

spherical_diffusion

for terms of the form r^-2*Dr(r^2*Dr(u(t, r)))

Based on https://web.mit.edu/braatzgroup/analysisoffinitedifferencediscretizationschemesfordiffusioninsphereswithvariablediffusivity.pdf

See scheme 1 in appendix A. The r = 0 case is treated in a later appendix

Substitution rules for independent variables that are not the derivative direction x. Values are taken from the argument-ordered index II via x2i, never from the hash-order-dependent s.x̄ permutation.

upwind_difference

Generate a finite difference expression in u using the upwind difference at point II::CartesianIndex in the direction of x

Patterns the staggered pointwise path cannot discretize either; falling back keeps this strategy's behaviour identical to the pointwise form for them.

Implements the WENO scheme of Jiang and Shu. Specified in https://repository.library.brown.edu/studio/item/bdr:297524/PDF/ (Page 8-9) Implementation heavily inspired by https://github.com/ranocha/HyperbolicDiffEq.jl/blob/84c2d882e0c8956457c7d662bf7f18e3c27cfa3d/src/finitevolumes/wenojiang_shu.jl by H. Ranocha.

The index range a tap slice takes across a wrapping seam, mirroring _wrapinterface: indices at or below the first point come from the lower-end destination, indices past the last point from the upper-end destination. A range that straddles the seam is not a slice.

Returns (dest, range) where dest is nothing (same array, unshifted interior), :self (same array, remapped across a periodic seam), or the partner named tuple of a two-domain interface. wrap_periodic_range is the self-periodic special case.

Replace the PDESystem with an equivalent PDESystem which is compatible with MethodOfLines, mutates boundarymap and v

Modified copilot explanation: