A simplified Cauchy-Kowalewskaya procedure for the implicit solution of generalized Riemann problems of hyperbolic balance laws

TL;DR

This paper introduces a simplified, recursive Cauchy-Kowalewskaya procedure for implicit Taylor series solutions of hyperbolic balance laws, improving computational efficiency and maintaining accuracy in generalized Riemann problem solvers.

Contribution

A new recursive formula for the CK procedure simplifies implementation and enhances efficiency in solving hyperbolic balance laws within the GRP framework.

Findings

The proposed method achieves expected theoretical orders of accuracy.

It improves computational efficiency in solving GRPs.

The scheme maintains stability while enhancing performance.

Abstract

The Cauchy-Kowalewskaya (CK) procedure is a key building block in the design of solvers for the Generalised Rieman Problem (GRP) based on Taylor series expansions in time. The CK procedure allows us to express time derivatives in terms of purely space derivatives. This is a very cumbersome procedure, which often requires the use of software manipulators. In this paper, a simplification of the CK procedure is proposed in the context of implicit Taylor series expansion for GRP, for hyperbolic balance laws in the framework of [Journal of Computational Physics 303 (2015) 146-172]. A recursive formula for the CK procedure, which is straightforwardly implemented in computational codes, is obtained. The proposed GRP solver is used in the context of the ADER approach and several one-dimensional problems are solved to demonstrate the applicability and efficiency of the present scheme. An enhancement in terms of efficiency, is obtained. Furthermore, the expected theoretical orders of accuracy are achieved, conciliating accuracy and stability.

Tables

Table 1. Table 1: Coefficients a l , k subscript 𝑎 𝑙 𝑘 a_{l,k} in expression ( 52 ).

$a_{l,\left(l-4\right)}$	$a_{l,\left(l-3\right)}$	$a_{(l,l-2)}$	$a_{(l,l-1)}$	$a_{(l,l)}$	$a_{(l,l+1)}$	$l$
0	0	0	0	1	1	1
0	0	0	1	2	1	2
0	0	1	3	3	1	3
0	1	4	6	4	1	4
1	5	10	10	5	1	5

Table 2. Table 2: Coefficients b l . k subscript 𝑏 formulae-sequence 𝑙 𝑘 b_{l.k} in expression ( 52 ).

$b_{l,\left(l-4\right)}$	$b_{l,\left(l-3\right)}$	$b_{(l,l-2)}$	$b_{(l,l-1)}$	$b_{(l,l)}$	$b_{(l,l+1)}$	$l$
0	0	0	0	0	1	1
0	0	0	0	1	1	2
0	0	0	1	2	1	3
0	0	1	3	3	1	4
0	1	4	6	4	1	5

Table 3. Table 3: Linear system. Output time t o ​ u ​ t = 1 subscript 𝑡 𝑜 𝑢 𝑡 1 t_{out}=1 with C c ​ f ​ l = 0.9 , subscript 𝐶 𝑐 𝑓 𝑙 0.9 C_{cfl}=0.9, β = − 1 𝛽 1 \beta=-1 , λ = 1 𝜆 1 \lambda=1 .

Mesh	$L_{\infty }$ - err	$L_{\infty }$- ord	$L_1$ - err	$L_1$ - ord	$L_2$ - err	$L_2$ - ord	CPU
8	-	$2.45e-02$	-	$1.33e-02$	-	$1.74e-02$	0.0064
16	2.35	$4.81e-03$	2.38	$2.55e-03$	2.55	$2.97e-03$	0.0114
32	1.67	$1.52e-03$	2.85	$3.54e-04$	2.34	$5.84e-04$	0.0287
64	2.95	$1.96e-04$	3.71	$2.71e-05$	3.47	$5.29e-05$	0.0853
128	1.49	$6.96e-05$	1.83	$7.62e-06$	1.91	$1.41e-05$	0.3395

Table 4. Table 4: Non-linear system. Output time t o ​ u ​ t = 0.1 subscript 𝑡 𝑜 𝑢 𝑡 0.1 t_{out}=0.1 with C c ​ f ​ l = 0.9 , subscript 𝐶 𝑐 𝑓 𝑙 0.9 C_{cfl}=0.9, β = − 1 𝛽 1 \beta=-1 .

Mesh	$L_{\infty }$ - err	$L_{\infty }$- ord	$L_1$ - err	$L_1$ - ord	$L_2$ - err	$L_2$ - ord	CPU
32	-	$1.95e-02$	-	$5.45e-03$	-	$7.81e-03$	0.0061
64	1.53	$6.75e-03$	2.12	$1.25e-03$	1.94	$2.04e-03$	0.0201
128	1.75	$2.01e-03$	2.30	$2.54e-04$	2.18	$4.49e-04$	0.0490
256	1.80	$5.79e-04$	2.18	$5.62e-05$	2.15	$1.01e-04$	0.1703
512	1.38	$2.23e-04$	2.13	$1.29e-05$	2.06	$2.44e-05$	0.7851

Table 5. Table 5: Euler equations.Output time t o ​ u ​ t = 1 subscript 𝑡 𝑜 𝑢 𝑡 1 t_{out}=1 with C c ​ f ​ l = 0.9 . subscript 𝐶 𝑐 𝑓 𝑙 0.9 C_{cfl}=0.9.

Mesh	$L_{\infty }$ - err	$L_{\infty }$- ord	$L_1$ - err	$L_1$ - ord	$L_2$ - err	$L_2$ - ord	CPU
8	-	$1.51e-01$	-	$1.05e-01$	-	$1.14e-01$	0.0119
16	1.22	$6.45e-02$	1.76	$3.09e-02$	1.59	$3.79e-02$	0.0401
32	1.43	$2.40e-02$	1.62	$1.01e-02$	1.65	$1.21e-02$	0.1569
64	1.49	$8.50e-03$	1.90	$2.71e-03$	1.79	$3.49e-03$	0.4242
128	1.53	$2.95e-03$	2.07	$6.46e-04$	1.85	$9.69e-04$	1.6852

Table 6. Table 6: Coefficients a M + 1 , k subscript 𝑎 𝑀 1 𝑘 a_{M+1,k} for interpolation in space, P ​ ( x ) = ∑ k = 0 M a M + 1 , k ⋅ ξ k 𝑃 𝑥 superscript subscript 𝑘 0 𝑀 ⋅ subscript 𝑎 𝑀 1 𝑘 superscript 𝜉 𝑘 P(x)=\sum_{k=0}^{M}a_{M+1,k}\cdot\xi^{k} .

Second order ($M=1$).
$a_{2,0}=f_2$
$a_{2,1}=\left(f_2-f_1\right)$
Third order ($M=2$).
$a_{3,0}=f_2$
$a_{3,0}=\left(f_3-f_1\right)$
$a_{3,0}=2\left(f_3-2f_2+f_1\right)$
Fourth order ($M=3$).
$a_{4,0}=-\left(f_4-9f_3-9f_2+f_1\right)/16$
$a_{4,1}=+\left(\left(-f_4+27f_3-27f_2+f_1\right)\right)/8$
$a_{4,2}=\left(\left(9f_4-9f_3-9f_2+9f_1\right){\xi }^2\right)/4$
$a_{4,3}=-\left(\left(-9f_4+27f_3-27f_2+9f_1\right){\xi }^3\right)/2$
Fifth order ($M=4$).
$a_{5,0}=f_3$
$a_{5,1}=+\left(\left(-f_5+8f_4-8f_2+f_1\right)\right)/3$
$a_{5,2}=-\left(\left(2f_5-32f_4+60f_3-32f_2+2f_1\right)\right)/3$
$a_{5,3}=-\left(\left(-16f_5+32f_4-32f_2+16f_1\right)\right)/3$
$a_{5,4}=\left(\left(32f_5-128f_4+192f_3-128f_2+32f_1\right)\right)/3$

Table 7. Table 7: Coefficients b M + 1 , k subscript 𝑏 𝑀 1 𝑘 b_{M+1,k} for interpolation in time, T ​ ( τ ) = ∑ k = 0 M − 1 b M + 1 , k ⋅ τ k 𝑇 𝜏 superscript subscript 𝑘 0 𝑀 1 ⋅ subscript 𝑏 𝑀 1 𝑘 superscript 𝜏 𝑘 T(\tau)=\sum_{k=0}^{M-1}b_{M+1,k}\cdot\tau^{k} .

Third order ($M=2$).
$b_{3,0}=(-\sqrt{3}{f}_2+f_2+\frac{(\sqrt{3}+1)f_1)}{2}$
$b_{3,1}=(\sqrt{3}{f}_2-\sqrt{3}{f}_1))$
Fourth order ($M=3$).
$b_{4,0}=\left(-\sqrt{15}{f}_3+5f_3-4f_2+\left(\sqrt{15}+5\right)f_1\right)/6$
$b_{4,1}=-\frac{\left(-\sqrt{15}{f}_3+10f_3-20f_2+\left(\sqrt{15}+10\right)f_1\right)}{3}$
$b_{4,2}=\frac{\left(10f_3-20f_2+10f_1\right)}{3}$
Fifth order ($M=4$).
$b_{5,0}=-0.1139171962819898f_4+0.4007615203116506f_3-0.8136324494869276f_2+1.526788125457266f_1$
$b_{5,1}=2.15592710364526f_4-7.41707042146264f_3+13.80716692568958f_2-8.546023607872199f_1$
$b_{5,2}=-7.935761849944949f_4+24.99812585921913f_3-31.38822236344606f_2+14.32585835417188f_1$
$b_{5,3}=7.420540068038946f_4\cdot -18.79544940755506f_3+18.79544940755506f_2-7.420540068038946f_1$

Table 8. Table 8: Numerical source term. Here, ω j subscript 𝜔 𝑗 \omega_{j} are the corresponding Gaussian weights.

Second order ($M=1$).
${𝐒}_i=\frac{1}{2}{\sum }_{j=1}^{n_{GP}}{\omega }_j(𝐒({𝐐}_i({\xi }_2,{\tau }_j)+𝐒\left({𝐐}_i({\xi }_1,{\tau }_j)\right)$
Third order ($M=2$).
${𝐒}_i=\frac{1}{6}{\sum }_{j=1}^{n_{GP}}{\omega }_j(𝐒({𝐐}_i({\xi }_3,{\tau }_j)+4𝐒({𝐐}_i({\xi }_2,{\tau }_j)+𝐒\left({𝐐}_i({\xi }_1,{\tau }_j)\right)$
Fourth order ($M=3$).
${𝐒}_i=\frac{1}{8}{\sum }_{j=1}^{n_{GP}}{\omega }_j(𝐒({𝐐}_i({\xi }_4,{\tau }_j)+3𝐒({𝐐}_i({\xi }_3,{\tau }_j)+3𝐒\left({𝐐}_i({\xi }_2,{\tau }_j)\right)+𝐒\left({𝐐}_i({\xi }_1,{\tau }_j)\right)$
Fifth order ($M=4$).
${𝐒}_i=\frac{1}{90}{\sum }_{j=1}^{n_{GP}}{\omega }_j(7𝐒({𝐐}_i({\xi }_5,{\tau }_j)+32𝐒({𝐐}_i({\xi }_4,{\tau }_j)+12𝐒({𝐐}_i({\xi }_3,{\tau }_j)+32𝐒\left({𝐐}_i({\xi }_2,{\tau }_j)\right)+7𝐒\left({𝐐}_i({\xi }_1,{\tau }_j)\right)$

Equations

$$\displaystyle\begin{array}[]{c}\partial_{t}\mathbf{Q}+\partial_{x}\mathbf{F}(\mathbf{Q})=\mathbf{S}(\mathbf{Q})\;,\\ \mathbf{Q}(x,0)=\mathbf{H}_{0}(x)\;,\end{array}$$

$$\displaystyle\begin{array}[]{c}\mathbf{Q}_{i}^{n+1}=\mathbf{Q}_{i}^{n}-\frac{\Delta t}{\Delta x}\biggl{(}\mathbf{F}_{i+\frac{1}{2}}-\mathbf{F}_{i-\frac{1}{2}}\biggr{)}+\Delta t\mathbf{S}_{i}\;,\end{array}$$

$$\displaystyle\begin{array}[]{c}\mathbf{Q}_{i}^{n}=\frac{1}{\Delta x}\displaystyle\int_{x-\frac{1}{2}}^{x+\frac{1}{2}}\mathbf{Q}(x,t^{n})dx\end{array}$$

$$\displaystyle\begin{array}[]{c}\mathbf{F}_{i+\frac{1}{2}}=\frac{1}{\Delta t}\displaystyle\int_{t^{n}}^{t^{n+1}}\mathbf{F}_{h}(\mathbf{Q}_{i}(x_{i+\frac{1}{2}},t),\mathbf{Q}_{i+1}(x_{i+\frac{1}{2}},t))dt\;,\par\\ \mathbf{S}_{i}=\frac{1}{\Delta t\Delta x}\displaystyle\int_{t^{n}}^{t^{n+1}}\int_{x-\frac{1}{2}}^{x+\frac{1}{2}}\mathbf{S}(\mathbf{Q}_{i}(x,t)dxdt\;,\end{array}$$

$$\displaystyle\begin{array}[]{c}\mathbf{Q}_{i}(x,\tau)=\mathbf{Q}_{i}(x,0_{+})-\sum_{k=1}^{M}\frac{(-\tau)^{k}}{k!}\partial_{t}^{(k)}\mathbf{Q}_{i}(x,\tau)\;,\end{array}$$

$$\displaystyle\begin{array}[]{c}\mathbf{Q}_{i}(x,\tau)=\mathbf{Q}_{i}(x,0_{+})-\sum_{k=1}^{M}\frac{(-\tau)^{k}}{k!}\mathbf{G}^{(k)}(\mathbf{Q}_{i}(x,\tau),\partial_{x}\mathbf{Q}_{i}(x,\tau),...,\partial_{x}^{(k)}\mathbf{Q}_{i}(x,\tau))\;,\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{t}\mathbf{Q}=-\mathbf{A}(\mathbf{Q})\partial_{x}\mathbf{Q}+\mathbf{S}(\mathbf{Q})\;,\end{array}$$

$$\displaystyle\begin{array}[]{c}\displaystyle\partial_{t}^{(2)}\mathbf{Q}_{i}=-\sum_{j=1}^{m}\sum_{l=1}^{m}\frac{\partial\mathbf{A}_{i,j}(\mathbf{Q})}{\partial\mathbf{Q}_{l}}\partial_{t}\mathbf{Q}_{l}\partial_{x}\mathbf{Q}_{j}-\sum_{j=1}^{m}\mathbf{A}_{i,j}(\mathbf{Q})\partial_{t}(\partial_{x}\mathbf{Q}_{j})\\ \displaystyle+\sum_{j=1}^{m}\mathbf{B}(\mathbf{Q})_{i,j}\partial_{t}\mathbf{Q}_{j}\;,\end{array}$$

$$\displaystyle\begin{array}[]{c}\displaystyle\partial_{x}(\partial_{t}\mathbf{Q}_{j})=-\sum_{j=1}^{m}\sum_{l=1}^{m}\frac{\partial\mathbf{A}_{i,j}(\mathbf{Q})}{\partial\mathbf{Q}_{l}}\partial_{x}\mathbf{Q}_{l}\partial_{x}\mathbf{Q}_{j}-\sum_{j=1}^{m}\mathbf{A}_{i,j}(\mathbf{Q})\partial_{x}^{(2)}\mathbf{Q}_{j}\\ \displaystyle+\sum_{j=1}^{m}\mathbf{B}(\mathbf{Q})_{i,j}\partial_{x}\mathbf{Q}_{j}\;,\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{t}\mathbf{Q}=-\mathbf{A}(\mathbf{Q})\partial_{x}\mathbf{Q}+\mathbf{S}(\mathbf{Q})\;.\par\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{t}\mathbf{Q}=-\mathbf{A}(x,t)\partial_{x}\mathbf{Q}+\mathbf{S}(\mathbf{Q})\;,\par\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{t}(\partial_{t}\mathbf{Q})=-\mathbf{A}_{t}\partial_{x}\mathbf{Q}-\mathbf{A}\partial_{t}(\partial_{x}\mathbf{Q})+\mathbf{B}\partial_{t}\mathbf{Q}\;,\par\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{t}^{(2)}\mathbf{Q}=-\mathbf{A}_{t}\partial_{x}\mathbf{Q}-\mathbf{A}\partial_{x}(\partial_{t}\mathbf{Q})+\mathbf{B}\partial_{t}\mathbf{Q}\;.\par\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{x}(\partial_{t}\mathbf{Q})=-\mathbf{A}_{x}\partial_{x}\mathbf{Q}-\mathbf{A}\partial_{x}^{(2)}\mathbf{Q}+\mathbf{B}\partial_{x}\mathbf{Q}\\ =-\mathbf{A}\partial_{x}^{(2)}\mathbf{Q}+(\mathbf{B}-\mathbf{A}_{x}^{(1)})\partial_{x}\mathbf{Q}\;.\par\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{t}^{(2)}\mathbf{Q}=-\mathbf{A}_{t}\partial_{x}\mathbf{Q}-\mathbf{A}(-\mathbf{A}\partial_{x}^{(2)}\mathbf{Q}+(\mathbf{B}-\mathbf{A}_{x}^{(1)})\partial_{x}\mathbf{Q})+\mathbf{B}\partial_{t}\mathbf{Q}\\ =\mathbf{A}^{2}\partial_{x}^{(2)}\mathbf{Q}+(-\mathbf{A}_{t}-\mathbf{A}(\mathbf{B}-\mathbf{A}_{x}^{(1)}))\partial_{x}\mathbf{Q}+\mathbf{B}\partial_{t}\mathbf{Q}\;.\par\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{t}^{(2)}\mathbf{Q}=\mathbf{C}(2,2)\partial_{x}^{(2)}\mathbf{Q}+\mathbf{C}(2,1)\partial_{x}\mathbf{Q}+\mathbf{B}\partial_{t}\mathbf{Q}\;,\par\end{array}$$

$$\displaystyle\begin{array}[]{c}\mathbf{C}(2,2)=\mathbf{A}^{2}\;,\mathbf{C}(2,1)=-\mathbf{A}_{t}-\mathbf{A}(\mathbf{B}-\mathbf{A}_{x}^{(1)})\;.\par\end{array}$$

$$\displaystyle\begin{array}[]{cl}\partial_{x}(\partial_{t}\mathbf{Q})=&-\mathbf{A}\partial_{x}^{(2)}\mathbf{Q}+(\mathbf{B}-\mathbf{A}_{x})\partial_{x}\mathbf{Q}\;,\\ \partial_{x}^{(2)}(\partial_{t}\mathbf{Q})=&-\mathbf{A}\partial_{x}^{(3)}\mathbf{Q}+(\mathbf{B}-2\mathbf{A}_{x})\partial_{x}^{(2)}\mathbf{Q}+(\mathbf{B}_{x}-\mathbf{A}_{x}^{(2)})\partial_{x}\mathbf{Q}\;,\\ \partial_{x}^{(3)}(\partial_{t}\mathbf{Q})=&-\mathbf{A}\partial_{x}^{(4)}\mathbf{Q}+(\mathbf{B}-3\mathbf{A}_{x})\partial_{x}^{(3)}\mathbf{Q}+(2\mathbf{B}_{x}-3\mathbf{A}_{x}^{(2)})\partial_{x}^{(2)}\mathbf{Q}\\ &+(\mathbf{B}_{x}^{(2)}-\mathbf{A}_{x}^{(3)})\partial_{x}\mathbf{Q}\;,\\ \partial_{x}^{(4)}(\partial_{t}\mathbf{Q})=&-\mathbf{A}\partial_{x}^{(5)}\mathbf{Q}+(\mathbf{B}-4\mathbf{A}_{x})\partial_{x}^{(4)}\mathbf{Q}+(3\mathbf{B}_{x}-6\mathbf{A}_{x}^{(2)})\partial_{x}^{(3)}\mathbf{Q}\\ &+(3\mathbf{B}_{x}^{(2)}-4\mathbf{A}_{x}^{(3)})\partial_{x}^{(2)}\mathbf{Q}+(\mathbf{B}_{x}^{(3)}-\mathbf{A}_{x}^{(4)})\partial_{x}\mathbf{Q}\;,\\ \partial_{x}^{(5)}(\partial_{t}\mathbf{Q})=&-\mathbf{A}\partial_{x}^{(6)}\mathbf{Q}+(\mathbf{B}-5\mathbf{A}_{x})\partial_{x}^{(5)}\mathbf{Q}+(4\mathbf{B}_{x}-10\mathbf{A}_{x}^{(2)})\partial_{x}^{(4)}\mathbf{Q}\\ &+(6\mathbf{B}_{x}^{(2)}-10\mathbf{A}_{x}^{(3)})\partial_{x}^{(3)}\mathbf{Q}+(4\mathbf{B}_{x}^{(3)}-5\mathbf{A}_{x}^{(4)})\partial_{x}^{(2)}\mathbf{Q}\\ &+(\mathbf{B}_{x}^{(4)}-\mathbf{A}_{x}^{(5)})\partial_{x}\mathbf{Q}\;,\\ \end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{x}^{(l)}(\partial_{t}\mathbf{Q})=\par\sum_{k=1}^{l+1}(b_{l,k}\mathbf{B}_{x}^{(l+1-k)}-a_{l,k}\mathbf{A}_{x}^{(l+2-k)})\partial_{x}^{(k)}\mathbf{Q}\;.\par\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{x}^{(l)}(\partial_{t}\mathbf{Q})=\displaystyle\sum_{k=1}^{l+1}\par\mathbf{D}(l+1,k)\par\partial_{x}^{(k)}\mathbf{Q}\;,\end{array}$$

$$\displaystyle\begin{array}[]{c}\mathbf{D}(l+1,k)=\biggl{(}\par\left(\begin{array}[]{c}l-1\\ l-k\end{array}\right)\par\mathbf{B}_{x}^{(l-k)}\par-\par\left(\begin{array}[]{c}l\\ l+1-k\end{array}\right)\mathbf{A}_{x}^{(l+1-k)}\par\par\biggr{)}\par\end{array}$$

$$\displaystyle\begin{array}[]{c}\left(\begin{array}[]{c}l\\ -1\end{array}\right)=0\;,\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{x}^{(l)}(\partial_{t}\mathbf{Q})=\displaystyle\sum_{k=1}^{l+1}\biggl{(}\par\left(\begin{array}[]{c}l-1\\ l-k\end{array}\right)\par\mathbf{B}_{x}^{(l-k)}\par-\par\left(\begin{array}[]{c}l\\ l+1-k\end{array}\right)\mathbf{A}_{x}^{(l+1-k)}\par\par\biggr{)}\partial_{x}^{(k)}\mathbf{Q}\;.\end{array}$$

$$\displaystyle\begin{array}[]{ll}\partial_{x}^{(l+1)}(\partial_{t}\mathbf{Q})&=\par\displaystyle\sum_{k=1}^{l+1}\par\biggl{(}\par\left(\begin{array}[]{c}l-1\\ l-k\end{array}\right)\par\mathbf{B}_{x}^{(l+1-k)}\par-\par\left(\begin{array}[]{c}l\\ l+1-k\end{array}\right)\mathbf{A}_{x}^{(l+2-k)}\par\biggr{)}\partial_{x}^{(k)}\mathbf{Q}\par\par\par\\ &\displaystyle+\sum_{k=1}^{l+1}\biggl{(}\par\left(\begin{array}[]{c}l-1\\ l-k\end{array}\right)\mathbf{B}_{x}^{(l+1-k)}\par-\par\left(\begin{array}[]{c}l\\ l+1-k\end{array}\right)\mathbf{A}_{x}^{(l+1-k)}\par\par\biggr{)}\partial_{x}^{(k+1)}\mathbf{Q}\par\\ &=\par\par\biggl{(}\par\left(\begin{array}[]{c}l-1\\ l-1\end{array}\right)\par\mathbf{B}_{x}^{(l+1)}\par-\par\left(\begin{array}[]{c}l\\ l\end{array}\right)\mathbf{A}_{x}^{(l+2)}\par\biggr{)}\partial_{x}\mathbf{Q}\par\\ &+\par\displaystyle\sum_{k=2}^{l+1}\par\bigg{[}\par\left(\begin{array}[]{c}l-1\\ l+1-k\end{array}\right)+\left(\begin{array}[]{c}l-1\\ l-k\end{array}\right)\par\biggr{]}\mathbf{B}_{x}^{(l+1-k)}\par\partial_{x}^{(k)}\mathbf{Q}\par\\ &-\displaystyle\sum_{k=2}^{l+1}\biggl{[}\left(\begin{array}[]{c}l\\ l+2-k\end{array}\right)+\left(\begin{array}[]{c}l\\ l+1-k\end{array}\right)\par\biggr{]}\mathbf{A}_{x}^{(l+2-k)}\par\par\partial_{x}^{(k)}\mathbf{Q}\par\\ &+\par\displaystyle\biggl{(}\par\left(\begin{array}[]{c}l-1\\ -1\end{array}\right)\par\mathbf{B}_{x}^{(-1)}\par-\par\left(\begin{array}[]{c}l\\ 0\end{array}\right)\mathbf{A}_{x}^{(0)}\par\biggr{)}\partial_{x}^{(l+2)}\mathbf{Q}\;.\end{array}$$

$$\displaystyle\begin{array}[]{c}\left(\begin{array}[]{c}l-1\\ l+1-k\end{array}\right)+\left(\begin{array}[]{c}l-1\\ l-k\end{array}\right)\par=\par\left(\begin{array}[]{c}l\\ l+1-k\end{array}\right)\;,\par\\ \\ \left(\begin{array}[]{c}l-1\\ l-1\end{array}\right)=\left(\begin{array}[]{c}l\\ l\end{array}\right)=\left(\begin{array}[]{c}l\\ 0\end{array}\right)=\left(\begin{array}[]{c}l-1\\ 0\end{array}\right)=1\end{array}$$

$$\displaystyle\begin{array}[]{c}\left(\begin{array}[]{c}m\\ -1\end{array}\right)=0\;,\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{x}^{(l+1)}(\partial_{t}\mathbf{Q})=\displaystyle\sum_{k=1}^{l+2}\biggl{(}\par\left(\begin{array}[]{c}l\\ l+1-k\end{array}\right)\par\mathbf{B}_{x}^{(l+1-k)}\par-\par\left(\begin{array}[]{c}l+1\\ l+2-k\end{array}\right)\mathbf{A}_{x}^{(l+2-k)}\par\par\biggr{)}\partial_{x}^{(k)}\mathbf{Q}\;.\end{array}$$

$$\displaystyle\begin{array}[]{c}\displaystyle\par\partial_{t}^{(k)}\mathbf{Q}=\sum_{l=1}^{k}\mathbf{C}(k,l)\partial_{x}^{(l)}\mathbf{Q}+\partial_{t}^{(k-2)}(\mathbf{B}\partial_{t}\mathbf{Q})\;,\par\end{array}$$

$$\displaystyle\begin{array}[]{c}\mathbf{C}(k,l)=\left\{\begin{array}[]{cc}\mathbf{C}(k-1,k-1)\mathbf{D}(k,k),&l=k\;,\\ \mathbf{C}(k-1,l)_{t}+\sum_{m=l-1}^{k-1}\mathbf{C}(k-1,m)\mathbf{D}(m+1,l),&l<k\;,\\ \end{array}\right.\end{array}$$

$$\displaystyle\begin{array}[]{c}\partial_{t}\mathbf{Q}=-\mathbf{A}\partial_{x}\mathbf{Q}+\mathbf{S}(\mathbf{Q})\;.\end{array}$$

Full text

††journal: Computers and Fluids

1 Introduction

This paper concerns the solution of Generalized Riemann Problems (GRP) in the context of high-order finite volume methods. The ADER (Arbitrary Accuracy DERivative Riemann problem method), first put forward by Toro et al. [20], is of particular interest in this work. The method in [20] was devoted to develop a procedure able to compute the numerical solution, of the one-dimensional linear advection problem, of arbitrary order of accuracy in both space and time. This method can be considered as a generalization of Godunov’s method, where the numerical fluxes can be obtained from the local solution of GRP where the initial condition consists of polynomial functions of suitable order. Subsequently, ADER was extended to solve linear systems of hyperbolic conservation laws in [26, 31]. The ADER philosophy was extended by Toro and Titarev in [34] to solve the non-linear systems; inhomogeneous Burgers equation and the nonlinear shallow-water equations with variable bed elevation. In [29], the ADER approach was extended, by Toro and Titarev, to nonlinear but homogeneous hyperbolic systems. Furthermore, the extension of ADER to scalar balance laws was investigated by Toro and Takakura, [28].

The original ADER scheme [20, 29], has two main steps, reconstruction and flux calculation. The marching in time generates cell averages, then the reconstruction procedure generates a special type of interpolation polynomial of the solution from cell averages. The flux calculation is carried out from the solution of the GRP, which is proposed in terms of a Taylor series expansion in time, where the time derivatives are completely expressed in terms of spatial derivatives by means of the Cauchy-Kowalewskaya or Lax-Wendroff procedure. The space derivatives are obtained from homogeneous linearized Riemann problems constructed from the governing equation and an initial condition given by the derivatives of interpolation polynomials. This was the summary of the pioneering ADER method, which in principle is able to generate approximations of arbitrary order of accuracy.

The accuracy of ADER methods depend on the number of terms in the Taylor series expansions. In the particular case of the first order, it recovers the Godunov method and for the second order, ADER recovers the second-order GRP method of Ben-Artzi and Falcoviz [5]. In [6], a re-interpretation carried out by Castro and Toro of the high-order numerical method proposed by Harten et al. [17] has allowed us to formulate GRP solutions and thus ADER schemes in a different way. In this new interpretation of ADER, a classical Riemann problem is built from the governing equation and a piece-wise initial condition which is formed from two constant states. In this approach, the constant states are local predictors of the solution within computational cells, that is, these are extrapolation values of the solution at both sides of the cell interfaces at a given time. The evolution of these extrapolated values is carried out by using Taylor series expansion in time, where the time derivatives are still expressed via Cauchy-Kowalewskaya functionals but filled with the spatial derivates of the reconstruction interpolation functions. In this form the predictors in two adjacent computational cells are interacted at the cell interface through the classical Riemann problem. In [6], the approach based on the Harten et al. is called the HEOC solver and the original GRP solver of Toro and Titarev scheme is referred to as the TT solver. So, the difference between the HEOC formulation and the TT formulation is that, in the HEOC case only one classical Riemann problem is required but it needs to be solved at each quadrature point, whereas, in the TT approach a sequence of classical Riemann problems are required, only once, one for the leading term of Taylor expansions and linearised Riemann problems for the spatial derivatives. So, after these terms are available the computation of GRP solution via TT is reduced jut to evaluate a polynomial in time. The similarities between both approaches are the use of Taylor series expansion and the use of the Cauchy-Kowalewskaya procedure. A detailed review of GRP solvers is done in [6, 21]. Similarly, further details of ADER schemes can be found in Chapters 19 and 20 of the textbook by Toro [30]. Notice that, the GRP solvers require the ability of solving classical Riemann problems. The number of hyperbolic systems where the exact solution of Riemann problems is available, is limited. In general, the exact solution of Riemann problems for several hyperbolic system can be very difficult to be obtained. Fortunately, the ADER approaches described above can use approximate Riemann solvers, see [16]. In this sense the search for approximate Riemann solvers for general hyperbolic system is very relevant area of research, where the ADER philosophy can benefit. Balsara [1, 2] and Balsara et al. [3] have extended multidimensional HLL and HLLC to Euler and MHD equations. Goetz et al. [15] have shown that approximate Riemann solver can be obtained from the well-known HLL solver. See also [13, 4] where universal approximate GRP solvers based on HLL method and the inclusion of intermediate waves, have been reported.

The ADER approach allows flexibility to incorporate the finite element approach into the finite volume framework. In an intermediate stage in the mixing between finite volume and finite element approach, Dumbser and Munz, [10, 11], have implemented the ADER approach for the Discontinuous Galerkin approach applied to the aeroacoustics and the Euler equations in two dimensions. In this approach, the ADER is used to evolve polynomials in space and time, through the evolution of their degrees of freedom by using Taylor series expansion and the Cauchy-Kowalewskaya procedure but instead of using reconstruction and the derivatives of reconstruction polynomials, the authors proposed to use the test functions of the finite element space. In this sense, the work of Dumbser et al. [9] in their pioneering work, has introduced the Galerkin framework for obtaining the predictor within cells. The difference between this approach and that of Dumbser and Munz is that this neither requires the use of any Taylor series expansion nor Cauchy-Kowalewskaya procedure. Some contributions to the developments of this class of solver can be found in [7, 8, 12, 14, 24], to mention but a few.

The methodologies based on Galerkin approaches require the inversion of matrices and the solution of non-linear algebraic equations which are very time-consuming processes. On the other hand, the methodologies based on Taylor series suffer from the Cauchy-Kowaleskaya procedure which becomes cumbersome when the accuracy increases. Furthermore, the Taylor series expansions without modifications cannot deal with stiff source terms. In this sense, in [22, 32] Montecinos and Toro have introduced the implicit Taylor series expansion and Cauchy-Kowakeskaya procedure to deal with hyperbolic balance laws with stiff source terms. In this approach, the Cauchy-Kowakeskaya procedure requires the spatial derivatives evolved in time, which cannot be obtained straightforwardly from Riemann problems as in conventional ADER methods discussed above. This is basically because conventional ADER methods based on TT solver uses linearised Riemann problems which are homogeneous. Thus, the influence of the source term is only involved in the Cauchy-Kowalewskaya procedure. In [32] the spatial derivatives are evolved by using two approaches, which differ in the number of terms in the Taylor expansion, so Complete Implicit Taylor Approach (CITA) and Reduced Implicit Taylor Approach (RITA), in both HEOC and TT approaches, are investigated there. A limitation of the approach in [32], is the computational cost, for high orders of accuracy, the CPU time increases dramatically. However, second order approaches work very well as reported in [35] where an extension to transport phenomena on unstructured meshes has been reported.

In this paper, we propose a strategy related to the implicit Taylor series expansion and the Cauchy-Kowalewskaya procedure. The difference between the conentional approaches described above and the proposed one is that spatial derivatives do not need to be evolved and the Cauchy-Kowalewskaya procedure is modified by expressing high-order time derivatives not only in terms of spatial derivatives of the data but also on space and time derivatives of the Jacobian matrices of the flux and source functions as well. The derivatives are obtained from an interpolation fashion on selected nodal points. This simplification allows us to provide a closed form for the Cauchy-Kowalewskaya functionals, in a recursive formula. Furthermore, this approach requires the solution of one algebraic equation but the number of variables is the same for all orders of accuracy. Closed forms of the Cauchy-Kowaleskaya functional are available for some partial differential equations as linear advection systems, linear systems with constant matrices [25, 18] and for the non-linear two-dimensional Euler equation, [11]. However, the expression obtained here is useful for all hyperbolic balance laws.

This paper is organized as follows. In section 2, the general framework is presented. In section 3, the new predictor step is introduced. In the section 4, numerical tests are concerned. Finally, in section 5 conclusions and remarks are drawn.

2 The framework

In this paper we present a strategy for solving a hyperbolic balance law in the conservative form

[TABLE]

where $\mathbf{H}_{0}(x)$ is a prescribed function in $\mathbb{R}^{m}$. Here $\mathbf{Q}(x,t)\in\mathbb{R}^{m}$ is the vector of unknowns, $\mathbf{F}(\mathbf{Q})\in\mathbb{R}^{m}$ is the physical flux function and $\mathbf{S}(\mathbf{Q})\in\mathbb{R}^{m}$ is the source term.

To compute a numerical solution of (3), we divide computational domains into $N$ uniform cells of the form $I_{i}^{n}:=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[t^{n},t^{n+1}]$ and then by integrating on $I_{i}^{n}$ we obtain the well-known one-step formula.

[TABLE]

where

[TABLE]

and the numerical flux $\mathbf{F}_{i+\frac{1}{2}}$ as well as the source term $\mathbf{S}_{i}$ are computed by adopting the ADER strategy, [20, 33, 34, 30]. In this paper

[TABLE]

where $\mathbf{F}_{h}(\mathbf{Q}_{L},\mathbf{Q}_{R})$ is a classical numerical flux function, which can be seen as a function of two states $\mathbf{Q}_{L}$ and $\mathbf{Q}_{R}$. This can be an approximate resolved flux in a Riemann problem. It is possible to design an approximate GRP solver out of the HLL, [15]. In this paper we will use the Rusanov solver, obtained from HLL by taking extreme left and right maximum waves speed to be the same but in opposite directions. Here $\mathbf{Q}_{i}(x,t)$ corresponds to a predictor within the computational cell $I_{i}^{n}$. In the next section, further details of the predictor step are provided.

3 The predictor step

In this section we provide the details to obtain the predictor $\mathbf{Q}_{i}(x,t)$. We adopt the strategy of the implicit Taylor series expansion and the Cauchy-Kowalewskaya procedure presented in [32]. However, instead of the conventional Cauchy-Kowalewskaya procedure we use a simplified version of this procedure. For the sake of completeness, we provide a brief review of the approach in [32]. The predictor is computed as

[TABLE]

which by means of the Cauchy-Kowalewskaya procedure, can be written as

[TABLE]

where $M$ is an integer which corresponds to the order of accuracy $M+1$, in both space and time. Here, $\mathbf{G}^{(k)}$ corresponds to the Cauchy-Kowalewskaya functional and it is the function which expresses the time derivatives in terms of spatial derivatives, $\partial_{t}^{(k)}\mathbf{Q}_{i}(x,\tau)=\mathbf{G}^{(k)}(\mathbf{Q}_{i}(x,\tau),\partial_{x}\mathbf{Q}_{i}(x,\tau),...,\partial_{x}^{(k)}\mathbf{Q}_{i}(x,\tau))$. Notice that this functional requires the information of spatial derivatives at $\tau$, which must be estimated. In [32] two strategies are proposed, where the time spatial derivatives are obtained from implicit Taylor series as well. These approaches require the solution of algebraic equations.

In the next section we provide a brief review of the conventional Cauchy-Kowalewskaya procedure and subsequently, a simplification of this procedure is presented.

3.1 A brief review of the Cauchy-Kowalewskaya procedure

Here, we briefly describe the Cauchy-Kowalewskaya procedure for obtaining the time derivatives of the data. For the sake of simplicity, in this section we omit the sub index $i$ in $\mathbf{Q}_{i}$, to indicate the approximation within the cell $[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}].$

The time derivatives are obtained from the governing equation (3). As for example, the first derivative is given by

[TABLE]

where $\mathbf{A}(\mathbf{Q})$ is the Jacobian matrix of $\mathbf{F}(\mathbf{Q})$ with respect to $\mathbf{Q}$. So, to obtain the second time derivative we differentiate in time the equation (16), so we have

[TABLE]

where $\mathbf{B}(\mathbf{Q})$ is the Jacobian matrix of $\mathbf{S}(\mathbf{Q})$ with respect to $\mathbf{Q}$, $\partial_{t}\mathbf{Q}_{i}$ and $\partial_{t}\mathbf{Q}_{j}$ are the $i$th and the $j$th components of the vector state $\partial_{t}\mathbf{Q}$, respectively and then by differentiating the expression (16) with respect to $x$, we obtain

[TABLE]

the same procedure is applied to obtain $\partial_{t}^{(3)}\mathbf{Q}$ and so on, in principle any high order time derivative can be obtained through this procedure. However, the procedure becomes very cumbersome for derivatives of orders higher than two, furthermore the complexity scales with the number of unknowns $m$ and the order of accuracy as well. This justifies the requirement of finding some efficient strategy to approximate temporal derivatives by following the Cauchy-Kowalewskaya ideas.

3.2 The simplified Cauchy-kowalewskaya procedure

In this section, we derive a simplified Cauchy-kowalewskaya procedure, to approximate time derivatives. For the sake of simplicity, in this section we also omit the sub index $i$ in $\mathbf{Q}_{i}$, to indicate the approximation within the cell $[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}].$ As seen in the previous section, the conventional Cauchy-Kowalewsky procedure provides the first time derivative as

[TABLE]

At this point we introduce the first simplification. Instead of considering the previous equation, we are going to use the approximation

[TABLE]

which means, the matrix $\mathbf{A}$ is considered just a space-time dependent matrix.

By taking into account the approximation (26), the second time derivative can be approximated in three steps described below.

Step I. We differentiate $\partial_{t}\mathbf{Q}$ with respect to $t$, thus

[TABLE]

where $\mathbf{B}$ is the Jacobian matrix of the source term with respect to $\mathbf{Q}$. For the remaining part of this paper, we use the notation $\partial_{x}\mathbf{A}=\mathbf{A}_{x}$ for any matrix $\mathbf{A}$. Similarly, we use the convention $\partial_{x}^{(l)}\mathbf{A}=\mathbf{A}_{x}^{(l)}$ for the $l$-th partial derivative of the matrix $\mathbf{A}$ with respect to $x$. Do not confuse with $\mathbf{A}^{l}$ which means matrix multiplication, $l$ times.

Step II. At this point we introduce the second simplification, which is to consider also the matrix $\mathbf{B}$ as a space-time dependent matrix rather than a state dependent matrix. We assume regularity enough such that the spatial and time derivatives can be interchanged. So

[TABLE]

Step III. Here, we differentiate in space the expression $\partial_{t}\mathbf{Q}$, taking into account the simplifications introduced above, to obtain

[TABLE]

This step regards the main difference with respect to the conventional Cauchy-Kowalewskaya procedure. By inserting the previous expression into (30), we obtain

[TABLE]

The last expression can be written as

[TABLE]

where $\mathbf{C}(2,i)$, $i=1,2$ represent the matrix coefficients

[TABLE]

The novel contribution of this paper is as follows. We are going to show that this procedure can be generalized. Before to give this main result, we need to express $\partial_{x}^{(l)}(\partial_{t}\mathbf{Q})$ only in terms of spatial derivatives of the data and the Jacobian matrices. Notice that from simplifications introduced above, we obtain

[TABLE]

so by inspection we observe that these derivatives, we can be arranged as

[TABLE]

Notice that, $\partial_{x}^{(l)}$ stands by the $l$-th spatial derivative, with the convention $\partial_{x}^{(0)}\mathbf{M}=\mathbf{M}$ for any function $\mathbf{M}$, which may be a scalar, vector or matrix function. This notation is also extended to temporal derivatives.

The table 1, shows the coefficient $a_{l,k}$. Similarly, the table 2 shows the coefficients $b_{l,k}$. We observe that they follow the structure of the Pascal triangle, in the combinatorial theory. In fact, the structure is given by the following.

Lemma 3.1**.**

[TABLE]

where

[TABLE]

and

[TABLE]

for all integer $l$.

Proof.

Let us prove it by induction.

We already know that this is true for $l=1$.

• 2.

Let us assume it is true for $l$, that is

[TABLE]

• 3.

Let us prove it is true for $l+1$. Indeed

[TABLE]

By considering the properties of the combinatorial factors

[TABLE]

and the assumption

[TABLE]

for all $m$, after grouping terms we obtain

[TABLE]

This completes the proof.

∎

Proposition 3.2**.**

The high-order time derivatives have the following recursive form

[TABLE]

where

[TABLE]

here, the matrix $\mathbf{D}$ is given by (60). We impose $\mathbf{C}(k,0)=\mathbf{0}$ $\forall k>0$, $\mathbf{C}(1,1)=-\mathbf{A}$ and $\partial_{t}^{(-1)}(\mathbf{B}\partial_{t}\mathbf{Q})=\mathbf{S}(\mathbf{Q})$.

Proof.

Let us prove this proposition by induction.

The result is true for $k=2$. In fact, we know that

[TABLE]

Since, $\partial_{t}^{(2)}\mathbf{Q}=\partial_{t}(\partial_{t}\mathbf{Q})$ from (137) we obtain

[TABLE]

Here, we have used the chain rule

[TABLE]

Therefore, from the expressions $\mathbf{C}(1,1)=-\mathbf{A}$, $\mathbf{D}(2,2)=-\mathbf{A}$, $\mathbf{D}(2,1)=\mathbf{B}-\mathbf{A}_{x}$ and by identifying terms, the induction hypothesis is valid for $k=2$.

• 2.

We assume the induction hypothesis is valid for $k=n$ and thus

[TABLE]

for all $k\leq n$.

• 3.

Let us prove this is valid for $k=n+1$. In fact

[TABLE]

with $\mathbf{C}(n,0)=\mathbf{0}$. So by collecting terms and defining

[TABLE]

we can write (150) as

[TABLE]

this proves the sought result.

Notice that the condition $\partial_{t}^{(-1)}(\mathbf{B}\partial_{t}\mathbf{Q})=\partial_{t}^{(-1)}(\partial_{t}\mathbf{S}(\mathbf{Q}))=\partial_{t}^{(0)}\mathbf{S}(\mathbf{Q})=\mathbf{S}(\mathbf{Q})$ is natural, which also justifies the expression (137), it also corresponds to $k=1$ in the formula (131).

∎

Notice that the expression (131) is only possible from the simplifications proposed in this work. This is not possible, in general, for the conventional Cauchy-Kowalewskaya procedure. Notice that (131) expresses the time derivatives in terms of spatial derivatives of $\mathbf{Q}$, space and time derivatives of both $\mathbf{A}$ and $\mathbf{B}$.

Corollary 3.3**.**

The expression (131) can be written as

[TABLE]

where

[TABLE]

Proof.

This result follows from the manipulation of (131) in Proposition 3.2. Particularly the term $\partial_{t}^{(k-2)}(\mathbf{B}\partial_{t}\mathbf{Q})$ can be expressed, by using the result in Proposition A.1, as

[TABLE]

By collecting terms and isolating for $l=k-1$, we obtain

[TABLE]

and thus the result holds. ∎

Proposition 3.4**.**

[TABLE]

where $\mathbf{M}_{r}$ are those in (162).

Proof.

This is a consequence of the corollary 3.3.

∎

As will be seen in next sections, the previous results provide the closed form for approximations to Cauchy Kowalewskaya functionals, $\mathbf{G}^{(k)}$, which will be important to design fixed-point iteration procedures. The appendix B.1 shows operational details for generating the matrix $\mathbf{C}(k,l)$ involved into the simplified Cauchy-Kowaleskaya procedure introduced in this section.

3.3 The predictor step based on a modified implicit Taylor series expansion

Notice that the predictor $\mathbf{Q}_{i}$, within $I_{i}^{n}$, is required for evaluating integrals in (10). On the other hand, the evaluation of these integrals is carried out by means of quadrature rules in space and time. So, for the temporal integration we use the Gaussian rule, which involves $\tau_{j}$, $j=1,...,n_{T}$ Gaussian points. Whereas, for the spatial integration, we use the Newton-Cotes rule, which involves $\xi_{m}$, $m=1,...,n_{S}$ equidistant quadrature points. See appendix B.2 for further details about the set up of these quadrature points through reference elements. This set of quadrature points allows us to build space-time nodal points $(\xi_{m},\tau_{j})$ within the space-time cell $I_{i}^{n}$, as illustrated in the figure 1. So, from the previous comment it is evident that for flux and source evaluations, we only need the information of $\mathbf{Q}_{i}$ at $(\xi_{m},\tau_{j})$.

To obtain approximation of the predictor at every space-time node $(\xi_{m},\tau_{j})$, we propose the following strategy.

Provide a starting guess for $\mathbf{Q}_{i}(\xi_{m},\tau_{j})$, $m=1,...,n_{S}$ and $j=1,...,n_{T}$.

This is done by using the formula

[TABLE]

which corresponds to the second order accurate expression in [32]. Here, $\mathbf{W}_{i}(\xi)$ represents the reconstruction polynomial obtained within the space-time cell $I_{i}^{n}$. Any reconstruction procedure can be implemented, however, in this work we use the Weighted Essentially Non-Oscillatory (WENO) reconstruction method described in [9]. 2. 2.

Compute the approximation of high-order derivatives in time and space as well of the state and Jacobian matrices. For this purpose, we use the following approach.

Let $\mathbf{M}$ be a function, which may represent the state function $\mathbf{Q}$ and the Jacobian matrices $\mathbf{A}$ and $\mathbf{B}$ as well.

• (a)

Then, for obtaining $\mathbf{M}_{x}^{(l)}(\xi,\tau_{j})$, we first interpolate the function $\mathbf{M}$ on the nodes $(\xi_{m},\tau_{j})$ with $j$ fix and varying $m=1,...,n_{S}$. So an interpolation function $\tilde{\mathbf{M}}(\xi,\tau_{j})$ is obtained. Then, we are able to provide approximations of spatial derivatives of $\mathbf{M}(\xi,\tau_{j})$ for any order $l$ by using the spatial derivatives of $\tilde{\mathbf{M}}(\xi,\tau_{j})$.

• (b)

Similarly, to obtain $\mathbf{M}_{t}^{(l)}(\xi_{m},\tau)$, we first interpolate the function $\mathbf{M}$ on the nodes $(\xi_{m},\tau_{j})$ with $m$ fix and varying $j=1,...,n_{T}$. So an interpolation function $\tilde{\mathbf{M}}(\xi_{m},\tau)$ is obtained. Then, we are able to provide approximations of temporal derivatives of $\mathbf{M}(\xi_{m},\tau)$ for any order $l$ by using the temporal derivatives of $\tilde{\mathbf{M}}(\xi_{m},\tau)$.

In the appendix B.2, is shown the form of these interpolation polynomials for the orders of accuracy considered in this paper. 3. 3.

Update $\mathbf{Q}_{i}$ at every $(\xi_{m},\tau_{j})$ by using

[TABLE]

where $\tilde{\mathbf{G}}^{(k)}=\tilde{\mathbf{G}}^{(k)}(\mathbf{Q}_{i},...,\partial_{x}^{(k)}\mathbf{Q}_{i},\mathbf{A}_{t}^{(l)},..,\mathbf{A}_{x}^{(l)},..,\mathbf{B}_{t}^{(l)},...,\mathbf{B}_{x}^{(l)},...)$ is given by (131). The derivatives of the state function and matrices are evaluated at $(\xi_{m},\tau_{j})$ and computed in the previous step. The equation (176) corresponds to the implicit Taylor series expansion in [32] with the difference that $\tilde{\mathbf{G}}^{(k)}$ is a simplification of the conventional Cauchy-Kowalewskaya functional.

For solving (176), we use the following nested Picard iteration procedure.

[TABLE]

where $s$ is an iteration index and $\mathbf{M}_{l}$ comes from the proposition 3.4. We have omitted the arguments of $\mathbf{Q}_{i}$.

To solve it, we build an algebraic system, which has the form

[TABLE]

So, the update of $\mathbf{Q}_{i}$ is carried out as $\mathbf{Q}_{i}^{s+1}=\mathbf{Q}_{i}^{s}-\delta$, where $\delta$ is the solution to $\mathcal{J}(\mathbf{Q}_{i}^{s})\delta=\mathcal{H}(\mathbf{Q}_{i}^{s})$, where

[TABLE]

is the Jacobian matrix of $\mathcal{H}(\mathbf{Y})$ with respect to $\mathbf{Y}$. The update is carried out $M$ times, where $M+1$ corresponds to the order of accuracy. Notice that the same algebraic equation for $\mathbf{Y}$ has to be solved for any order of accuracy. So, the size of $\mathbf{Y}$ does not depend on the accuracy. 4. 4.

Go to step 2. Finish if the global loop has been done by $M$ times. By virtue of the efficiency we use a limited number of iterations. From experiments, not shown here, the result does not vary in terms of accuracy if a stop criterion, based on the tolerance for the relative error between subsequent approximations, is implemented.

Once $\mathbf{Q}_{i}$ is computed for each cell $I_{i}^{n}$, the numerical flux and source terms can be easily evaluated. In the appendix B.3, is shown the form in which the integrals in (10) are evaluated.

This completes the description of the proposed strategy for obtaining the predictor within the computational cell $I_{i}^{n}$ by using the implicit GRP approach.

4 Numerical results

In this section we shall consider numerical test aimed at assessing the accuracy and performance of the present scheme. The time step, $\Delta t$, will be computed by using the well-known CFL condition

[TABLE]

where $\lambda_{abs}=\max_{i}(\max_{j}(|\lambda_{j}(\mathbf{Q}_{i}^{n})|))$, here $\lambda_{j}$, $j=1,...,m$ are the eigenvalues of the Jacobian matrix of $\mathbf{F}$ evaluated at $\mathbf{Q}_{i}^{n}$, the data at each cell $I_{i}^{n}$ and the maximum is taken over all cells $I_{i}^{n}$. On the hand, to assess the empirically the convergence rate, we are going to use the norm

[TABLE]

where $\mathbf{Q}^{e}(x,t)$ is the exact solution, $\mathbf{W}_{i}(x)$ is the reconstruction polynomial within the interval $[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$ obtained at the output time of the global simulation. In this paper we are going to use (190) with $p=1$, $p=2$ and the maximum norm given by

[TABLE]

Let us remark that, in the case of hyperbolic systems in which the solution vector $\mathbf{Q}$ contains more than one variables, we carry out the computation of errors as indicated above but for some particular component of the solution vector.

4.1 A linear system of hyperbolic balance laws

Here, we consider the linear system in [22], given by

[TABLE]

where

[TABLE]

The problem is endowed with periodic boundary conditions. This system has the exact solution

[TABLE]

where

[TABLE]

Here we consider $\lambda=1$ and $\beta=-1$. This is a simple test aimed at evaluating the accuracy of the present scheme. As can be seen in the Table 3, the expected theoretical orders of accuracy are achieved.

4.2 A system of non-linear hyperbolic balance laws

Here we assess the present methods, applied to the non-linear system

[TABLE]

where $\mathbf{F}(\mathbf{Q})$ and $\mathbf{S}(\mathbf{Q})$ are given by

[TABLE]

where $\beta\leq 0$ is a constant value, see [32]. The exact solution is given by

[TABLE]

where $w_{1}$ and $w_{2}$ are the solutions to

[TABLE]

where the initial condition for each equation is

[TABLE]

Notice that system (225) requires the solution of the Burgers equation with a non linear source term, in [32] this solution is reported. Table 4, shows the empirical orders of accuracy and the CPU times. Comparing with CPU times of the implicit Taylor series expansion and conventional Cauchy-Koealewskaya procedure in [32], we observe that the present scheme depicts important improvements in the performance. An improvement of one order of magnitude compared with the strategies in [32], is obtained.

4.3 The LeVeque and Yee test

Here, we apply our schemes to the well-known and challenging scalar test problem proposed by LeVeque and Yee [19], given by

[TABLE]

We solve this PDE on the computational domain $[0,1]$ with transmissive boundary conditions and the initial condition given by

[TABLE]

The solution on the characteristic curves satisfies de ordinary differential equation $\frac{d(x(t),t)}{dt}=\beta q(x(t),t)(q(x(t),t)-1)(q(x(t),t)-\frac{1}{2})$, which has two stable solutions $q\equiv 0$ and $q\equiv 1$ and one unstable solution in $q\equiv\frac{1}{2}$ where any solution trays to away from this. Similarly, any solution associated to characteristic curves necessarily must converge to one of the two stable solutions. On the other hand, a numerical scheme which is not able to solve stiff source terms, may introduce an excessive numerical diffusion and so the numerical solution, following characteristic curves, converges to the wrong stable solution. This penalizes the right propagation. Figure 2 shows the comparison between the exact solution and the numerical approximations provided by the present scheme of second, third, fourth and fifth orders of accuracy. The figure shows a good agreement for $\beta=-10000$ at $t_{out}=0.3$, which correspond to the stiff regime. We have used $C_{cfl}=0.2$ and $300$ cells. This test illustrates the ability of the present scheme for solving hyperbolic balance laws with stiff source terms.

4.4 The Euler equations

Now let us consider the Euler equations, given by

[TABLE]

where the pressure $p$ is related with the conserved variables through the equation

[TABLE]

for an ideal gas $\gamma=1.4$. Notice that, the choice of the initial condition given by the functions

[TABLE]

provides the exact solution for the system (239), which corresponds to the set of functions

[TABLE]

Notice that, the variables $\rho,u,p$ correspond to the non-conservative variables, the corresponding translation to conserved variables needs to be done. This test has a complex eigenstructure, which is a challenge for numerical methods. Table 5, shows the results of the empirical convergence rate assessment for the density variable $\rho$, at $t_{out}=1$ and $C_{cfl}=0.9$, we observe that the scheme achieves the expected theoretical orders of accuracy.

4.5 The Shu and Osher test

Here, we consider the test problem, first time proposed by Shu and Osher in [27], which is given by (239) and the initial condition given, in terms of non-conserved variables $\mathbf{W}=[\rho,u,p]^{T}$, as

[TABLE]

The problem is solved on $[-1,1]$ up to the output time $t_{out}$, see [14] for further details. Figure 3 shows a comparison between a reference solution and numerical approximations. The reference solution has been obtained with the scheme of third order using $2000$ cells. The numerical results correspond to second and third orders of accuracy for which we have used $300$ cells, $C_{cfl}=0.5$ and $t_{out}=0.47$. This test illustrates the ability of the present scheme for solving complex fluids, a good agreement is observed for the scheme of second and third order of accuracy on the smooth region, whereas, the third order scheme present a better performance in both the smooth region and the high frequency region as well.

5 Conclusions

In this paper, a simplified Cauchy-Kowalewkaya procedure has been proposed. The strategy uses not only the spatial derivatives of the data but also de derivatives of the Jacobian matrices in both space and time. The simplification allows us to propose a recursive formula which requires the ability of obtaining time and space derivatives of the data as well as matrices. This is achieved by using interpolations within a suitable arrangement of nodal points, which allows us to extract the information for flux and source term evaluations in a straightforward manner. This method is implemented in the context of GRP’s solvers based on implicit Taylor series expansions. The solver in [32] uses the same elements, that is, implicit Taylor series expansions and Cauchy-Kowaleskaya procedure. The Taylor expansion is used to obtain the data and the evolution of spatial derivatives required for the scheme. However, in the present approach the Taylor series expansion is used only once and the evolution of the space derivatives is not required. Despite, in both approaches, that is in [32] and the present one, the solution of an algebraic equation is required, that in [32] increases the number of unknowns as the accuracy increases, whereas, the number of unknowns in the present approach remains constant when the accuracy increases. We have implemented the GRP solver in the context of ADER methods and several tests reported in the literature have been solved. An empirical convergence rate assessment has been done for some of them. We have observed that the performance of the present scheme is at least one order of magnitude cheaper than schemes in [32]. Furthermore, the expected theoretical orders of accuracy have been achieved up to fifth order of accuracy. The extension to hyperbolic systems in 2D and 3D is the subject of ongoing research.

Acknowledgements

G.M thanks to the National chilean Fund for Scientific and Technological Development, FONDECYT, in the frame of the research project for Initiation in Research, number 11180926.

Appendix A Matrix-vector multiplication

In this appendix, we provide the algebraic details required for the main results in section 3.2.

Proposition A.1**.**

[TABLE]

with the convention $\mathbf{A}_{t}^{(0)}=\mathbf{A}$, here $\mathbf{A}\cdot\mathbf{B}$ can be any matrix-matrix or matrix-vector multiplication.

Proof.

The proof will be carried out by mathematical induction. Let us simplify the notation by considering matrices and vectors as elements which only depend on the variable $t$. Furthermore, we are going to assume that both $\mathbf{A}$ and $\mathbf{B}$ are regular enough in the variable $t$.

Let us verify for $l=1$. In fact, it is asy to verify that

[TABLE]

• 2.

Let us assume the formula is valid up to $l=n$, that is

[TABLE]

and we are going to prove, in the next step, that this is valid also for $l=n+1$.

• 3.

Since $\mathbf{A}$ and $\mathbf{B}$ are regular enough, we have

[TABLE]

by expanding this expression we have

[TABLE]

By collecting terms, and using elemental combinatorial algebra we have

[TABLE]

Thus the proof is completed.

∎

Notice that the proof is carried out for the derivative with respect to the variable $t$, however, this result applies to derivatives with respect to any variable.

Appendix B Operational details

In this appendix, Fortran 90 codes of the recursive formula for obtaining time derivatives are reported. Furthermore, the set up of quadrature points through reference elements and how these are used to obtain polynomials for approximating the spatial and temporal derivatives of high-order, are provided. The numerical flux and the source term evaluation strategy is also reported.

B.1 Fortran codes for computing the matrix coefficients $\mathbf{D}$ and $\mathbf{C}$

The codes in Fortran 90 for the main expressions obtained in the section 3.2, are reported.

The subroutines $MATRIX\_D$ and $MATRIX\_C$, provide the matrix $D$ and $C$ corresponding to the formulas (60) and (135), respectively.

! ------------------------------------------------------------------------ SUBROUTINE MATRIX_D ( l, k, i, j, A, B, DxA, DxB, D) ! ---------------------------------------------------------------------- ! This subroutine computes the matrix D for the recursive simplified ! Cauchy-Kowalewskaya procedure in Proposition 3.2 and Lemma 3.1. ! The inputs of this subroutine are: ! ! l and k which are the indices in the equation (21). ! i, j are the index for the space and time nodal location. ! DxA, DxB, the spatial derivatives of matrix A and B. ! ! This subroutine has been implemented with the following global variables: ! "Accuracy" is the accuracy of the scheme, which coincides with the ! number of quadrature points in space. ! "nGP" the number of quadrature points for time, here nGP = Accuracy-1. ! "NVAR" number of unknowns. ! ---------------------------------------------------------------------- DOUBLE PRECISION DxA ( NVAR, NVAR, Accuracy, nGP, Accuracy-1), & & DxB ( NVAR, NVAR, Accuracy, nGP, Accuracy-1), & & A ( NVAR, NVAR, Accuracy, nGP), & & B ( NVAR, NVAR, Accuracy, nGP), & & D ( NVAR, NVAR)

INTEGER i, j, k, l
! ----------------------------------------------------------------------
! Variables passing this point are local.
! ----------------------------------------------------------------------
DOUBLE PRECISION comB, comA
DOUBLE PRECISION DA ( NVAR, NVAR), DB ( NVAR, NVAR)
! ----------------------------------------------------------------------

if (l - k - 1 < 0) then
   DB = 0.0
elseif (l - k - 1 == 0) then
   DB = B ( :, :, i, j)
else
   DB = DxB ( :, :, i, j , l - k - 1)
end if

!--------------------------------------------------------------------
! Notice that the zero-derivative of a matrix is the same matrix.
! Due to the form in which matrices DxA and DxB where constructed,
! the last entry represents the order of the space derivative.
!--------------------------------------------------------------------
if (l - k  == 0) then
   DA = A ( :, :, i, j)
else
   DA = DxA ( :, :, i, j , l - k)
end if
!--------------------------------------------------------------------
! "FUN_COMBINATORY(n,k)" is a function which returns the combinatorial
! function of n over k.
!--------------------------------------------------------------------
comB = FUN_COMBINATORY ( l - 2, l - k - 1)
comA = FUN_COMBINATORY ( l - 1, l - k)
!
D = comB * DB - comA * DA

! ----------------------------------------------------------------------

END SUBROUTINE MATRIX_D ! ------------------------------------------------------------------------

! ------------------------------------------------------------------------ SUBROUTINE MATRIX_C ( A, B, DxA, DxB, CM) ! ---------------------------------------------------------------------- ! This subroutine computes the matrix "C" for the recursive simplified ! Cauchy-Kowalewskaya procedure in Proposition 3.2. The inputs of this ! matrix are: ! "DxA", "DxB", the spatial derivatives of matrix A and B. These are ! derivatives from first order to Accuracy-th order. ! ! This subroutine has been implemented with the following global variables: ! "Accuracy" is the accuracy of the scheme, which coincides with the ! number of quadrature points in space. ! "nGP" the number of quadrature points for time, here nGP = Accuracy-1. ! "NVAR" number of unknowns. ! ---------------------------------------------------------------------- DOUBLE PRECISION DxA ( NVAR, NVAR, Accuracy, nGP, Accuracy-1), & & DxB ( NVAR, NVAR, Accuracy, nGP, Accuracy-1), & & A ( NVAR, NVAR, Accuracy, nGP), & & B ( NVAR, NVAR, Accuracy, nGP), & & CM ( NVAR, NVAR, Accuracy, nGP, Accuracy, Accuracy) ! ---------------------------------------------------------------------- ! Variables passing this point are local. ! ---------------------------------------------------------------------- INTEGER i, j, k1, l1, m, NTIME, NSPACE

DOUBLE PRECISION DtCM ( NVAR, NVAR, Accuracy, nGP, nGP-1, Accuracy-1, &
     & Accuracy), Dmat ( NVAR, NVAR)
! ----------------------------------------------------------------------

! ----------------------------------------------------------------------
NTime  = nGP         ! Rename the number of quadrature points in time.
NSpace = Accuracy    ! Rename the number of quadrature points in space.

! ----------------------------------------------------------------------
! Initialize the coefficient-matrix at zero.
! ----------------------------------------------------------------------
CM = 0.0

!--------------------------------------------------------------------
! Compute the matrices contributing to the Cauchy-Kowalewskaya functional
! of the simplified approach. These are the matrices which are the
! matrix coefficients of spatial derivatives. See formula (29).
!--------------------------------------------------------------------
! Set "C(1,1) = -A". See Proposition 3.2.
!--------------------------------------------------------------------
do j = 1, NTime
   do i = 1, NSpace

      CM ( :, :, i, j, 1, 1) = - A ( :, :, i, j)

   end do
end do
!----------------------------------------------------------------
! Initialize the recursive step to generate matrix coefficient C.
!----------------------------------------------------------------
do k1 = 2, Accuracy - 1
   !
   do j = 1, NTime
      do i = 1, NSpace
         !-------------------------------------------------------------
         ! The "MATRIX_D" implements the formula (21), providing "Dmat".
         !-------------------------------------------------------------
         call MATRIX_D ( k1, k1, i, j, A, B, DxA, DxB, Dmat)
         !
         CM (  :, :, i, j, k1, k1) =  matmul (                   &
              & CM ( :, :, i, j, k1-1, k1-1), Dmat)
      end do
   end do
   !-------------------------------------------------------------
   do l1 = 1,  k1 - 1
      !
      !----------------------------------------------------
      ! Compute the time derivatives from values at space
      ! time-nodes.
      !----------------------------------------------------
      CALL MATRIX_TIME_GRADIENT (                          &
           & CM ( :, :, :, :, k1-1, l1),                   &
           & DtCM ( :, :, :, :, :, k1-1, l1)  )
      !----------------------------------------------------
      ! Implement the formula (29) of the paper.
      !----------------------------------------------------
      do j = 1, NTime
         do i = 1, NSpace
            !-------------------------------------------------------
            !-------------------------------------------------------
            !
            CM (  :, :, i, j, k1, l1) = DtCM( :, :, i, j, 1, k1-1, l1)
            !
            do m = l1-1, k1-1

               if( m > 0) then
                  call MATRIX_D ( m+1, l1, i, j, A, B, DxA, DxB, Dmat)

                  CM (  :, :, i, j, k1, l1)  = CM (  :, :, i, j, k1, l1) + &
                       &   matmul (                    &
                       & CM ( :, :, i, j, k1-1, m), Dmat)
               end if
            end do
            !
            !----------------------------------------------------
         end do
      end do
      !-------------------------------------------------------------
   end do !  End loop: "do l1 = 1, k1 - 1"
   !----------------------------------------------------------------
end do ! End loop: "do k1 = 2, Accuracy - 1".

! ----------------------------------------------------------------------

END SUBROUTINE MATRIX_C ! ------------------------------------------------------------------------

B.2 Approximation of space and time derivatives

The interpolation polynomials obtained here can be straightforward generalized to scalar, vector and matrix functions. So, we provide the polynomials for a generic function $f(\xi)$. To carry out interpolations, we have used the following strategy.

For interpolation in space, we first cast an interval $[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$ into $[-\frac{1}{2},\frac{1}{2}]$ by the change of variable $x=x_{i-\frac{1}{2}}+(\xi-\frac{1}{2})\Delta x$. Second, we consider $\xi_{j}=-\frac{1}{2}+\frac{(j-1)}{M}$ with $j=1,...,M+1$. So, the polynomial interpolation in space for $M+1$ order has the form $P(x)=\sum_{k=0}^{M}a_{M+1,k}\cdot\xi^{k}$. The coefficients $a_{M+1,k}$ are given in the table 6 for $M=1,2,3,4$, we have used the convention $f_{j}=f(\xi_{j})$. Because of the change of variable, each space derivative of order $l$ needs to be scaled by $\Delta x^{-l}$.

To obtain interpolation in time, we cast an interval $[t^{n},t^{n+1}]$ into the reference element $[0,1]$ by the change of variable $t=t^{n}+\tau\Delta t$. Then we construct an interpolation polynomial from the Gaussian points, $\tau_{j}$ in $[0,1]$ with $j=1,...,n_{GP}$. Interpolation in time are only needed for accuracy higher than 3. In such a case. The interpolation polynomial has the form $T(\tau)=\sum_{k=0}^{M-1}b_{M+1,k}\cdot\tau^{k}$. The coefficients $b_{M+1,k}$ are shown in the table 7 for $M=2,3,4$. We have used the convention $f_{j}=f(\tau_{j})$. Because of the change of variable, each time derivative of order $l$ needs to be scaled by $\Delta t^{-l}$.

B.3 Evaluation of the numerical flux and source term

Regarding the evaluation of (10). The numerical flux can be easily evaluated as $\mathbf{F}_{i+\frac{1}{2}}=\sum_{j=1}^{n_{GP}}\omega_{j}\mathbf{F}_{h}(\mathbf{Q}_{i}(\xi_{M+1},\tau_{j}),\mathbf{Q}_{i+1}(\xi_{1},\tau_{j}))$, where $\omega_{j}$ corresponds to the $j$th Gaussian weight in $[0,1]$. Whereas, the source terms, can be easily obtained from the interpolation in space and quadrature points in time. Table 8 shows the form in which the source is computed for several orders of accuracy.