implicit

Definitions of implicit time integration schemes for conservative methods.

class BDF2(mesh, root=None, **default)

Class responsible for implementing an implicit 2nd-order backward differentiation formula that updates the current solution (\(t = t^{n}\)) to the next time step (\(t = t^{n+1}\)), using also the previous solution (\(t = t^{n-1}\)). Assuming an HDG formulation,

\[\begin{split}\widetilde{\bm{M}} \bm{U}^{n+1} + \bm{f}(\bm{U}^{n+1}, \hat{\bm{U}}) &= \widetilde{\bm{M}} \Big( \frac{4}{3} \bm{U}^{n} - \frac{1}{3} \bm{U}^{n-1}\Big),\\ \bm{g}(\bm{U}^{n+1}, \hat{\bm{U}}) &= \bm{0},\end{split}\]

where \(\widetilde{\bm{M}} = 3 \bm{M} / (2\delta t)\) is the weighted mass matrix and \(\bm{M}\) is the mass matrix and \(\bm{f}\) and \(\bm{g}\) arise from the spatial discretization of the PDE on the physical elements and the AE on the facets, respectively.

class DIRK34_LDD(mesh, root=None, **default): Updates the solution via a 3-stage 4th-order diagonally-implicit Runge-Kutta (DIRK) with a low-dispersion and dissipation. Taken from Table A.1 in [6]. Its corresponding Butcher tableau is:

\[\begin{split}\begin{array}{c|ccc} 0.257820901066211 & 0.377847764031163 & \phantom{-}0 & 0 \\ 0.434296446908075 & 0.385232756462588 & \phantom{-}0.461548399939329 & 0 \\ 0.758519768667167 & 0.675724855841358 & -0.061710969841169 & 0.241480233100410 \\ \hline & 0.750869573741408 & -0.362218781852651 & 0.611349208111243 \end{array}\end{split}\]

class DIRK43_WSO2(mesh, root=None, **default)

Updates the solution via a 4-stage 3rd-order (stiffly-accurate) diagonally-implicit Runge-Kutta (DIRK) with a weak stage order (WSO) of 3. Taken from Section 3 in [7]. Its corresponding Butcher tableau is:

\[\begin{split}\begin{array}{c|cccc} 0.01900072890 & 0.01900072890 & \phantom{-}0 & 0 & 0 \\ 0.78870323114 & 0.40434605601 & \phantom{-}0.38435717512 & 0 & 0 \\ 0.41643499339 & 0.06487908412 & -0.16389640295 & 0.51545231222 & 0 \\ 1 & 0.02343549374 & -0.41207877888 & 0.96661161281 & 0.42203167233\\ \hline & 0.02343549374 & -0.41207877888 & 0.96661161281 & 0.42203167233 \end{array}\end{split}\]

Note:: No need to explicitly form the solution at the next time step, since this is a stiffly-accurate method, i.e. \(\bm{U}^{n+1} = \bm{y}_{4}\).

class DIRKSchemes(mesh, root=None, **default)

Assuming an HDG formulation, we are solving

\[\begin{split}\bm{M} \partial_t \bm{U} + \bm{f}\big(\bm{U}, \hat{\bm{U}} \big) &= \bm{0},\\ \bm{g}\big(\bm{U}, \hat{\bm{U}} \big) &= \bm{0},\end{split}\]

where \(\bm{M}\) is the mass matrix, and \(\bm{f}\) and \(\bm{g}\) arise from the spatial discretization of the PDE on the physical elements and the AE on the facets, respectively.

To obtain the solution at \(t = t^{n+1}\), an s-stage DIRK update formula is

\[\bm{U}^{n+1} = \bm{U}^{n} - \delta t \sum_{i=1}^{s} b_{i} \bm{M}^{-1} \bm{f}(\bm{z}_{i}),\]

where \(\bm{z}_i = \big( \bm{y}_i, \hat{\bm{y}}_i \big)^T\) denotes the physical \(\bm{y}\) and facet \(\hat{\bm{y}}\) solution at the ith stage, obtained from

\[\begin{split}\bm{y}_{i} &= \bm{U}^{n} - \delta t \sum_{j=1}^{s} a_{ij} \bm{M}^{-1} \bm{f}(\bm{z}_j),\\ \bm{g}(\bm{z}_i) &= \bm{0}.\end{split}\]

The above can be used to define a residual for solving the nonlinear system iteratively

\[\begin{split}\bm{r}_i &= \widetilde{\bm{M}}_i \bm{y}_i - \widetilde{\bm{M}}_i \bm{U}^{n} + \frac{1}{a_{ii}} \sum_{j=1}^{i-1} a_{ij} \bm{f}(\bm{z}_j) + \bm{f}(\bm{z}_i),\\ \hat{\bm{r}}_i &= \bm{g}(\bm{z}_i),\end{split}\]

where \(\widetilde{\bm{M}}_i = \bm{M}/(a_{ii} \delta t)\) is the weighted mass matrix associated with the ith stage.

Therefore, a Newton-Rhapson update formula for the kth iteration can be written as

\[\bm{N}_{i}^{k} \delta \bm{z}_{i}^{k} = - \bm{R}_{i}^{k},\]

where the overall residual is \(\bm{R}_{i}^{k} = \Big( \bm{r}(\bm{z}_{i}^{k}), \hat{\bm{r}}(\bm{z}_{i}^{k}) \Big)^T\) and the iteration matrix is defined as

\[\begin{split}\bm{N}_{i}^{k} &= \partial \bm{R} (\bm{z}_{i}^{k} ) / \partial \bm{z}_i, \\[2ex] &= \begin{pmatrix} \widetilde{\bm{M}}_{i} + \partial \bm{f}^{k} / \partial \bm{y}_i & \partial \bm{f}^{k} / \partial \hat{\bm{y}}_i\\ \partial \bm{g}^{k} / \partial \bm{y}_i & \partial \bm{g}^{k} / \partial \hat{\bm{y}}_i \end{pmatrix}.\end{split}\]

Note:: In terms of implementation, this is based on two bilinear forms: blf and blfs

\[\begin{split}\bf{blf} = \begin{pmatrix} \widetilde{\bm{M}}_i \bm{y}_i + \bm{f}(\bm{z}_i)\\ \bm{g}(\bm{z}_i) \end{pmatrix}, \qquad \bf{blfs} = \begin{pmatrix} \bm{f}(\bm{z}_i)\\ \bm{0} \end{pmatrix}.\end{split}\]

class ImplicitEuler(mesh, root=None, **default)

Class responsible for implementing an implicit (backwards-)Euler time-marching scheme that updates the current solution (\(t = t^{n}\)) to the next time step (\(t = t^{n+1}\)). Assuming an HDG formulation,

\[\begin{split}\widetilde{\bm{M}} \bm{U}^{n+1} + \bm{f}(\bm{U}^{n+1}, \hat{\bm{U}}) &= \widetilde{\bm{M}} \bm{U}^{n},\\ \bm{g}(\bm{U}^{n+1}, \hat{\bm{U}}) &= \bm{0},\end{split}\]

where \(\widetilde{\bm{M}} = \bm{M} / \delta t\) is the weighted mass matrix, \(\bm{M}\) is the mass matrix and \(\bm{f}\) and \(\bm{g}\) arise from the spatial discretization of the PDE on the physical elements and the AE on the facets, respectively.

class ImplicitSchemes(mesh, root=None, **default)

class SDIRK22(mesh, root=None, **default)

Updates the solution via a 2-stage 2nd-order (stiffly-accurate) singly diagonally-implicit Runge-Kutta (SDIRK). Taken from Section 2.6 in [2]. Its corresponding Butcher tableau is:

\[\begin{split}\begin{array}{c|cc} \alpha & \alpha & 0 \\ 1 & 1 - \alpha & \alpha \\ \hline & 1 - \alpha & \alpha \end{array}\end{split}\]

where \(\alpha = (2 - \sqrt{2})/2\).

Note:: No need to explicitly form the solution at the next time step, since this is a stiffly-accurate method, i.e. \(\bm{U}^{n+1} = \bm{y}_{2}\).

class SDIRK33(mesh, root=None, **default)

Updates the solution via a 3-stage 3rd-order (stiffly-accurate) singly diagonally-implicit Runge-Kutta (SDIRK). Taken from Section 2.7 in [2]. Its corresponding Butcher tableau is:

\[\begin{split}\begin{array}{c|ccc} 0.4358665215 & 0.4358665215 & 0 & 0 \\ 0.7179332608 & 0.2820667392 & \phantom{-}0.4358665215 & 0 \\ 1 & 1.2084966490 & -0.6443631710 & 0.4358665215 \\ \hline & 1.2084966490 & -0.6443631710 & 0.4358665215 \end{array}\end{split}\]

Note:: No need to explicitly form the solution at the next time step, since this is a stiffly-accurate method, i.e. \(\bm{U}^{n+1} = \bm{y}_{3}\).

class SDIRK54(mesh, root=None, **default)

Updates the solution via a 5-stage 4th-order (stiffly-accurate) singly diagonally-implicit Runge-Kutta (SDIRK). Taken from Table 6.5 in [8]. Its corresponding Butcher tableau is:

\[\begin{split}\begin{array}{c|ccccc} \frac{1}{4} & \frac{1}{4} & \phantom{-}0 & 0 & \phantom{-}0 & 0 \\ \frac{3}{4} & \frac{1}{2} & \phantom{-}\frac{1}{4} & 0 & \phantom{-}0 & 0 \\ \frac{11}{20} & \frac{17}{50} & -\frac{1}{25} & \frac{1}{4} & \phantom{-}0 & 0 \\ \frac{1}{2} & \frac{371}{1360} & -\frac{137}{2720} & \frac{15}{544} & \phantom{-}\frac{1}{4} & 0 \\ 1 & \frac{25}{24} & -\frac{49}{48} & \frac{125}{16} & -\frac{85}{12} & \frac{1}{4}\\ \hline & \frac{25}{24} & -\frac{49}{48} & \frac{125}{16} & -\frac{85}{12} & \frac{1}{4} \end{array}\end{split}\]

Note:: No need to explicitly form the solution at the next time step, since this is a stiffly-accurate method, i.e. \(\bm{U}^{n+1} = \bm{y}_{5}\).