Advanced Tutorial 1 : Zero-to-Nodes (Z2N) and Node-to-Node (N2N) formulations for SDC

📜 If you already know about SDC from theoriginal paperof Dutt, Greengard & Rokhlin, you may notice that their description is very different from the one givenin Step 4…

Indeed, this tutorial introduced SDC using a Zero-to-Nodes formulation (Z2N), which describes the SDC node updates from the initial step solution (zero) to the node. This approach is identical as the one used to describe Runge-Kutta methods with Butcher tables in the literature.

The SDC authors however used a different formulation, namely the Node-to-Node formulation (N2N), which describes the node update from one node to the next. While both formulations can produce identical SDC schemes, they have some fundamental differences from an implementation perspective, and leads to different generalizations of SDC.

Deriving the N2N formulation

We start from the Z2N update on the Dahlquist problem :

\[u^{k+1} - \lambda\Delta{t}Q_\Delta u^{k+1} = u_n + \lambda\Delta{t}(Q-Q_\Delta)u^k,\]

we can expend it to an update for each node solution \(u_{m} \simeq u(t_n+\tau_m\Delta{t}),\; m \in \{1,\dots,M\}\):

\[u^{k+1}_{m+1} - \lambda\Delta{t}\sum_{j=1}^{m+1}q^\Delta_{m+1,j} u^{k+1}_{j} = u_n + \lambda\Delta{t}\sum_{j=1}^{M}q_{m+1,j}u^{k}_{j} - \lambda\Delta{t}\sum_{j=1}^{m+1}q^\Delta_{m+1,j}u^{k}_{j},\]

where \(u_n\) is the initial solution for the time-step (scalar, abusing notation again …), and we note \((q^\Delta)_{i,j} := Q_\Delta\) and \((q)_{i,j} := Q\). Rearranging and regrouping terms, we can write it like this :

\[u^{k+1}_{m+1} = u_n + \lambda\Delta{t}\sum_{j=1}^{m+1}q^\Delta_{m+1,j} (u^{k+1}_{j} - u^{k}_{j}) + \lambda\Delta{t}\sum_{j=1}^{M}q_{m+1,j}u^{k}_{j}.\]

Now subtracting the update formula for \(u^{k+1}_m\) from \(u^{k+1}_{m+1}\), we get for \(m > 0\) (starting from the second node) :

\[u^{k+1}_{m+1} = u^{k+1}_m + \lambda\Delta{t}\sum_{j=1}^{m+1}\left(q^\Delta_{m+1,j} - q^\Delta_{m,j}\right)\left(u^{k+1}_{j} - u^{k}_{j}\right) + \lambda\Delta{t}\sum_{j=1}^{M}\left(q_{m+1,j}-q_{m,j}\right)u^{k}_{j},\]

and for the first node (no difference) :

\[u^{k+1}_{1} = u_n + \lambda\Delta{t}q^\Delta_{1,1} (u^{k+1}_{1} - u^{k}_{1}) + \lambda\Delta{t}\sum_{j=1}^{M}q_{1,j}u^{k}_{j}\]

We call this the Node-to-node (N2N) update, since each node update for \(u^{k+1}_{m+1}\) explicitly depends on the previously computed update \(u^{k+1}_{m}\) (and the correction terms \(\sum ...\)), except for the first node where it depends on the initial step solution \(u_n\).

🔔 Note that \(q^\Delta_{m,m+1}=0\) because of the lower triangular nature of \(Q_\Delta\), so we can add this coefficient in the generic N2N sweep formula to simplify it.

Defining \(s_{m+1,j} = q_{m+1,j}-q_{m,j} \; \forall m \in \{1, M-1\}\) and \(s_{1,j} = q_{1,j}\), we note \(S\) the matrix built with the \((s)_{i,j}\) coefficients, and write the generic N2N sweep update into matrix formulation :

\[\begin{split}\begin{pmatrix} 1 \\ -1 & 1 \\ & \ddots & \ddots \\ & & -1 & 1 \end{pmatrix} (u^{k+1} - u^{k}) = \begin{pmatrix} u0 \\ 0 \\ \vdots \\ 0 \end{pmatrix} + \lambda\Delta{t}S_\Delta (u^{k+1}-u^{k}) + \lambda\Delta{t}S u^k,\end{split}\]

where \(S_\Delta\) is built from the \(Q_\Delta\) matrix the same way as \(S\) from \(Q\). Now you may remark that the N2N formula is actually more complex than the Z2N formula given above. Thing is : some specific \(Q_\Delta\) coefficients allows to simplify it a lot …

Backward-Euler based sweep

Consider now one of the original SDC form proposed by Dutt, Greengard & Rokhlin : the Backward-Euler based sweep. The associated \(Q_\Delta\) coefficients are implemented in qmat, and considering a simple illustrative node distribution we obtain (using the default N2N formulation) :

from qmat import genQDeltaCoeffs
genQDeltaCoeffs("BE", nodes=[0.1, 0.3, 0.7, 1.0])

array([[0.1, 0. , 0. , 0. ],
       [0.1, 0.2, 0. , 0. ],
       [0.1, 0.2, 0.4, 0. ],
       [0.1, 0.2, 0.4, 0.3]])

This is a triangular matrix with non-zero diagonal (implicit sweep), but with all non-zero coefficients in each columns being identical. This implies that for \(m \in \{1, M\}\) :

\[\lambda\Delta{t}\sum_{j=1}^{m} s^\Delta_{m+1,j}\left(u^{k+1}_{j} - u^{k}_{j}\right) = \lambda\Delta{t}\sum_{j=1}^{m}\left(q^\Delta_{m+1,j} - q^\Delta_{m,j}\right)\left(u^{k+1}_{j} - u^{k}_{j}\right) = 0\]

thus simplifies the N2N formula into :

\[u^{k+1}_{m+1} = u^{k+1}_m + \lambda\Delta{t}s^\Delta_{m+1,m+1}\left(u^{k+1}_{m+1} - u^{k}_{m+1}\right) + \lambda\Delta{t}\sum_{j=1}^{M}s_{m+1,j}u^{k}_{j},\]

where we note \(\forall k\;u^{k}_0 := u_n\), such that the formula can be applied on all nodes. In fact, the \(S_\Delta\) matrix is diagonal, as we can see using qmat setting form="N2N" for the genQDeltaCoeffs function :

genQDeltaCoeffs("BE", form="N2N", nodes=[0.1, 0.3, 0.7, 1.0])

array([[0.1, 0. , 0. , 0. ],
       [0. , 0.2, 0. , 0. ],
       [0. , 0. , 0.4, 0. ],
       [0. , 0. , 0. , 0.3]])

In this case in particular, we have \(s^\Delta_{m+1,m+1} = q^\Delta_{m+1,m+1} = \tau_{m+1}-\tau_{m}\), where \((\tau_m)_{1\leq m \leq M} \in [0, 1]\) are the nodes positions and \(\tau_0=0\). This is then exactly the SDC formula for Backward Euler introduced by Dutt, Greengard & Rokhlin.

🔔 Note that the correction term \(\lambda\Delta{t}s^\Delta_{m+1,m+1}\left(u^{k+1}_{m+1} - u^{k}_{m+1}\right)\) depends on the result of the sweep update \(u^{k+1}_{m+1}\), requiring then an implicit solver.

Forward Euler based sweep

An other sweep update was proposed by Dutt, Greengard & Rokhlin in the same spirit as for Backward Euler, but this time using a Forward Euler based correction.

Consider the \(Q_\Delta\) coefficients implemented in qmat for the same node distribution as before :

genQDeltaCoeffs("FE", nodes=[0.1, 0.3, 0.7, 1.0])

array([[0. , 0. , 0. , 0. ],
       [0.2, 0. , 0. , 0. ],
       [0.2, 0.4, 0. , 0. ],
       [0.2, 0.4, 0.3, 0. ]])

Those are similar to the ones obtained for Backward Euler, except they are “shifted” to the left and the diagonal contains only zeros. Looking at the corresponding \(S_\Delta\) coefficients, we obtain :

\[u^{k+1}_{m+1} = u^{k+1}_m + \lambda\Delta{t}s^\Delta_{m,m}\left(u^{k+1}_{m} - u^{k}_{m}\right) + \lambda\Delta{t}\sum_{j=1}^{M}s_{m+1,j}u^{k}_{j},\]

with \(s^\Delta_{m,m}=q^\Delta_{m,m}=\tau_{m+1}-\tau_{m}\) where \((\tau_m)_{1\leq m \leq M} \in [0, 1]\) are the nodes positions and \(\tau_0=0\) (as before), as we can see looking at the \(S_\Delta\) coefficient generated by qmat :

genQDeltaCoeffs("FE", form="N2N", nodes=[0.1, 0.3, 0.7, 1.0])

array([[0. , 0. , 0. , 0. ],
       [0.2, 0. , 0. , 0. ],
       [0. , 0.4, 0. , 0. ],
       [0. , 0. , 0.3, 0. ]])

This is similar as for Backward Euler, except the correction term does not depend on \(u^{k+1}_{m+1}\) this time, hence making this update fully explicit.

Additional notes

Other sweep types have a simplified form in N2N formulation, e.g using the trapezoid rule (Crank-Nicholson) :

genQDeltaCoeffs("TRAP", form="N2N", nodes=[0.1, 0.3, 0.7, 1.0])

array([[0.05, 0.  , 0.  , 0.  ],
       [0.1 , 0.1 , 0.  , 0.  ],
       [0.  , 0.2 , 0.2 , 0.  ],
       [0.  , 0.  , 0.15, 0.15]])

which produce a sweep update of this form :

\[u^{k+1}_{m+1} = u^{k+1}_m + \lambda\Delta{t}s^\Delta_{m+1,m+1}\left(u^{k+1}_{m+1} - u^{k}_{m+1}\right) + \lambda\Delta{t}s^\Delta_{m,m}\left(u^{k+1}_{m} - u^{k}_{m}\right) + \lambda\Delta{t}\sum_{j=1}^{M}s_{m+1,j}u^{k}_{j}.\]

From an algorithmic perspective, implementing SDC into N2N form for Backward / Forward Euler or the trapezoid rule is usually more efficient than the Z2N form, since it requires less floating point operations during sweep (correction terms use all node solution in Z2N). However, the N2N formulation has two major issues when considering generic SDC methods :

1. Only a few type of \(Q_\Delta\) coefficients allows a simplified N2N formulation, while some other (like LU for instance), don’t have a simplified N2N formulation, hence making the Z2N implementation more efficient

2. The N2N formulation is inherently sequential, and do not allow the design of parallel SDC variant that have a diagonal \(Q_\Delta\) matrix, allowing to perform all node update in parallel across the nodes

Furthermore, while the N2N formulation was historically used to describe SDC, many later SDC-related publications have been using the Z2N formulation since it’s more convenient for the analysis, when looking to SDC as a preconditioned fixed point iteration.