10. Higher-Order Temporal Integration and Splitting#

This chapter explains how the baseline first-order pressure-based finite-volume algorithm is extended to higher-order time accuracy, and how that temporal update interacts with source splitting and Strang splitting. The finite-volume transport chapter focused on how the spatial terms are assembled into algebraic rows. The present chapter asks a different question: once those spatial rows are known, how are they advanced from time level \(n\) to time level \(n+1\) in a way that is formally higher-order in time?

The discussion is organized around three ideas:

  1. one-step second-order time integration, represented here by Crank-Nicolson,

  2. multistep second-order time integration, represented by BDF2,

  3. operator/source splitting, including first-order source splitting and second-order Strang splitting for reacting flows.

Throughout the chapter it is useful to distinguish time levels from outer-iteration states.

  • Superscripts \(n-1\), \(n\), and \(n+1\) denote physical time levels.

  • Superscripts \(k\) and \(k+1\) denote nonlinear or segregated outer iterations at the same new time level.

  • A superscript \(*\) denotes a provisional state available during the predictor stage of the pressure-based algorithm.

  • A prime denotes a correction variable.

That distinction is important. A scheme may be one-step in time and still require several outer iterations at time level \(n+1\). Those outer iterations do not introduce extra physical time levels; they are only the nonlinear solution procedure used to converge the equations at the current time step.

10.1. Semi-Discrete Starting Point#

After spatial discretization, the control-volume equation for one transported quantity may be written schematically as

(10.1.1)#\[\frac{d Q_P}{dt} + \mathcal{R}_P(Q) = 0,\]

where \(Q_P\) is the cell-integrated unknown and \(\mathcal{R}_P\) is the spatial residual assembled from the convective, diffusive, pressure, and source contributions described in Finite-Volume Transport Assembly on Unstructured Grids. For example,

(10.1.2)#\[Q_P = (\rho \phi)_P\,\Omega_P\]

for a generic transported scalar, while for continuity one may take

(10.1.3)#\[Q_P = \rho_P\,\Omega_P.\]

Equation (10.1.1) is the common starting point for the time integrators below. The difference among the schemes is therefore not in the spatial operator itself, but in how the time derivative is approximated over the interval \([t^n,t^{n+1}]\).

10.2. Crank-Nicolson: One-Step Second-Order Time Integration#

Crank-Nicolson is the classical trapezoidal-rule update applied to the semi-discrete equation. For a generic ODE,

(10.2.1)#\[\frac{dy}{dt}=f(y,t),\]

the Crank-Nicolson update is

(10.2.2)#\[y^{n+1}=y^n+\frac{\Delta t}{2}\left[f^{n+1}+f^n\right].\]

The important point for readers is that this is a one-step method. It uses the two endpoints of the current time interval, \(t^n\) and \(t^{n+1}\), and does not require the older value at \(t^{n-1}\). Its second-order accuracy comes from trapezoidal averaging over the current interval, not from using a longer time history.

Applied to the semi-discrete finite-volume equation in Eq. (10.1.1), the Crank-Nicolson update becomes

(10.2.3)#\[\frac{Q_P^{n+1}-Q_P^n}{\Delta t} + \frac{1}{2}\left[\mathcal{R}_P(Q^{n+1})+\mathcal{R}_P(Q^n)\right] =0.\]

Equivalently,

(10.2.4)#\[\frac{Q_P^{n+1}-Q_P^n}{\Delta t} + \frac{1}{2}\mathcal{R}_P(Q^{n+1}) = -\frac{1}{2}\mathcal{R}_P(Q^n).\]

This form makes the solver interpretation transparent:

  • the new-time half of the spatial residual is treated implicitly,

  • the old-time half of the spatial residual is known and therefore goes to the right-hand side.

If the spatial operator for the new-time transported-variable row is written in the form derived in Finite-Volume Transport Assembly on Unstructured Grids,

(10.2.5)#\[\mathcal{R}_P(Q^{n+1}) \equiv a_{P,\mathrm{sp}}^{n+1}\,\phi_P^{n+1} - \sum_{\text{nbr}} a_{\text{nbr},\mathrm{sp}}^{n+1}\,\phi_{\text{nbr}}^{n+1} - b_{P,\mathrm{sp}}^{n+1},\]

then Crank-Nicolson leads schematically to the algebraic row

(10.2.6)#\[\left[ \frac{\rho_P^n\Omega_P}{\Delta t} + \frac{1}{2}a_{P,\mathrm{sp}}^{n+1} \right]\phi_P^{n+1} = \sum_{\text{nbr}} \frac{1}{2}a_{\text{nbr},\mathrm{sp}}^{n+1}\,\phi_{\text{nbr}}^{n+1} + b_{P,\mathrm{CN}},\]

where \(b_{P,\mathrm{CN}}\) collects the known old-time transient state together with the old-time half of the spatial operator.

10.2.1. Why Crank-Nicolson Does Not Need \(n-1\)#

The reason no \(n-1\) value appears is simply that Crank-Nicolson is a one-step endpoint average. It approximates the time integral over the single interval \([t^n,t^{n+1}]\) by using the residual at the beginning and the end of that interval. Nothing from an earlier interval is required.

This is worth stating explicitly because the notation can otherwise be misleading. In a segregated pressure-based solver, one still solves for intermediate iterates such as

(10.2.1.1)#\[Q^{n+1,k} \longrightarrow Q^{n+1,k+1},\]

but those are not extra time levels. They are only nonlinear iterations used to converge the Crank-Nicolson update at the same physical time \(t^{n+1}\).

10.2.2. Practical Interpretation#

Crank-Nicolson is attractive because it is second-order accurate and less temporally dissipative than first-order BDF. The trade-off is that, because it averages old and new residuals, it can be less strongly damping than BDF-based schemes. In stiff problems or in flows with sharp transients, many solvers therefore blend Crank-Nicolson with a more dissipative first-order form for robustness. That blending does not change the main conceptual point made here: Crank-Nicolson is second-order and one-step, so it advances from \(n\) to \(n+1\) without needing \(n-1\).

10.3. Second-Order Backward Differencing (BDF2)#

BDF2 is also second-order accurate, but it is a multistep scheme rather than a one-step scheme. For a generic transported scalar \(\phi\), the exact BDF2 approximation of the time derivative is

(10.3.1)#\[\left.\frac{\partial(\rho\phi)}{\partial t}\right|^{n+1} \approx \frac{3(\rho\phi)^{n+1}-4(\rho\phi)^n+(\rho\phi)^{n-1}}{2\Delta t} .\]

For continuity,

(10.3.2)#\[\left.\frac{\partial\rho}{\partial t}\right|^{n+1} \approx \frac{3\rho^{n+1}-4\rho^n+\rho^{n-1}}{2\Delta t} .\]

So, unlike Crank-Nicolson, BDF2 does require the older time level \(n-1\). That is the defining difference between a one-step and a multistep second-order scheme in this chapter.

Applied to Eq. (10.1.1), the BDF2 update is

(10.3.3)#\[\frac{3Q_P^{n+1}-4Q_P^n+Q_P^{n-1}}{2\Delta t} + \mathcal{R}_P(Q^{n+1}) =0.\]

If the new-time spatial operator is assembled in the standard transported-row form of Finite-Volume Transport Assembly on Unstructured Grids, the corresponding algebraic row is schematically

(10.3.4)#\[\left[ \frac{3\rho_P^n\Omega_P}{2\Delta t} + a_{P,\mathrm{sp}}^{n+1} \right]\phi_P^{n+1} = \sum_{\text{nbr}} a_{\text{nbr},\mathrm{sp}}^{n+1}\,\phi_{\text{nbr}}^{n+1} + b_{P,\mathrm{BDF2}},\]

where \(b_{P,\mathrm{BDF2}}\) contains the history contribution built from the \(n\) and \(n-1\) states together with the explicit spatial/source terms.

10.3.1. Solver-Friendly BDF2 Split#

In segregated implementations it is often convenient to view BDF2 as a first-order backbone plus a second-order correction. The exact identity is

(10.3.1.1)#\[\frac{3X^{n+1}-4X^n+X^{n-1}}{2\Delta t} = \frac{X^{n+1}-X^n}{\Delta t} + \frac{X^{n+1}-2X^n+X^{n-1}}{2\Delta t},\]

where \(X\) may stand for \(\rho\phi\), \(p\), or another transported quantity.

The first term on the right is the BDF1 backbone. The second term is the second-order correction. During one segregated outer iteration, that correction may be lagged using the current iterate at the new time level, for example

(10.3.1.2)#\[\frac{X^{n+1}-2X^n+X^{n-1}}{2\Delta t} \;\longrightarrow\; \frac{X^{n+1,*}-2X^n+X^{n-1}}{2\Delta t}.\]

This does not define a different time integrator. It is only a solver rearrangement. At outer-iteration convergence, \(X^{n+1,*}\rightarrow X^{n+1}\), and the exact BDF2 formula is recovered.

10.3.2. Split Transient Derivatives and Time-Integration Factors#

The BDF2 split above is useful beyond ordinary transported quantities. It also applies when the transient part of an equation contains a difference of two cell-integrated quantities, but only one of them is kept as the primary solved unknown. A simple model problem is

(10.3.2.1)#\[\frac{d}{dt}\left(Y-X\right) + \mathcal{R}(Y) = 0,\]

where \(Y\) is the quantity represented in the algebraic row and \(X\) is evaluated from the current lagged thermodynamic or pressure state. Rearranging Eq. (10.3.2.1) gives

(10.3.2.2)#\[\frac{dY}{dt} + \mathcal{R}(Y) = \frac{dX}{dt}.\]

This rearrangement is the whole reason that a time-derivative term can appear on the right-hand side. It is not a different physical model; it is the result of choosing \(Y\) as the solved variable while the conservative transient quantity is \(Y-X\).

For first-order BDF over a time step \(\Delta t_n\), Eq. (10.3.2.2) gives

(10.3.2.3)#\[\frac{Y^{n+1}-Y^n}{\Delta t_n} + \mathcal{R}(Y^{n+1}) = \frac{X^{n+1,*}-X^n}{\Delta t_n},\]

where the \(*\) indicates the current available new-time estimate of \(X\) during the nonlinear or segregated iteration. The factor multiplying the split right-hand-side difference is therefore simply

(10.3.2.4)#\[\alpha_X^{\mathrm{BDF1}} = \frac{1}{\Delta t_n}.\]

For BDF2 with variable time step, let

(10.3.2.5)#\[h = \Delta t_n = t^{n+1}-t^n, \qquad k = \Delta t_{n-1} = t^n-t^{n-1}.\]

The three-point backward derivative of \(X\) at \(t^{n+1}\) is

(10.3.2.6)#\[D_2 X^{n+1} = \frac{2h+k}{h(h+k)}X^{n+1} - \frac{h+k}{hk}X^n + \frac{h}{k(h+k)}X^{n-1}.\]

This can be split into a BDF1 backbone plus a correction:

(10.3.2.7)#\[\begin{aligned} D_2 X^{n+1} &= \frac{X^{n+1}-X^n}{h} + \left[ \frac{1}{h+k}X^{n+1} - \frac{1}{k}X^n + \frac{h}{k(h+k)}X^{n-1} \right]. \end{aligned}\]

Thus the right-hand-side derivative can be represented by one backbone factor and three correction factors:

(10.3.2.8)#\[\alpha_X^{\mathrm{BDF2,base}} = \frac{1}{h}, \qquad c_{F2} = \frac{1}{h+k}, \qquad c_{F1} = \frac{1}{k}, \qquad c_{F0} = \frac{h}{k(h+k)}.\]

For constant time step, \(h=k=\Delta t\), this reduces to

(10.3.2.9)#\[\frac{X^{n+1}-X^n}{\Delta t} + \frac{X^{n+1}-2X^n+X^{n-1}}{2\Delta t} = \frac{3X^{n+1}-4X^n+X^{n-1}}{2\Delta t}.\]

The important point is that the BDF2 right-hand-side derivative is not obtained by multiplying the BDF1 backbone by the full BDF2 current-time storage coefficient. The BDF1 backbone keeps the factor \(1/\Delta t_n\), and the second-order accuracy comes from the separate correction term.

The blended Crank-Nicolson/BDF form gives a different set of factors because it is a one-step residual average rather than a multistep history formula. Let \(\beta\) be the blend parameter, with \(\beta=0\) corresponding to BDF1 and \(\beta=1\) corresponding to pure Crank-Nicolson. Blending the BDF1 residual with the Crank-Nicolson residual gives

(10.3.2.10)#\[\frac{(Y-X)^{n+1}-(Y-X)^n}{\Delta t_n} + \left(1-\frac{\beta}{2}\right)\mathcal{R}^{n+1} + \frac{\beta}{2}\mathcal{R}^{n} =0.\]

It is convenient in the assembled row to normalize the new-time residual so that \(\mathcal{R}^{n+1}\) has the same coefficient as in the first-order row. Dividing by \(1-\beta/2\) gives

(10.3.2.11)#\[F_1 \frac{(Y-X)^{n+1}-(Y-X)^n}{\Delta t_n} + \mathcal{R}^{n+1} = - F_0\mathcal{R}^{n},\]

with

(10.3.2.12)#\[F_1 = \frac{1}{1-\beta/2} = \frac{2}{2-\beta}, \qquad F_0 = \frac{\beta/2}{1-\beta/2} = \frac{\beta}{2-\beta}.\]

After the split transient quantity is moved to the right-hand side, the corresponding \(X\) contribution is therefore

(10.3.2.13)#\[F_1 \frac{X^{n+1,*}-X^n}{\Delta t_n}.\]

For pure BDF1, \(\beta=0\) and \(F_1=1\). For pure Crank-Nicolson, \(\beta=1\) and \(F_1=2\). This is why the factor multiplying a split right-hand-side derivative is scheme dependent, even though the underlying model problem in Eq. (10.3.2.1) is the same.

10.3.3. Continuity-Added View for Scalar Rows#

For scalar rows assembled in the continuity-added form, the same continuum identity used in Finite-Volume Transport Assembly on Unstructured Grids remains helpful:

(10.3.3.1)#\[\frac{\partial(\rho\phi)}{\partial t} - \phi\,\frac{\partial\rho}{\partial t} = \rho\,\frac{\partial\phi}{\partial t}.\]

The role of this identity is unchanged in the higher-order setting. It allows the transient part of the transported-scalar equation to be written in a form that produces a clean scalar coefficient multiplying the new unknown while the remaining history-dependent pieces are carried on the right-hand side.

Note

At the first time step of a BDF2 run, the solver falls back to first-order BDF because the value at \(n-1\) is not yet available.

10.3.4. Pressure Time Derivative in Enthalpy/Temperature Forms#

The split-transient discussion above is the generic time-integration result. The total-enthalpy equation supplies one important specialization. As shown in the implemented-equations appendix, the conservative transient quantity in the total-enthalpy form is \(\rho H-p\). After integration over a control volume, the split quantity has the form

(10.3.4.1)#\[Y_P = (\rho H)_P\Omega_P, \qquad X_P = p_P\Omega_P.\]

Keeping \(H\) as the solved unknown moves the \(X_P=p_P\Omega_P\) time derivative to the right-hand side. For BDF2, Eq. (10.3.2.7) gives

(10.3.4.2)#\[\begin{split}\begin{aligned} D_2(p\Omega)_P^{n+1} &= \frac{(p\Omega)_P^{n+1}-(p\Omega)_P^n}{h} \\ &\quad + \left[ \frac{1}{h+k}(p\Omega)_P^{n+1} - \frac{1}{k}(p\Omega)_P^n + \frac{h}{k(h+k)}(p\Omega)_P^{n-1} \right]. \end{aligned}\end{split}\]

On a fixed grid with constant \(\Omega_P\), this reduces to the familiar pressure derivative multiplied by cell volume. On a moving or deforming grid, the cell-integrated quantity \(p\Omega\) should be kept intact.

For constant time step, the BDF2 pressure-volume derivative is

(10.3.4.3)#\[\frac{(p\Omega)_P^{n+1}-(p\Omega)_P^n}{\Delta t} + \frac{(p\Omega)_P^{n+1}-2(p\Omega)_P^n+(p\Omega)_P^{n-1}}{2\Delta t} = \frac{3(p\Omega)_P^{n+1}-4(p\Omega)_P^n+(p\Omega)_P^{n-1}}{2\Delta t}.\]

For blended Crank-Nicolson/BDF, Eq. (10.3.2.13) gives the right-hand-side pressure-volume contribution

(10.3.4.4)#\[F_1 \frac{(p\Omega)_P^{n+1,*}-(p\Omega)_P^n}{\Delta t_n}.\]

Thus the pressure-time contribution in total-enthalpy or temperature forms is not special because it is pressure. It is special because it is a split transient quantity placed on the right-hand side. BDF2 uses a first-order backbone plus a separate history correction. Blended Crank-Nicolson uses the one-step normalized factor \(F_1\).

10.4. Source Splitting versus Unsplit Treatment#

The higher-order time schemes above determine how one advances a semi-discrete equation in time once the operator has been chosen. A different question is whether that operator is advanced all at once or in parts.

Suppose the full evolution law is written schematically as

(10.4.1)#\[\frac{dU}{dt}=F_T(U)+F_S(U),\]

where \(F_T\) denotes the transport part and \(F_S\) denotes a stiff source operator such as finite-rate chemistry, particle coupling, or another local production/destruction mechanism.

Two broad strategies are then possible.

First, in an unsplit treatment, both parts remain in the same time-discrete equation and are advanced together within the same solve or nonlinear iteration loop.

Second, in a source-split treatment, the source operator is peeled off and advanced separately. The simplest sequential split is the first-order Lie or source-splitting composition,

(10.4.2)#\[U^{(a)}=\mathcal{S}_{\Delta t}U^n, \qquad U^{n+1}=\mathcal{T}_{\Delta t}U^{(a)},\]

or the reverse ordering,

(10.4.3)#\[U^{(a)}=\mathcal{T}_{\Delta t}U^n, \qquad U^{n+1}=\mathcal{S}_{\Delta t}U^{(a)}.\]

Here \(\mathcal{T}_{\Delta t}\) is the transport propagator and \(\mathcal{S}_{\Delta t}\) is the source propagator over a full time step. In general, this sequential splitting is only first-order accurate in time because the composition is not symmetric.

So, in the terminology of this chapter:

  • higher-order time integration answers the question “how is a chosen operator advanced in time?”

  • source splitting / operator splitting answers the different question “is that operator advanced as one block or as several sub-operators?”

Those two ideas are related but not identical.

10.5. Second-Order Strang Splitting for Reacting Flows#

Strang splitting is the symmetric second-order version of source splitting. It is especially useful when the source term is much stiffer than the transport operator, as is common in reacting flows.

Let the full flow state be denoted schematically by

(10.5.1)#\[U=(\rho,\rho \vec{u},\rho E,\rho Y_1,\ldots,\rho Y_{N_s}), \qquad \frac{dU}{dt}=F_T(U)+F_R(U),\]

where \(F_T\) is the transport operator and \(F_R\) is the reaction operator.

For the species equations, the reaction subproblem is

(10.5.2)#\[\frac{d(\rho Y_k)}{dt}=\dot{\omega}_k, \qquad k=1,\ldots,N_s, \qquad \sum_{k=1}^{N_s}\dot{\omega}_k=0,\]

and the transport subproblem is the same species transport equation with the chemistry source removed:

(10.5.3)#\[\frac{\partial(\rho Y_k)}{\partial t} +\nabla\cdot(\rho \vec{u}Y_k) = \nabla\cdot \vec{J}_k .\]

The chemistry subproblem is cell-local: there are no face fluxes in Eq. (10.5.2). The pressure-based momentum, pressure, and scalar transport solves all remain in the transport operator \(F_T\).

10.5.1. Symmetric Strang Composition#

Two common second-order symmetric compositions are

(10.5.1.1)#\[\mathcal{G}_{\Delta t}^{(TRT)} = \mathcal{T}_{\Delta t/2}\,\mathcal{R}_{\Delta t}\,\mathcal{T}_{\Delta t/2}, \qquad \mathcal{G}_{\Delta t}^{(RTR)} = \mathcal{R}_{\Delta t/2}\,\mathcal{T}_{\Delta t}\,\mathcal{R}_{\Delta t/2},\]

where \(\mathcal{T}\) and \(\mathcal{R}\) denote the transport and reaction propagators, respectively. Both are second-order accurate in time because the composition is symmetric[Str68, VHB98].

For reacting-flow calculations in stream, the reaction-centered view is usually the more useful one conceptually, because the chemistry solve produces an intermediate reacted state that then seeds the transport solve. Written in that form,

(10.5.1.2)#\[U^{(a)}=\mathcal{R}_{\Delta t/2}U^n, \qquad U^{(b)}=\mathcal{T}_{\Delta t}U^{(a)}, \qquad U^{n+1}=\mathcal{R}_{\Delta t/2}U^{(b)} .\]

The first half-reaction step produces a reacted intermediate state, the transport solve advances that reacted state over the flow time step, and the final half-reaction step restores the symmetry needed for second-order accuracy.

10.5.2. Why One-Step Transport Schemes Fit Naturally Inside Strang Splitting#

A one-step transport scheme such as Crank-Nicolson fits naturally into a split algorithm because each transport leg begins from a fully defined intermediate state and only needs that state together with the end-of-substep state. A multistep scheme such as BDF2, by contrast, also needs a history value from an earlier time level. That extra time history can still be managed, but it makes the construction of split substeps less natural and less transparent.

That is why Strang-split workflows are usually paired with a one-step second-order transport scheme rather than with a multistep BDF2 transport update.

10.5.3. Chemistry Substep in stream#

In the implementation, the split chemistry solve is carried out with a cell-local ODE integrator (CVODE)[CH94]. The primary unknowns advanced in the chemistry substep are the species densities. Denoting the reacted intermediate state by

(10.5.3.1)#\[(\rho Y_k)^{\mathrm{ODE}} = \mathcal{R}_{\Delta t_\chi}\left[(\rho Y_k)^{\mathrm{in}}\right],\]

the corresponding intermediate mass fractions are

(10.5.3.2)#\[Y_k^{\mathrm{ODE}}=\frac{(\rho Y_k)^{\mathrm{ODE}}}{\rho^{\mathrm{in}}}.\]

Here \(\Delta t_\chi\) denotes the chemistry interval used to construct the reacted intermediate state. In the formal Strang composition, each chemistry leg uses \(\Delta t_\chi=\Delta t/2\).

The reaction rates require thermodynamic closure. In stream, this closure is supplied by the equation of state in one of two conceptually distinct ways:

  1. Frozen thermodynamic state. The incoming pressure and temperature are held fixed during the chemistry ODE solve.

  2. EOS-refreshed thermodynamic state. As the composition changes during the chemistry solve, the thermodynamic state is recomputed from the EOS using the current composition together with the selected transported energy variable (\(e\) or \(h\)) held fixed for the chemistry substep.

The second form is the real-gas-compatible version, because the temperature and pressure seen by the chemistry integrator come from the EOS rather than from an ideal-gas specialization.

10.5.4. Transport Solve After the Chemistry Update#

Once the reacted intermediate composition has been constructed, the transport solve advances the usual pressure-based system. For the species equations, the important change is that the transient term is now written relative to the reacted intermediate state rather than the unsplit old-time composition. In first-order BDF form, this is

(10.5.4.1)#\[\frac{\rho^n\left(Y_k^{n+1}-Y_k^{\mathrm{ODE},n}\right)}{\Delta t} +\mathcal{T}_k^{n+1} = 0,\]

where \(\mathcal{T}_k\) denotes the spatial transport operator for species \(k\) (advection plus diffusion, with chemistry omitted because it has already been handled by the split chemistry solve).

For a one-step second-order transport discretization such as Crank-Nicolson, the same idea appears as

(10.5.4.2)#\[\frac{\rho^n\left(Y_k^{n+1}-Y_k^{\mathrm{ODE},n}\right)}{\Delta t} + \frac{1}{2}\left[\mathcal{T}_k^{n+1}+\mathcal{T}_k^{(a)}\right] = 0 .\]

Here \(\mathcal{T}_k^{(a)}\) denotes the transport operator evaluated on the reacted intermediate state produced by the first half-reaction step. In words: chemistry is advanced first over a half-step, the transport equation is then advanced from that reacted intermediate state using the selected one-step transport integrator, and the second chemistry half-step completes the symmetric update.

The momentum and pressure-correction equations are not solved during the chemistry ODE substep. They remain part of the transport solve described in Anatomy of a Segregated Pressure-Based Solver and Finite-Volume Transport Assembly on Unstructured Grids. The chemistry step influences them only indirectly through the updated composition and the EOS-consistent thermodynamic state used to evaluate density, temperature, and transport properties.

Note

In practical stream usage, the Strang-splitting reacting-flow workflow is paired with a one-step second-order transport update, typically Crank-Nicolson or a blended CN/BDF form, rather than with a multistep BDF2 transport update.

The main distinction between unsplit, source-split, and Strang-split treatments is therefore:

  • in an unsplit formulation, transport and source terms remain in the same time-discrete solve,

  • in a source-split formulation, the source operator is advanced separately in a one-sided sequential update,

  • in a Strang-split formulation, the operator split is made symmetric so that second-order time accuracy is recovered.

10.6. References#

[CH94]

S.D. Cohen and A.C. Hindmarsh. Cvode user guide. Technical Report, Lawrence Livermore National Laboratory, 1994.

[Str68]

G. Strang. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis, 5:506–517, 1968.

[VHB98]

J.G. Verwer, W.H. Hundsdorfer, and J.G. Blom. Numerical time integration for air pollution models. Technical Report, CWI, 1998.