Advection–diffusion operator

AdvectionDiffusion{T,C}

Flux-form advection-diffusion operator on a KineticGrid,

\[\begin{aligned} \frac{\partial f}{\partial t} &= -\frac{1}{\mathcal{V}'} \left[\partial_p (\mathcal{V}' \Gamma^p) + \partial_\xi (\mathcal{V}' \Gamma^\xi) + \partial_r (\mathcal{V}' \Gamma^r)\right], \\ \Gamma^p &= A^p f - D^{pp} \partial_p f - D^{p\xi} \partial_\xi f, \\ \Gamma^\xi &= A^\xi f - D^{\xi\xi} \partial_\xi f - D^{\xi p} \partial_p f, \\ \Gamma^r &= A^r f - D^{rr} \partial_r f, \end{aligned}\]

with coefficients given on the faces of each direction (advection_p is $A^p$, diffusion_p_xi is $D^{p\xi}$, and so on). With p in units of mₑc, the coefficients in p and ξ are in s⁻¹, advection_r in m/s and diffusion_r_r in m²/s. Particle conservation is exact up to the fluxes through the outer boundaries.

The pitch grid may differ between radii, with the same cell counts and the same p faces at all radii. The pitch cells of a radial face are then the overlap cells of the pitch grids of its two radii (radial_face_overlap), and the radial coefficients are given per overlap cell: advection_r[i, s, k] for momentum cell i and overlap cell s of face k, at the pitch radial_face_overlap(grid, k).xi_centers[s]. Their second dimension is the largest count over the faces, and entries beyond the count of a face are not used. The radial flux is computed once per overlap cell, from the values of the two cells it lies in, and given to both with its share of each cell, so no particle is lost or created between different pitch grids. Both values count as constant over their cells. Where the pitch grids of two radii differ, the difference of the two values therefore contains a pitch difference of the order of the cell width, which acts as a diffusion in the pitch of the order Dʳʳ (Δξ/Δr)² times the fraction of the cells that do not match. With one pitch grid for all radii the overlap cells are the pitch cells and nothing changes.

Boundary conditions:

• p = 0 (p_zero_boundary): :outflow, the flux (𝒱'/p²)|₀ · (p²Aᵖ)|₀ · f₁, with the limit p²Aᵖ given in op.advection_p_times_p2_at_zero[j, ir] (nonzero for collisional friction, which removes particles through p = 0 as DREAM's hot-tail grid does). With :quick and :upwind_second_order the scheme's own face value at p = 0 replaces f₁, as in DREAM's AdvectionTerm; the limiters give f₁ there (mirrored ghost cells). Or :zero_flux, no flux through p = 0, as on DREAM's runaway grid: DREAM multiplies the limit (src/Equations/Kinetic/SlowingDownTerm.cpp:54-56) by the interpolation coefficients of the face p = 0 (fvm/Equation/AdvectionTerm.set.cpp:98-118), which its default lower condition AD_BC_DIRICHLET (fvm/Equation/AdvectionTerm.cpp:36) sets to zero (fvm/Equation/AdvectionInterpolationCoefficient.cpp:438-442). Only a distribution with an internal boundary condition switches to AD_BC_MIRRORED (src/Settings/Equations/distribution.cpp:234-237): f_hot has one (src/Settings/Equations/f_hot.cpp:96), f_re none (src/Settings/Equations/f_re.cpp:87). The ghost cells of the four-point schemes below p = 0 then hold zero, as DREAM's Dirichlet condition leaves out the cells beyond the grid. p_zero_boundary applies to a grid from p = 0 only, and :zero_flux needs p_min = 0; a grid with p_min > 0 takes p_lower_boundary.
• First radial face: no flux.
• p = p_min > 0 (p_lower_boundary): :zero_flux, or :prescribed with the flux Γᵖ through the face given in op.p_flux_at_lower_boundary[j, ir] [m⁻³ s⁻¹] (positive into the domain, independent of f), or :interface: the face is shared with the grid of another kinetic unknown below (a runaway grid above a hot-tail grid), the operator adds no flux through it, and a MomentumInterfaceTerm in the equation system carries the flux between the two grids. p_lower_boundary applies to a grid with p_min > 0 only: on a grid from p = 0, :prescribed and :interface are refused, and the boundary there is p_zero_boundary.
• ξ = ±1: no flux.
• p = p_max (p_upper_boundary): :outflow, advective outflow only (Aᵖ > 0); or :f_zero, outflow into a ghost cell with f = 0 one cell spacing beyond the last centre, adding diffusive loss Dᵖᵖ f_N / (p_N - p_{N-1}) (DREAM BC_F_0; with one momentum cell there is no spacing and the boundary is :outflow, which the pitch-only tests on one momentum cell use); or :constant_flux, 𝒱'Γᵖ at the boundary face held equal to 𝒱'Γᵖ at the last interior face (DREAM BC_PHI_CONST, its default; needs at least 2 momentum cells); or :extrapolated_flux, 𝒱'Γᵖ at the boundary linearly extrapolated in p from the last two interior faces (DREAM BC_DPHI_CONST, "known to cause instabilities near p_max"; needs at least 3 momentum cells, no :exponential_fitting counterpart). Ported literally from DREAM PXiExternalLoss.cpp, including an apparent index inconsistency: both terms share the last interior face's 𝒱' and are multiplied by that face's own Aᵖ, while the second term's interpolation (upstream direction and weights) and diffusion coefficient are the second-to-last face's own; or :extrapolated_flux_consistent, 𝒱'Γᵖ linearly extrapolated with every face's own 𝒱' and Γᵖ throughout, matching DREAM's stated definition of BC_DPHI_CONST ("dΦ/dp constant", Φ = 𝒱'Γᵖ) rather than its code, which besides the Aᵖ index above also shares one 𝒱' between the two terms (needs at least 3 momentum cells). DREAM comparisons keep :extrapolated_flux, so they stay exact to round-off; this boundary follows the stated definition instead (see the docs page Comparison with DREAM). Or :interface, the counterpart of the lower boundary: the face is shared with the grid of a kinetic unknown above, and the MomentumInterfaceTerm carries the flux. The ghost cells of the four-point schemes beyond an interface hold zero, as beyond any other boundary in p, so the last interior face on each side sees the other grid as empty.
• r = a (r_boundary): :zero_flux, or :absorbing (f = 0 on the edge face: the diffusive flux is Dʳʳ f_N over the distance from the last cell centre to the face), or :absorbing_dream, DREAM's BC_F_0 of its transport terms, with the distance between the last two cell centres, which places f = 0 half a cell outside the edge (1 with one radius, as in DREAM). On a uniform grid the edge flux of :absorbing is twice that of :absorbing_dream. The lowest decay rate of a uniform Dʳʳ in a cylinder converges at second order with :absorbing and at first order with DREAM's form: 0.997 and 0.911 of the exact rate at 11 radii. The advective part is the upwind outflow in all of them.

Schemes (scheme), for the advection and diffusion along each direction between two cell centres; the cross terms Dᵖᵠ, Dᵠᵖ are always central:

• :central: face value interpolated linearly between the centres, second order. Loses positivity where the cell Péclet number |A|Δ/D exceeds 2 (the matrix gets negative off-diagonal couplings).
• :upwind: face value of the upstream cell, first order, positivity preserving.
• :exponential_fitting: the Scharfetter-Gummel (Chang-Cooper) flux Γ = (D/Δ) [B(-w) f_L - B(w) f_R], w = AΔ/D, B(x) = x/(eˣ - 1), with Δ the distance between the centres. Exact for steady solutions of the one-dimensional face problem with constant A/D between the centres (so the collisional Maxwellian is steady up to the variation of A/D over a cell), positivity preserving for every Péclet number (the off-diagonal couplings are (D/Δ) B(±w) ≥ 0), second order where the Péclet number is small, where it reduces to the central flux with midpoint weights, and upwind where it is large. Not in DREAM (Runes2P departure); used in CQL3D (Chang-Cooper). Falls back to upwind advection with central diffusion where D ≤ 0. At the :f_zero boundary the ghost cell enters the same flux.
• Four-point schemes (in p and ξ, diffusion central): DREAM's :quick and :upwind_second_order (linear in f) and flux limiters :tcdf, :muscl, :smart, :ospre, plus :van_albada (Runes2P), which are nonlinear in f. The face value uses two upstream cells and one downstream cell (src/operators/flux_limiters.jl), with ghost cells as DREAM sets them: mirrored at p = 0 with p_zero_boundary = :outflow (DREAM's hot-tail distribution), zero with p_zero_boundary = :zero_flux (DREAM's runaway distribution) and at a lower boundary pmin > 0, zero beyond pmax, mirrored at ξ = ±1. The fluxes through the boundary faces keep the rules listed above (the scheme's face value at p = 0, upwind at pmax as DREAM's `BCF_0`). With fewer than three cells in a direction the face value is central, as in DREAM.

The scheme can differ per direction (p_scheme, xi_scheme, r_scheme). scheme sets p_scheme and xi_scheme, and r_scheme when it is a two-point scheme. With a four-point scheme the radial scheme is :central unless r_scheme names one.

A four-point scheme in r gives the face value of the radial advection. The radial diffusion is central. Ghost cells as in DREAM: mirrored at the first face, zero beyond the edge of the plasma. The flux through the edge is the upwind outflow. It needs one pitch grid at all radii, and with fewer than three radii the face value is central.

The operator declares which coefficients beyond those of momentum space it can carry, and its Jacobian pattern (stencil_pattern) holds only the cells they couple: radial_transport (advection_r, diffusion_r_r) the radial neighbours, cross_diffusion (diffusion_p_xi, diffusion_xi_p) the corners of the 3×3 block in p and ξ. Without them the radii of a Fokker-Planck operator are not coupled, and a direct factorization treats the radii as separate blocks. A nonzero coefficient that the operator does not declare is refused where the coefficients are checked: the operator's place in an EquationSystem, a rebuild by a FokkerPlanckTerm, a RechesterRosenbluthUpdate, operator_matrix and kinetic_conductivity. In the coloring of a Jacobian with that pattern its entries would otherwise be added to declared entries without an error. A write during a run by an IntervalUpdate is not checked unless the update checks it. A FokkerPlanckTerm with a plasma source clears both declarations of its operator, since the rebuild of its coefficients writes neither. A declaration may be cleared after the operator is built but never set: the pattern of an EquationSystem is fixed when the system is built, and a declaration set later would not reach it.

AdvectionDiffusion(grid; scheme=:upwind, p_scheme=scheme, xi_scheme=scheme,
                   r_scheme=(scheme is four-point ? :central : scheme),
                   p_lower_boundary=:zero_flux, p_upper_boundary=:outflow,
                   p_zero_boundary=:outflow, r_boundary=:zero_flux,
                   radial_transport=true, cross_diffusion=true)

Operator with all coefficients zero; set them in place, e.g. op.diffusion_p_p .= 1. radial_transport = false declares that the operator carries no advection_r and diffusion_r_r, and cross_diffusion = false that it carries no diffusion_p_xi and diffusion_xi_p: its stencil then leaves out the radial neighbours or the corners of the 3×3 block in p and ξ (stencil_pattern). Both are declared by default, so that an operator filled by hand may hold every coefficient.

apply!(df, op, f; scale=1, limiter_state=nothing, electric_field=nothing)

Add scale times the operator applied to f to df (in place, allocation free). With electric_field (one value per radius, in V/m) the advection of every face is that of the operator plus the field times the coefficients of a unit field (set_unit_field_coefficients!): the field as an argument, for an electric field that is an unknown of the system. The operator is linear in f except for a prescribed lower-p flux, which adds a constant source, and for the flux limiters, whose face values depend on f. With limiter_state, the limiter weights are evaluated at that state and applied linearly to f: the limiter frozen at limiter_state (DREAM's linearly implicit solver takes them from the previous time level; the :frozen_limiter Jacobian of FokkerPlanckTerm).

momentum_flux(op::AdvectionDiffusion, f, face, ir=1; electric_field=nothing)

Particle flux per unit volume [m⁻³ s⁻¹] through the momentum face face (between cells face - 1 and face) at radius index ir, (1/V') Σ_j Δξ_j 𝒱'ᵖ Γᵖ, positive towards larger momentum. This is the rate at which particles cross the surface p = p_faces[face], counted in the same normalization as density, so a runaway rate measured through a surface above the critical momentum is directly comparable with connor_hastie_rate. Face 1 is the flux through p = p_min (the p = 0 limit or the prescribed boundary flux) and face momentum_cell_count + 1 the flux out of the grid. On a grid from p = 0, face 1 under p_zero_boundary = :zero_flux reads zero: this operator takes no flux through p = 0 there. A face with an :interface boundary (face 1 under p_lower_boundary = :interface, the last face under p_upper_boundary = :interface) reads zero: this operator carries no flux there, and interface_flux gives the flux between the grids. electric_field (one value per radius, in V/m) is the field of an operator that takes it as an argument (apply!).

with_scheme(op::AdvectionDiffusion, scheme; p_scheme=scheme, xi_scheme=scheme)

The operator op with other interpolation schemes in p and ξ, sharing every coefficient array, boundary value, boundary condition and the declarations of radial transport and cross diffusion with op (a rebuild of op's coefficients, set_fokker_planck_coefficients!, is seen by both). The radial scheme is unchanged. Used by the Jacobian modes of FokkerPlanckTerm: the residual of one scheme with the stored Jacobian of another.

_flux_eltype(op::AdvectionDiffusion, f, field=nothing)

The element type of the fluxes of op applied to f with the electric field field (one value per radius, nothing for none): that of f promoted with those of the coefficients and of the field, so that every branch of a flux returns one type. With dual coefficients (the derivatives with respect to the plasma values) and a Float64 f, the element type of f alone would make the fluxes return Union{Float64, Dual}, and apply! would then return wrong partials under exponential fitting at the default optimization level, though not at -O0, with 23 partials on Julia 1.13.1.

apply_radius!(df, op, f, ir, scale)

Operator contribution to all momentum cells at radius index ir. Writes only to df[:, :, ir], so different radii can be processed concurrently. radial_mode is the computation of the radial fluxes (_radial_flux_mode), a function barrier: the loop over the cells is compiled for it.

Each p face of a pitch row is evaluated once: the flux through the face above a cell is carried to the next cell as the flux through its face below, also past a cell with no phase-space volume. The exception is the last one or two interior p faces, which the boundary at pmax evaluates again when it continues or extrapolates their flux (`:constantflux,:extrapolatedflux,:extrapolatedflux_consistent`). The ξ faces are evaluated from both of their cells, since a buffer of them would need storage per call or in the operator. Each r face is evaluated from the radius on each side of it: with per-radius pitch grids the two sides give different values, and a buffer across radii would tie a radius to the one before it.

similar_coefficients(op::AdvectionDiffusion, ::Type{C})

An operator with the grid, the schemes, the boundary conditions and the declarations of radial transport and cross diffusion of op, and zero coefficients of the number type C: coefficients that hold dual numbers, whose application to a distribution gives the derivative of the residual with respect to the plasma values the coefficients are built from, through the scheme's own face weights.

Schemes whose advective face value uses four cells: DREAM's :quick and :upwind_second_order (linear in f) and the flux limiters :tcdf, :muscl, :smart, :ospre (DREAM) and :van_albada (Runes2P), which are nonlinear in f.

_four_point_face_value(scheme, x_face, x_d, x_u, x_uu, y_d, y_u, y_uu)

Face value from the three cells. For the limiters the weights are evaluated on the values alone and the derivative is carried as Σ c_k (y_k - value(y_k)) with the exact coefficients of _four_point_weights: differentiating the ratio of slopes directly would overflow in the far tail of a Maxwellian, where the slopes are subnormal, although the derivative is bounded (DREAM writes its full Jacobian the same way).

_four_point_weights(scheme, x_face, x_d, x_u, x_uu, v_d, v_u, v_uu)

Interpolation weights (w_d, w_u, w_uu) of the face value at x_face on the downstream cell (x_d), the upstream cell (x_u) and the next upstream cell (x_uu), in DREAM's form for nonuniform grids (DREAM's deltas). For the limiters they depend on the cell values v_d, v_u, v_uu (plain numbers) through the slope ratio r; also returned are the derivatives (c_d, c_u, c_uu) of the face value with respect to the three values (ψ' of the limiter; zero where r is clamped). For the linear schemes both triples are equal.