and differentiate with respect to time to obtain

(5.42) $\displaystyle\frac{d^{2}\widehat{\boldsymbol{v}}}{dt^{2}}$
$\displaystyle=-\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}})\frac{d\tilde{\boldsymbol{w}}}{dt}$
(5.43) $\displaystyle\frac{d^{2}\widehat{\boldsymbol{e}}}{dt^{2}}$
$\displaystyle=\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\mathbf{\Phi}_{e}^{T}\mathcal{P}_{\mathsf{F}^{tv}}\boldsymbol{J}_{w}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}})\frac{d\tilde{\boldsymbol{w}}}{dt}$
(5.44) $\displaystyle\frac{d^{2}\widehat{\boldsymbol{x}}}{dt^{2}}$
$\displaystyle=\mathbf{\Phi}_{x}^{T}\frac{d\tilde{\boldsymbol{v}}}{dt}.$

Substituting back to (5.37), we observe that

(5.45) $\displaystyle\widehat{\boldsymbol{v}}(t_{n+1})$
$\displaystyle=\widehat{\boldsymbol{v}}(t_{n})-\Delta
t_{n}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\left(\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}(t_{n}))+\dfrac{\Delta
t_{n}}{2}\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}(t_{n}))\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}(t_{n}))\right)+(\Delta
t_{n})^{3}\boldsymbol{\mathcal{R}}_{v}^{(1)}$
$\displaystyle\widehat{\boldsymbol{e}}(t_{n+1})$
$\displaystyle=\widehat{\boldsymbol{e}}(t_{n})+\Delta
t_{n}\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\mathbf{\Phi}_{e}^{T}\mathcal{P}_{\mathsf{F}^{tv}}\left(\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}(t_{n}))+\dfrac{\Delta
t_{n}}{2}\boldsymbol{J}_{w}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}(t_{n}))\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}(t_{n}))\right)+(\Delta
t_{n})^{3}\boldsymbol{\mathcal{R}}_{e}^{(1)}$
$\displaystyle\widehat{\boldsymbol{x}}(t_{n+1})$
$\displaystyle=\widehat{\boldsymbol{x}}(t_{n})+\Delta
t_{n}\mathbf{\Phi}_{x}^{T}\left(\tilde{\boldsymbol{v}}(t_{n})-\dfrac{\Delta
t_{n}}{2}\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}(t_{n}))\right)+(\Delta
t_{n})^{3}\boldsymbol{\mathcal{R}}_{x}^{(1)},$

where
$\frac{d\tilde{\boldsymbol{w}}}{dt}(t_{n})=\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}(t_{n}))$
is used through Eq. (5.34). Next, we are going to find a similar relation for
the time-discrete reduced-order coefficients. First, we rewrite the first
Runge-Kutta stage of (3.22) to obtain

(5.46) $\displaystyle\widehat{\boldsymbol{v}}_{n+\frac{1}{2}}$
$\displaystyle=\widehat{\boldsymbol{v}}_{n}-\frac{\Delta
t_{n}}{2}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})$
(5.47) $\displaystyle\widehat{\boldsymbol{e}}_{n+\frac{1}{2}}$
$\displaystyle=\widehat{\boldsymbol{e}}_{n}+\frac{\Delta
t_{n}}{2}\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\mathbf{\Phi}_{e}^{T}\mathcal{P}_{\mathsf{F}^{tv}}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n})+(\Delta
t_{n})^{2}\boldsymbol{\mathcal{R}}_{e}^{(2)}$ (5.48)
$\displaystyle\widehat{\boldsymbol{x}}_{n+\frac{1}{2}}$
$\displaystyle=\widehat{\boldsymbol{x}}_{n}+\frac{\Delta
t_{n}}{2}\mathbf{\Phi}_{x}^{T}\tilde{\boldsymbol{v}}_{n}+(\Delta
t_{n})^{2}\boldsymbol{\mathcal{R}}_{x}^{(2)},$

where $\mathsf{F}^{\mathfrak{1}}$ and $\mathsf{F}^{tv}$ are defined in Eq.
(2.6). The first Runge-Kutta stage correctors are given by

(5.49) $\displaystyle\boldsymbol{\mathcal{R}}_{e}^{(2)}$
$\displaystyle=-\dfrac{1}{4}\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\mathbf{\Phi}_{e}^{T}\mathcal{P}_{\mathsf{F}^{tv}}\left(\boldsymbol{F}(\tilde{\boldsymbol{w}}_{n})\right)^{T}\cdot\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})$
$\displaystyle\boldsymbol{\mathcal{R}}_{x}^{(2)}$
$\displaystyle=-\dfrac{1}{4}\mathbf{\Phi}_{x}^{T}\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n}).$

By denoting
$\boldsymbol{\mathcal{R}}_{w}^{(2)}=[\boldsymbol{0}_{N_{\mathcal{V}}};\mathbf{\Phi}_{e}\boldsymbol{\mathcal{R}}_{e}^{(2)};\mathbf{\Phi}_{x}\boldsymbol{\mathcal{R}}_{x}^{(2)}]^{T}$
and multiplying the basis matrices to (5.46), we obtain

(5.50)
$\tilde{\boldsymbol{w}}_{n+\frac{1}{2}}=\tilde{\boldsymbol{w}}_{n}+\frac{\Delta
t_{n}}{2}\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+(\Delta
t_{n})^{2}\boldsymbol{\mathcal{R}}_{w}^{(2)},$

where $\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})$ is defined in (5.34).
Multiplying the basis matrix $\mathbf{\Phi}_{v}$ to
$\widehat{\boldsymbol{v}}_{n+\frac{1}{2}}$ in (5.46) and
$\widehat{\boldsymbol{v}}_{n+1}$ in (3.22), we have

(5.51) $\displaystyle\tilde{\boldsymbol{v}}_{n+\frac{1}{2}}$
$\displaystyle=\tilde{\boldsymbol{v}}_{n}-\dfrac{\Delta
t_{n}}{2}\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})$
$\displaystyle\tilde{\boldsymbol{v}}_{n+1}$
$\displaystyle=\tilde{\boldsymbol{v}}_{n}-\Delta
t_{n}\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}}),$

By the definition of $\bar{\tilde{\boldsymbol{v}}}_{n+\frac{1}{2}}$ below
(3.22) and the second equation in (5.51), we observe that

(5.52) $\displaystyle\bar{\tilde{\boldsymbol{v}}}_{n+\frac{1}{2}}$
$\displaystyle=\tilde{\boldsymbol{v}}_{n}-\dfrac{\Delta
t_{n}}{2}\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}}).$

By subtracting the first equation in (5.51) from (5.52), we obtain

(5.53) $\displaystyle\bar{\tilde{\boldsymbol{v}}}_{n+\frac{1}{2}}$
$\displaystyle=\tilde{\boldsymbol{v}}_{n+\frac{1}{2}}+\dfrac{\Delta
t_{n}}{2}\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\left(\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})-\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})\right),$

which allows us to rewrite the second Runge-Kutta stage of (3.22) as

(5.54) $\displaystyle\widehat{\boldsymbol{v}}_{n+1}$
$\displaystyle=\widehat{\boldsymbol{v}}_{n}-\Delta
t_{n}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})$
(5.55) $\displaystyle\widehat{\boldsymbol{e}}_{n+1}$
$\displaystyle=\widehat{\boldsymbol{e}}_{n}+\Delta
t_{n}\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\mathbf{\Phi}_{e}^{T}\mathcal{P}_{\mathsf{F}^{tv}}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})+(\Delta
t_{n})^{3}\boldsymbol{\mathcal{R}}_{e}^{(3)}$ (5.56)
$\displaystyle\widehat{\boldsymbol{x}}_{n+1}$
$\displaystyle=\widehat{\boldsymbol{x}}_{n}+\Delta
t_{n}\mathbf{\Phi}_{x}^{T}\tilde{\boldsymbol{v}}_{n+\frac{1}{2}}+(\Delta
t_{n})^{3}\boldsymbol{\mathcal{R}}_{x}^{(3)},$

where the second Runge-Kutta correctors are given by

(5.57) $\displaystyle\boldsymbol{\mathcal{R}}_{e}^{(3)}$
$\displaystyle=\dfrac{1}{2\Delta
t_{n}}\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\mathbf{\Phi}_{e}^{T}\mathcal{P}_{\mathsf{F}^{tv}}\left(\boldsymbol{F}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})\right)^{T}\cdot\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\left(\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})-\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})\right)$
$\displaystyle\boldsymbol{\mathcal{R}}_{x}^{(3)}$
$\displaystyle=\dfrac{1}{2\Delta
t_{n}}\mathbf{\Phi}_{x}^{T}\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\left(\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})-\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})\right).$

Using Taylor’s remainder theorem again, we have

(5.58)
$\displaystyle\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})$
$\displaystyle=\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})+\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\left(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}}-\tilde{\boldsymbol{w}}_{n}\right)+(\Delta
t_{n})^{2}\boldsymbol{\mathcal{R}}_{\mathsf{F}^{\mathfrak{1}}}^{(4)}$
$\displaystyle\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})$
$\displaystyle=\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n})+\boldsymbol{J}_{w}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n})\left(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}}-\tilde{\boldsymbol{w}}_{n}\right)+(\Delta
t_{n})^{2}\boldsymbol{\mathcal{R}}_{\mathsf{F}^{tv}}^{(4)},$

where the remainders are given by

(5.59)
$\displaystyle\boldsymbol{\mathcal{R}}_{\mathsf{F}^{\mathfrak{1}}}^{(4)}$
$\displaystyle=\dfrac{2}{(\Delta
t_{n})^{2}}\sum_{|\boldsymbol{I}|=2}\dfrac{(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}}-\tilde{\boldsymbol{w}}_{n})^{\boldsymbol{I}}}{\boldsymbol{I}!}\int_{0}^{1}(1-s)\partial^{\boldsymbol{I}}\mathsf{F}^{\mathfrak{1}}\left((1-s)\tilde{\boldsymbol{w}}_{n}+s\tilde{\boldsymbol{w}}_{n+\frac{1}{2}}\right)ds$
$\displaystyle\boldsymbol{\mathcal{R}}_{\mathsf{F}^{tv}}^{(4)}$
$\displaystyle=\dfrac{2}{(\Delta
t_{n})^{2}}\sum_{|\boldsymbol{I}|=2}\dfrac{(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}}-\tilde{\boldsymbol{w}}_{n})^{\boldsymbol{I}}}{\boldsymbol{I}!}\int_{0}^{1}(1-s)\partial^{\boldsymbol{I}}\mathsf{F}^{tv}\left((1-s)\tilde{\boldsymbol{w}}_{n}+s\tilde{\boldsymbol{w}}_{n+\frac{1}{2}}\right)ds.$

Here $\boldsymbol{I}=(I_{1},I_{2},\ldots,I_{N})$ is a multi-index with order
$|\boldsymbol{I}|=I_{1}+I_{2}+\ldots+I_{N}$, and the notations
$\boldsymbol{w}^{\boldsymbol{I}}$, $\boldsymbol{I}!$ and
$\partial^{\boldsymbol{I}}$ are formally defined as

(5.60) $\displaystyle\boldsymbol{w}^{\boldsymbol{I}}$
$\displaystyle=w_{1}^{I_{1}}w_{2}^{I_{2}}\ldots w_{N}^{I_{N}}$
$\displaystyle\boldsymbol{I}!$ $\displaystyle=I_{1}!I_{2}!\ldots I_{N}!$
$\displaystyle\partial^{\boldsymbol{I}}$
$\displaystyle=\dfrac{\partial^{|\boldsymbol{I}|}}{\partial
w_{1}^{I_{1}}\partial w_{2}^{I_{2}}\cdots\partial w_{N}^{I_{N}}}.$

Substituting (5.50) and (5.58) into (5.54), we have

(5.61) $\displaystyle\widehat{\boldsymbol{v}}_{n+1}$
$\displaystyle=\widehat{\boldsymbol{v}}_{n}-\Delta
t_{n}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\left[\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})+\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\left(\frac{\Delta
t_{n}}{2}\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+(\Delta
t_{n})^{2}\boldsymbol{\mathcal{R}}_{w}^{(2)}\right)+(\Delta
t_{n})^{2}\boldsymbol{\mathcal{R}}_{\mathsf{F}^{\mathfrak{1}}}^{(4)}\right]$
(5.62) $\displaystyle\widehat{\boldsymbol{e}}_{n+1}$
$\displaystyle=\widehat{\boldsymbol{e}}_{n}+\Delta
t_{n}\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\mathbf{\Phi}_{e}^{T}\mathcal{P}_{\mathsf{F}^{tv}}\left[\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n})+\boldsymbol{J}_{w}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n})\left(\frac{\Delta
t_{n}}{2}\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+(\Delta
t_{n})^{2}\boldsymbol{\mathcal{R}}_{w}^{(2)}\right)+(\Delta
t_{n})^{2}\boldsymbol{\mathcal{R}}_{\mathsf{F}^{tv}}^{(4)}\right]+(\Delta
t_{n})^{3}\boldsymbol{\mathcal{R}}_{e}^{(3)}$ (5.63)
$\displaystyle\widehat{\boldsymbol{x}}_{n+1}$
$\displaystyle=\widehat{\boldsymbol{x}}_{n}+\Delta
t_{n}\mathbf{\Phi}_{x}^{T}\left(\tilde{\boldsymbol{v}}_{n}-\frac{\Delta
t_{n}}{2}\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\right)+(\Delta
t_{n})^{3}\boldsymbol{\mathcal{R}}_{x}^{(3)}.$

Combining (5.45), (5.61) and (5.36) yields the error identities

(5.64) $\displaystyle\boldsymbol{\delta}_{v,n+1}$
$\displaystyle=\boldsymbol{\delta}_{v,n}+\Delta
t_{n}\mathbf{a}^{(1)}_{v,n}+(\Delta t_{n})^{2}\mathbf{a}^{(2)}_{v,n}+(\Delta
t_{n})^{3}\mathbf{a}^{(3)}_{v,n}$ $\displaystyle\boldsymbol{\delta}_{e,n+1}$
$\displaystyle=\boldsymbol{\delta}_{e,n}+\Delta
t_{n}\mathbf{a}^{(1)}_{e,n}+(\Delta t_{n})^{2}\mathbf{a}^{(2)}_{e,n}+(\Delta
t_{n})^{3}\mathbf{a}^{(3)}_{e,n}$ $\displaystyle\boldsymbol{\delta}_{x,n+1}$
$\displaystyle=\boldsymbol{\delta}_{x,n}+\Delta
t_{n}\mathbf{a}^{(1)}_{x,n}+(\Delta t_{n})^{2}\mathbf{a}^{(2)}_{x,n}+(\Delta
t_{n})^{3}\mathbf{a}^{(3)}_{x,n},$

where the coefficient vectors are defined as

(5.65) $\displaystyle\mathbf{a}^{(1)}_{v,n}$
$\displaystyle=-\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\left(\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}(t_{n}))-\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\right)$
$\displaystyle\mathbf{a}^{(1)}_{e,n}$
$\displaystyle=\mathbf{\Phi}_{e}\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\mathbf{\Phi}_{e}^{T}\mathcal{P}_{\mathsf{F}^{tv}}\left(\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}(t_{n}))-\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n})\right)$
$\displaystyle\mathbf{a}^{(1)}_{x,n}$
$\displaystyle=\mathbf{\Phi}_{x}\mathbf{\Phi}_{x}^{T}\left(\tilde{\boldsymbol{v}}(t_{n})-\tilde{\boldsymbol{v}}_{n}\right)$
$\displaystyle\mathbf{a}^{(2)}_{v,n}$
$\displaystyle=-\frac{1}{2}\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\left(\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}(t_{n}))\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}(t_{n}))-\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})\right)$
$\displaystyle\mathbf{a}^{(2)}_{e,n}$
$\displaystyle=\frac{1}{2}\mathbf{\Phi}_{e}\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\mathbf{\Phi}_{e}^{T}\mathcal{P}_{\mathsf{F}^{tv}}\left(\boldsymbol{J}_{w}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}(t_{n}))\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}(t_{n}))-\boldsymbol{J}_{w}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n})\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})\right)$
$\displaystyle\mathbf{a}^{(2)}_{x,n}$
$\displaystyle=-\frac{1}{2}\mathbf{\Phi}_{x}\mathbf{\Phi}_{x}^{T}\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\left(\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}(t_{n}))-\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\right)$
$\displaystyle\mathbf{a}^{(3)}_{v,n}$
$\displaystyle=\mathbf{\Phi}_{v}(\boldsymbol{\mathcal{R}}_{v}^{(1)}+\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}(\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\boldsymbol{\mathcal{R}}_{w}^{(2)}+\boldsymbol{\mathcal{R}}_{\mathsf{F}^{\mathfrak{1}}}^{(4)}))$
$\displaystyle\mathbf{a}^{(3)}_{e,n}$
$\displaystyle=\mathbf{\Phi}_{e}(\boldsymbol{\mathcal{R}}_{e}^{(1)}-\boldsymbol{\mathcal{R}}_{e}^{(3)}-\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\mathbf{\Phi}_{e}^{T}\mathcal{P}_{\mathsf{F}^{tv}}(\boldsymbol{J}_{w}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n})\boldsymbol{\mathcal{R}}_{w}^{(2)}+\boldsymbol{\mathcal{R}}_{\mathsf{F}^{tv}}^{(4)}))$
$\displaystyle\mathbf{a}^{(3)}_{x,n}$
$\displaystyle=\mathbf{\Phi}_{x}(\boldsymbol{\mathcal{R}}_{x}^{(1)}-\boldsymbol{\mathcal{R}}_{x}^{(3)}).$

We remark that all the residual vectors only appear in
$(\mathbf{a}^{(3)}_{v,n},\mathbf{a}^{(3)}_{e,n},\mathbf{a}^{(3)}_{x,n})$ and
therefore scale with $(\Delta t_{n})^{3}$. We estimate the first order terms
using a similar argument as (5.23) in Theorem 1, i.e.

(5.66)
$\displaystyle\left\lVert\boldsymbol{L}_{\mathcal{V}}^{-1}\mathbf{a}^{(1)}_{v,n}\right\rVert_{2}$
$\displaystyle\leq\left\lVert\boldsymbol{M}_{\mathcal{V}}^{-1/2}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}(\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}(t_{n}))-\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n}))\right\rVert_{2}$
$\displaystyle\left\lVert\boldsymbol{L}_{\mathcal{E}}^{-1}\mathbf{a}^{(1)}_{e,n}\right\rVert_{2}$
$\displaystyle\leq\left\lVert\boldsymbol{M}_{\mathcal{E}}^{-1/2}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}(\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}(t_{n}))-\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n}))\right\rVert_{2}$
$\displaystyle\left\lVert\mathbf{a}^{(1)}_{x,n}\right\rVert_{2}$
$\displaystyle\leq\left\lVert\boldsymbol{M}_{\mathcal{V}}^{-1/2}\right\rVert_{2}\left\lVert\tilde{\boldsymbol{v}}(t_{n})-\tilde{\boldsymbol{v}}_{n}\right\rVert_{\boldsymbol{M}_{\mathcal{V}}}.$

With the assumption (5.6), we have

(5.67)
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\mathbf{a}^{(1)}_{v,n},\mathbf{a}^{(1)}_{e,n},\mathbf{a}^{(1)}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq
C_{1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,n},\boldsymbol{\delta}_{e,n},\boldsymbol{\delta}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}.$

where $C_{1}$ is the constant defined in (5.25). Similarly, with the
assumptions (5.33) and (5.6), we estimate the second order terms by

(5.68)
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\mathbf{a}^{(2)}_{v,n},\mathbf{a}^{(2)}_{e,n},\mathbf{a}^{(2)}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq
C_{3}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,n},\boldsymbol{\delta}_{e,n},\boldsymbol{\delta}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|},$

with
$C_{3}=\frac{1}{2}\max\left\\{\left\lVert\boldsymbol{M}_{\mathcal{V}}^{-1/2}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}K_{3},\left\lVert\boldsymbol{M}_{\mathcal{E}}^{-1/2}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}K_{4},\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}K_{1}\right\\}$.
Using the error identities (5.64) and applying triangle inequality, we have

(5.69)
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,n+1},\boldsymbol{\delta}_{e,n+1},\boldsymbol{\delta}_{x,n+1})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq(1+C_{1}\Delta
t_{n}+C_{3}(\Delta
t_{n})^{2}){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,n},\boldsymbol{\delta}_{e,n},\boldsymbol{\delta}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+(\Delta
t_{n})^{3}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\mathbf{a}^{(3)}_{v,n},\mathbf{a}^{(3)}_{e,n},\mathbf{a}^{(3)}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}.$

Next, we find a bound for the third order terms. By direct computation, we
observe that

(5.70)
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\mathbf{a}^{(3)}_{v,n},\mathbf{a}^{(3)}_{e,n},\mathbf{a}^{(3)}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
$\displaystyle\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\mathbf{\Phi}_{v}\boldsymbol{\mathcal{R}}_{v}^{(1)},\mathbf{\Phi}_{e}\boldsymbol{\mathcal{R}}_{e}^{(1)},\mathbf{\Phi}_{x}\boldsymbol{\mathcal{R}}_{x}^{(1)})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+$
$\displaystyle\quad\quad
C_{4}C_{5}\left\lVert\left(\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\boldsymbol{\mathcal{R}}_{w}^{(2)},\boldsymbol{J}_{w}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n})\boldsymbol{\mathcal{R}}_{w}^{(2)},\boldsymbol{0}_{N_{\mathcal{V}}}\right)\right\rVert_{2}+$
$\displaystyle\quad\quad
C_{4}\left\lVert(\boldsymbol{0}_{N_{\mathcal{V}}},\mathbf{\Phi}_{e}\boldsymbol{\mathcal{R}}_{e}^{(3)},\mathbf{\Phi}_{x}\boldsymbol{\mathcal{R}}_{x}^{(3)})\right\rVert_{2}+C_{4}C_{5}\left\lVert\left(\boldsymbol{\mathcal{R}}_{\mathsf{F}^{\mathfrak{1}}}^{(4)},\boldsymbol{\mathcal{R}}_{\mathsf{F}^{tv}}^{(4)},\boldsymbol{0}_{N_{\mathcal{V}}}\right)\right\rVert_{2},$

where
$C_{4}=\max\left\\{\left\lVert\boldsymbol{M}_{\mathcal{V}}^{1/2}\right\rVert_{2},\left\lVert\boldsymbol{M}_{\mathcal{E}}^{1/2}\right\rVert_{2},1\right\\}$
and
$C_{5}=\max\left\\{\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2},\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}\right\\}$.
We will estimate each of terms on the right hand side of (5.70). We start with
the last term on the right hand side of (5.70), where we invoke (5.50) to
estimate (5.59) by

(5.71)
$\displaystyle\left\lVert\boldsymbol{\mathcal{R}}_{\mathsf{F}^{\mathfrak{1}}}^{(4)}\right\rVert_{2}$
$\displaystyle\leq\dfrac{1}{8}NM^{(4)}_{\mathsf{F}^{\mathfrak{1}}}\left\lVert\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+2\Delta
t_{n}\boldsymbol{\mathcal{R}}_{w}^{(2)}\right\rVert_{2}^{2}$
$\displaystyle\left\lVert\boldsymbol{\mathcal{R}}_{\mathsf{F}^{tv}}^{(4)}\right\rVert_{2}$
$\displaystyle\leq\dfrac{1}{8}NM^{(4)}_{\mathsf{F}^{tv}}\left\lVert\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+2\Delta
t_{n}\boldsymbol{\mathcal{R}}_{w}^{(2)}\right\rVert_{2}^{2},$

where the constants $M^{(2)}_{\mathsf{F}^{\mathfrak{1}}}$ and
$M^{(2)}_{\mathsf{F}^{tv}}$ are defined as

(5.72) $\begin{split}M^{(2)}_{\mathsf{F}^{\mathfrak{1}}}&=\max_{0\leq
j\leq\widetilde{N}_{t}-1}\max_{|\boldsymbol{I}|=2}\max_{0\leq s\leq
1}\left\lVert\partial^{\boldsymbol{I}}\mathsf{F}^{\mathfrak{1}}\left((1-s)\tilde{\boldsymbol{w}}_{j}+s\tilde{\boldsymbol{w}}_{j+\frac{1}{2}}\right)\right\rVert_{2}\\\
M^{(2)}_{\mathsf{F}^{tv}}&=\max_{0\leq
j\leq\widetilde{N}_{t}-1}\max_{|\boldsymbol{I}|=2}\max_{0\leq s\leq
1}\left\lVert\partial^{\boldsymbol{I}}\mathsf{F}^{tv}\left((1-s)\tilde{\boldsymbol{w}}_{j}+s\tilde{\boldsymbol{w}}_{j+\frac{1}{2}}\right)\right\rVert_{2}.\end{split}$

Therefore, the last term on the right hand side of (5.70) can be estimated by

(5.73)
$\left\lVert\left(\boldsymbol{\mathcal{R}}_{\mathsf{F}^{\mathfrak{1}}}^{(4)},\boldsymbol{\mathcal{R}}_{\mathsf{F}^{tv}}^{(4)},\boldsymbol{0}_{N_{\mathcal{V}}}\right)\right\rVert_{2}\leq\dfrac{1}{8}N(M_{\mathsf{F}^{\mathfrak{1}}}^{(2)}+M_{\mathsf{F}^{tv}}^{(2)})\left\lVert\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+2\Delta
t_{n}\boldsymbol{\mathcal{R}}_{w}^{(2)}\right\rVert_{2}^{2}.$

For the third term on the right hand side of (5.70), substituting (5.58) and
(5.50) into (5.57) gives

(5.74) $\displaystyle\boldsymbol{\mathcal{R}}_{e}^{(3)}$
$\displaystyle=-\dfrac{1}{4}\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\mathbf{\Phi}_{e}^{T}\mathcal{P}_{\mathsf{F}^{tv}}\left(\boldsymbol{F}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})\right)^{T}\cdot\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\left(\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+2\Delta
t_{n}(\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\boldsymbol{\mathcal{R}}_{w}^{(2)}+\boldsymbol{\mathcal{R}}_{\mathsf{F}^{\mathfrak{1}}}^{(4)})\right)$
$\displaystyle\boldsymbol{\mathcal{R}}_{x}^{(3)}$
$\displaystyle=-\dfrac{1}{4}\mathbf{\Phi}_{x}^{T}\mathbf{\Phi}_{v}\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\mathbf{\Phi}_{v}^{T}\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\left(\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+2\Delta
t_{n}(\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\boldsymbol{\mathcal{R}}_{w}^{(2)}+\boldsymbol{\mathcal{R}}_{\mathsf{F}^{\mathfrak{1}}}^{(4)})\right).$

(5.74) and (5.71) together imply

(5.75)
$\begin{split}\left\lVert(\boldsymbol{0}_{N_{\mathcal{V}}},\mathbf{\Phi}_{e}\boldsymbol{\mathcal{R}}_{e}^{(3)},\mathbf{\Phi}_{x}\boldsymbol{\mathcal{R}}_{x}^{(3)})\right\rVert_{2}&\leq\dfrac{1}{4}\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}\left\lVert\boldsymbol{F}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})\right\rVert_{2}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}\\\
&\quad\quad\left[\left\lVert\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\right\rVert_{2}\left\lVert\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+2\Delta
t_{n}\boldsymbol{\mathcal{R}}_{w}^{(2)}\right\rVert_{2}+\dfrac{1}{4}NM^{(2)}_{\mathsf{F}^{\mathfrak{1}}}\Delta
t_{n}\left\lVert\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+2\Delta
t_{n}\boldsymbol{\mathcal{R}}_{w}^{(2)}\right\rVert_{2}^{2}\right].\end{split}$

For the second term on the right hand side of (5.70), we have

(5.76)
$\displaystyle\left\lVert\left(\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\boldsymbol{\mathcal{R}}_{w}^{(2)},\boldsymbol{J}_{w}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n})\boldsymbol{\mathcal{R}}_{w}^{(2)},\boldsymbol{0}_{N_{\mathcal{V}}}\right)\right\rVert_{2}$
$\displaystyle\leq\left(\left\lVert\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{n})\right\rVert_{2}^{2}+\left\lVert\boldsymbol{J}_{w}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{n})\right\rVert_{2}^{2}\right)^{\frac{1}{2}}\left\lVert\boldsymbol{\mathcal{R}}_{w}^{(2)}\right\rVert_{2}.$

Using (5.73), (5.75) and (5.76), we can rewrite (5.70) as

(5.77)
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\mathbf{a}^{(3)}_{v,n},\mathbf{a}^{(3)}_{e,n},\mathbf{a}^{(3)}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
$\displaystyle\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\mathbf{\Phi}_{v}\boldsymbol{\mathcal{R}}_{v}^{(1)},\mathbf{\Phi}_{e}\boldsymbol{\mathcal{R}}_{e}^{(1)},\mathbf{\Phi}_{x}\boldsymbol{\mathcal{R}}_{x}^{(1)})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+C_{4}C_{5}\left(M^{(1)}_{\mathsf{F}^{\mathfrak{1}}}+M^{(1)}_{\mathsf{F}^{tv}}\right)\left\lVert\boldsymbol{\mathcal{R}}_{w}^{(2)}\right\rVert_{2}+$
$\displaystyle\quad\quad\dfrac{C_{4}}{4}\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}\left\lVert\boldsymbol{F}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})\right\rVert_{2}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}M^{(1)}_{\mathsf{F}^{\mathfrak{1}}}\left\lVert\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+2\Delta
t_{n}\boldsymbol{\mathcal{R}}_{w}^{(2)}\right\rVert_{2}+$
$\displaystyle\quad\quad\dfrac{C_{4}}{16}\left[2C_{5}N(M_{\mathsf{F}^{\mathfrak{1}}}^{(2)}+M_{\mathsf{F}^{tv}}^{(2)})+\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}\left\lVert\boldsymbol{F}(\tilde{\boldsymbol{w}}_{n+\frac{1}{2}})\right\rVert_{2}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}NM_{\mathsf{F}^{\mathfrak{1}}}^{(2)}\Delta
t_{n}\right]$
$\displaystyle\quad\quad\left\lVert\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+2\Delta
t_{n}\boldsymbol{\mathcal{R}}_{w}^{(2)}\right\rVert_{2}^{2},$

where the constants $M^{(1)}_{\mathsf{F}^{\mathfrak{1}}}$ and
$M^{(1)}_{\mathsf{F}^{tv}}$ are defined as

(5.78) $\begin{split}M^{(1)}_{\mathsf{F}^{\mathfrak{1}}}&=\max_{0\leq
j\leq\widetilde{N}_{t}-1}\left\lVert\boldsymbol{J}_{w}\mathsf{F}^{\mathfrak{1}}(\tilde{\boldsymbol{w}}_{j})\right\rVert_{2}\\\
M^{(1)}_{\mathsf{F}^{tv}}&=\max_{0\leq
j\leq\widetilde{N}_{t}-1}\left\lVert\boldsymbol{J}_{w}\mathsf{F}^{tv}(\tilde{\boldsymbol{w}}_{j})\right\rVert_{2}.\end{split}$

By the definition of $\boldsymbol{\mathcal{R}}_{w}^{(2)}$ before (5.50), we
obtain

(5.79)
$\displaystyle\left\lVert\boldsymbol{\mathcal{R}}_{w}^{(2)}\right\rVert_{2}$
$\displaystyle\leq\dfrac{1}{4}\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}\left\lVert\boldsymbol{F}(\tilde{\boldsymbol{w}}_{n})\right\rVert_{2}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}\left\lVert\boldsymbol{F}(\tilde{\boldsymbol{w}}_{n})\right\rVert_{2}\sqrt{N_{\mathcal{V}}},$

which further implies

(5.80)
$\displaystyle\left\lVert\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})+2\Delta
t_{n}\boldsymbol{\mathcal{R}}_{w}^{(2)}\right\rVert_{2}$
$\displaystyle\leq\left\lVert\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{n})\right\rVert_{2}+\dfrac{1}{2}\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}\left\lVert\boldsymbol{F}(\tilde{\boldsymbol{w}}_{n})\right\rVert_{2}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}\left\lVert\boldsymbol{F}(\tilde{\boldsymbol{w}}_{n})\right\rVert_{2}\sqrt{N_{\mathcal{V}}}\Delta
t_{n}.$

There, (5.77) can be further rewritten as

(5.81)
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\mathbf{a}^{(3)}_{v,n},\mathbf{a}^{(3)}_{e,n},\mathbf{a}^{(3)}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
$\displaystyle\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\mathbf{\Phi}_{v}\boldsymbol{\mathcal{R}}_{v}^{(1)},\mathbf{\Phi}_{e}\boldsymbol{\mathcal{R}}_{e}^{(1)},\mathbf{\Phi}_{x}\boldsymbol{\mathcal{R}}_{x}^{(1)})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+$
$\displaystyle\quad\quad\dfrac{C_{4}C_{5}}{4}\left(M^{(1)}_{\mathsf{F}^{\mathfrak{1}}}+M^{(1)}_{\mathsf{F}^{tv}}\right)\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}\sqrt{N_{\mathcal{V}}}+$
$\displaystyle\quad\quad\dfrac{C_{4}}{8}\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}M^{(1)}_{\mathsf{F}^{\mathfrak{1}}}$
$\displaystyle\quad\quad\left[2\tilde{M}^{(0)}_{\tilde{\mathsf{F}}}+\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}\sqrt{N_{\mathcal{V}}}\Delta
t_{n}\right]+$
$\displaystyle\quad\quad\dfrac{C_{4}}{64}\left[2C_{5}N(M_{\mathsf{F}^{\mathfrak{1}}}^{(2)}+M_{\mathsf{F}^{tv}}^{(2)})+\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}NM_{\mathsf{F}^{\mathfrak{1}}}^{(2)}\Delta
t_{n}\right]$
$\displaystyle\quad\quad\left[2\tilde{M}^{(0)}_{\tilde{\mathsf{F}}}+\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}\sqrt{N_{\mathcal{V}}}\Delta
t_{n}\right]^{2},$

where the constants $M^{(0)}_{\boldsymbol{F}}$ and
$\tilde{M}^{(0)}_{\mathsf{F}}$ are defined as

(5.82) $\begin{split}M^{(0)}_{\boldsymbol{F}}&=\max_{0\leq
j\leq\widetilde{N}_{t}-1}\left\\{\left\lVert\boldsymbol{F}(\tilde{\boldsymbol{w}}_{j})\right\rVert_{2},\left\lVert\boldsymbol{F}(\tilde{\boldsymbol{w}}_{j+\frac{1}{2}})\right\rVert_{2}\right\\}\\\
\tilde{M}^{(0)}_{\tilde{\mathsf{F}}}&=\max_{0\leq
j\leq\widetilde{N}_{t}-1}\left\lVert\tilde{\mathsf{F}}(\tilde{\boldsymbol{w}}_{j})\right\rVert_{2}.\end{split}$

Finally, for the first term on the right hand side of (5.70), from (5.38), we
see that

(5.83)
${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\mathbf{\Phi}_{v}\boldsymbol{\mathcal{R}}_{v}^{(1)},\mathbf{\Phi}_{e}\boldsymbol{\mathcal{R}}_{e}^{(1)},\mathbf{\Phi}_{x}\boldsymbol{\mathcal{R}}_{x}^{(1)})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq\dfrac{1}{6}\max_{0\leq
s\leq
t_{f}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\left(\dfrac{d^{3}\tilde{\boldsymbol{v}}}{dt^{3}}(s),\dfrac{d^{3}\tilde{\boldsymbol{e}}}{dt^{3}}(s),\dfrac{d^{3}\tilde{\boldsymbol{x}}}{dt^{3}}(s)\right)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|},$

where we have used that fact that $\tilde{\boldsymbol{w}}$ is $C^{3}$ on
$[0,t_{f}]$, since $\boldsymbol{F}$ is of class $C^{2}$ on
$\mathcal{V}\times\mathcal{E}\times\mathcal{V}$ and so are
$\mathsf{F}^{\mathfrak{1}}$ and $\mathsf{F}^{tv}$. Combining (5.81) and
(5.83), we conclude that,

(5.84)
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\mathbf{a}^{(3)}_{v,n},\mathbf{a}^{(3)}_{e,n},\mathbf{a}^{(3)}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
$\displaystyle\leq A_{0}+A_{1}\Delta t+A_{2}(\Delta t)^{2}+A_{3}(\Delta
t)^{3},$

where the coefficients $A_{0},A_{1},A_{2},A_{3}$ are given by

(5.85) $\displaystyle A_{0}$ $\displaystyle=\dfrac{1}{6}\max_{0\leq s\leq
t_{f}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\left(\dfrac{d^{3}\tilde{\boldsymbol{v}}}{dt^{3}}(s),\dfrac{d^{3}\tilde{\boldsymbol{e}}}{dt^{3}}(s),\dfrac{d^{3}\tilde{\boldsymbol{x}}}{dt^{3}}(s)\right)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+$
$\displaystyle\quad\quad\dfrac{C_{4}C_{5}}{4}\left(M^{(1)}_{\mathsf{F}^{\mathfrak{1}}}+M^{(1)}_{\mathsf{F}^{tv}}\right)\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}\sqrt{N_{\mathcal{V}}}+$
$\displaystyle\quad\quad\dfrac{C_{4}}{4}\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}M^{(1)}_{\mathsf{F}^{\mathfrak{1}}}\tilde{M}^{(0)}_{\tilde{\mathsf{F}}}+\dfrac{C_{4}C_{5}}{8}N(M_{\mathsf{F}^{\mathfrak{1}}}^{(2)}+M_{\mathsf{F}^{tv}}^{(2)})\left(\tilde{M}^{(0)}_{\tilde{\mathsf{F}}}\right)^{2}$
$\displaystyle A_{1}$
$\displaystyle=\dfrac{C_{4}}{8}\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)^{2}\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}^{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}^{2}M^{(1)}_{\mathsf{F}^{\mathfrak{1}}}M^{(0)}_{\boldsymbol{F}}\sqrt{N_{\mathcal{V}}}+$
$\displaystyle\quad\quad\dfrac{C_{4}}{32}\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}NM_{\mathsf{F}^{\mathfrak{1}}}^{(2)}\tilde{M}^{(0)}_{\tilde{\mathsf{F}}}+$
$\displaystyle\quad\quad\dfrac{C_{4}C_{5}}{8}N(M_{\mathsf{F}^{\mathfrak{1}}}^{(2)}+M_{\mathsf{F}^{tv}}^{(2)})\tilde{M}^{(0)}_{\tilde{\mathsf{F}}}\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}\sqrt{N_{\mathcal{V}}}$
$\displaystyle A_{2}$
$\displaystyle=\dfrac{C_{4}C_{5}}{32}N(M_{\mathsf{F}^{\mathfrak{1}}}^{(2)}+M_{\mathsf{F}^{tv}}^{(2)})\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)^{2}\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}^{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}^{2}\left(M^{(0)}_{\boldsymbol{F}}\right)^{2}N_{\mathcal{V}}+$
$\displaystyle\quad\quad\dfrac{C_{4}}{16}\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)^{2}\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}^{2}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}^{2}NM_{\mathsf{F}^{\mathfrak{1}}}^{(2)}\tilde{M}^{(0)}_{\tilde{\mathsf{F}}}M^{(0)}_{\boldsymbol{F}}\sqrt{N_{\mathcal{V}}}$
$\displaystyle A_{3}$
$\displaystyle=\dfrac{C_{4}}{64}\left(\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{E}}^{-1}\right\rVert_{2}\left\lVert\mathcal{P}_{\mathsf{F}^{tv}}\right\rVert_{2}M^{(0)}_{\boldsymbol{F}}+1\right)^{3}\left\lVert\widehat{\boldsymbol{M}}_{\mathcal{V}}^{-1}\right\rVert_{2}^{3}\left\lVert\mathcal{P}_{\mathsf{F}^{\mathfrak{1}}}\right\rVert_{2}^{3}NM_{\mathsf{F}^{\mathfrak{1}}}^{(2)}\Delta
t_{n}\left(M^{(0)}_{\boldsymbol{F}}\right)^{2}N_{\mathcal{V}}.$

Substituting (5.84) back to (5.69), we see that

(5.86)
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,n+1},\boldsymbol{\delta}_{e,n+1},\boldsymbol{\delta}_{x,n+1})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+A(\Delta
t)^{2}$ $\displaystyle\leq(1+C_{1}\Delta t_{n}+C_{3}(\Delta
t_{n})^{2})\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,n},\boldsymbol{\delta}_{e,n},\boldsymbol{\delta}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+A(\Delta
t)^{2}\right)$ $\displaystyle\leq e^{(C_{1}+C_{3}\Delta t)\Delta
t_{n}}\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,n},\boldsymbol{\delta}_{e,n},\boldsymbol{\delta}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+A(\Delta
t)^{2}\right),$

where $A=C_{1}^{-1}(A_{0}+A_{1}\Delta t+A_{2}(\Delta t)^{2}+A_{3}(\Delta
t)^{3})$. By induction, we have

(5.87)
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,n},\boldsymbol{\delta}_{e,n},\boldsymbol{\delta}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+A(\Delta
t)^{2}\leq e^{(C_{1}+C_{3}\Delta
t)t_{n}}\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,0},\boldsymbol{\delta}_{e,0},\boldsymbol{\delta}_{x,0})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+A(\Delta
t)^{2}\right),$

which implies

(5.88)
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,n},\boldsymbol{\delta}_{e,n},\boldsymbol{\delta}_{x,n})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
$\displaystyle\leq e^{(C_{1}+C_{3}\Delta
t)t_{n}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,0},\boldsymbol{\delta}_{e,0},\boldsymbol{\delta}_{x,0})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+\left(e^{(C_{1}+C_{3}\Delta
t)t_{n}}-1\right)A(\Delta t)^{2}$ $\displaystyle\leq e^{(C_{1}+C_{3}\Delta
t)t_{n}}\left({\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\boldsymbol{\delta}_{v,0},\boldsymbol{\delta}_{e,0},\boldsymbol{\delta}_{x,0})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+A(\Delta
t)^{2}(C_{1}+C_{3}\Delta t)t_{n}\right),$

which yields the desired result. ∎

##  6 Numerical experiments

In this section, we present some numerical results to test the performance of
our proposed method. Our ROM is applied to several Lagrangian hydrodynamics
problems that can be simulated with Laghos444GitHub page,
https://github.com/CEED/Laghos/tree/rom. and libROM555GitHub page,
https://github.com/LLNL/libROM.. Simple command-line options for simulating
these problems are provided in Appendix A for the purpose of reproducible
research. The relative error for each ROM field is measured against the
corresponding FOM solution at the final time $t_{f}$, which is defined as:

(6.1) $\displaystyle\mathcal{\epsilon}_{v,t_{f}}$
$\displaystyle=\frac{\|\boldsymbol{v}_{N_{t}}-\tilde{\boldsymbol{v}}_{\widetilde{N}_{t}}\|_{2}}{\|\boldsymbol{v}_{N_{t}}\|_{2}},$
$\displaystyle\mathcal{\epsilon}_{e,t_{f}}$
$\displaystyle=\frac{\|\boldsymbol{e}_{N_{t}}-\tilde{\boldsymbol{e}}_{\widetilde{N}_{t}}\|_{2}}{\|\boldsymbol{e}_{N_{t}}\|_{2}},$
$\displaystyle\mathcal{\epsilon}_{x,t_{f}}$
$\displaystyle=\frac{\|\boldsymbol{x}_{N_{t}}-\tilde{\boldsymbol{x}}_{\widetilde{N}_{t}}\|_{2}}{\|\boldsymbol{x}_{N_{t}}\|_{2}}.$

The speed-up of each ROM simulation is measured by dividing the wall-clock
time for the ROM time loop by the wall-clock time for the corresponding FOM
time loop. A visualization tool, VisIt [84], is used to plot the meshes and
solution fields both for FOM and ROM. All the simulations in this numerical
section use Quartz in Livermore Computing Center666High performance computing
at LLNL, https://hpc.llnl.gov/hardware/platforms/Quartz, on Intel Xeon CPUs
with 128 GB memory, peak TFLOPS of 3251.4, and peak single CPU memory
bandwidth of 77 GB/s. For simplicity, we only report serial tests, but we have
also observed good speed-up for parallel FOM simulation and serial ROM
simulation.

###  6.1 Problem setting

In this subsection, we introduce the settings of several benchmark problems,
including the Gresho vortex problem, the Sedov Blast problem, the Taylor–Green
vortex problem, and the triple-point problem.

####  6.1.1 Gresho vortex problem

The Gresho vortex problem is a two-dimensional standard benchmark test for the
incompressible inviscid Navier–Stokes equations [85]. In this problem, we
consider a manufactured smooth solution from extending the steady state Gresho
vortex solution to the compressible Euler equations. The computational domain
is the unit square $\widetilde{\Omega}=[-0.5,0.5]^{2}$ with wall boundary
conditions on all surfaces, i.e. $v\cdot n=0$. Let $(r,\phi)$ denote the polar
coordinates of a particle $\widetilde{x}\in\widetilde{\Omega}$. The initial
angular velocity is given by

(6.2) $v_{\phi}=\begin{cases}5r&\text{ for }0\leq r<0.2\\\ 2-5r&\text{ for
}0.2\leq r<0.4\\\ 0&\text{ for }r\geq 0.4.\end{cases}$

The initial density is given by $\rho=1$. The initial thermodynamic pressure
is given by

(6.3) $p=\begin{cases}5+\frac{25}{2}r^{2}&\text{ for }0\leq r<0.2\\\
9-4\log(0.2)+\frac{25}{2}-20r+4\log(r)&\text{ for }0.2\leq r<0.4\\\
3+4\log(2)&\text{ for }r\geq 0.4.\end{cases}$

The initial energy is related to the pressure and the density by (2.2). The
adiabatic index in the ideal gas equations of state is set to be a constant
$\gamma=5/3$. The initial mesh is a uniform Cartesian hexahedral mesh, which
deforms over time. No artificial viscosity stress is added. We compute the
resultant source terms for driving the time-dependent simulation up to some
point in time, and perform normed error analysis on the final computational
mesh. The visualized solution of the Gresho vortex problem is given in the
first column of Fig. 1.

####  6.1.2 Sedov Blast problem

The Sedov blast problem is a three-dimensional standard shock hydrodynamic
benchmark test [86]. In this problem, we consider a delta source of internal
energy initially deposited at the origin of a three-dimensional cube,
analogous to the two-dimensional test in [71], to which we refer for further
details. The computational domain is the unit cube
$\widetilde{\Omega}=[0,1]^{3}$ with wall boundary conditions on all surfaces,
i.e. $v\cdot n=0$. The initial velocity is given by $v=0$. The initial density
is given by $\rho=1$. The initial energy is given by a delta function at the
origin. In our implementation, the delta function energy source is
approximated by setting the internal energy to zero in all degrees of freedom
except at the origin. The default value of the initial internal energy at the
origin is $e(0,0,0)=0.25$. The adiabatic index in the ideal gas equations of
state is set to be a constant $\gamma=1.4$. The initial mesh is a uniform
Cartesian hexahedral mesh, which deforms over time. Artificial viscosity
stress is added. We compute the resultant source terms for driving the time-
dependent simulation up to some point in time, and perform normed error
analysis on the final computational mesh. The visualized solution of the Sedov
blast problem is given in the second column of Fig. 1. It can be seen that the
radial symmetry is maintained in the shock wave propagation, thanks to a high-
order artificial viscosity formulation.

####  6.1.3 Taylor–Green vortex problem

The Taylor–Green vortex problem is a three-dimensional standard flow benchmark
test for the incompressible inviscid Navier–Stokes equations [87]. In this
problem, we consider a manufactured smooth solution from extending the steady
state Taylor–Green vortex solution to the compressible Euler equations,
analogous to the two-dimensional test in [71], to which we refer for further
details. The computational domain is the unit cube
$\widetilde{\Omega}=[0,1]^{3}$ with wall boundary conditions on all surfaces,
i.e. $v\cdot n=0$. The initial velocity is given by

(6.4) $v=(\sin(\pi x)\cos(\pi y)\cos(\pi z),-\cos(\pi x)\sin(\pi y)\cos(\pi
z),0).$

The initial density is given by $\rho=1$. The initial thermodynamic pressure
is given by

(6.5) $p=100+\dfrac{(\cos(2\pi x)+\cos(2\pi y))(\cos(2\pi z)+2)-2}{16}.$

The initial energy is related to the pressure and the density by (2.2). The
adiabatic index in the ideal gas equations of state is set to be a constant
$\gamma=5/3$. The initial mesh is a uniform Cartesian hexahedral mesh, which
deforms over time. Artificial viscosity stress is added. We compute the
resultant source terms for driving the time-dependent simulation up to some
point in time, and perform normed error analysis on the final computational
mesh. The visualized solution of the Taylor–Green vortex problem is given in
the third column of Fig. 1.

####  6.1.4 Triple–point problem

The triple-point problem is a three-dimensional shock test with two materials
in three states [88]. In this problem, we consider a three-state, two-
material, 2D Riemann problem which generates vorticity, analogous to the two-
dimensional test in [71], to which we refer for further details. The
computational domain is $\widetilde{\Omega}=[0,7]\times[0,3]\times[0,1.5]$
with wall boundary conditions on all surfaces, i.e. $v\cdot n=0$. The initial
velocity is given by $v=0$. The initial density is given by

(6.6) $\rho=\begin{cases}1&\text{ if }x\leq 1\text{ or }y\leq 1.5,\\\
1/8&\text{ if }x>1\text{ and }y>1.5.\end{cases}$

The initial thermodynamic pressure is given by

(6.7) $p=\begin{cases}1&\text{ if }x\leq 1,\\\ 0.1&\text{ if }x>1.\end{cases}$

The initial energy is related to the pressure and the density by (2.2). The
adiabatic index in the ideal gas equations of state is set to be

(6.8) $\gamma=\begin{cases}1.5&\text{ if }x\leq 1\text{ or }y>1.5,\\\
1.4&\text{ if }x>1\text{ and }y\leq 1.5.\end{cases}$

The initial mesh is a uniform Cartesian hexahedral mesh, which deforms over
time. Artificial viscosity stress is added. We compute the resultant source
terms for driving the time-dependent simulation up to some point in time, and
perform normed error analysis on the final computational mesh. The visualized
solution of the triple point problem is given in the fourth column of Fig. 1.

We remark that the FOM simulation for all tests is characterized by several
FOM user-defined input values, namely the CFL constant $\alpha$ in (2.15), the
mesh refinement level $m$ which controls the mesh size $h=2^{-m}h_{0}$ where
$h_{0}$ is the coarsest mesh size, and the polynomial order $k$ used in the
$Q_{k}$–$Q_{k-1}$ pair for finite element discretization in [71]. In
particular, the dimensions $N_{\mathcal{V}}$ and $N_{\mathcal{E}}$ of the FOM
state space increase with the mesh size $h$ and the polynomial order $k$.
Table 1 summarizes the FOM user-defined input values and discretization
parameters in various benchmark problems, which characterizes the FOM
simulation in the later sections. The resultant system of semi-discrete
Lagrangian hydrodynamics equations is numerically solved by the RK2-average
scheme introduced in Section 2.2.1.

Problem | Gresho vortex | Sedov Blast | Taylor–Green vortex | Triple–point
---|---|---|---|---
$\alpha$ | 0.5 | 0.5 | 0.1 | 0.5
$m$ | 4 | 2 | 2 | 2
$k$ | 3 | 2 | 2 | 2
$N_{\mathcal{V}}$ | 18818 | 14739 | 14739 | 38475
$N_{\mathcal{E}}$ | 9216 | 4096 | 4096 | 10752
Table 1: FOM user-defined input values and discretization parameters in
various benchmark problems.

At the end of this section, we demonstrate some representative simulation
results on the long-time simulation in various benchmark problems, which will
be further discussed in Section 6.3. The final time $t_{f}$ of the simulation
is taken to be 0.62 for the Gresho vortex problem, 0.8 for the Sedov Blast
problem, 0.25 for the Taylor–Green vortex problem and 0.8 for the triple–point
problem. Figure 1 shows the visualization of final-time velocity and energy on
the Lagrangian frame for long-time FOM simulation in various benchmark
problems. It can be seen that the Lagrangian mesh is extensively distorted and
compressed in the Gresho vortex problem and the Sedov Blast problem, which
suggests that sharp gradients are developed in these problems. Figure 2 shows
the visualization of final-time velocity and energy on the Lagrangian frame
for long-time ROM simulation in various benchmark problems using the time-
windowing approach. It can be seen that the ROM simulation results are in good
agreement with the FOM simulation results in Figure 1. The ROM simulation
results successfully capture the extreme mesh distortion.

Fig. 1: Visualization of final-time velocity (top) and energy (bottom) on the
Lagrangian frame for long-time FOM simulation in various benchmark problems:
the Gresho vortex problem (first column), the Sedov Blast problem (second
column), the Taylor–Green vortex problem (third column) and the triple–point
problem (fourth column).

Fig. 2: Visualization of final-time velocity (top) and energy (bottom) on the
Lagrangian frame for long-time ROM simulation in various benchmark problems:
the Gresho vortex problem (first column), the Sedov Blast problem (second
column), the Taylor–Green vortex problem (third column) and the triple–point
problem (fourth column).

###  6.2 Short-time ROM simulation

In this section, we use the spatial ROM introduced in Section 3 to perform
short-time simulation in various benchmark problems. Only reproductive cases
are considered, i.e. the problem setting in the ROM is identical to that in
the FOM. First, in the offline phase, the fully discrete FOM scheme in Section
2.2 is used to compute snapshot solution for performing POD discussed in
Section 3.5. The reduced basis matrices for the solution variables are then
used to formulate the reduced order model (3.8) in the online phase. In the
case of hyper-reduction, we also obtain the reduced basis matrices for the
nonlinear terms using SVD or the SNS relation, and then follow the DEIM
procedure discussed in Section 3.7 to obtain sampling matrices and formulate
the hyper-reduced system (3.17). Using a time integrator again, the fully
discrete (hyper-)reduced order model is used to compute an approximate
solution. The ROM solution is compared with the FOM solution at the final time
$t_{f}$ with respect to the relative errors (6.1). Also, the wall time of ROM
simulation is compared with that of the FOM simulation. We will study the
dependence of the accuracy and speed-up of ROM on the ratios of dimension
reduction

(6.9)
$\zeta_{v}=\dfrac{n_{v}}{N_{\mathcal{V}}},\quad\zeta_{e}=\dfrac{n_{e}}{N_{\mathcal{E}}},\quad\zeta_{x}=\dfrac{n_{x}}{N_{\mathcal{V}}},\quad\zeta_{\breve{\mathsf{F}^{\mathfrak{1}}}}=\dfrac{n_{\breve{\mathsf{F}^{\mathfrak{1}}}}}{N_{\mathcal{V}}},\quad\zeta_{\breve{\mathsf{F}^{tv}}}=\dfrac{n_{\breve{\mathsf{F}^{tv}}}}{N_{\mathcal{E}}}.$

On the other hand, for hyper-reduction, the numbers of sampling indices
$\\{n_{\breve{\mathsf{F}^{\mathfrak{1}}}},n_{\breve{\mathsf{F}^{tv}}}\\}$ are
controlled by the product of the over-sampling factors
$\\{\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}},\lambda_{\breve{\mathsf{F}^{tv}}}\\}\in[1,\infty)$
and the reduced basis dimensions
$\\{n_{\mathsf{F}^{\mathfrak{1}}},n_{\mathsf{F}^{tv}}\\}$, i.e.

(6.10)
$n_{\breve{\mathsf{F}^{\mathfrak{1}}}}=\min\left\\{N_{\mathcal{V}},\,\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}}n_{\mathsf{F}^{\mathfrak{1}}}\right\\},\quad
n_{\breve{\mathsf{F}^{tv}}}=\min\left\\{N_{\mathcal{E}},\,\lambda_{\breve{\mathsf{F}^{tv}}}n_{\mathsf{F}^{tv}}\right\\}.$

The resultant reduced system of semi-discrete Lagrangian hydrodynamics
equations is numerically solved by the RK2-average scheme introduced in
Section 3.4.

In Figure 3, we compare the ROM performance on the short-time simulation in
the Gresho vortex problem against reduced basis dimensions. This comparison
examines how the solution representability of the reduced order model depends
on the dimensions of the POD reduced bases. To this end, we fix all other ROM
parameters and do not perform hyper-reduction, to avoid complicating the
results. Table 2 shows the reduced basis dimensions being tested and their
ratios to the dimensions of the corresponding FOM finite element spaces. From
Figure 3, we see that the relative errors are overall decreasing with the
reduced basis dimensions, and the accuracy is very outstanding even when low
dimensional subspaces are used. However, since the nonlinear terms scale with
the FOM size, the speed-up achieved is only around 1.4 for all the cases
tested. This highlights the need of hyper-reduction for the sake of remarkable
speed-up.

Fig. 3: ROM performance comparison for short-time simulation in the Gresho vortex problem with varying ratio of reduced basis dimensions: relative $L^{2}$ error for velocity (top-left), relative $L^{2}$ error for energy (top-right), relative $L^{2}$ error for position (bottom-left) and speed-up (bottom-right). Table 2 provides test points.  Test point | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$n_{v}$ | 6 | 11 | 17 | 23 | 28
$n_{e}$ | 9 | 18 | 27 | 36 | 44
$n_{x}$ | 2 | 4 | 6 | 8 | 10
$\zeta_{v}$ | 0.000318844 | 0.000584547 | 0.00090339 | 0.00122223 | 0.00148794
$\zeta_{e}$ | 0.000976562 | 0.00195312 | 0.00292969 | 0.00390625 | 0.00477431
$\zeta_{x}$ | 0.000106281 | 0.000212562 | 0.000318844 | 0.000425125 | 0.000531406
Table 2: List of reduced basis dimensions being tested in Figure 3 and their
ratio to the dimensions of the corresponding FOM finite element spaces for
short-time simulation in the Gresho vortex problem.

In Figure 4, we compare the ROM performance on the short-time simulation in
the Gresho vortex problem against the number of sampling indices. This
comparison examines how the solution accuracy depends on the projection error
of the nonlinear terms in the DEIM nonlinear model reduction. We compare
different approaches for obtaining the nonlinear term bases, i.e. directly
applying snapshot SVD and using the SNS relation [13]. To avoid complicating
the results with respect to different ROM parameters, we fix the reduced basis
dimensions as
$(n_{v},n_{e},n_{x},n_{\mathsf{F}^{\mathfrak{1}}},n_{\mathsf{F}^{tv}})=(28,44,10,28,44)$.
We also compare different algorithms for obtaining the sampling indices,
namely the over-sampling DEIM (see Algorithm 3 of [81] and Algorithm 5 of
[82]) and QDEIM (see [79]). Table 3 shows the number of sampling indices being
tested and their ratio to the dimensions of the corresponding FOM finite
element spaces. From Figure 4, we see that the relative errors and speed-up
are overall decreasing with the number of sampling indices, while the accuracy
is still remarkable for each of the nonlinear model reduction technique. We
remark that the first test point with DEIM has a numerical instability as the
time step size eventually vanishes. In this experiment, we observe that SNS-
DEIM has the highest accuracy and best speed-up compared with the other
approaches. In the rest of the paper, unless otherwise specified, we will
mainly focus on the SNS-DEIM approach as suggested by this experiment.

Fig. 4: ROM performance comparison for short-time simulation in the Gresho vortex problem with varying ratio of sampling indices: relative $L^{2}$ error for velocity (top-left), relative $L^{2}$ error for energy (top-right), relative $L^{2}$ error for position (bottom-left) and speed-up (bottom-right). Table 3 provides test points.  Test point | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
$\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}}$ | 8 | 15 | 23 | 30 | 37 | 45
$\lambda_{\breve{\mathsf{F}^{tv}}}$ | 7 | 14 | 21 | 28 | 35 | 41
$n_{\breve{\mathsf{F}^{\mathfrak{1}}}}$ | 616 | 1232 | 1848 | 2464 | 3080 | 3696
$n_{\breve{\mathsf{F}^{tv}}}$ | 396 | 792 | 1188 | 1584 | 1980 | 2376
$\zeta_{\breve{\mathsf{F}^{\mathfrak{1}}}}$ | 0.0327346 | 0.0654692 | 0.0982038 | 0.130938 | 0.163673 | 0.196408
$\zeta_{\breve{\mathsf{F}^{tv}}}$ | 0.0429688 | 0.0859375 | 0.128906 | 0.171875 | 0.214844 | 0.257812
Table 3: List of number of sampling indices being tested in Figure 4 and their
ratio to the dimensions of the corresponding FOM finite element spaces for
short-time simulation in the Gresho vortex problem.

Table 4 shows the short-time numerical results using SNS-DEIM spatial ROM for
various benchmark problems. The quantities reported in Table 4 include

  * •

the final time $t_{f}$ and the number of time steps $N_{t}$ in the FOM
simulations,

  * •

the ROM user-defined input values, namely the threshold $\epsilon_{\sigma}$
and the over-sampling factors
$\\{\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}},\lambda_{\breve{\mathsf{F}^{tv}}}\\}$,

  * •

the ROM parameters, namely the reduced dimensions $\\{n_{v},n_{e},n_{x}\\}$
(note that with our SNS setting, we have $n_{\mathsf{F}^{\mathfrak{1}}}=n_{v}$
and $n_{\mathsf{F}^{tv}}=n_{e}$) and the number of sampling indices
$\\{n_{\breve{\mathsf{F}^{\mathfrak{1}}}},n_{\breve{\mathsf{F}^{tv}}}\\}$, and

  * •

the ROM results, namely the number of time steps $\widetilde{N}_{t}$ in the
ROM simulations, the relative errors
$\\{\mathcal{\epsilon}_{v,t_{f}},\mathcal{\epsilon}_{e,t_{f}},\mathcal{\epsilon}_{x,t_{f}}\\}$
at the final time, and the relative speed-up of the ROM simulations to FOM
simulations.

In each of the benchmark problems, it can be observed that the number of time
steps $\widetilde{N}_{t}$ in the ROM simulations is close to but slightly
differ from $N_{t}$ in the FOM simulations. The ROM simulation benefits from
the relatively cheap reduced dimension computation in each time step to
achieve an overall relative speed-up. The relative error is small for each
variable, indicating that the ROM solutions provide accurate approximations to
the FOM solutions.

Problem | Gresho vortex | Sedov Blast | Taylor–Green vortex | Triple–point
---|---|---|---|---
$t_{f}$ | 0.1 | 0.1 | 0.05 | 0.2
$N_{t}$ | 87 | 245 | 122 | 28
$\epsilon_{\sigma}$ | 0.9999 | 0.9999 | 0.9999 | 0.9999
$\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}}$ | 8 | 2 | 2 | 2
$\lambda_{\breve{\mathsf{F}^{tv}}}$ | 8 | 2 | 2 | 2
$n_{v}$ | 32 | 53 | 9 | 10
$n_{e}$ | 72 | 13 | 16 | 10
$n_{x}$ | 10 | 17 | 3 | 7
$n_{\breve{\mathsf{F}^{\mathfrak{1}}}}$ | 256 | 106 | 18 | 20
$n_{\breve{\mathsf{F}^{tv}}}$ | 576 | 26 | 32 | 20
$\widetilde{N}_{t}$ | 87 | 225 | 122 | 29
$\mathcal{\epsilon}_{v,t_{f}}$ | 1.385366e-02 | 2.74144e-03 | 9.05004e-05 | 9.41694e-03
$\mathcal{\epsilon}_{e,t_{f}}$ | 2.02558e-03 | 8.39586e-04 | 3.69690e-06 | 4.91346e-03
$\mathcal{\epsilon}_{x,t_{f}}$ | 1.76761e-04 | 2.04953e-05 | 6.81088e-07 | 4.03934e-05
speed-up | 8.27130 | 26.55553 | 14.75786 | 53.02794
Table 4: FOM results, ROM user-defined input values, parameters and results
for short-time simulation in various benchmark problems.

###  6.3 Long-time ROM simulation

In this section, we use the spatial ROM introduced in Section 3 and the time-
windowing ROM introduced in Section 4 to perform long-time simulation in
various benchmark problems. Again, only reproductive cases are considered. For
the spatial ROM, we use SNS-DEIM for hyper-reduction and the offline-online
procedure splitting is the same as in Section 6.2. Table 5 shows the long-time
numerical results using SNS-DEIM spatial ROM for various benchmark problems.
It can be observed that the relative errors are still reasonable, but the
relative speed-up is lower than that in the short-time simulations. In these
simulations, large over-sampling factors
$\\{\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}},\lambda_{\breve{\mathsf{F}^{tv}}}\\}$
are used to maintain sufficient numbers of sampling indices
$\\{n_{\breve{\mathsf{F}^{\mathfrak{1}}}},n_{\breve{\mathsf{F}^{tv}}}\\}$, for
otherwise the time step size would eventually vanish due to the adaptive time
stepping control and the CFL constraints. However, larger numbers of sampling
indices lead to more expensive overdetermined systems (3.13) and (3.15) for
the hyper-reduction and hurt the overall speed-up in ROM simulations. It is
especially important to point out that the spatial ROM does not achieve a
speed-up in the Gresho vortex problem since more reduced basis vectors are
required to ensure the solution representability of the reduced linear space,
and almost full sampling has to be used for the sake of numerical stability.

Problem | Gresho vortex | Sedov Blast | Taylor–Green vortex | Triple–point
---|---|---|---|---
$t_{f}$ | 0.62 | 0.8 | 0.25 | 0.8
$N_{t}$ | 1672 | 719 | 897 | 193
$\epsilon_{\sigma}$ | 0.9999 | 0.9999 | 0.9999 | 0.9999
$\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}}$ | 100 | 6 | 80 | 6
$\lambda_{\breve{\mathsf{F}^{tv}}}$ | 100 | 6 | 40 | 6
$n_{v}$ | 165 | 225 | 40 | 20
$n_{e}$ | 331 | 27 | 62 | 16
$n_{x}$ | 34 | 29 | 6 | 10
$n_{\breve{\mathsf{F}^{\mathfrak{1}}}}$ | 16500 | 1350 | 3200 | 120
$n_{\breve{\mathsf{F}^{tv}}}$ | 9216 | 162 | 2480 | 96
$\widetilde{N}_{t}$ | 1670 | 742 | 898 | 193
$\mathcal{\epsilon}_{v,t_{f}}$ | 4.29468e-02 | 9.72044e-03 | 3.60272e-02 | 2.52978e-03
$\mathcal{\epsilon}_{e,t_{f}}$ | 6.82899e-03 | 5.95163e-04 | 3.28395e-03 | 2.36835e-03
$\mathcal{\epsilon}_{x,t_{f}}$ | 7.1045e-05 | 3.67286e-04 | 2.77386e-04 | 4.72467e-05
speed-up | 0.79646 | 3.35523 | 1.96808 | 15.29550
Table 5: FOM results, ROM user-defined input values, parameters and results
for long-time simulation in various benchmark problems.

Next, we examine the performance of the time windowing ROM approach. In our
implementation, the end points of the time windows $\\{T_{w}\\}_{w=1}^{N_{w}}$
can either be determined by physical time as in Section 4.4.1, or by the
number of samples as in Section 4.4.2. We focus our discussion on the SNS-DEIM
time windowing ROM approach with the end points of the time windows determined
by the number of samples. In the offline phase, the fully discrete RK2-average
FOM scheme is first used to compute solution snapshots. The time windows are
determined according to (4.21). Following Section 4.5, the collected solution
data are clustered into different windows for performing POD, and the
resultant reduced solution bases for velocity and energy are left-multiplied
by the corresponding mass matrix to obtain the nonlinear term bases, which in
turn determine the sampling indices by DEIM. In the online phase, the reduced
basis matrices and the sampling matrices are then used to formulate the hyper-
reduced system (4.10) in each time window. The fully discrete hyper-reduced
system follows from applying the RK2-average scheme as discussed in Section
3.4.

Table 6 shows the long-time numerical results for various benchmark problems
with time windowing ROM. In the case of decomposing time windows by number of
samples, the number of time windows $N_{w}$ is an additional ROM parameter
controlled by the number of samples per window $N_{\text{sample}}$ as an
additional user-defined input. We use $N_{\text{sample}}=10$ in each of the
benchmark problems. Moreover, in our implementation, the over-sampling factors
$\\{\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}},\lambda_{\breve{\mathsf{F}^{tv}}}\\}$
are user-defined constants over all time windows, which control the number of
sampling indices in the window $w$ by

(6.11)
$n_{\breve{\mathsf{F}^{\mathfrak{1}}}}^{w}=\min\left\\{N_{\mathcal{V}},\,\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}}n_{\mathsf{F}^{\mathfrak{1}}}^{w}\right\\},\quad
n_{\breve{\mathsf{F}^{tv}}}^{w}=\min\left\\{N_{\mathcal{E}},\,\lambda_{\breve{\mathsf{F}^{tv}}}n_{\mathsf{F}^{tv}}^{w}\right\\}.$

In Table 6, the reported values of reduced basis dimensions and the number of
sampling indices are taken as the greatest among all windows, i.e.

(6.12) $n_{v}=\max_{1\leq w\leq N_{w}}n_{v}^{w},\quad n_{e}=\max_{1\leq w\leq
N_{w}}n_{e}^{w},\quad n_{x}=\max_{1\leq w\leq N_{w}}n_{x}^{w},\quad
n_{\breve{\mathsf{F}^{\mathfrak{1}}}}=\max_{1\leq w\leq
N_{w}}n_{\breve{\mathsf{F}^{\mathfrak{1}}}}^{w},\quad
n_{\breve{\mathsf{F}^{tv}}}=\max_{1\leq w\leq
N_{w}}n_{\breve{\mathsf{F}^{tv}}}^{w}.$

It can be observed from comparing Table 6 with Table 5 that it is possible to
achieve a higher relative speed-up and improve solution accuracy using time
windowing ROM. In particular, for the Gresho vortex problem, when the spatial
ROM fails to achieve any relative speed-up, the time windowing approach yields
relative error or order $10^{-6}$ in all solution variables and a $25$ times
speed-up. This suggests that the time windowing approach is more useful for
the advection-dominated cases, in which the intrinsic dimensions of reduced
subspaces in the spatial ROM are large.

Problem | Gresho vortex | Sedov Blast | Taylor–Green vortex | Triple–point
---|---|---|---|---
$t_{f}$ | 0.62 | 0.8 | 0.25 | 0.8
$N_{t}$ | 1672 | 719 | 897 | 193
$N_{\text{sample}}$ | 10 | 10 | 10 | 10
$\epsilon_{\sigma}$ | 0.9999 | 0.9999 | 0.9999 | 0.9999
$\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}}$ | 2 | 2 | 2 | 2
$\lambda_{\breve{\mathsf{F}^{tv}}}$ | 2 | 2 | 2 | 2
$N_{w}$ | 335 | 144 | 180 | 39
$n_{v}$ | 8 | 10 | 3 | 7
$n_{e}$ | 10 | 7 | 5 | 8
$n_{x}$ | 6 | 7 | 2 | 5
$n_{\breve{\mathsf{F}^{\mathfrak{1}}}}$ | 16 | 20 | 6 | 14
$n_{\breve{\mathsf{F}^{tv}}}$ | 20 | 14 | 10 | 16
$\widetilde{N}_{t}$ | 1599 | 697 | 897 | 155
$\mathcal{\epsilon}_{v,t_{f}}$ | 3.50163e-06 | 4.58111e-04 | 1.22769e-06 | 9.31072e-04
$\mathcal{\epsilon}_{e,t_{f}}$ | 3.51595e-07 | 2.95541e-05 | 5.35795e-08 | 6.99656e-04
$\mathcal{\epsilon}_{x,t_{f}}$ | 4.28729e-07 | 3.86835e-05 | 6.95479e-09 | 3.82897e-05
speed-up | 24.59098 | 26.47649 | 23.46642 | 75.03020
Table 6: FOM results, time windowing ROM user-defined input values, parameters
and results for long-time simulation in various benchmark problems.

At the end of this subsection, we compare the performance of the ROM approach
with SNS-DEIM across different window size on long-time simulation in the
Gresho vortex problem and the Sedov Blast problem. We will consider the
spatial ROM in Section 3 and the time windowing ROM in Section 4. The
endpoints of the time windows $\\{T_{w}\\}_{w=1}^{N_{w}}$ are determined by
(4.21) with the number of samples $N_{\text{sample}}$. More precisely, the
cases under consideration in the Gresho vortex problem are:

  * •

$N_{w}=1$ corresponds to serial ROM,

  * •

$N_{\text{sample}}=10$ corresponds to time windowing ROM which generates
$N_{w}=335$ time windows, and

  * •

$N_{\text{sample}}=20$ corresponds to time windowing ROM which generates
$N_{w}=168$ time windows, and

  * •

$N_{\text{sample}}=40$ corresponds to time windowing ROM which generates
$N_{w}=84$ time windows.

Meanwhile, the cases under consideration in the Sedov Blast problem are:

  * •

$N_{w}=1$ corresponds to serial ROM,

  * •

$N_{\text{sample}}=10$ corresponds to time windowing ROM which generates
$N_{w}=144$ time windows, and

  * •

$N_{\text{sample}}=20$ corresponds to time windowing ROM which generates
$N_{w}=72$ time windows, and

  * •

$N_{\text{sample}}=40$ corresponds to time windowing ROM which generates
$N_{w}=36$ time windows.

Depending on the reduced dimensions
$\\{n_{\mathsf{F}^{\mathfrak{1}}},n_{\mathsf{F}^{tv}}\\}$ of the nonlinear
term bases, each case is tested with various combinations of over-sampling
factors
$(\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}},\lambda_{\breve{\mathsf{F}^{tv}}})$
to obtain adjustable solution accuracy and speed-up.

For each of the problems, we construct a Pareto front in Figure 5 for each of
the cases, which is characterized by the ROM user-defined input values that
minimize the competing objectives of relative $L^{2}$ error for velocity and
relative wall time. An overall Pareto front that selects the ROM parameters
that are Pareto-optimal across all groups is also illustrated for each
problem. Table 7 reports the ROM user-defined input values that yielded the
results on the overall Pareto front. The results also suggest that the time
windowing approaches can produce a remarkable speed-up as well as accurate
approximations with the appropriate use of hyper-reduction.

Fig. 5: ROM performance comparison for long-time simulation in the Gresho vortex problem (top) and the Sedov Blast problem (bottom). Relative $L^{2}$ error for velocity versus relative wall time for varying ROM parameters.  Label | 1-(1) | 1-(2) | 1-(3) | 2-(1) | 2-(2) | 2-(3) | 2-(4)
---|---|---|---|---|---|---|---
$N_{w}$ | – | – | – | – | – | – | –
$N_{\text{sample}}$ | 20 | 20 | 10 | 10 | 20 | 20 | 40
$\epsilon_{\sigma}$ | 0.9999 | 0.9999 | 0.9999 | 0.9999 | 0.9999 | 0.9999 | 0.9999
$\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}}$ | 2 | 2 | 4 | 2 | 2 | 4 | 8
$\lambda_{\breve{\mathsf{F}^{tv}}}$ | 2 | 2 | 4 | 2 | 2 | 4 | 8
Table 7: ROM user-defined input values yielding Pareto-optimal performance for
long-time simulation in the Gresho vortex problem and the Sedov Blast problem.
Figure 5 provides labels.

###  6.4 Parametric ROM simulation

In this section, we use the time-windowing ROM introduced in Section 4 to
perform long-time simulation in the Sedov Blast problem in a parametric
problem setting. In our numerical experiments, the initial internal energy
deposited at the origin is parametrized by the problem parameter
$\boldsymbol{\mu}\in\mathsf{D}=\mathbb{R}^{+}$, where we set
$e(0,0,0;\boldsymbol{\mu})=0.25\boldsymbol{\mu}$. First, in the offline phase,
the fully discrete FOM scheme is used to compute snapshot solution with
several problem parameters $\\{\boldsymbol{\mu}_{k}\\}_{k=1}^{n_{\mu}}$. Then,
we consider the SNS-DEIM time windowing ROM approach with the end points of
the time windows determined by number of samples as in Section 4.4.2.

The FOM solutions with $n_{\mu}=3$ problem parameters, namely
$\\{\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2},\boldsymbol{\mu}_{3}\\}=\\{0.8,1.0,1.2\\}$,
are computed and collected as snapshots for determining the time windows. Then
POD and SNS-DEIM are performed in each window using each choice of the offset
vectors as discussed in Section 4.6. The time windows are determined according
to (4.21) with $N_{\text{sample}}=20$. Then the ROM is tested on a generic
problem parameter $\boldsymbol{\mu}$, where the offset vectors
$\boldsymbol{v}_{\text{os}}^{w}(\boldsymbol{\mu}),\boldsymbol{e}_{\text{os}}^{w}(\boldsymbol{\mu}),\boldsymbol{x}_{\text{os}}^{w}(\boldsymbol{\mu})$
are computed in each of the time windows $w$ according to the choice of offset
vectors used in constructing the solution bases. In Figure 6, we compare the
ROM performance on the long-time simulation in the Sedov Blast problem against
the generic problem parameter $\boldsymbol{\mu}\in[0.7,1.3]$. This comparison
examines how the solution representability of the reduced order model depends
on the problem parameter and the choice of offset vectors. To this end, we fix
all other ROM parameters and do not perform hyper-reduction to avoid
complicating the results. The comparison is also made against the spatial ROM.
As shown in Figure 6, the time windowing approach provides a much better
speed-up than the serial ROM. We also remark that the serial ROM has a
numerical instability in $\boldsymbol{\mu}\in\\{1.2,1.3\\}$ as the time step
size eventually vanishes, which we do not see in the time windowing ROM. We
can observe from Figure 6 that using the interpolation scheme in Section 4.6.3
outperforms the other approaches of computing offset vectors in terms of the
relative errors in the reproductive cases
$\boldsymbol{\mu}\in\\{0.8,1.0,1.2\\}$, at which the offset vectors are
exactly the interpolating values for each time window. However, this approach
produces larger error at the predictive cases. This suggests that the error of
offset determination by the interpolation scheme plays a crucial part in the
final-time error. On the other hand, using the the initial state as offset
vectors, which is described in Section 4.6.1, provides uniformly improved
accuracy when compared with serial ROM. Moreover, it provides a speed-up of
around 15 times at all the tested parameter values, which outperforms the
interpolation scheme.

Fig. 6: Time-windowing ROM performance comparison for long-time simulation in
the parametric Sedov Blast problem with varying problem parameter: relative
$L^{2}$ error for velocity (top-left), relative $L^{2}$ error for energy (top-
right), relative $L^{2}$ error for position (bottom-left) and speed-up
(bottom-right).

At the end of this subsection, we compare the performance of various
parametric time windowing ROM approaches across a wide variation of all ROM
parameters on long-time simulation in the Sedov Blast problem in the setting
of varying the energy factor. The FOM solutions with $n_{\mu}=3$ problem
parameters, namely
$\\{\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2},\boldsymbol{\mu}_{3}\\}=\\{0.8,1.0,1.2\\}$,
are computed and collected as snapshots for determining the time windows, and
then performing POD and SNS-DEIM in each window using each choice of the
offset vectors as discussed in Section 4.6. The time windows are determined
according to (4.21) with $N_{\text{sample}}=20$. The ROM performance is tested
on the problem parameters $\boldsymbol{\mu}\in\\{0.7,0.8,0.9\\}$. Again, we
compare different choices of the offset vectors as discussed in Section 4.6,
and the comparison is also made against the spatial ROM. Depending on the
reduced dimensions $\\{n_{\mathsf{F}^{\mathfrak{1}}},n_{\mathsf{F}^{tv}}\\}$
of the nonlinear term bases, each case is tested with various combinations of
over-sampling factors
$(\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}},\lambda_{\breve{\mathsf{F}^{tv}}})$
to obtain adjustable solution accuracy and speed-up.

We remark that the tested parameter $\boldsymbol{\mu}=0.8$ is a reproductive
case, the tested parameter $\boldsymbol{\mu}=0.9$ is a predictive case in the
interpolation regime, while the case $n_{\mu}=0.7$ is predictive in the
extrapolation regime. For each of the tested problem parameters
$\boldsymbol{\mu}\in\\{0.7,0.8,0.9\\}$, we construct a Pareto front in Figure
7 for each of the cases, which is characterized by the ROM user-defined input
values that minimize the competing objectives of relative $L^{2}$ error for
velocity and relative wall time. An overall Pareto front that selects the ROM
parameters that are Pareto-optimal across all groups is also illustrated for
each problem. Table 8 reports the ROM user-defined input values that yielded
the results on the overall Pareto front. The numerical results show that for
the predictive cases $\boldsymbol{\mu}\in\\{0.7,0.9\\}$, using the initial
state as offset vectors described in Section 4.6.1 is Pareto optimal in any
region. It is also Pareto optimal for smaller relative wall time for the
reproductive case $\boldsymbol{\mu}=0.8$. In larger relative wall time, using
the interpolation scheme in Section 4.6.3 is Pareto optimal for
$\boldsymbol{\mu}=0.8$.

Fig. 7: ROM performance comparison for long-time simulation in the parametric Sedov Blast problem with the tested problem parameters $\boldsymbol{\mu}=0.7$ (top), $\boldsymbol{\mu}=0.8$ (mid) and $\boldsymbol{\mu}=0.9$ (bottom). Relative $L^{2}$ error for velocity versus relative wall time for varying ROM parameters.  Label | 3-(1) | 3-(2) | 3-(3) | 3-(4) | 4-(1) | 4-(2) | 4-(3) | 5-(1) | 5-(2) | 5-(3)
---|---|---|---|---|---|---|---|---|---|---
$N_{w}$ | – | – | – | 1 | – | – | – | – | – | –
$N_{\text{sample}}$ | 20 | 40 | 80 | – | 20 | 80 | 20 | 20 | 40 | 80
$\epsilon_{\sigma}$ | 0.9999 | 0.9999 | 0.9999 | 0.9999 | 0.9999 | 0.9999 | 0.9999 | 0.9999 | 0.9999 | 0.9999
Offset | Initial | Initial | Initial | – | Initial | Initial | Interpolate | Initial | Initial | Initial
$\lambda_{\breve{\mathsf{F}^{\mathfrak{1}}}}$ | 4 | 4 | 4 | 32 | 4 | 4 | 8 | 4 | 4 | 32
$\lambda_{\breve{\mathsf{F}^{tv}}}$ | 4 | 4 | 4 | 32 | 4 | 4 | 8 | 4 | 4 | 32
Table 8: ROM user-defined input values yielding Pareto-optimal performance for
long-time simulation in the parametric Sedov Blast problem. Figure 7 provides
labels.

##  7 Conclusion

In this paper, we develop an efficient reduced order modeling approach for
Lagrangian hydrodynamics simulations. A time-windowing approach is introduced
to construct temporally-local ROM spaces which are small but accurate within a
short period in advection-dominated problems. Different techniques of window
division and construction of offset variables are discussed. Over-sampling
hyper-reduction on the nonlinear terms are introduced to ensure adequate
speed-up. Error bounds for the fully discrete reduced order model with
adaptive time-step control are developed using the continuous-in-time full
order model solution as the reference. The error bounds are controlled by
several quantities, including the oblique projection error of the solution and
the nonlinear terms onto reduced subspaces, the mismatch of initial condition,
an exponential factor of the time of evaluation, and the maximum time step
size. Numerical examples are shown to verify the capability of our method. For
the long-time simulation in benchmark advection-dominated problems, the time
windowing reduced order model successfully captures the extreme mesh
distortion exhibited in the Lagrangian full order model, and attains
remarkable speed-up of 25 to 75 times and outstanding solution accuracy with
the order of -9 to -4 in the reproductive cases. As suggested by the numerical
results, the subspace relation is useful for extracting quality nonlinear term
bases and improving the solution accuracy and overall speed-up. Various offset
approaches in the time windowing reduced order model are compared. In the
reproductive case, using the interpolation scheme is Pareto optimal in the
region of higher solution accuracy, while using the initial states is Pareto
optimal in the region of smaller relative clock time. In the predictive cases,
using the initial states provides uniformly improved accuracy when compared
with serial ROM and a speed-up of around 15 times at all the tested parameter
values.

## Acknowledgments

This work was performed at Lawrence Livermore National Laboratory and
partially funded by LDRD (21-FS-042). Lawrence Livermore National Laboratory
is operated by Lawrence Livermore National Security, LLC, for the U.S.
Department of Energy, National Nuclear Security Administration under Contract
DE-AC52-07NA27344 and LLNL-JRNL-820660.

## Disclaimer

This document was prepared as an account of work sponsored by an agency of the
United States government. Neither the United States government nor Lawrence
Livermore National Security, LLC, nor any of their employees makes any
warranty, expressed or implied, or assumes any legal liability or
responsibility for the accuracy, completeness, or usefulness of any
information, apparatus, product, or process disclosed, or represents that its
use would not infringe privately owned rights. Reference herein to any
specific commercial product, process, or service by trade name, trademark,
manufacturer, or otherwise does not necessarily constitute or imply its
endorsement, recommendation, or favoring by the United States government or
Lawrence Livermore National Security, LLC. The views and opinions of authors
expressed herein do not necessarily state or reflect those of the United
States government or Lawrence Livermore National Security, LLC, and shall not
be used for advertising or product endorsement purposes.

##  Appendix A Command line options of Laghos

In this section, we present some examples of the command line options of
Laghos simulation for the purpose of reproducible research. Due to rapid
software development in the repository, the commands lines are subject to
change. However, we try our best to present the command lines compatible with
recent versions of different dependent softwares and maintain a simple usage
of the program. The following command lines are compatible with the recent
commits of the master branch of MFEM 777GitHub page,
https://github.com/mfem/mfem, commit 4fb1a54, the master branch of libROM
888GitHub page, https://github.com/LLNL/libROM, commit 143d4f6, and the rom
branch of Laghos 999GitHub page, https://github.com/CEED/Laghos/tree/rom,
commit 934badf. In order to use the ROM capability of Laghos, a user has to
navigate to the rom subdirectory.

###  A.1 Problem specification

First, we present the commands lines for the FOM user-defined input values
reported in Table 1 using the executable laghos. These command lines can also
be used along with command line options for ROM user-defined input values. We
use the long-time simulation for illustration. The command line options for
specifying the Gresho vortex problem are

    ./laghos -p 4 -m data/square_gresho.mesh -rs 4 -ok 3 -ot 2 -tf 0.62 -s 7

The command line options for specifying the Sedov Blast problem are

    ./laghos -p 1 -m data/cube01_hex.mesh -pt 211 -tf 0.8 -s 7

The command line options for specifying the Taylor–Green vortex problem are

    ./laghos -p 0 -m data/cube01_hex.mesh -cfl 0.1 -tf 0.25 -s 7

The command line options for specifying the triple–point problem are

    ./laghos -p 3 -m data/box01_hex.mesh -tf 0.8 -s 7

###  A.2 Long-time serial ROM simulation

Next, we present the command line options for the long-time serial ROM
simulation in the Sedov Blast problem reported in the second column of Table
5. The options are appended to the command lines for specifying the Sedov
Blast problem in Section A.1. The other problems in Table 5 can be simulated
using similar options. For computing FOM reference solution, snapshot sampling
and basis generation in the offline phase, one appends

    -offline -writesol -romsns -ef 0.9999 -sample-stages

For hyper-reduction preprocessing in the online phase, one appends

    -online -romhrprep -romsns -rdimv 225 -rdime 27 -rdimx 29 -sfacv 6 -sface 6

For ROM simulation in the online phase, one appends

    -online -romhr -romsns -rdimv 225 -rdime 27 -rdimx 29 -sfacv 6 -sface 6

Finally, for solution postprocessing and calculating the relative error, one
appends

    -restore -soldiff -romsns -rdimv 225 -rdime 27 -rdimx 29

###  A.3 Long-time time windowing ROM simulation

Next, we present the command line options for the long-time time windowing ROM
simulation in the Sedov Blast problem reported in the second column of Table
6. The options are appended to the command lines for specifying the Sedov
Blast problem in Section A.1. The other problems in Table 6 can be simulated
using similar options. For computing FOM reference solution, snapshot sampling
and basis generation in the offline phase, one appends

    -offline -writesol -romsns -ef 0.9999 -sample-stages -nwinsamp 10 -rostype load

For hyper-reduction preprocessing in the online phase, one appends

    -online -romhrprep -romsns -sfacv 2 -sface 2 -nwin 144 -rostype load

For ROM simulation in the online phase, one appends

    -online -romhr -romsns -sfacv 2 -sface 2 -nwin 144 -rostype load

Finally, for solution postprocessing and calculating the relative error, one
appends

    -restore -soldiff -romsns -nwin 144 -rostype load

###  A.4 Parametric time windowing ROM simulation

Finally, we present the command line options for the long-time time windowing
ROM simulation in the parametric Sedov Blast problem reported in the Label
4-(1) of Table 8. The options are appended to the command lines for specifying
the Sedov Blast problem in Section A.1, except that for the basis generation a
separate executable merge is used. The other labels in Table 8 can be
simulated using similar options. For computing FOM reference solution with
$\boldsymbol{\mu}=0.8$ and snapshot sampling in the offline phase, one appends

    -offline -romsns -bef 1.0 -rpar 0 -sample-stages -rostype initial
    -offline -romsns -bef 1.2 -rpar 1 -sample-stages -rostype initial
    -offline -romsns -bef 0.8 -rpar 2 -sample-stages -rostype initial -writesol

For basis generation in the offline phase of the parametric problem, one uses
a separate executable merge

    ./merge -nset 3 -romsns -ef 0.9999 -nwinsamp 20 -rostype initial

For hyper-reduction preprocessing in the online phase, one appends

    -online -romhrprep -romsns -bef 0.8 -sfacv 4 -sface 4 -nwin 80 -rostype initial

For ROM simulation in the online phase, one appends

    -online -romhr -romsns -bef 0.8 -sfacv 4 -sface 4 -nwin 80 -rostype initial

Finally, for solution postprocessing and calculating the relative error, one
appends

    -restore -soldiff -romsns -nwin 80 -rostype initial

## References

  * [1] Shun Wang, Eric de Sturler, and Glaucio H Paulino.  Large-scale topology optimization using preconditioned Krylov subspace methods with recycling.  International journal for numerical methods in engineering, 69(12):2441–2468, 2007.
  * [2] Miguel A Salazar de Troya and Daniel A Tortorelli.  Three-dimensional adaptive mesh refinement in stress-constrained topology optimization.  Structural and Multidisciplinary Optimization, 62(5):2467–2479, 2020.
  * [3] Miguel A Salazar De Troya and Daniel A Tortorelli.  Adaptive mesh refinement in stress-constrained topology optimization.  Structural and Multidisciplinary Optimization, 58(6):2369–2386, 2018.
  * [4] Daniel A White, Youngsoo Choi, and Jun Kudo.  A dual mesh method with adaptivity for stress-constrained topology optimization.  Structural and Multidisciplinary Optimization, 61(2):749–762, 2020\.
  * [5] Youngsoo Choi, Charbel Farhat, Walter Murray, and Michael Saunders.  A practical factorization of a Schur complement for PDE-constrained distributed optimal control.  Journal of Scientific Computing, 65(2):576–597, 2015.
  * [6] Youngsoo Choi.  Simultaneous analysis and design in PDE-constrained optimization.  PhD thesis, Stanford University, 2012.
  * [7] Ralph C Smith.  Uncertainty quantification: theory, implementation, and applications, volume 12.  Siam, 2013.
  * [8] Lorenz Biegler, George Biros, Omar Ghattas, Matthias Heinkenschloss, David Keyes, Bani Mallick, Luis Tenorio, Bart van Bloemen Waanders, Karen Willcox, and Youssef Marzouk.  Large-scale inverse problems and quantification of uncertainty, volume 712.  John Wiley & Sons, 2011.
  * [9] David Galbally, Krzysztof Fidkowski, Karen Willcox, and Omar Ghattas.  Non-linear model reduction for uncertainty quantification in large-scale inverse problems.  International journal for numerical methods in engineering, 81(12):1581–1608, 2010.
  * [10] Chi Hoang, Youngsoo Choi, and Kevin Carlberg.  Domain-decomposition least-squares Petrov-Galerkin (DD-LSPG) nonlinear model reduction.  arXiv preprint arXiv:2007.11835, 2020.
  * [11] Felix Fritzen, Bernard Haasdonk, David Ryckelynck, and Sebastian Schöps.  An algorithmic comparison of the hyper-reduction and the discrete empirical interpolation method for a nonlinear thermal problem.  Mathematical and computational applications, 23(1):8, 2018.
  * [12] Youngsoo Choi and Kevin Carlberg.  Space–time least-squares Petrov–Galerkin projection for nonlinear model reduction.  SIAM Journal on Scientific Computing, 41(1):A26–A58, 2019.
  * [13] Youngsoo Choi, Deshawn Coombs, and Robert Anderson.  SNS: a solution-based nonlinear subspace method for time-dependent model order reduction.  SIAM Journal on Scientific Computing, 42(2):A1116–A1146, 2020.
  * [14] Kevin Carlberg, Youngsoo Choi, and Syuzanna Sargsyan.  Conservative model reduction for finite-volume models.  Journal of Computational Physics, 371:280–314, 2018.
  * [15] Rambod Mojgani and Maciej Balajewicz.  Lagrangian basis method for dimensionality reduction of convection dominated nonlinear flows.  arXiv preprint arXiv:1701.04343, 2017.
  * [16] Youngkyu Kim, Karen May Wang, and Youngsoo Choi.  Efficient space-time reduced order model for linear dynamical systems in Python using less than 120 lines of code.  arXiv preprint arXiv:2011.10648, 2020.
  * [17] Dunhui Xiao, Fangxin Fang, Andrew G Buchan, Christopher C Pain, Ionel Michael Navon, Juan Du, and G Hu.  Non-linear model reduction for the navier–Stokes equations using residual deim method.  Journal of Computational Physics, 263:1–18, 2014.
  * [18] John Burkardt, Max Gunzburger, and Hyung-Chun Lee.  POD and CVT-based reduced-order modeling of Navier–Stokes flows.  Computer methods in applied mechanics and engineering, 196(1-3):337–355, 2006.
  * [19] David Amsallem, Matthew Zahr, Youngsoo Choi, and Charbel Farhat.  Design optimization using hyper-reduced-order models.  Structural and Multidisciplinary Optimization, 51(4):919–940, 2015\.
  * [20] Youngsoo Choi, Gabriele Boncoraglio, Spenser Anderson, David Amsallem, and Charbel Farhat.  Gradient-based constrained optimization using a database of linear reduced-order models.  Journal of Computational Physics, 423:109787, 2020.
  * [21] Youngsoo Choi, Geoffrey Oxberry, Daniel White, and Trenton Kirchdoerfer.  Accelerating design optimization using reduced order models.  arXiv preprint arXiv:1909.11320, 2019.
  * [22] Sean McBane and Youngsoo Choi.  Component-wise reduced order model lattice-type structure design.  arXiv preprint arXiv:2010.10770, 2020.
  * [23] Mohamadreza Ghasemi and Eduardo Gildin.  Localized model reduction in porous media flow.  IFAC-PapersOnLine, 48(6):242–247, 2015.
  * [24] Rui Jiang and Louis J Durlofsky.  Implementation and detailed assessment of a GNAT reduced-order model for subsurface flow simulation.  Journal of Computational Physics, 379:192–213, 2019.
  * [25] Yanfang Yang, Mohammadreza Ghasemi, Eduardo Gildin, Yalchin Efendiev, Victor Calo, et al.  Fast multiscale reservoir simulations with pod-deim model reduction.  SPE Journal, 21(06):2–141, 2016.
  * [26] Min Wang, Siu Wun Cheung, Eric T. Chung, Maria Vasilyeva, and Yuhe Wang.  Generalized multiscale multicontinuum model for fractured vuggy carbonate reservoirs.  Journal of Computational and Applied Mathematics, 366:112370, 2020\.
  * [27] Huanhuan Yang and Alessandro Veneziani.  Efficient estimation of cardiac conductivities via POD-DEIM model order reduction.  Applied Numerical Mathematics, 115:180–199, 2017.
  * [28] Hongfei Fu, Hong Wang, and Zhu Wang.  POD/DEIM reduced-order modeling of time-fractional partial differential equations with applications in parameter identification.  Journal of Scientific Computing, 74(1):220–243, 2018.
  * [29] Pengfei Zhao, Cai Liu, and Xuan Feng.  POD-DEIM based model order reduction for the spherical shallow water equations with Turkel-Zwas finite difference discretization.  Journal of Applied Mathematics, 2014, 2014.
  * [30] R Ştefănescu and Ionel Michael Navon.  POD/DEIM nonlinear model order reduction of an ADI implicit shallow water equations model.  Journal of Computational Physics, 237:95–114, 2013.
  * [31] Youngsoo Choi, Peter Brown, Bill Arrighi, Roberti Anderson, and Kevin Huynh.  Space-time reduced order model for large-scale linear dynamical systems with application to Boltzmann transport problems.  Journal of Computational Physics, 424:109845, 2021.
  * [32] M Mordhorst, Timm Strecker, D Wirtz, Thomas Heidlauf, and Oliver Röhrle.  POD-DEIM reduction of computational EMG models.  Journal of Computational Science, 19:86–96, 2017.
  * [33] Gabriel Dimitriu, Ionel M Navon, and Răzvan Ştefănescu.  Application of POD-DEIM approach for dimension reduction of a diffusive predator-prey system with allee effect.  In International conference on large-scale scientific computing, pages 373–381. Springer, 2013.
  * [34] Harbir Antil, Matthias Heinkenschloss, Ronald HW Hoppe, Christopher Linsenmann, and Achim Wixforth.  Reduced order modeling based shape optimization of surface acoustic wave driven microfluidic biochips.  Mathematics and Computers in Simulation, 82(10):1986–2003, 2012\.
  * [35] Ming-C Cheng.  A reduced-order representation of the Schrödinger equation.  AIP Advances, 6(9):095121, 2016.
  * [36] Serkan Gugercin and Athanasios C Antoulas.  A survey of model reduction by balanced truncation and some new results.  International Journal of Control, 77(8):748–766, 2004.
  * [37] Peter Benner, Serkan Gugercin, and Karen Willcox.  A survey of projection-based model reduction methods for parametric dynamical systems.  SIAM review, 57(4):483–531, 2015.
  * [38] Changhong Mou, Birgul Koc, Omer San, Leo Rebholz, and Traian Iliescu.  Data-driven variational multiscale reduced order models.  arXiv preprint arXiv:2002.06457, 2020.
  * [39] Eric J Parish and Karthik Duraisamy.  Non-Markovian closure models for large eddy simulations using the Mori-Zwanzig formalism.  Physical Review Fluids, 2(1):014604, 2017.
  * [40] Mahmoud Gadalla, Marta Cianferra, Marco Tezzele, Giovanni Stabile, Andrea Mola, and Gianluigi Rozza.  On the comparison of LES data-driven reduced order approaches for hydroacoustic analysis.  Computers & Fluids, page 104819, 2020.
  * [41] Michel Bergmann, C-H Bruneau, and Angelo Iollo.  Enablers for robust POD models.  Journal of Computational Physics, 228(2):516–538, 2009.
  * [42] Jan Östh, Bernd R Noack, Siniša Krajnović, Diogo Barros, and Jacques Borée.  On the need for a nonlinear subscale turbulence term in POD models as exemplified for a high-Reynolds-number flow over an Ahmed body.  Journal of Fluid Mechanics, 747:518–544, 2014.
  * [43] Joan Baiges, Ramon Codina, and Sergio Idelsohn.  Reduced-order subscales for POD models.  Computer Methods in Applied Mechanics and Engineering, 291:173–196, 2015.
  * [44] Omer San and Romit Maulik.  Extreme learning machine for reduced order modeling of turbulent geophysical flows.  Physical Review E, 97(4):042322, 2018.
  * [45] Hannah Lu and Daniel M Tartakovsky.  Lagrangian dynamic mode decomposition for construction of reduced-order models of advection-dominated phenomena.  Journal of Computational Physics, 407:109229, 2020.
  * [46] Rémi Abgrall, David Amsallem, and Roxana Crisovan.  Robust model reduction by ${L}^{1}$-norm minimization and approximation via dictionaries: application to nonlinear hyperbolic problems.  Advanced Modeling and Simulation in Engineering Sciences, 3(1):1–16, 2016.
  * [47] Kevin Carlberg.  Adaptive $h$-refinement for reduced-order models.  International Journal for Numerical Methods in Engineering, 102(5):1192–1210, 2015.
  * [48] Eric J Parish and Kevin T Carlberg.  Windowed least-squares model reduction for dynamical systems.  arXiv preprint arXiv:1910.11388, 2019.
  * [49] Yukiko S Shimizu and Eric J Parish.  Windowed space-time least-squares Petrov-Galerkin method for nonlinear model order reduction.  arXiv preprint arXiv:2012.06073, 2020.
  * [50] Benjamin Peherstorfer.  Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling.  arXiv preprint arXiv:1812.02094, 2018.
  * [51] PG Constantine and G Iaccarino.  Reduced order models for parameterized hyperbolic conservations laws with shock reconstruction.  Center for Turbulence Research Annual Brief, 2012.
  * [52] Tommaso Taddei and Lei Zhang.  Space-time registration-based model reduction of parameterized one-dimensional hyperbolic PDEs.  arXiv preprint arXiv:2004.06693, 2020.
  * [53] Julius Reiss, Philipp Schulze, Jörn Sesterhenn, and Volker Mehrmann.  The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena.  SIAM Journal on Scientific Computing, 40(3):A1322–A1344, 2018.
  * [54] Donsub Rim, Scott Moe, and Randall J LeVeque.  Transport reversal for model reduction of hyperbolic partial differential equations.  SIAM/ASA Journal on Uncertainty Quantification, 6(1):118–150, 2018\.
  * [55] G Welper.  Transformed snapshot interpolation with high resolution transforms.  SIAM Journal on Scientific Computing, 42(4):A2037–A2061, 2020.
  * [56] Michael Kirby and Dieter Armbruster.  Reconstructing phase space from PDE simulations.  Zeitschrift für angewandte Mathematik und Physik ZAMP, 43(6):999–1022, 1992.
  * [57] Kookjin Lee and Kevin T Carlberg.  Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.  Journal of Computational Physics, 404:108973, 2020.
  * [58] Kookjin Lee and Kevin Carlberg.  Deep conservation: A latent dynamics model for exact satisfaction of physical conservation laws.  arXiv preprint arXiv:1909.09754, 2019.
  * [59] Youngkyu Kim, Youngsoo Choi, David Widemann, and Tarek Zohdi.  A fast and accurate physics-informed neural network reduced order model with shallow masked autoencoder.  arXiv preprint arXiv:2009.11990, 2020.
  * [60] Youngkyu Kim, Youngsoo Choi, David Widemann, and Tarek Zohdi.  Efficient nonlinear manifold reduced order model.  arXiv preprint arXiv:2011.07727, 2020.
  * [61] Donsub Rim, Benjamin Peherstorfer, and Kyle T Mandli.  Manifold approximations via transported subspaces: Model reduction for transport-dominated problems.  arXiv preprint arXiv:1912.13024, 2019.
  * [62] Donsub Rim, Luca Venturi, Joan Bruna, and Benjamin Peherstorfer.  Depth separation for reduced deep networks in nonlinear model reduction: Distilling shock waves in nonlinear hyperbolic problems.  arXiv preprint arXiv:2007.13977, 2020.
  * [63] C. S. Peskin.  The immersed boundary method.  Acta Numerica, 11:479–517, 2002.
  * [64] D. Boffi and L. Gastaldi.  Discrete models for fluid-structure interactions: the finite element immersed boundary method.  Discrete and Continuous Dynamical Systems - Series S, 9:89–107, 2016.
  * [65] Siu Wun Cheung, Eric T. Chung, and Hyea Hyun Kim.  A mass conservative scheme for fluid–structure interaction problems by the staggered discontinuous Galerkin method.  Journal of Scientific Computing, 74:1423–1456, 2018.
  * [66] P. K. Yeung.  Lagrangian investigations of turbulence.  Annual Review of Fluid Mechanics, 34:115–142, 2002.
  * [67] B. Lüthi, A. Tsinober, and W. Kinzelbach.  Lagrangian measurement of vorticity dynamics in turbulent flow.  Journal of Fluid Mechanics, 528:87–118, 2005.
  * [68] C. W. Hirt, A. A. Amsden, and J. L. Cook.  An arbitrary Lagrangian-Eulerian computing method for all flow speeds.  Journal of Computational Physics, 14:227–253, 1974.
  * [69] David J. Benson.  An efficient, accurate, simple ALE method for nonlinear finite element programs.  Computer Methods in Applied Mechanics and Engineering, 72:305–350, 1989.
  * [70] J. L. Guermond, B. Popov, L. Saavedra, and Y. Yang.  Invariant domains preserving arbitrary Lagrangian Eulerian approximation of hyperbolic systems with continuous finite elements.  SIAM Journal on Scientific Computing, 39:A385–A414, 2017.
  * [71] Veselin A Dobrev, Tzanio V Kolev, and Robert N Rieben.  High-order curvilinear finite element methods for Lagrangian hydrodynamics.  SIAM Journal on Scientific Computing, 34(5):B606–B641, 2012.
  * [72] F.H. Harlow and A.A. Amsfen.  Fluid Dynamics: A LASL Monograph.  Tech. rep. LA-4700. Los Alamos Scientific Laboratory, 1971.
  * [73] Saifon Chaturantabut and Danny C Sorensen.  Nonlinear model reduction via discrete empirical interpolation.  SIAM Journal on Scientific Computing, 32(5):2737–2764, 2010.
  * [74] Gal Berkooz, Philip Holmes, and John L Lumley.  The proper orthogonal decomposition in the analysis of turbulent flows.  Annual review of fluid mechanics, 25(1):539–575, 1993.
  * [75] Harold Hotelling.  Analysis of a complex of statistical variables into principal components.  Journal of educational psychology, 24(6):417, 1933.
  * [76] Michel Loeve.  Probability Theory.  D. Van Nostrand, New York, 1955.
  * [77] Michael Hinze and Stefan Volkwein.  Proper orthogonal decomposition surrogate models for nonlinear dynamical systems: Error estimates and suboptimal control.  In Dimension reduction of large-scale systems, pages 261–306. Springer, 2005.
  * [78] Karl Kunisch and Stefan Volkwein.  Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics.  SIAM Journal on Numerical analysis, 40(2):492–515, 2002.
  * [79] Zlatko Drmac and Serkan Gugercin.  A new selection operator for the discrete empirical interpolation method—improved a priori error bound and extensions.  SIAM Journal on Scientific Computing, 38(2):A631–A648, 2016.
  * [80] Zlatko Drmac and Arvind Krishna Saibaba.  The discrete empirical interpolation method: Canonical structure and formulation in weighted inner product spaces.  SIAM Journal on Matrix Analysis and Applications, 39(3):1152–1180, 2018.
  * [81] Kevin Carlberg, Charbel Farhat, Julien Cortial, and David Amsallem.  The GNAT method for nonlinear model reduction: Effective implementation and application to computational fluid dynamics and turbulent flows.  Journal of Computational Physics, 242:623–647, 2013.
  * [82] Kevin Carlberg, Charbel Bou‐Mosleh, and Charbel Farhat.  Efficient non-linear model reduction via a least–squares Petrov–Galerkin projection and compressive tensor approximations.  International Journal for Numerical Methods in Engineering, 86:155–181, 2011.
  * [83] Ming Gu and Stanley C Eisenstat.  Efficient algorithms for computing a strong rank-revealing QR factorization.  SIAM Journal on Scientific Computing, 17(4):848–869, 1996.
  * [84] Sean Ahern, Eric Brugger, Brad Whitlock, Jeremy S Meredith, Kathleen Biagas, Mark C Miller, and Hank Childs.  Visit: Experiences with sustainable software.  arXiv preprint arXiv:1309.1796, 2013.
  * [85] Philip M. Gresho and Stevens T. Chan.  On the theory of semi‐implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. part 2: Implementation.  International Journal for Numerical Methods in Fluids, 11:621–659, 1990.
  * [86] L.I. Sedov.  Similarity and Dimensional Methods in Mechanics, Tenth Edition.  Taylor & Francis, 1993.
  * [87] Taylor Geoffrey Ingram and Green Albert Edward.  Mechanism of the production of small eddies from large ones.  Proceedings of the Royal Society of London A, 158(895):499–521, 1937.
  * [88] S. Galera, P.-H. Maire, and J. Breil.  A two-dimensional unstructured cell-centered multi-material ALE scheme using VOF interface reconstruction.  Journal of Computational Physics, 229:621–659, 2010.