Interpolatory HDG Method for Parabolic Semilinear PDEs

Abstract

We propose the interpolatory hybridizable discontinuous Galerkin (Interpolatory HDG) method for a class of scalar parabolic semilinear PDEs. The Interpolatory HDG method uses an interpolation procedure to efficiently and accurately approximate the nonlinear term. This procedure avoids the numerical quadrature typically required for the assembly of the global matrix at each iteration in each time step, which is a computationally costly component of the standard HDG method for nonlinear PDEs. Furthermore, the Interpolatory HDG interpolation procedure yields simple explicit expressions for the nonlinear term and Jacobian matrix, which leads to a simple unified implementation for a variety of nonlinear PDEs. For a globally-Lipschitz nonlinearity, we prove that the Interpolatory HDG method does not result in a reduction of the order of convergence. We display 2D and 3D numerical experiments to demonstrate the performance of the method.

Implementation Details for General Nonlinearities

1.1 The Interpolatory HDG Formulation

The full Interpolatory HDG discretization is to find $(\varvec{q}^n_h,u^n_h,\widehat{u}^n_h)\in \varvec{V}_h\times W_h\times M_h$ such that

$$\begin{aligned} \begin{aligned}&(\varvec{q}^n_h,\varvec{r})_{\mathcal {T}_h}-(u^n_h,\nabla \cdot \varvec{r})_{\mathcal {T}_h}+\left\langle \widehat{u}^n_h,\varvec{r} \cdot \varvec{n} \right\rangle _{\partial {\mathcal {T}_h}} = 0, \\&(\partial ^+_tu^n_h,w)_{{\mathcal {T}}_h}+(\nabla \cdot \varvec{q}^n_h, w)_{\mathcal {T}_h}+\langle \tau (u_h^n - \widehat{u}_h^n),w\rangle _{\partial {\mathcal {T}_h}} + ( {\mathcal {I}}_h F(-\varvec{q}_h^n, u_h^n),w)_{\mathcal {T}_h}{=} (f^n,w)_{\mathcal {T}_h},\\&\langle {\varvec{q}}^n_h\cdot \varvec{n} + \tau (u_h^n - \widehat{u}_h^n), \mu \rangle _{\partial {\mathcal {T}_h}\backslash \varepsilon ^{\partial }_h} =0,\\&u^0_h =\varPi _W u_0, \end{aligned} \end{aligned}$$

(45)

for all $(\varvec{r},w,\mu )\in \varvec{V}_h\times W_h\times M_h$ and $n=1,2,\ldots ,N$. Similar to Sect. 3.2, we have

$$\begin{aligned} ( {\mathcal {I}}_h F(-\varvec{q}_h^n, u_h^n),w)_{\mathcal {T}_h} = A_1 {\mathcal {F}}(\varvec{\alpha }^n, \varvec{\beta }^{n},\varvec{\gamma }^{n}), \end{aligned}$$

(46)

where

$$\begin{aligned} {\mathcal {F}}(\varvec{\alpha }^{n},\varvec{\beta }^{n},\varvec{\gamma }^{n}) = [F(\alpha _1^{n},\beta _1^{n}, \gamma _1^{n}),\ldots ,F(\alpha _{N_1}^{n},\beta _{N_1}^{n}, \gamma _{N_1}^{n})]^T. \end{aligned}$$

(47)

Then the system (45) can be rewritten as

$$\begin{aligned} \underbrace{\begin{bmatrix} A_1&\quad 0&\quad -\,A_2&\quad A_4 \\ 0&\quad A_1&\quad -\,A_3&\quad A_5 \\ A_2^T&\quad A_3^T&\quad A_6 +{\varDelta t}^{-1}A_1&\quad -\,A_7\\ A_4^T&\quad A_5^T&\quad A_7^T&\quad -\,A_8 \end{bmatrix}}_{M} \underbrace{\left[ {\begin{array}{*{20}{l}} \varvec{\alpha }^{n}\\ \varvec{\beta }^{n}\\ \varvec{\gamma }^{n}\\ \varvec{\zeta }^{n} \end{array}} \right] }_{\varvec{x}_{n}}+ \underbrace{\left[ {\begin{array}{*{20}{l}} 0\\ 0\\ A_1 {\mathcal {F}}(\varvec{\alpha }^{n},\varvec{\beta }^{n},\varvec{\gamma }^{n})\\ 0 \end{array}} \right] }_{{\mathscr {F}}(\varvec{x}_{n})} =\underbrace{\left[ {\begin{array}{*{20}{l}} 0\\ 0\\ b_1^n+{\varDelta t}^{-1}A_1\varvec{\gamma }^{n-1} \\ 0 \end{array}} \right] }_{\varvec{b}_n}, \end{aligned}$$

(48)

i.e., $M\varvec{x}_n + {\mathscr {F}}(\varvec{x}_n) = \varvec{b}_n$.

Newton’s method proceeds as in Sect. 3.2, but the Jacobian matrix $G'(\varvec{x}_n^{(m-1)})$ is now given by

$$\begin{aligned} G'(\varvec{x}_n^{(m-1)}) = M+{\mathscr {F}}'(\varvec{x}_n^{(m-1)}), \quad {\mathscr {F}}'(\varvec{x}_n^{(m-1)}) = \begin{bmatrix} 0&\quad 0&\quad 0&\quad 0 \\ 0&\quad 0&\quad 0&\quad 0 \\ A_{11}^{n,(m)}&\quad A_{12}^{n,(m)}&\quad A_{13}^{n,(m)}&\quad 0\\ 0&\quad 0&\quad 0&\quad 0 \end{bmatrix}, \end{aligned}$$

where for $k = 1, 2, 3,$ we define

$$\begin{aligned}&A_{1k}^{n,(m)} = A_1\text {diag}\big ({\mathcal {F}}_k'(\varvec{\alpha }^{n,(m-1)},\varvec{\beta }^{n,(m-1)},\varvec{\gamma }^{n,(m-1)})\big ),\\&{\mathcal {F}}_k'(\varvec{\alpha }^{n,(m-1)},\varvec{\beta }^{n,(m-1)},\varvec{\gamma }^{n,(m-1)}) \\&\quad = \big [F_k'(\alpha _1^{n,(m-1)},\beta _1^{n,(m-1)}, \gamma _1^{n,(m-1)}),\cdots ,F_k'(\alpha _{N_1}^{n,(m-1)},\beta _{N_1}^{n,(m-1)}, \gamma _{N_1}^{n,(m-1)})\big ]^T, \end{aligned}$$

and $F_k'$ denotes the partial derivative of F with respect to the kth variable. Therefore, the linear system that must be solved is now given by

$$\begin{aligned} \begin{bmatrix} A_1&\quad 0&\quad -\,A_2&\quad A_4 \\ 0&\quad A_1&\quad -\,A_3&\quad A_5 \\ A_2^T+ A_{11}^{n,(m)}&\quad A_3^T+ A_{12}^{n,(m)}&\quad A_6 +{\varDelta t}^{-1}A_1+ A_{13}^{n,(m)}&\quad -\,A_7\\ A_4^T&\quad A_5^T&\quad A_7^T&\quad -\,A_8 \end{bmatrix} \left[ {\begin{array}{*{20}{c}} \varvec{\alpha }^{n,(m)}\\ \varvec{\beta }^{n,(m)}\\ \varvec{\gamma }^{n,(m)}\\ \varvec{\zeta }^{n,(m)} \end{array}} \right] ={\widetilde{\varvec{b}}}, \end{aligned}$$

(49)

where

$$\begin{aligned} {\widetilde{\varvec{b}}} = G'(\varvec{x}_n^{(m-1)}) \varvec{x}_n^{(m-1)} - G(\varvec{x}_n^{(m-1)}). \end{aligned}$$

(50)

1.2 Local Solver

The system (49) can be rewritten as

$$\begin{aligned} \begin{bmatrix} B_1&\quad B_2&\quad B_3\\ B_4&\quad B_5&\quad -\,B_6\\ B_3^T&\quad B_6^T&\quad B_7\\ \end{bmatrix} \left[ {\begin{array}{*{20}{c}} \varvec{x}\\ \varvec{y}\\ \varvec{z} \end{array}} \right] =\left[ {\begin{array}{*{20}{c}} b_1\\ b_2\\ b_3 \end{array}} \right] , \end{aligned}$$

(51)

where $\varvec{x}=[\varvec{\alpha ^{n,(m)}};\varvec{\beta ^{n,(m)}}]$, $\varvec{y}=\varvec{\gamma }^{n,(m)}$, $\varvec{z}=\varvec{\zeta }^{n,(m)}$, ${\widetilde{\varvec{b}}} = [ b_1;b_2;b_3]$, and $\{B_i\}_{i=1}^7$ are the corresponding blocks of the coefficient matrix in (49). The system (51) is equivalent with following equations:

$$\begin{aligned} B_1 \varvec{x} + B_2\varvec{y} +B_3\varvec{z}&= b_1, \end{aligned}$$

(52a)

$$\begin{aligned} B_4 \varvec{x} +B_5\varvec{y} -B_6\varvec{z}&= b_2, \end{aligned}$$

(52b)

$$\begin{aligned} B_3^T\varvec{x}+ B_6^T\varvec{y} + B_7 \varvec{z}&=b_3. \end{aligned}$$

(52c)

Similar to before, the matrices $B_1$ and $B_5$ are block diagonal with small blocks and they can be easily inverted. Use (52a) and (52b) to express $\varvec{x}$ and $\varvec{y}$ in terms of $\varvec{z}$ as follows:

$$\begin{aligned} \varvec{x}&= B_1^{-1}B_2\left( B_4B_1^{-1}B_2+B_5\right) ^{-1}\left( (B_6+B_4B_1^{-1}B_3)\varvec{z}+ b_2-B_4B_1^{-1}b_1\right) \nonumber \\&\quad -B_1^{-1}B_3\varvec{z} + B_1^{-1}b_1\nonumber \\&=:{\tilde{B}}_1 \varvec{z} +{\tilde{b}}_1, \end{aligned}$$

(53)

$$\begin{aligned} \varvec{y}&=\left( B_4B_1^{-1}B_2+B_5\right) ^{-1} \left( (B_6+B_4B_1^{-1}B_3)\varvec{z} + b_2-B_4B_1^{-1}b_1\right) \nonumber \\&=:{\tilde{B}}_2 \varvec{\gamma }^n +{\tilde{b}}_2, \end{aligned}$$

(54)

where

$$\begin{aligned} Q = B_4B_1^{-1}B_2+B_5 = B_4B_1^{-1}B_2 +A_6 +{\varDelta t}^{-1}A_1+ A_{13}^{n,(m)}. \end{aligned}$$

As in Sect. 3.3, the matrix Q is block diagonal with small blocks. Since $A_1$ is positive definite, if $\varDelta t$ is small enough then Q is easily inverted. Then we insert $\varvec{x}$ and $\varvec{y}$ into (26c) and obtain the final system only involving $\varvec{z}$:

$$\begin{aligned} (B_3^T {\tilde{B}}_1 + B_5^T {\tilde{B}}_2 + B_6) \varvec{z} = b_3 -B_3^T{\tilde{b}}_1 -B_5^T {\tilde{b}}_2. \end{aligned}$$

(55)