Jacobians

SparseJacobian(residual!, grid; backend=sparse_ad_backend(pattern=stencil_pattern(grid)))

Prepared sparse Jacobian of an in-place residual residual!(dF, f) acting on arrays shaped like grid. Sparsity pattern and coloring are computed once; jacobian!(J, sparse_jacobian, f) then re-evaluates the values.

residual! must contain only stencil-local terms. Terms with global momentum coupling, such as the knock-on operator, make the Jacobian dense per radius and are applied matrix-free in the Krylov solver instead. The default pattern has width 1; for an operator with four-point schemes pass backend=sparse_ad_backend(pattern=stencil_pattern(op)), or the coloring folds the second neighbours silently onto the pattern.

check_stencil_pattern(op)

Return true if every entry of the Jacobian of op lies inside stencil_pattern(op), which holds the radial neighbours and the corners of the 3×3 block only when op declares radial transport and cross diffusion: an operator that holds a coefficient it does not declare gives false. A linear operator is traced on its grid; one with a scheme that is nonlinear in f gets its dense Jacobian by ForwardDiff at a random positive state (tracing cannot follow the branches of a limiter). Slow; intended for tests of local operators on small grids.

is_linear(op)
is_linear(term)

Whether the operator is linear in f (up to a constant source); false for flux-limiter schemes, whose face values depend on f. For an EquationTerm: whether the term is linear in the state up to a constant (false by default; terms opt in), which lets NewtonSolver accept a step of a linear system on the update test alone.

is_local(op)

Whether op couples each cell only to its stencil_pattern neighbours. Only local operators may be assembled with the stencil pattern: with a known pattern, column coloring silently adds entries outside the pattern onto entries inside it. Operators default to non-local; local operator types opt in.

jacobian!(J, sparse_jacobian, f)

Overwrite the values of J (with the pattern from jacobian(sparse_jacobian, f)) at f.

jacobian(sparse_jacobian, f)

Sparse Jacobian of sparse_jacobian.residual! at f (new matrix).

operator_matrix(op; backend=sparse_ad_backend(pattern=stencil_pattern(op)), check_local=true)

Sparse matrix M with vec(apply(op, f)) = M * vec(f) + c for an operator op that is linear in f up to a constant source c. The default pattern follows the operator's schemes, boundary conditions and declarations (stencil_pattern), and an operator that holds a coefficient it does not declare is refused with it (AdvectionDiffusion). check_local = true refuses operators that are not local (is_local); pass false only with a backend that detects the pattern itself.

sparse_ad_backend(; pattern=nothing, tag=ForwardDiff.Tag{OperatorJacobianTag,Float64}(),
                  chunk_limit=jacobian_chunk_limit())

ForwardDiff with greedy column coloring. The sparsity pattern is taken from pattern when given, otherwise detected by operator tracing (slow for large grids). The number of Jacobian-vector products equals the number of colors, which stays bounded as the grid is refined for local stencils. tag is the tag of the dual numbers, one type for every function that is differentiated (OperatorJacobianTag). tag = nothing takes ForwardDiff's own tag, which is needed for a function that itself takes a Jacobian with this backend of data that depend on its argument. chunk_limit caps the chunk size of the dual numbers for a given pattern (jacobian_chunk_limit).

stencil_pattern(grid; p_width=1, xi_width=1, r_width=1, p_upper_boundary=:f_zero,
                radial_transport=true, cross_diffusion=true)
stencil_pattern(op::AdvectionDiffusion)

Sparsity pattern of the local transport operators: the 3×3 (p, ξ) neighbourhood, the cells up to p_width away in p and xi_width away in ξ along the axes (2 for the four-point schemes), and with radial_transport the radial neighbours of each cell (the cell of the same pitch at the two neighbouring radii, or with different pitch grids the cells that overlap it) and the cells up to r_width radii away. The corners of the 3×3 block, which only the p-ξ cross diffusion couples, are part of it with cross_diffusion. With p_upper_boundary = :extrapolated_flux or :extrapolated_flux_consistent the last momentum row is additionally coupled up to 3 cells away in p along j (the extrapolated flux uses the last two interior faces), and, for :extrapolated_flux_consistent with cross_diffusion, diagonally to (np - 2, j ± 1) through its second term's p-ξ cross contribution (:extrapolated_flux's second term has none).

For an operator, the widths follow its schemes, and the boundary at p_max and the declarations of radial transport and cross diffusion are its own (AdvectionDiffusion).

JacobianTag

Tag of the dual numbers of the automatic Jacobian of an EquationSystem. ForwardDiff's default tag holds the type of the differentiated function, which holds the terms of the system. With it every system would have dual numbers of a type of its own, and the right-hand side of every term would be compiled again for each system: about 120 s per system in the tests. With this tag the dual numbers of all systems have one type, and the code of a term is compiled once.

OperatorJacobianTag

Tag of the dual numbers of the Jacobians of single operators and residual functions (sparse_ad_backend, SparseJacobian, operator_matrix, check_stencil_pattern), for the reason of JacobianTag. It differs from the tag of the systems, so a Jacobian of an operator that a term of a system takes inside the automatic Jacobian of that system keeps its dual numbers apart.

jacobian_chunk_limit()

The upper limit of the chunk size of the dual numbers of the sparse Jacobians, from the environment variable RUNES2P_JACOBIAN_CHUNK_SIZE, or nothing without it.

The chunk size is the number of colors that one evaluation of the residual carries. ForwardDiff takes up to 12 colors in one chunk, so systems with 9, 10 and 11 colors have dual numbers of three types, and the code of the terms is compiled once for each. A limit of 8, the setting for the test suite, gives every pattern with 8 colors or more the chunk size 8, one type of dual numbers and one compiled copy of the code, at the price of a second evaluation for a pattern of 9 to 12 colors (ForwardDiff splits 13 colors or more into two batches in any case). The values of the Jacobian do not depend on the chunk size. A run without the variable keeps ForwardDiff's choice, which needs the fewest evaluations.