Symmetry-preserving finite-difference discretizations of arbitrary order on structured curvilinear staggered grids

TL;DR

This paper introduces a new finite-difference discretization method that preserves symmetry and conservation properties on structured curvilinear grids, supporting arbitrary order and complex operators, demonstrated through shallow-water equations.

Contribution

It presents a novel symmetry-preserving finite-difference scheme compatible with arbitrary order, various grid types, and complex operators, expanding the applicability of mimetic discretizations.

Findings

Exact conservation of physical quantities demonstrated

Convergence matches theoretical order

Applicable to complex operators like chain rules and nonlinear advection

Abstract

Symmetry-preserving (mimetic) discretization aims to preserve certain properties of a continuous differential operator in its discrete counterpart. For these discretizations, stability and (discrete) conservation of mass, momentum and energy are proven in the same way as for the original continuous model. This paper presents a new finite-difference symmetry-preserving space discretization. Boundary conditions and time integration are not addressed. The novelty is that it combines arbitrary order of convergence, orthogonal and non-orthogonal structured curvilinear staggered meshes, and the applicability to a wide variety of continuous operators, involving chain rules and nonlinear advection, as illustrated by the shallow-water equations. Experiments show exact conservation and convergence corresponding to expected order.

Figures

Tables

Table 1. Table 1: Propagation speeds and Riemann invariants used to construct exact one-dimensional solutions for the three models [ 20 ] . The one-dimensional velocity is given by v 𝑣 v .

	Linear-wave eq.	Compressible-wave eq.	Shallow-water eq.
Propagation speed $V$
forward $V_{+}$	$c$	$\frac{1}{2}\left(v+\sqrt{v^2+4c^2}\right)$	$v+\sqrt{g\rho }$
backward $V_{-}$	$-c$	$\frac{1}{2}\left(v-\sqrt{v^2+4c^2}\right)$	$v-\sqrt{g\rho }$
Riemann invariant $F$
forward $F_{+}$	${\rho }_{0}v+c\rho$	$\ln \left(\rho V_{+}\right)-\frac{v{V}_{-}}{2c^2}$	$v+2\sqrt{g\rho }$
backward $F_{-}$	${\rho }_{0}v-c\rho$	$\ln \left(\frac{\rho }{V_{+}}\right)-\frac{v{V}_{+}}{2c^2}$	$v-2\sqrt{g\rho }$

Table 2. Table 2: Several properties of the operators used in the models. The left null space properties are used to prove conservation of mass and momentum. For energy conservation, all properties are required.

Left null space properties
Continuous version	Discrete equivalent
$\int \nabla \cdot \overrightarrow{f}\text{d}V={\oint }_{\delta V}\overrightarrow{f}\cdot \overrightarrow{n}\text{d}S=0$	${\mathrm{𝖣𝖨𝖵}}^{\ast }\mathrm{𝚌𝟷}=0$
$\int \nabla f\text{d}V={\oint }_{\delta V}f\overrightarrow{n}\text{d}S=\overrightarrow{0}$	${\mathrm{𝖦𝖱𝖠𝖣}}^{\ast }\mathrm{𝚌𝟷𝟶𝟶}={\mathrm{𝖦𝖱𝖠𝖣}}^{\ast }\mathrm{𝚌𝟶𝟷𝟶}={\mathrm{𝖦𝖱𝖠𝖣}}^{\ast }\mathrm{𝚌𝟶𝟶𝟷}=0$
${\int }_V\nabla \cdot \left(\rho \overrightarrow{v}\otimes \overrightarrow{f}\right)\text{d}V={\oint }_{\delta V}\rho \overrightarrow{f}\overrightarrow{v}\cdot \overrightarrow{n}\text{d}S=\overrightarrow{0}$	${\mathrm{𝖠𝖣𝖵𝖤𝖢}}^{\ast }\mathrm{𝚌𝟷}=0$
${\int }_V\nabla \cdot \rho \overrightarrow{f}\text{d}V={\oint }_{\delta V}\rho \overrightarrow{f}\cdot \overrightarrow{n}\text{d}S=0$	$\mathrm{𝖣𝖨𝖵}{\widetilde{𝗋}}^{\ast }\mathrm{𝚌𝟷}=0$
Zero derivative for constants (right null space properties)
Continuous version	Discrete equivalent
$\nabla \cdot (1,0,0)=\nabla \cdot (0,1,0)=\nabla \cdot (0,0,1)=0$	$\mathrm{𝖣𝖨𝖵}\mathrm{𝚌𝟷𝟶𝟶}=\mathrm{𝖣𝖨𝖵}\mathrm{𝚌𝟶𝟷𝟶}=\mathrm{𝖣𝖨𝖵}\mathrm{𝚌𝟶𝟶𝟷}=0$
$\nabla 1=\overrightarrow{0}$	$\mathrm{𝖦𝖱𝖠𝖣}\mathrm{𝚌𝟷}=0$
Chain-rule properties
Continuous version	Discrete equivalent
$R\left(p\right)\nabla S\left(p\right)=\nabla Q\left(p\right)$	$\widetilde{𝗋}\mathrm{𝖦𝖱𝖠𝖣}S\left(𝚙\right)=\mathrm{𝖦𝖱𝖠𝖣}Q\left(𝚙\right)$
$R\left(p\right)\nabla Q\left(p\right)=\nabla p$	$\mathrm{𝗋𝖦𝖱𝖠𝖣}Q\left(𝚙\right)=\mathrm{𝖦𝖱𝖠𝖣}𝚙$
with $Q\left(p\right):={\int }^p\frac{1}{R\left(q\right)}\text{d}q$ and $S\left(p\right):={\int }^p\frac{1}{R^2\left(q\right)}\text{d}q$
Symmetry properties
Continuous version	Discrete equivalent
${\int }_{V}f\nabla \cdot \overrightarrow{g}\text{d}V+{\int }_V\overrightarrow{g}\cdot \nabla f\text{d}V={\oint }_{\delta V}f\overrightarrow{g}\cdot \overrightarrow{n}\text{d}S=0$	$\mathrm{𝖣𝖨𝖵}+{\mathrm{𝖦𝖱𝖠𝖣}}^{\ast }=0$
${\int }_{V}f\nabla \cdot \rho \overrightarrow{g}\text{d}V+{\int }_V\overrightarrow{g}\cdot \rho \nabla f\text{d}V={\oint }_{\delta V}\rho f\overrightarrow{g}\cdot \overrightarrow{n}\text{d}S=0$	$\mathrm{𝖣𝖨𝖵}\widetilde{𝗋}+\widetilde{𝗋}{\mathrm{𝖦𝖱𝖠𝖣}}^{\ast }=0$
${\int }_{V}f\nabla \cdot \rho \overrightarrow{g}\text{d}V+{\int }_V\overrightarrow{g}\cdot \rho \nabla f\text{d}V={\oint }_{\delta V}\rho f\overrightarrow{g}\cdot \overrightarrow{n}\text{d}S=0$	$\mathrm{𝖣𝖨𝖵𝗋}+{\mathrm{𝗋𝖦𝖱𝖠𝖣}}^{\ast }=0$
${\int }_V\overrightarrow{g}\cdot \nabla \cdot \rho \overrightarrow{v}\otimes f\text{d}V+{\int }_V\overrightarrow{f}\cdot \nabla \cdot \rho \overrightarrow{v}\otimes g\text{d}V$	$\mathrm{𝖠𝖣𝖵𝖤𝖢}+{\mathrm{𝖠𝖣𝖵𝖤𝖢}}^{\ast }=\text{diag}\left({\mathrm{𝖨𝗇𝗍𝖾𝗋𝗉}}_{v\leftarrow c}\mathrm{𝖣𝖨𝖵𝗋}𝚟\right)$
$={\int }_V\overrightarrow{f}\cdot \overrightarrow{g}\nabla \cdot \rho \overrightarrow{v}\text{d}V+{\oint }_{\delta V}\rho \left(\overrightarrow{f}\cdot \overrightarrow{g}\right)\left(\overrightarrow{v}\cdot \overrightarrow{n}\right)\text{d}S={\int }_V\overrightarrow{f}\cdot \overrightarrow{g}\nabla \cdot \rho \overrightarrow{v}\text{d}V$

Table 3. Table 3: Expressions for the conserved quantities mass, momentum and energy in the three models. Here, ρ = R ​ ( p ) 𝜌 𝑅 𝑝 \rho=R(p) for the compressible-wave model and isentropic Euler model.

	Linear-wave	Compressible-wave	Isentropic Euler
Mass
continuous $M$	${\int }_V\rho \text{d}V$	${\int }_V\rho \text{d}V$	${\int }_V\rho \text{d}V$
discrete M	${⟨\mathrm{𝚌𝟷},\mathrm{𝚛𝚑𝚘}⟩}_c$	${⟨\mathrm{𝚌𝟷},\mathrm{𝚛𝚑𝚘}⟩}_c$	${⟨\mathrm{𝚌𝟷},\mathrm{𝚛𝚑𝚘}⟩}_c$
Momentum
continuous $\overrightarrow{M}$	${\rho }_0{\int }_V\overrightarrow{v}\text{d}V$	${\rho }_0{\int }_V\overrightarrow{v}\text{d}V$	${\int }_V\rho \overrightarrow{v}\text{d}V$
discrete $\overrightarrow{𝙼}$	${\rho }_0\left(\begin{array}{c}{⟨\mathrm{𝚌𝟷𝟶𝟶},𝚟⟩}_v\\ {⟨\mathrm{𝚌𝟶𝟷𝟶},𝚟⟩}_v\\ {⟨\mathrm{𝚌𝟶𝟶𝟷},𝚟⟩}_v\end{array}\right)$	${\rho }_0\left(\begin{array}{c}{⟨\mathrm{𝚌𝟷𝟶𝟶},𝚟⟩}_v\\ {⟨\mathrm{𝚌𝟶𝟷𝟶},𝚟⟩}_v\\ {⟨\mathrm{𝚌𝟶𝟶𝟷},𝚟⟩}_v\end{array}\right)$	$\left(\begin{array}{c}{⟨\mathrm{𝚌𝟷𝟶𝟶},\mathrm{𝚛𝚟}⟩}_v\\ {⟨\mathrm{𝚌𝟶𝟷𝟶},\mathrm{𝚛𝚟}⟩}_v\\ {⟨\mathrm{𝚌𝟶𝟶𝟷},\mathrm{𝚛𝚟}⟩}_v\end{array}\right)$
Energy
continuous $E$	${\int }_V\left(e_{int}+e_{kin}\right)\text{d}V$	${\int }_V\left(e_{int}+e_{kin}\right)\text{d}V$	${\int }_V\left(e_{int}+e_{kin}\right)\text{d}V$
$e_{int}$	$\frac{c^2{\rho }^2}{2{\rho }_0}$	${\rho }_0{\int }^p\frac{R\left(p\right)-R\left(q\right)}{R^2\left(q\right)}\text{d}q$	${\int }^p\frac{R\left(p\right)-R\left(q\right)}{R\left(q\right)}\text{d}q$
$e_{kin}$	$\frac{{\rho }_0}{2}{\left|\overrightarrow{v}\right|}^2$	$\frac{{\rho }_0}{2}{\left|\overrightarrow{v}\right|}^2$	$\frac{\rho }{2}{\left|\overrightarrow{v}\right|}^2$
discrete $𝙴$	${⟨\mathrm{𝚌𝟷},e_{int}\left(𝚙\right)⟩}_c$	${⟨\mathrm{𝚌𝟷},e_{int}\left(𝚙\right)⟩}_c$	${⟨\mathrm{𝚌𝟷},e_{int}\left(𝚙\right)⟩}_c$
	$+\frac{{\rho }_0}{2}{⟨𝚟,𝚟⟩}_v$	$+\frac{{\rho }_0}{2}{⟨𝚟,𝚟⟩}_v$	$+\frac{1}{2}{⟨\mathrm{𝚛𝚟},𝚟⟩}_v$

Table 4. Table 4: Relative errors (in 2-norm) for pressure in the linear-wave equations (Galerkin-type approach) at time T = 10 𝑇 10 T=10 , for various interpolation orders, using reltol equal to 1e-11 in lsode . For 160 × 160 160 160 160\times 160 grid points on the uniform grid and 8th order, reltol= 1e-13 is used.

Uniform grid
	order=2		order=4		order=6		order=8
	error	order	error	order	error	order	error	order
$10\times 10$	1.12e+00		7.64e-01		1.68e-01		4.71e-01
$20\times 20$	1.43e+00	-0.35	4.10e-01	0.90	1.43e-01	0.23	4.49e-02	3.39
$40\times 40$	4.95e-01	1.53	7.85e-02	2.38	6.86e-03	4.38	7.82e-04	5.84
$80\times 80$	4.02e-01	0.30	6.34e-03	3.63	1.22e-04	5.82	3.55e-06	7.78
$160\times 160$	1.46e-01	1.46	4.01e-04	3.98	1.93e-06	5.98	1.43e-08	7.96

Table 5. Table 5: Relative losses for linear-wave equations (Galerkin-type approach) at time T = 10 𝑇 10 T=10 , for varying tolerances reltol in the time-integration method lsode . Fourth-order discretization on a 20 × 20 20 20 20\times 20 grid.

	Uniform grid			Curvilinear grid
reltol	mass loss	mom. loss	energy loss	mass loss	mom. loss	energy loss
1e-07	5.29e-14	1.20e-13	9.94e-07	4.14e-15	6.33e-15	1.60e-06
1e-08	5.27e-14	1.20e-13	6.97e-08	2.83e-15	4.78e-15	1.56e-08
1e-09	5.16e-14	1.21e-13	6.39e-08	3.92e-15	6.57e-15	3.25e-09
1e-10	5.29e-14	1.20e-13	1.15e-09	3.92e-15	5.58e-15	1.50e-11
1e-11	5.36e-14	1.17e-13	5.54e-11	2.61e-15	7.73e-15	1.87e-11
1e-12	5.23e-14	1.18e-13	4.74e-11	1.96e-15	4.43e-15	5.84e-12

Table 6. Table 6: Relative errors (in 2-norm) for pressure in the linear-wave equations (finite-difference approach) at time T = 10 𝑇 10 T=10 , for various (same) orders for derivatives and interpolations, using reltol equal to 1e-11 in lsode . For 160 × 160 160 160 160\times 160 grid points and 8th order, reltol= 1e-13 is used.

Uniform grid
	order=2		order=4		order=6		order=8
	error	order	error	order	error	order	error	order
$10\times 10$	9.29e-01		4.52e-01		1.33e-01		1.62e-01
$20\times 20$	6.22e-01	0.58	8.90e-02	2.34	1.88e-02	2.82	5.54e-03	4.87
$40\times 40$	2.58e-01	1.27	7.68e-03	3.53	4.14e-04	5.51	3.59e-05	7.27
$80\times 80$	7.67e-02	1.75	4.91e-04	3.97	6.77e-06	5.94	1.49e-07	7.92
$160\times 160$	1.96e-02	1.97	3.08e-05	4.00	9.79e-08	6.11	6.00e-10	7.96

Table 7. Table 7: Relative losses for linear-wave equations (finite-difference approach) at time T = 10 𝑇 10 T=10 , for varying tolerances reltol in the time-integration method lsode . Fourth-order discretization on a 20 × 20 20 20 20\times 20 grid.

	Uniform grid			Curvilinear grid
reltol	mass loss	mom. loss	energy loss	mass loss	mom. loss	energy loss
1e-07	2.18e-16	2.39e-15	1.94e-06	2.18e-16	1.36e-14	3.29e-08
1e-08	1.09e-15	1.57e-15	1.22e-06	2.18e-16	1.43e-14	4.61e-09
1e-09	4.35e-16	1.12e-15	1.78e-09	4.35e-16	1.15e-14	4.33e-10
1e-10	1.09e-15	1.66e-15	2.98e-09	1.52e-15	1.10e-14	7.85e-11
1e-11	0.00e+00	9.86e-16	1.35e-09	6.53e-16	1.52e-14	1.28e-11
1e-12	2.18e-16	1.66e-15	3.16e-11	1.31e-15	1.27e-14	8.80e-13

Table 8. Table 8: Relative errors (in 2-norm) for the densities in the compressible-wave equations, for various (same) orders for derivatives and interpolations, using reltol equal to 1e-11 in lsode .

Uniform grid, $t=t_N/2$
	order=2		order=4		order=6		order=8
	error	order	error	order	error	order	error	order
$10\times 10$	3.03e-01		6.65e-02		3.52e-02		2.48e-02
$20\times 20$	1.04e-01	1.54	1.24e-02	2.42	5.29e-03	2.74	3.49e-03	2.83
$40\times 40$	3.05e-02	1.77	1.32e-03	3.24	3.13e-04	4.08	1.53e-04	4.51
$80\times 80$	7.98e-03	1.93	9.50e-05	3.79	8.55e-06	5.19	2.29e-06	6.06
$160\times 160$	2.01e-03	1.99	6.10e-06	3.96	1.59e-07	5.75	1.56e-08	7.20

Table 9. Table 9: Relative losses for the compressible-wave equations, for varying tolerances reltol in the time-integration method lsode , at the time when the solution becomes discontinuous. Fourth-order discretization on a 20 × 20 20 20 20\times 20 grid.

	Uniform grid			Curvilinear grid
reltol	mass loss	mom. loss	energy loss	mass loss	mom. loss	energy loss
1e-05	6.82e-16	1.85e-14	6.53e-09	1.02e-15	1.67e-13	6.41e-10
1e-06	3.41e-16	1.94e-14	7.99e-10	1.19e-15	1.69e-13	1.57e-10
1e-07	1.02e-15	2.14e-14	4.18e-11	8.52e-16	1.64e-13	6.64e-12
1e-08	3.41e-16	2.39e-14	1.68e-13	0.00e+00	1.68e-13	6.49e-13
1e-09	6.82e-16	1.75e-14	3.30e-13	8.52e-16	1.69e-13	5.20e-14
1e-10	1.70e-16	1.81e-14	3.33e-14	5.11e-16	1.66e-13	1.09e-15
1e-11	3.41e-16	1.97e-14	2.90e-15	2.56e-15	1.66e-13	1.45e-15
1e-12	1.02e-15	2.10e-14	1.81e-16	1.88e-15	1.66e-13	1.81e-16

Table 10. Table 10: Relative errors (in 2-norm) for the densities in the shallow-water equations at time t = t N / 2 𝑡 subscript 𝑡 𝑁 2 t=t_{N}/2 , for various (same) orders for derivatives and interpolations, using reltol equal to 1e-11 in lsode .

Uniform grid
	order=2		order=4		order=6		order=8
	error	order	error	order	error	order	error	order
$10\times 10$	7.07e-02		2.23e-02		1.33e-02		1.10e-02
$20\times 20$	2.62e-02	1.43	4.60e-03	2.28	2.31e-03	2.53	1.65e-03	2.74
$40\times 40$	8.16e-03	1.68	6.00e-04	2.94	1.89e-04	3.61	1.05e-04	3.97
$80\times 80$	2.21e-03	1.89	4.81e-05	3.64	6.40e-06	4.89	2.03e-06	5.69
$160\times 160$	5.61e-04	1.98	3.15e-06	3.93	1.24e-07	5.69	1.46e-08	7.12

Table 11. Table 11: Relative losses for the shallow-water equations, for varying tolerances reltol in the time-integration method lsode , at the time when the solution becomes discontinuous. Fourth-order discretization on a 20 × 20 20 20 20\times 20 grid.

	Uniform grid			Curvilinear grid
reltol	mass loss	mom. loss	energy loss	mass loss	mom. loss	energy loss
1e-05	8.07e-16	1.49e-15	1.16e-05	2.02e-16	5.40e-14	3.35e-06
1e-06	1.41e-15	9.39e-16	8.07e-07	4.04e-16	5.55e-14	4.85e-08
1e-07	0.00e+00	7.04e-16	2.07e-07	8.07e-16	5.56e-14	1.59e-09
1e-08	0.00e+00	1.37e-15	5.05e-09	2.02e-16	5.55e-14	7.83e-10
1e-09	1.01e-15	5.25e-16	1.36e-10	4.04e-16	5.44e-14	1.05e-11
1e-10	0.00e+00	1.11e-15	3.91e-11	0.00e+00	5.68e-14	3.12e-12
1e-11	6.06e-16	9.39e-16	2.12e-12	1.01e-15	5.48e-14	5.79e-13
1e-12	4.04e-16	8.30e-16	5.61e-13	6.06e-16	5.53e-14	4.02e-14

Equations

$$\frac{\partial \rho }{\partial t}+{\rho }_0\nabla \cdot \vec{v}=0,$$

$$\frac{\text{d}\mathtt{rho}}{\text{d}t}+{\rho }_0\mathsf{DIV}\mathtt{v}=0,$$

$$\frac{\partial e}{\partial t}+\nabla \cdot p\vec{v}=0,$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \rho \vec{v}=0,$$

$$\frac{\text{d}\mathtt{rho}}{\text{d}t}+\mathsf{DIV}\tilde{\mathsf{r}}\mathtt{v}=0,\frac{\text{d}\mathtt{v}}{\text{d}t}+\mathsf{GRAD}Q(\mathtt{p})=0,\mathtt{rho}=R(\mathtt{p}),$$

$$\frac{\partial e}{\partial t}+\nabla \cdot \left(e_{int}+{\rho }_0{\int }^p\frac{1}{R(q)}\text{mbox}dq\right)\vec{v}=0.$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \rho \vec{v}=0,$$

$$\frac{\text{d}\mathtt{rho}}{\text{d}t}+\mathsf{DIVr}\mathtt{v}=0,\frac{\text{d}\mathtt{rv}}{\text{d}t}+\mathsf{ADVEC}\mathtt{v}+\mathsf{GRAD}\mathtt{p}=0,\mathtt{rho}=R(\mathtt{p}).$$

$$\frac{\partial e}{\partial t}+\nabla \cdot (e+p)\vec{v}=0.$$

$${\vec{\xi }}_{c,i+m_{x}j+m_x{m}_{y}k}=(i\Delta \xi ,j\Delta \eta ,k\Delta \zeta {)}^{\top },$$

$${\vec{x}}_{c,i}=\vec{X}({\vec{\xi }}_{c,i}).$$

$$\displaystyle\vec{x}_{e,i}:=\vec{X}\left(\vec{\xi}_{e,i}\right):=\vec{X}\left(\vec{\xi}_{c,i}+\frac{1}{2}{\mbox{\scriptsize$\left(\!\!\begin{array}[]{l}\Delta\xi\\ 0\\ 0\end{array}\!\!\right)$}}\right),$$

$$\displaystyle\vec{x}_{t,i}:=\vec{X}\left(\vec{\xi}_{t,i}\right):=\vec{X}\left(\vec{\xi}_{c,i}+\frac{1}{2}{\mbox{\scriptsize$\left(\!\!\begin{array}[]{l}0\\ 0\\ \Delta\zeta\end{array}\!\!\right)$}}\right),$$

$$\vec{v}(\vec{x})=v_x(\vec{x}){\vec{r}}_x(\vec{x})+v_y(\vec{x}){\vec{r}}_y(\vec{x})+v_z(\vec{x}){\vec{r}}_z(\vec{x}),$$

$$v_x(\vec{x})={\vec{r}}_x(\vec{x})\cdot \vec{v}(\vec{x}),$$

$$\displaystyle\left(\begin{array}[]{lll}\displaystyle\frac{\partial\vec{X}}{\partial\xi}&\displaystyle\frac{\partial\vec{X}}{\partial\eta}&\displaystyle\frac{\partial\vec{X}}{\partial\zeta}\end{array}\right)={\bf P}~{}\mbox{diag}({\bf h})~{}{\bf Q}.$$

$$[{\vec{r}}_x,{\vec{r}}_y,{\vec{r}}_z]:=\mathbf{P}\mathbf{Q}.$$

$${\mathtt{vx}}_i={\vec{r}}_x({\vec{x}}_{e,i})\cdot \vec{v}({\vec{x}}_{e,i}),{\mathtt{vy}}_i={\vec{r}}_y({\vec{x}}_{n,i})\cdot \vec{v}({\vec{x}}_{n,i}),{\mathtt{vz}}_i={\vec{r}}_z({\vec{x}}_{t,i})\cdot \vec{v}({\vec{x}}_{t,i}).$$

$$\displaystyle\vec{r}_{x}(\vec{x}_{c,i})=\left(\begin{array}[]{c}{\tt rxx\_at\_c}_{i}\\ {\tt rxy\_at\_c}_{i}\\ {\tt rxz\_at\_c}_{i}\end{array}\right),$$

$$\vec{v}(\vec{x})={\vec{c}}_{(1,0,0)}(\vec{x})=(1,0,0),\vec{v}(\vec{x})={\vec{c}}_{(0,1,0)}(\vec{x})=(0,1,0),\vec{v}(\vec{x})={\vec{c}}_{(0,0,1)}(\vec{x})=(0,0,1),$$

$$c100=\left(\begin{array}{c}rxx\_at\_e\\ ryx\_at\_n\\ rzx\_at\_t\end{array}\right),c010=\left(\begin{array}{c}rxy\_at\_e\\ ryy\_at\_n\\ rzy\_at\_t\end{array}\right),c001=\left(\begin{array}{c}rxz\_at\_e\\ ryz\_at\_n\\ rzz\_at\_t\end{array}\right).$$

$$⟨\mathtt{a},\mathtt{b}{⟩}_c:={\mathtt{a}}^{\top }\text{mbox}diag(\mathtt{dVc})\mathtt{b},$$

$$⟨\mathtt{v},\mathtt{w}{⟩}_v$$

$$:={\mathtt{vx}}^{\top }\text{mbox}diag(\mathtt{dVe})\mathtt{wx}+{\mathtt{vy}}^{\top }\text{mbox}diag(\mathtt{dVn})\mathtt{wy}+{\mathtt{vz}}^{\top }\text{mbox}diag(\mathtt{dVt})\mathtt{wz},$$

$$\vec{\mathtt{M}}:={\left(⟨c100,\mathtt{rv}{⟩}_v,⟨c010,\mathtt{rv}{⟩}_v,⟨c001,\mathtt{rv}{⟩}_v\right)}^{\top }.$$

$${\mathsf{A}}^{\ast }=\text{mbox}diag(\mathtt{dVc}{)}^{-1}{\mathsf{A}}^{\top }\text{mbox}diag(\mathtt{dVe}),$$

$$⟨\mathtt{x},\mathsf{A}\mathtt{y}{⟩}_e={\mathtt{x}}^{\top }\text{mbox}diag(\mathtt{dVe})\mathsf{A}\mathtt{y}={\mathtt{x}}^{\top }({\mathsf{A}}^{\ast }{)}^{\top }\text{mbox}diag(\mathtt{dVc})\mathtt{y}=⟨{\mathsf{A}}^{\ast }\mathtt{x},\mathtt{y}{⟩}_c.$$

$$\frac{\partial M}{\partial t}=\frac{\partial }{\partial t}{\int }_V\rho \text{mbox}dV=-{\int }_V\nabla \cdot \rho \vec{v}\text{mbox}dV=-{\oint }_{\delta V}\rho \vec{v}\cdot \vec{n}\text{mbox}dS.$$

$${\mathsf{DIV}}^{\ast }c1=0$$

$$⟺{\mathtt{f}}^{\top }\text{mbox}diag(\mathtt{dVc})\left\{\text{mbox}diag(\mathtt{dVc}{)}^{-1}{\mathsf{DIV}}^{\top }\text{mbox}diag(\mathtt{dVc})\right\}c1=0\forall \mathtt{f}$$

$$⟺{\mathtt{f}}^{\top }{\mathsf{DIV}}^{\top }\text{mbox}diag(\mathtt{dVc})c1=0\forall \mathtt{f}$$

$$⟺{c1}^{\top }\text{mbox}diag(\mathtt{dVc})\mathsf{DIV}\mathtt{f}=0\forall \mathtt{f}$$

$$⟺\text{The discrete integral of }\mathsf{DIV}\mathtt{f}\text{ equals zero}\forall \mathtt{f},$$

Full text

This paper presents a new finite-difference symmetry-preserving space discretization. Boundary conditions and time integration are not addressed. The novelty is that it combines arbitrary order of convergence, orthogonal and non-orthogonal structured curvilinear staggered meshes, and the applicability to a wide variety of continuous operators, involving chain rules and nonlinear advection, as illustrated by the shallow-water equations. Experiments show exact conservation and convergence corresponding to expected order.

Symmetry-preserving discretizations, Mimetic methods, Finite-difference methods, Mass, momentum and energy conservation, Curvilinear staggered grid

1 Introduction and motivation

Computer simulations of a flow phenomenon require the discretization of the flow properties, reducing the number of values needed to represent the flow state from infinite to some large finite number. In the resulting discrete model of the flow phenomenon, the (continuous) differential operators have been replaced by (discrete) difference operators. Unfortunately, not all properties of the differential operators are automatically inherited by their discrete approximations. The chain and product rules needed in the manipulation of nonlinear equations, for example, do not always work in discrete cases. Moreover, symmetry and positiveness may be lost in the discretization process, mass, momentum, and energy may not be conserved, aliasing errors can occur, and duality and self-adjointness of the differential operators may be violated [2, 17, 21].

Symmetry-preserving methods, or mimetic methods, aim to preserve certain properties of a continuous operator in its discrete counterpart [17]. The use of symmetry-preserving discretizations makes it possible to construct discrete models which allow all the manipulations needed to prove stability and (discrete) conservation in the same way they were proven in the original continuous model. In this paper, we will present a new finite-difference symmetry-preserving discretization of arbitrary order on structured curvilinear staggered grids.

There are a variety of symmetry-preserving discretizations available in the literature. In [35], an exhaustive overview is given of different techniques to obtain mass- or energy-conserving methods. Typically, symmetry properties of differential operators are only automatically preserved in central-difference approximations on uniform, rectilinear grids [17]. Finite-volume methods can be used to construct conservative discretizations for mass and momentum, but it is in general not possible to also obtain energy conservation [35].

In [42, 43], a fourth-order symmetry-preserving finite-volume method is constructed using Richardson extrapolation of a second-order symmetry-preserving method [40]. The extension to unstructured collocated meshes is presented in [32], and an application can be seen in [41]. The extension to upwind discretizations was made in [39], and a discretization for the advection operator for curvilinear collocated meshes was found in [16]. In [17], the method is extended to non-uniform structured curvilinear collocated grids by deriving a discrete product rule. Furthermore, a symmetry-preserving method that conserves mass and energy for compressible-flow equations with a state equation is described in [35]. For rectilinear grids, this method works well, but it is challenging to let this method work for unstructured grids [35]. Finally, in [26], a symmetry-preserving discretization for curvilinear collocated and rectilinear staggered meshes is found exploiting the skew-symmetric nature of the advection operator on square-root variables.

Another option to preserve symmetry is to use discrete filters to regularize the convective terms of the equation [33, 19]. The combination of a symmetry-preserving discretization and regularization for compressible flows is studied in [27].

Mimetic finite-difference methods also mimic the important properties of differential operators. An interesting review is given in [21], and recently, a second-order mimetic discretization of the Navier-Stokes equations conserving mass, momentum, and kinetic energy was presented in [23].

Castillo et al. have provided a framework for mimetic operators in [7, 8]. A second- and fourth-order mimetic approach is constructed for non-uniform rectilinear staggered meshes in [1, 3, 9]. In [10] their method is extended to curvilinear staggered meshes, but discrete conservation of mass, momentum and energy is not shown. A second-order mimetic finite-difference method for rectilinear staggered meshes is also constructed in [28].

Other mimetic finite-difference methods use algebraic topology to design and analyze compatible discrete operators corresponding to a continuous formulation [4, 18, 25]. In order to construct a discrete de Rham complex, certain conditions on reconstruction and reduction operators are imposed: they should be conforming, which means that the reconstruction is a right inverse of the reduction [4], they should be constant preserving [5], and the interpolation operator should commute with the differential operator [5]. In [6] a nice overview of mimetic methods is given. Discrete exterior calculus (DEC) is also related to these mimetic approaches [14]. However, in the papers about these mimetic methods, the discrete conservation of mass, momentum and energy is not studied. In [24], a mass, energy, enstrophy and vorticity conserving method is given for unstructured finite-element meshes, using also a conserving time integrator. Note that the method is not applied to staggered meshes, conservation of momentum is not mentioned, and that only incompressible models are used (chain rules are not needed and the advection operator is easier to process). More about mimetic time integration is found in [29].

Another class of symmetry-preserving methods uses the DG method [47, 12, 11]. In these papers, mass, momentum and total energy are conserved on curvilinear meshes, but staggering is not applied. The method is extended to shock capturing and positivity preservation in [48].

All the existing models in the literature have their own advantages and disadvantages. In general, they are not at the same time applicable for arbitrary discretization orders, sophisticated operators such as the advection operator, or the chain rule, and curvilinear staggered grids. The current paper presents a new discretization method that can handle these requirements simultaneously.

In this work, we first introduce some concepts by using a Galerkin-type approach, which is closely related to existing mimetic methods, and was also studied in [37]. Then, we present our new symmetry-preserving finite-difference technique for discretization in space. The novelty of this work is that it combines several important requirements for discretizations: the symmetry-preserving discretization is made for arbitrary order of accuracy; the method works for orthogonal and non-orthogonal structured curvilinear staggered meshes; and the method can be applied to a wide variety of continuous operators, involving chain rules and nonlinear advection, as will be illustrated by the shallow-water equations. The experiments show exact conservation of mass, momentum and energy, and convergence of the approximations corresponding to the expected order. The approach is very similar to the Richardson extrapolation scheme used in [16] and [17], but the current approach leads to smaller stencils, especially for high orders of accuracy. Apart from uniform, orthogonal grids, the experiments also use nonorthogonal curvilinear grids, in which the angles between grid lines are as small as 15${}^{\mbox{\scriptsize o}}$.

The energy equation, derived from the continuity, momentum and state equations, does not have its own discretization. Instead, discrete energy conservation is derived using the symmetries of the discrete operators [38]. One of the symmetries of interest in this paper is that the adjoint of the gradient is minus the divergence. See Table 2, for more details.

The subject of this paper is symmetry-preserving space discretization. In order to focus on this topic, issues concerning boundary conditions and time integration are not addressed. Instead, a periodic domain is used, and a standard time-integration method is used with such a small time step, that the time-integration errors are negligible.

Since the shallow-water equations are the main area of interest for the authors, all examples presented in this paper will be 2D, although the method can also be applied in 1D or 3D.

The outline of this paper is as follows: in Section 2, we present the three different models that will be used in the rest of the paper, and in Section 3, we give some information about the curvilinear grids we use. Section 4 contains the desired discretization properties, and the new symmetry-preserving discretization is explained in Section 5. The effectivity of this new method is investigated for the different models in Section 6. We conclude with a discussion of our method and future work in Section 7.

2 Models

In this paper, three different models are used to explain and test the symmetry-preserving ideas. Each of these models is presented in this section with respect to the following aspects:

• •

The model equations in terms of continuity, momentum, and state equations;

• •

A compact representation of the discrete model, including discrete operators like ${\sf DIV}$, the discrete divergence, and ${\sf GRAD}$, the discrete gradient;

• •

A consistent energy equation, derived from the continuity, momentum and state equations [38].

The energy density is given by the sum of the kinetic and internal energy densities: $e=e_{kin}+e_{int}$ (see Table 3 for the definition of energy density in each model). The energy equation will be used to show energy conservation of the models.

Section 4 discusses the properties that the discrete operators ${\sf DIV}$, ${\sf GRAD}$ and others are expected to have, for the discrete models to conserve mass, momentum and energy, after which Section 5 discusses how such operators can be constructed.

• •

Linear-wave equations

The simplest equation that can be used to discuss symmetry preservation, is the linear-wave equation, in which the evolution of the pressure $p$ in a domain $V$, the flow velocity $\vec{v}$ and the density $\rho$, are given by the continuity, momentum, and state equations:

[TABLE]

where $\rho_{0}$ is a constant ’reference’ density, and $c$ is the speed of sound (i.e. the propagation speed of waves). Initial conditions are specified at initial time $t=0$.

The discrete linear-wave equations are given by

[TABLE]

where p, rho and v are the vectors with the discrete pressures, densities and velocities, and where DIV is the discrete divergence and GRAD is the discrete gradient.

The continuity, momentum and state equations can be combined into the following energy equation:

[TABLE]

where $e$ is the energy density.

This model gives us the opportunity to introduce the relation between symmetry preservation and conservation, as well as the curvilinear staggered grid and staggered velocity components. Two approaches are used for the construction of a discretization. First, we use a Galerkin-type approach similar to the one used in [37], and secondly, we apply a finite-difference approach.

• •

Compressible-wave equations

We introduce a non-linearity into the system by including density variations in the continuity equation, and so obtain compressible-wave equations (without an advection term):

[TABLE]

where the pressure-dependent density function $R$ is assumed positive, continuous and monotonically non-decreasing, and where $Q$ is given by $Q(p):=\int^{p}1/R(y)~{}\mbox{d}y$. The discrete compressible-wave equations are given by

[TABLE]

where the operator ${\sf DIV}\tilde{\sf r}$ is a discretization of the operator $(\nabla\cdot\rho)$: this means that $({\sf DIV}\tilde{\sf r}\ {\tt v})_{i}$ approximates $(\nabla\cdot\rho\vec{v})(\vec{x}_{i})$. In Section 5.2.2, the explicit construction of the operator ${\sf DIV}\tilde{\sf r}$ will be discussed.

The corresponding energy equation is

[TABLE]

The model is analyzed for an arbitrary state equation $\rho=R(p)$, but tests are only conducted for the state equation $p=c^{2}\rho$, because then an exact solution is available to verify the results.

• •

Isentropic compressible Euler equations

Finally, after introducing an advection term, a symmetry-preserving discretization is formed for the following equations of isentropic compressible Euler gas dynamics [20]:

[TABLE]

The discrete system is given by

[TABLE]

The vector rv contains local momentum values, and is defined by ${\tt rv}:=\mbox{diag}({\sf Interp}_{v\leftarrow c}{\tt rho}){\tt v}$. ${\sf Interp}_{v\leftarrow c}$ indicates interpolation from the pressure grid points to the velocity grid points in the staggered grid. The operator ${\sf DIVr}$ is similar to the operator ${\sf DIV}\tilde{\sf r}$: it also approximates the operator $(\nabla\cdot\rho)$, but is constructed in a different way (Section 5.4). Section 5.4 also shows the construction of the discrete advection operator ADVEC.

The following energy equation can be derived from the continuity, momentum and state equations of the isentropic compressible Euler equations:

[TABLE]

Choosing the state equation $\rho=R(p)=\sqrt{2p/g}$ changes these equations into the shallow-water equations [20], which are investigated in this paper.

For each model, an exact solution is created using the propagation speeds and Riemann invariants as given in Table 1. For the compressible-wave equations and the isentropic compressible Euler equations, the exact solutions develop shocks after some time. After that moment, the (weak) solution becomes discontinuous, and although mass and momentum are still conserved, energy is not. In Section 6, we show that the symmetry-preserving method conserves mass, momentum and energy until the moment the shock occurs, and that the method converges with the expected order at a time before the shock occurs. At the moment the shock occurs, the convergence is slower. After the appearance of the shock, the discrete solution, which still conserves energy, cannot possibly be an accurate approximation of the correct (weak) solution, in which energy is lost. The comparison of the solutions therefore stops when the shock appears. The method could be extended with artificial viscosity (for instance in the form of flux limiters), which dissipate some of the energy in the system, and allow correct approximation of the discontinuous solution, but these techniques are beyond the scope of this article.

An example with a developing shock was chosen to show the properties of the method, because it shows that the method is accurate when it can be expected to, that it always conserves mass, momentum and energy, and that there are cases (discontinuous solutions) where the method needs an extension before it can be used to calculate meaningful solutions. All these aspects of the method would not be illustrated by an example where shocks do not develop.

3 Curvilinear grid

In this section, we introduce a parametrization for curvilinear grids. This is first done for a collocated mesh (used for scalar fields) and then for a staggered grid (used in vector fields).

3.1 Collocated grids and scalar fields

We introduce a uniform grid in computational-grid space, $\vec{\xi}_{c,i}$, that satisfies

[TABLE]

where $(m_{x},m_{y},m_{z})$ relate to the number of cells in each direction, and $(\Delta\xi,\Delta\eta,\Delta\zeta)$ are the corresponding cell widths. The subscript $c$ is chosen because the grid points can be seen as the cell centers of control volumes with vertices $\vec{\xi}_{c,i}+(\pm\Delta\xi,\ \pm\Delta\eta,\ \pm\Delta\zeta)^{\top}/2$.

The map ${\vec{X}}$ relates this uniform grid $\{\vec{\xi}_{c,0},\ldots,\vec{\xi}_{c,m_{x}m_{y}m_{z}-1}\}$ in computational-grid space to a curvilinear grid $\{\vec{x}_{c,0},\ldots,\vec{x}_{c,m_{x}m_{y}m_{z}-1}\}$ in physical space:

[TABLE]

The grid formed by the points $\xi_{c,i}$ is called the pressure grid, because this is the grid used to sample scalar fields, such as the pressure.

The properties of the mapping $\vec{X}$, along with those of the solution, determine the convergence rate of the discretization used. For clarity of presentation, we have only used smooth ($C^{\infty}$) mappings.

3.2 Staggered grids and vector fields

In a staggered grid, not only the pressure grid is used, but also a velocity grid. The velocity grid consists of grid points that are shifted by half a grid space, and therefore correspond to cell-face centers:

[TABLE]

where the subscripts $e$, $n$ and $t$ stand for ’east’, ’north’ and ’top’, respectively. The numbering of the grid points is visualized in Figure 1.

To find vector-field discretizations in curvilinear spaces, we should first introduce the notation of grid-aligned velocity components. A vector field $\vec{v}$ is represented by three scalar functions $v_{x}$, $v_{y}$ and $v_{z}$, using a local grid orientation $(\vec{r}_{x},\vec{r}_{y},\vec{r}_{z})$,

[TABLE]

where $\vec{r}_{x}$, $\vec{r}_{y}$ and $\vec{r}_{z}$ are orthonormal. The scalar functions $v_{x}$, $v_{y}$ and $v_{z}$ can be calculated from the vector field $\vec{v}$ using inner products:

[TABLE]

Typically [34, 46], a combination of covariant directions (i.e. parallel to grid lines) and contravariant directions (i.e. perpendicular to two of the grid lines) is used to construct a local grid orientation. Often [34], this leads to a different grid orientation in $\xi$-cell faces than in $\eta$- or $\zeta$-cell faces. For the method in this paper, a different kind of grid orientation is needed, because a simple representation of the kinetic energy requires a local grid orientation that is

• •

available in all points in space, not just in grid points;

• •

is the same for all quantities and equations;

• •

consists of three orthogonal directions.

In uniform Cartesian grids, the covariant and contravariant directions are the same, and the local grid orientation is given by $\vec{r}_{x}=(1,0,0)$, $\vec{r}_{y}=(0,1,0)$ and $\vec{r}_{z}=(0,0,1)$. In the general case, the contravariant and covariant directions are not the same, and an orthogonal grid orientation cannot be made from them.

An orthogonal, grid-aligned basis for vectors is based on the singular-value decomposition (SVD) of the grid directions, given by the vector ${\bf h}$ with the singular values and unitary matrices ${\bf P}$ and ${\bf Q}$, such that

[TABLE]

The orientation ($\vec{r}_{x},\vec{r}_{y},\vec{r}_{z}$) is now found by using only the unitary rotation matrices of the singular-value decomposition:

[TABLE]

This local-orientation matrix is exactly orthogonal, and forms a compromise between the covariant and contravariant directions. A 2D example of the computation of such local orientation is given in Figure 2.

A discrete vector field is represented by a vector with values for each direction: ${\tt v}=({\tt vx}^{\top},{\tt vy}^{\top},{\tt vz}^{\top})^{\top}$ contains components ${\tt vx}$, vy, and ${\tt vz}$ that are located at the grid points $\vec{x}_{e}$, $\vec{x}_{n}$, and $\vec{x}_{t}$, respectively. The components are given by

[TABLE]

Discrete samplings of the local grid orientation are stored in vectors ${\tt r**\_at\_*}$. As an example, we give the expressions for $\vec{r}_{x}$ at the cell centers or east cell-face centers, and for $\vec{r}_{y}$ at the north cell-face centers:

[TABLE]

Unlike in a uniform Cartesian grid, the discrete representation of the constant vector field $(1,0,0)$ is not a constant vector, and neither does it have zeros in the $n$-points and $t$-points of the grid. To compute the discrete representations c100, c010 and c001 of the constant fields

[TABLE]

we use equation (19) to find:

[TABLE]

4 Desired discretization properties

In this section, we investigate the properties that a symmetry-preserving discretization should satisfy. Therefore, we first introduce the different scalar products used.

4.1 Scalar products

Quantities such as densities and pressures form scalar fields, which are discretely approximated by vectors with a value for each pressure grid point. Let $a$ and $b$ be two continuous scalar fields, that correspond to the discrete scalar fields ${\tt a}$ and ${\tt b}$. The integral $\int_{V}ab\ \mbox{d}V$ can be approximated by the following scalar product:

[TABLE]

where only real numbers are used, and the vector ${\tt dVc}$ holds the integration weights for each pressure grid point. For instance, the discrete mass, M, is calculated from the vector ${\tt rho}$ of discrete densities according to ${\tt M}:=\langle{\tt c1},{\tt rho}\rangle_{c}$, where ${\tt c1}$ is the (constant) vector of only ones.

The situation is more complicated when vector fields are involved. Since the grid orientations have been chosen orthogonal, the scalar product $\langle{\tt v},{\tt w}\rangle_{v}$ of two discrete vector fields ${\tt v}$ and ${\tt w}$ can be written as the sum of the scalar products of the components:

[TABLE]

where ${\tt dVe}$, ${\tt dVn}$ and ${\tt dVt}$ are the vectors holding the integration weights for the grid points on the three directions of the velocity grid.

As an example, the total momentum $\vec{\tt M}$ is calculated for the compressible-wave model, using a discrete vector field rv of local momentum values:

[TABLE]

The scalar products introduced in this section determine the corresponding adjoints. For example, for a matrix ${\sf A}$ that maps values on pressure grid points to values on $e$-points, the adjoint ${\sf A}^{*}$ is given by

[TABLE]

since in that case, we have for all vectors ${\tt x}$ on $e$-points and ${\tt y}$ on pressure points:

[TABLE]

The scalar products and adjoints of this section are used to define the properties that discrete operators should satisfy.

4.2 Operator properties

In this section, we discuss Table 2, that shows some properties of the operators used to define the models from Section 2, along with their discrete equivalents. In [26], symmetry properties are used instead of the null space properties that are given in the table. These symmetry properties are applied to curvilinear collocated grids, and therefore this approach could not be applied to curvilinear staggered grids without the modifications presented in this paper.

The continuous versions of the properties often contain a boundary integral. For example, the change in total mass in the compressible-wave equations equals

[TABLE]

In this paper, we do not focus on boundary conditions, and therefore, we take periodic domains, such that equation (33) equals zero. We could also use closed walls, which satisfy $\vec{v}\cdot\vec{n}=0$. For other boundary conditions, the total mass, momentum and/or energy may change over time, due to an in- or outflux across the domain boundaries. This is beyond the scope of this paper.

Many discrete equalities in the table are of the same form, containing an adjoint operator and a vector consisting of ones. As an example, the equality for ${\sf DIV}^{*}{\tt c1}$ should be read as follows:

[TABLE]

which relates to its continuous equivalent. In Section 5, we will construct the discrete operators that satisfy these properties.

4.3 Conserved quantities

In each model, the conserved quantities (total mass, (total momentum) and total energy) can be identified, which are computed using integrals over the simulation domain $V$. For instance, for the linear-wave equations, the total mass $M:=\int_{V}\rho~{}\mbox{d}V$ is a conserved quantity. Table 3 displays the expressions for each of the conserved quantities in the models used in this paper. The symmetry-preserving discretizations should conserve these quantities in a discrete manner. In this section, we explain how discrete conservation can be proved, thereby using scalar products, and the properties from Table 2. Conservation of discrete mass in the linear-wave equations, for instance, is shown in the following very short proof:

[TABLE]

Note the transition from an inner product on the pressure grid to the velocity grid in the second last equality. The velocity $v$ is defined at the cell faces, and ${\sf DIV}\ {\tt v}$ is located at the cell centers. Therefore, the inner product $\left\langle{\tt c1},\rho_{0}{\sf DIV}~{}{\tt v}\right\rangle_{c}$ is computed in the pressure grid. Note that $\rho_{0}$ is constant, and therefore available both on the pressure and on the velocity grid. All proofs for the conservation of mass and momentum follow these same lines, using the left null space properties of the discretizations.

The energy-conservation proofs also use the symmetry properties. As an example, we show energy conservation for the isentropic Euler equations. The derivation of the continuous and discrete energy equations for all the models in this paper are presented in [38]. In the isentropic Euler equations, energy conservation is shown by calculating the time derivative of the discrete total energy E, which is the discrete approximation of the (continuous) total energy $E:=\int_{V}e\mbox{d}V$:

[TABLE]

If we apply the chain rule twice, we find that the first term equals

[TABLE]

For the second term, we use the fact that we only consider real numbers, such that for fields ${\tt a},{\tt b},{\tt c}$ we have

[TABLE]

since they all satisfy definition (31). If relation (34) and the product rule are applied several times, we find

[TABLE]

such that

[TABLE]

We can move the factor $e^{\prime}_{int}({\tt p})/R^{\prime}({\tt p})$ to the left, and use the continuity and momentum equations (3) to replace the time derivatives with spatial derivatives:

[TABLE]

The fraction $e^{\prime}_{int}/R^{\prime}$ is expanded using the definition of $e_{int}$ for the isentropic Euler equations shown in Table 3:

[TABLE]

Furthermore, the chain rule for ${\sf rGRAD}$ is applied (see Table 2), and the two advective terms are factorized:

[TABLE]

Finally, the symmetry property of ${\sf DIVr}$ and ${\sf rGRAD}$ is applied to the first terms:

[TABLE]

and also the symmetry property of ${\sf ADVEC}$ (see Table 2):

[TABLE]

since it holds that $\langle{\tt x},({\sf A}-{\sf A}^{*}){\tt x}\rangle=\langle{\tt x},{\sf A}{\tt x}\rangle-\langle{\tt x},{\sf A}^{*}{\tt x}\rangle=\langle{\tt x},{\sf A}{\tt x}\rangle-\langle{\sf A}{\tt x},{\tt x}\rangle=0$ for any matrix A. This means that total energy is indeed conserved: $\frac{\text{d}\tt E}{\text{d}t}=0$.

Now that we have investigated the properties that a discretization should satisfy to be symmetry preserving, we only need to construct such discretizations. This will cover the rest of this paper.

5 Numerical approaches

In this section, we describe two different techniques to spatially discretize the models: a Galerkin-type approach (which turns out to be very costly), and a new finite-difference technique. The Galerkin-type approach is only used to introduce some of the concepts that are needed for the finite-difference technique. Once the spatial discretization is found, time integration is needed for a simulation. We use the ODEPACK solver lsode [13], where it is verified that the time integration is accurate enough such that it has no influence on the results. Following this approach of highly-accurate time integration allows us to study the spatial discretization separately from the time-integration method, which in actual applications should be symmetry-preserving itself to guarantee the conservation of mass, momentum and/or total energy.

5.1 Galerkin-type approach

This section explains an approach that resembles the Galerkin method that is used in finite-element techniques. This technique will be used for the linear-wave model. The method is easy to understand and contains some concepts that help to understand the new finite-difference technique that is explained in the next section.

The traditional Galerkin method samples the solution to obtain a discrete vector of approximations. This discrete approximation is interpolated using a test function, and then the continuous differential operator is applied. The solution is then approximated by the discrete vector for which the inner product of the residual with all test functions equals zero [45].

In order to construct a symmetry-preserving discretization, we follow the same lines. We interpolate the discrete vector of approximations to a continuous function, apply the continuous differential operator, and sample the function back to a discrete vector. If the interpolation and the sampling operator are mutually adjoint, and the interpolation of constant functions is exact, then all properties in Table 2 (and thereby all conservation laws) are satisfied. The difference between the original Galerkin method and our approach is two-fold: (i) we focus on specific sampling and interpolation operators that are mutually adjoint, and (ii) in the original Galerkin method, the complete residual is sampled (multiplied by a test function), while we keep the original time derivative and only sample the space derivative and right-hand side. This resembles the idea of a lumped mass matrix. Note that the description below is closely related to mimetic discretizations that are defined using de Rham maps [31].

Let ${\cal J}_{c}$ be an interpolation operator that maps from the discrete field on the pressure grid to the continuous field, and ${\cal S}_{c}$ be the sampling operator that produces discrete values from a continuous function [44]. The interpolated fields are written using italic letters and sampled fields with truetype letters, for example, $f:={\cal J}_{c}~{}{\tt f}$, and ${\tt g}={\cal S}_{c}\ g$. The continuous differential operator ${\cal A}$ is applied to the continuous field obtained from interpolation of the discrete field using ${\cal J}_{c}$, and the result is mapped back using ${\cal S}_{c}$, which leads to the discrete operator

[TABLE]

The sampling operator and the interpolation operator are called mutually adjoint if ${\cal S}_{c}={\cal J}^{*}_{c}$. In the rest of this section, we will see that mutually-adjoint interpolation and sampling, combined with the exact interpolation of the constant functions, is in general enough to obtain the properties of the discrete operators in Table 2. For collocated grids, exact interpolation for constant functions means that ${\cal J}_{c}~{}{\tt c1}=1$. For staggered grids, we require that the interpolation operator $\vec{\mathcal{J}}_{v}$ satisfies

[TABLE]

5.1.1 Linear-wave equations: operators GRAD and DIV

In this section, we first define the mutually-adjoint interpolation and sampling operators that are used. Then, we present the symmetry-preserving operators DIV and GRAD used for the linear-wave equation. Earlier work on the Galerkin-type discretization was presented in [37].

For a standard uniform one-dimensional grid with $x_{i}=i$, standard Lagrange interpolation polynomials, $w$, are used. The corresponding curvilinear interpolation functions are constructed by applying the interpolation function $w$ in each of the three dimensions:

[TABLE]

The choice of the initial interpolation function $w$ determines the accuracy of the interpolation, as well as the sparseness of the discrete operators.

To see what the mutual adjointness actually means for the sampling and interpolation operators, they are both written in a more explicit form [44], in terms of the interpolation functions $w_{c,i}$ ($c$ for cell centers, with point index $i$). As usual, we use

[TABLE]

Mutual adjointness is found using the scalar product on the continuous space (defined by an integral) and the corresponding discrete scalar product on the pressure grid:

[TABLE]

such that the mutually-adjoint sampling functions should satisfy

[TABLE]

Here, we have chosen the integration weights ${\tt dVc}_{i}$ such that the integral of the sampling function equals 1: in that case, the sampling of a constant field is exact.

As an example, the standard linear interpolation of a one-dimensional function ${\tt f}$ between ${\tt f}_{i}$ and ${\tt f}_{i+1}$ is given by

[TABLE]

On the whole domain, the interpolated function satisfies equation (37a), with

[TABLE]

The corresponding sampling operator follows from equation (38):

[TABLE]

Similarly, higher-order interpolation and sampling operators can be constructed.

Next, we investigate the Galerkin-type approach for vector fields, in order to define the GRAD and DIV operators for the linear-wave equations. The interpolation of the discrete vector field ${\tt v}$ to the continuous vector field ${\vec{v}}=\vec{\mathcal{J}}_{v}{\tt v}$ is very similar to the interpolation for scalar fields (equation (37a)):

[TABLE]

with

[TABLE]

where ${w}_{e,i}$, ${w}_{n,i}$ and ${w}_{t,i}$ are interpolation functions similar to equation (36).

Sampling of a continuous vector field ${\vec{v}}$ to its discrete representation ${\tt v}=\mathcal{S}_{v}{\vec{v}}$ is similar to the way a continuous scalar field $f$ is sampled to its discrete representation ${\tt f}={\mathcal{S}}_{c}f$ (see equations (37b) and (38)), and is given by

[TABLE]

and analogously for ${\tt vy}$ and ${\tt vz}$.

The Galerkin-type approach requires that exact results are obtained when interpolating the discrete representations c100, c010 and c001 of the constant vector fields:

[TABLE]

Expanding the interpolation operator in ${\cal J}_{v}{\tt c100}$ (equation (39)) and taking the scalar product with $\vec{r}_{x}(\vec{x})$, an equation for the interpolation weights $w_{e,i}(\vec{x})$ is found. Here, we also use the definition of ${\tt c100}$ (equation (30)), the fact that the grid orientation is orthonormal, and definition (40):

[TABLE]

Combining this with the expected result ${\vec{r}}_{x}(\vec{x})\cdot(\vec{\mathcal{J}}_{v}{\tt c100})({\vec{x}})={\vec{r}}_{x}({\vec{x}})\cdot(1,0,0)$, and repeating the process for the constant fields c010 and c001, the following three equations are found for the weights $w_{e,i}(\vec{x})$:

[TABLE]

This can be written as a matrix-vector equation with three equations and involving one unknown for each point in the $e$-grid:

[TABLE]

This is an underdetermined system. We solve for small corrections (in the weighted least-squares sense) to the standard Lagrange interpolation polynomials while making sure that the support of the weight function is no larger than that of the Lagrange polynomials. Similar underdetermined systems have to be solved for $w_{n}$ and $w_{t}$.

The GRAD and DIV operators in the Galerkin-type approach are given by

[TABLE]

In practice, the divergence operator in $\vec{x}_{c,i}$ is computed as follows:

[TABLE]

with

[TABLE]

and similarly for ${\sf DIV}_{i,(n,j)}$ and ${\sf DIV}_{i,(t,j)}$, and where $s_{c,i}$ are the sampling functions of equation (38).

5.1.2 Computational costs

In the preceding sections, we have presented the Galerkin-type approach to compute discrete operators. In this section, the computational costs are investigated.

At the start of a simulation, the discrete operators are computed using volume integrals. To guarantee symmetry preservation, these integrals need to be calculated up to machine precision. We achieve high accuracy by using Richardson extrapolation of the composite trapezoidal rule. Although this is a time-consuming process, it is only done once to initialize a simulation. The result is stored in a matrix that corresponds to the discrete operator.

After initialization, the evaluation of the discrete operators is also quite expensive. As an example, we investigate the computational work needed to evaluate the divergence operator, which is related to the number of nonzeros in the DIV matrix. The interpolation function $w_{c,\tt ix,iy,iz}$ of order ${\tt N}$ is zero outside the support $\left[{\tt ix}-{\tt N/2}:{\tt ix}+{\tt N/2},{\tt iy}-{\tt N/2}:{\tt iy}+{\tt N/2},{\tt iz}-{\tt N/2}:{\tt iz}+{\tt N/2}\right]$. Similarly, the interpolation function $w_{e,\tt jx+1/2,jy,jz}$ is zero outside the support

${\tt\left[jx+1/2-N/2:jx+1/2+N/2,jy-N/2:jy+N/2,jz-N/2:jz+N/2\right]}$.

Therefore, interpolation functions $w_{c,\tt ix,iy,iz}$ and $w_{e,\tt jx+1/2,jy,jz}$ overlap for all indices jx,jy,jz for which

[TABLE]

which is ${\tt 2*N*(2*N-1)*(2*N-1)}$ functions. The divergence matrix has three times as many nonzeros (for $v_{x}$, $v_{y}$ and $v_{z}$), which is ${\tt 6*N*(2*N-1)*(2*N-1)}$ ${\tt\approx 24*N^{3}}$ nonzeros per matrix row. In 2D, we find ${\tt 4*N*(2*N-1)\approx 8*N^{2}}$ nonzeros per matrix row, and in 1D ${\tt 2*N}$.

A cheaper approach to arrive at symmetry-preserving discretizations is presented in the next section.

5.2 Finite-difference approach

The Galerkin-type approach described in the previous section has certain disadvantages. For example, the calculation of the discrete operators requires highly-accurate integration, and the resulting matrices have many nonzeros. Therefore, this method leads to relatively long calculation times. Furthermore, the application of the Galerkin-type approach to difficult operators such as the advection operator is not straightforward.

A more efficient and flexible approach is found using a finite-difference strategy. In this section, we explain the finite-difference technique for the different operators needed in the linear-wave, compressible-wave and isentropic Euler models.

5.2.1 Linear-wave equations: operators DIV and GRAD

A higher-order finite-difference approximation DIV for the divergence is constructed from an exact formulation (44) of the divergence in a cell center. Subsequently, standard finite-difference stencils are applied to this formulation, resulting in a higher-order, conservative approximation of the divergence. The resulting discrete operator DIV combines a series of steps such as interpolation and differentiation with respect to $\vec{\xi}$. All of these steps have continuous equivalents, and therefore its negative adjoint ${\sf GRAD}:=-{\sf DIV}^{*}$ is an accurate approximation of the gradient. The approach is very similar to the Richardson extrapolation scheme used in [16] and [17], but the current approach leads to smaller stencils, especially for high orders of accuracy.

The resulting finite-difference discretizations are locally conservative for the mass balance, because the change in a control volume can be written as the net effect of fluxes, each of which describes the transport between two ’neighbor’ grid points. In rectilinear, orthogonal, possibly stretched, grids the momentum equation is locally conservative in the same sense, with momentum control volumes drawn around the cell faces of the mass control volumes. In curvilinear or nonorthogonal grids, the momentum equation cannot be seen as locally conservative, because every sampled velocity contains momentum in a different direction, making it impossible to transfer momentum from one grid point/control volume to another. For similar reasons, the energy balance cannot be seen as locally conservative, because internal energy is sampled in the mass control volumes, and kinetic energy in the cell faces.

For the derivation of the finite-difference technique, a control volume $V_{i}$ is used, centered around $\vec{x}_{c,i}$. The control volume is an interval (1D), quadrilateral (2D) or hexahedron (3D). The method is described in 2D for simplicity, but it works the same way in 1D and 3D. In 2D, the control volume $V_{i}$ is given by

[TABLE]

with vertices given by

[TABLE]

see Figure 3(a).

The scaled flux of the vector field $\vec{v}$ out of the east cell face is denoted $F_{e,i}$ and given by

[TABLE]

and the scaled fluxes $F_{w,i}=-F_{e,i}(-\varepsilon)$, $F_{n,i}$ and $F_{s,i}=-F_{n,i}(-\varepsilon)$ out of the west, north and south cell faces analogously.

Gauss’s theorem equates the integrated divergence to the net outflux:

[TABLE]

For very small $\varepsilon$, the volume integral may be approximated by the midpoint rule (note the occurrence of the Jacobian determinant of the transformation $\vec{X}$):

[TABLE]

and this leads to the following exact representation of the divergence $\nabla\cdot\vec{v}$ (when taking $\varepsilon\rightarrow 0$):

[TABLE]

because the integration weights will be chosen equal to

[TABLE]

Note that the scaled fluxes $F_{e,i}$ and $F_{n,i}$ are given in terms of boundary integrals. They can be replaced by approximations in which only velocity values $\vec{v}$ occur, and no integration is needed. To do this, the cell-face centers $\vec{x}_{e,i}$ and $\vec{x}_{n,i}$ are defined by

[TABLE]

The flux $F_{e,i}$ is approximated by the approximate flux $\tilde{F}_{e,i}$, obtained using the midpoint rule:

[TABLE]

where the definition of the outward-directed normal on the east cell face is used, and the perp operator $\perp$ is defined as $(x,y)^{\perp}:=(-y,x)$. The location of $\tilde{F}_{e,i}(\varepsilon)$ for two different values of $\varepsilon$ is visualized in Figure 3(a).

Similarly, we define

[TABLE]

Since the midpoint rule is at least second-order accurate, the approximate fluxes and their derivatives are exact at $\varepsilon=0$. Hence, an exact formulation for the divergence $\nabla\cdot\vec{v}$ is given by

[TABLE]

which indeed does not include any integrals.

The discrete divergence ${\sf DIV}$ is now found by executing the following steps:

Interpolation of the velocities.

To compute $\tilde{F}_{e,i}(\varepsilon)$ (equation (42)), we need to evaluate the discrete value of

[TABLE]

see equation (15). Since $\vec{x}_{e,i}(\varepsilon)$ is a point in the $e$-grid, the discrete $y$-component of $\vec{v}$ (the sampled velocity components vy) is not available yet. The component ${\tt vy}$ is first interpolated to the pressure grid (using the destaggering matrix N2C) and then to the $e$-grid (with the staggering matrix C2E) using an interpolation procedure similar to the one in Section 5.1. The result will be stored in the vector ${\tt vy\_at\_e}={\sf C2E}~{}{\sf N2C}~{}{\tt vy}=:{\sf N2E}~{}{\tt vy}$.

Similarly, the components vx are interpolated to ${\tt vx\_at\_n}={\sf C2N}~{}{\sf E2C}~{}{\tt vx}=:{\sf E2N}~{}{\tt vx}$, so that the complete velocity vector is available at $e$- and $n$-points.

In matrix-vector notation, this interpolation step can be written as

[TABLE] 2. 2.

Computation of the fluxes.

All the approximate fluxes use the velocity vector in one point of the staggered grid: every $\tilde{F}_{e,i}$ uses the velocity vector at $e$-points, and $\tilde{F}_{n,i}$ at $n$-points. Therefore, the flux values are located at these grid points, for instance, the flux at the $e$-point $\vec{x}_{e,i}$ is

[TABLE]

and similar for $\overline{\tt F}_{n,k}$ for the $n$-point $\vec{x}_{n,k}$. In matrix-vector form this equals

[TABLE]

Note that the use of the interpolated velocity vectors is not the only difference between $\tilde{F}_{e,i}$ and $\overline{\tt F}_{e,i}$. Also the $x$-coordinate of

[TABLE]

used in $\tilde{F}_{e,i}$ (at distance $\varepsilon\Delta\xi$, see Figure 3(a)), differs from the $x$-coordinate of

[TABLE]

used in $\overline{\tt F}_{e,i}$ (fixed at cell faces, see Figure 3(b)). For constant vector fields, the interpolation is chosen to be exact, and we find the following relation:

[TABLE]

For non-constant vectors, this only approximately holds. Similarly, the following relation for constant vector fields holds:

[TABLE] 3. 3.

Finite-difference step.

The flux vectors $\overline{\tt F}_{e}=(\overline{\tt F}_{e,j})_{j}$ and $\overline{\tt F}_{n}=(\overline{\tt F}_{n,k})_{k}$ are used to calculate the discrete operator DIV corresponding to equation (44). Standard differentiation stencils are applied to calculate DIV (which lives on the pressure points). Therefore, we use the matrices ${\sf DIFFX}$ and ${\sf DIFFY}$, that make sure that central differences are used. In internal pressure points $\vec{x}_{c,i}$, row $i$ of DIFFX contains a $-1$ at position $(i,i-(\varepsilon+1/2))$ and a $1$ at $(i,i+(\varepsilon-1/2))$, and DIFFY is defined in a similar way, so that the $i$-th entry is given by

[TABLE]

The divergence (on the whole domain) can be computed as follows

[TABLE]

where $\alpha$ contains the coefficients of the differentiation stencil. For example, for second-order differentiation, $\alpha(1/2)=1$, and for fourth order, we use [15, 22]

[TABLE]

We need to show that this choice for the divergence operator indeed satisfies the null space properties (Table 2). The left null space property holds if ${\sf DIV}^{*}{\tt c1}=0$. Indeed, the integral of a divergence is zero:

[TABLE]

because ${\tt c1}^{\top}{\sf DIFFX}(\varepsilon)={\tt c1}^{\top}{\sf DIFFY}(\varepsilon)=0$. Next, we prove that the divergence of a constant vector field is zero. The interpolations used to calculate vx_at_n and vy_at_e are constructed such that they are exact for constant vector fields. Therefore, for a discrete vector field v representing a constant vector field ($\vec{v}(\vec{x})=\vec{c}$), a zero divergence is found, because at grid point $\vec{x}_{c,i}$, we have (note the exact relation between $\overline{\tt F}$ and $\tilde{F}$)

[TABLE]

The gradient operator ${\sf GRAD}$ is obtained from the divergence’s adjoint and is given by

[TABLE]

Using the notations ${\tt Nx}_{e}$, ${\tt Nx}_{n}$, ${\tt Ny}_{e}$, ${\tt Ny}_{n}$, E2N and N2E, introduced in this section, and the definition of the adjoint, the gradient operator GRAD is given by

[TABLE]

where the operator GRAD1 is given by

[TABLE]

In general, one cannot expect the (discrete) adjoint ${\sf GRAD}=-{\sf DIV}^{*}$ of an accurate discretization of the divergence to be an accurate discretization of the gradient. However, the divergence DIV can be seen as an operator in computational space, discretized on a uniform $(\xi,\eta)$-grid. On uniform grids, the adjoint of an interpolation operator is itself an accurate interpolation (in the opposite direction). Similarly, the adjoint of a discrete derivative operator is itself an accurate discrete derivative operator on uniform grids. Therefore, an accurate discretization of the continuous gradient operator is found if adjoint interpolations and adjoint derivative operators are combined to a (discrete) GRAD operator.

The finite-difference approach is computationally more efficient than the Galerkin-type approach. Similar to Section 5.1.2, we investigate the computational work needed to evaluate the divergence operator. Here, we assume interpolation and differentiation of order N. We start with the interpolation of $v_{x}$, $v_{y}$ and $v_{z}$ to the pressure grid points. Each interpolation involves N grid points, so this step is a matrix-vector multiplication with 3 * N nonzeros per grid point. Next, $v_{y}$ and $v_{z}$ are interpolated to the points of the $e$-grid, $v_{x}$ and $v_{z}$ to the points of the $n$-grid, and $v_{x}$ and $v_{y}$ to the points of the $t$-grid. This is three matrix-vector multiplications, each involving ${\tt 2*N}$ nonzeros. Finally, the evaluation of the divergence involves an $x$-, $y$- and $z$-derivative, each of which involves ${\tt N}$ nonzeros. Therefore, the evaluation of the divergence requires ${\tt 3*N+3*2*N+3*N}={\tt 12*N}$ nonzeros per grid point. In 2D, the calculation involves ${\tt 2*N}+{\tt 2*N+2*N}={\tt 6*N}$ nonzeros per grid point, and in 1D N nonzeros. Indeed, the calculation of the divergence operator is much cheaper than when using a Galerkin-type approach. This observation generally holds for all the operators that are designed in this paper.

5.2.2 Compressible-wave equations: operators $\tilde{\sf r}{\sf GRAD}$ and ${\sf DIV}\tilde{\sf r}$

The discretization of the compressible-wave equations requires the operator $\tilde{\sf r}{\sf GRAD}$, which is an approximation of the operator $(\tilde{\sf r}{\sf GRAD}~{}{\tt f})_{i}\approx(\rho\nabla f)(\vec{x}_{c,i})$. The operator $\tilde{\sf r}{\sf GRAD}$ is obtained by applying the same interpolations used in the GRAD-operator,

[TABLE]

where the intermediate densities $\widetilde{\tt rho}_{e}$ and $\widetilde{\tt rho}_{n}$ are given by

[TABLE]

Here, the functions $Q$ and $S$ that were introduced in Table 2 are applied. Obviously, the chain rule for $\tilde{\sf r}{\sf GRAD}$ holds: $\tilde{\sf r}{\sf GRAD}\ S({\tt p})={\sf GRAD}\ Q({\tt p})$. Of course, special care is necessary for the numerically-stable calculation of the intermediate densities, especially when elements of the diagonal matrix $\mbox{diag}({\sf DIFF}~{}S({\tt p}))$ are very small.

5.3 Isentropic Euler equations: operators rGRAD and DIVr

The techniques presented in the preceding sections are almost sufficient for a symmetry-preserving discretization of the isentropic Euler equations. Only the advection operator is not yet available. The advection operator is very similar to the divergence operator, on which it will be based. Here, we do not use the operators ${\sf DIV}\tilde{\sf r}$ and $\tilde{\sf r}{\sf GRAD}$, but the very similar operators ${\sf DIVr}$ and ${\sf rGRAD}$, given by

[TABLE]

where the intermediate densities ${\tt rho}_{e}$ and ${\tt rho}_{n}$ are given by

[TABLE]

5.4 Isentropic Euler equations: operator ADVEC

A symmetry-preserving approximation for the advection term, which also preserves momentum, is constructed using a ’splitting’ strategy very similar to the one presented in [48], [47] and [35], and classical results of Tadmor [30]. For the construction of the advection operator ADVEC, we start with a discrete advection operator ${\sf ADVEC}_{s}$ that works on a scalar field ${\tt f}_{v}$. Unlike a ’normal’ scalar field f, which is sampled at the cell centers, the scalar field ${\tt f}_{v}=({\tt f}_{e}^{\top},{\tt f}_{n}^{\top})^{\top}$ is sampled at the cell faces. The scalar advection operator ${\sf ADVEC}_{s}$ approximates the continuous advection operator: $({\sf ADVEC}_{s}{\tt f}_{v})_{c,i}\approx(\nabla\cdot\rho\vec{v}f)(\vec{x}_{c,i})$. This operator ${\sf ADVEC}_{s}$ has the following properties compared to ADVEC:

• •

Input field of the operator:

${\sf ADVEC}_{s}$ works on scalar fields sampled at the cell faces, and not, like ADVEC, on a staggered representation of a vector field;

• •

Output field of the operator:

The result ${\sf ADVEC}_{s}~{}{\tt f}_{v}$ lives on cell centers, and is not, like the result ${\sf ADVEC}~{}{\tt v}$, a staggered representation of a vector field;

• •

Applying the operator to a constant field:

When applied to a constant field, the operators ${\sf ADVEC}_{s}$ and ADVEC should both return the mass flux divergence ${\sf DIVr}~{}{\tt v}$.

• •

Left null space:

All constant fields should be in the null space of the adjoint operators ${\sf ADVEC}^{*}_{s}$ and ${\sf ADVEC}^{*}$.

• •

Meaningful adjoint:

${\sf ADVEC}_{s}$’s adjoint is an accurate approximation of the adjoint of the continuous advection operator, like ADVEC’s adjoint.

Because of the differences between ${\sf ADVEC}_{s}$ and ${\sf ADVEC}$ in the input and output spaces, certain interpolations and rotations are necessary which are not necessary in [48], [47] and [35].

The scalar advection operator ${\sf ADVEC}_{s}$ is defined as (cf. equation (47))

[TABLE]

In this computation, $\overline{\tt F}_{e}$ and $\overline{\tt F}_{n}$ are computed using $\vec{v}$ and equations (45) and (46). This advection operator ${\sf ADVEC}_{s}$ returns values at the pressure points of the grid. It satisfies the left null space property and the ’advection of a constant function is the divergence’-property

[TABLE]

The ’advection of a constant function is the divergence’-property is closely related to ADVEC’s symmetry property in Table 2.

In combination with the interpolation matrices C2E and C2N, the scalar advection operator can be used in models in which a scalar quantity is advected or convected, like temperature, or the concentration of a dissolved substance. The models discussed in this paper do not have such an advected/convected quantity.

The advection needed in the momentum equation of the isentropic Euler equations is not applied to a scalar field ${\tt f}_{v}$, but to a vector field ${\tt w}$. This is done by applying the scalar advection operator to the vector field’s components separately. To make sure the operator is momentum preserving, this should not be done in the local grid orientation, but in the Cartesian orientation. The discretization process below is described in detail for 2D for simplicity, but all steps have also been worked out for 3D.

cell-face interpolation and transformation to Cartesian orientation.

The interpolation operators E2N and N2E of Section 5.2.1 are needed to calculate the complete vectors, represented by wx, wy_at_e, wx_at_n and wy. Next, the vector field is transformed into the Cartesian grid. As an example, the vector field $\vec{w}_{e,i}=\vec{r}_{x}(\vec{x}_{e,i}){\tt wx}_{i}+\vec{r}_{y}(\vec{x}_{e,i}){\tt wy\_at\_e}_{i}$ is written as follows:

[TABLE]

This means that the component of $\vec{w}_{e,i}$ in the Cartesian direction $(1,0)$ equals

[TABLE]

such that

[TABLE]

Similarly, the representation of $\vec{w}_{n,i}$ in the Cartesian basis vector $(1,0)$ can be computed:

[TABLE]

The scalar field ${\tt w}_{v,(1,0)}=({\tt w}_{e,(1,0)}^{\top},{\tt w}_{n,(1,0)}^{\top})^{\top}$ combines the two results. This scalar field has values in the complete velocity grid, like the scalar field ${\tt f}_{v}$ to which the advection operator ${\sf ADVEC}_{s}$ was applied earlier. It contains the horizontal component of the vector field ${\tt w}$.

Using the same procedure, the vertical component is found:

[TABLE] 2. 2.

${\sf ADVEC}_{c}:$ application of ${\sf ADVEC}_{s}$ to each component.

The advection operator ${\sf ADVEC}_{s}$ can be applied to ${\tt w}_{v,(1,0)}$ and ${\tt w}_{v,(0,1)}$, and the result will be the operator ${\sf ADVEC}_{c}$ (advection at pressure points):

[TABLE]

This operator approximates the continuous advection operator: $\vec{a}(\vec{x}_{c,i})\approx(({\tt a}_{c,(1,0)})_{i},({\tt a}_{c,(0,1)})_{i})$. This operator satisfies the left null space property given by

[TABLE]

and the ’advection of a constant function is the divergence’-property given by

[TABLE] 3. 3.

${\sf ADVEC}_{a}$: interpolation to the velocity grid and transformation to the curvilinear grid.

The result of the advection operator ${\sf ADVEC}_{c}$ consists of standard Cartesian components at the pressure points of the grid, whereas the advection operator is needed for the components in local grid orientation at the velocity grid points. For this, we need to interpolate the components and transform the results. Interpolation is done by applying the adjoint operators ${\sf E2C}^{*}$ and ${\sf N2C}^{*}$, such that for each grid point $\vec{x}_{c,i}$, we compute

[TABLE]

Note that the interpolation from the cell centers to the cell faces is done with the adjoint interpolations ${\sf E2C}^{*}$ and ${\sf N2C}^{*}$ and not with the forward interpolations ${\sf C2E}$ and ${\sf C2N}$, as might have been expected. This is needed in order to satisfy the left null space property (55). In general, the adjoint of an accurate interpolation cannot be expected to be an accurate interpolation, but in the current case, the interpolations ${\sf E2C}$ and ${\sf N2C}$ are standard destaggering operations working on the computational grid, and their adjoints are accurate interpolations.

Finally, the results are transformed to the basis spanned by the local grid orientation vectors: inner products with $r_{x}(\vec{x}_{e,i})$ and $r_{y}(\vec{x}_{n,i})$ are taken. Doing so defines the advection operator ${\sf ADVEC}_{a}$, which has almost all the properties listed in Table 2:

[TABLE]

In Section 5.2.1, the interpolations E2C and N2C were required to be exact for the constant vector fields c10 and c01. The construction of the advection operator ${\sf ADVEC}_{a}$ also requires products of such fields to be exact, so (see equation (30) for the definition of c10 and c01)

[TABLE]

and similar for N2C.

This advection operator ${\sf ADVEC}_{a}$ also satisfies the left null space property. This can be seen by computing

[TABLE]

Due to orthonormality of the local grid orientation, we have

[TABLE]

such that

[TABLE]

Therefore,

[TABLE]

and similar for $\langle{\tt c01},{\sf ADVEC}_{a}~{}{\tt w}\rangle_{v}$.

The advection operator ${\sf ADVEC}_{a}$ also has the following ’advection of a constant function is the divergence’-property:

[TABLE]

for any discrete vector field c that is the discrete representation of a constant vector field. 4. 4.

ADVEC: construction of the symmetry-preserving advection operator.

The only property from Table 2 that the advection operator ${\sf ADVEC}_{a}$ does not have, is the symmetry property. We obtain the advection operator ADVEC that has the symmetry property (as well as the left null space property) by combining the operator ${\sf ADVEC}_{a}$ and its adjoint [35]:

[TABLE]

The advection operator ADVEC has precisely the symmetry property for the advection operator, provided that the interpolation operator ${\sf Interp}_{v\leftarrow c}$ used in the table is given by

[TABLE]

Negative densities may be the result of this interpolation, in cases where the density rapidly approaches zero, for instance in cases of wetting and drying, and this may affect the accuracy and stability of the numerical scheme. In such cases, the order of interpolation and discretization may be locally reduced (in a symmetry-preserving way). Alternatively, techniques like artificial porosity may be used to regularize the solution [36]. Such situations did not occur in the examples shown in this paper, and the techniques mentioned (local order reduction and artificial porosity) are left outside the scope of this paper.

All required symmetry-preserving discrete operators are now defined, and we can investigate the performance of these operators for the three models. This will be done in the next section.

6 Numerical results

In this section, the numerical results for the three different models are studied. The calculations for the linear-wave equations are done both with the Galerkin-type approach (Section 5.1) and with the finite-difference technique (Section 5.2). This makes it possible to compare the results and computation times. The more advanced models for compressible waves and shallow waves are only investigated with the finite-difference technique.

All calculations are done both on a uniform and on a periodic, non-orthogonal, curvilinear mesh. The periodicity of the curvilinear mesh gives us the opportunity to study symmetry preservation without focusing on boundary conditions (Section 4.2). An example of the uniform and curvilinear grids used for $20\times 20$ grid points is given in Figure 4.

6.1 Linear-wave equations

Choosing $p=\rho$ ($c=1$ in equation (1)), an exact solution for the linear-wave equation is discussed in Section 2. Since the propagation speed $V_{+}=1$ only has one value, the solution is a traveling wave, which does not change its shape, but only its position, traveling in ’southeast’ direction through the infinite periodic domain, given by

[TABLE]

In fact, the solution on the uniform mesh is a 1D function of $x-y$, so the results at time $T=10$ can be represented in a 1D plot (Figure 5). The numerical approximation, shown in red in the figure, is close to the exact solution, shown in blue. The relative errors of the Galerkin-type approximations are shown in Table 4 for different grid sizes and interpolation orders, and for the finite-difference approach they are shown in Table 6. It is seen that the solution becomes more accurate as more grid points are used and/or higher-order interpolation are used. The order of the approximation is equal to the interpolation order. Note that the results for our curvilinear grid are about the same as for the uniform grid. This is expected, since the exact solution, as well as the grid transformation $\vec{X}$ are infinitely smooth.

The effect of symmetry preservation is, that, assuming exact time integration, the total discrete mass M, momentum $\vec{\tt M}$ and total energy E are constant up to machine precision, even in the most inaccurate solutions shown in Table 4. In Tables 5 and 7, the mass, momentum and energy losses are given for a fourth-order discretization on a $20\times 20$ grid. Independently of the time-integration tolerance used, the mass is always conserved up to a very high precision, which depends on how accurately the matrix entries were calculated (using Richardson extrapolation of the composite trapezoidal rule). For the finite-difference models, where the discrete operators require no numerical integration, the error for mass and momentum is really in the order of machine accuracy. Considering energy conservation, the error depends on the accuracy of the time integration.

For linear models, the implicit midpoint rule and its higher-order Gauss-Runge-Kutta generalizations will conserve energy (in fact all quadratic and linear invariants) up to machine precision. For the nonlinear models, dedicated time-integration methods are needed. The current paper leaves time-integration methods outside the scope, because a standard time-integration method is used with such a small time step, that the time-integration errors are negligible.

Discrete mass, momentum and total energy are conserved up to machine precision, even in the most inaccurate solutions. The losses for different time-integration tolerances for the Galerkin-type approach are shown in Table 5. Similar results obtained in the finite-difference approach are shown in Table 7.

Because both the Galerkin-type approach and the finite-difference technique are available for the linear-wave equations, their relative computational performance can be compared. The Fortran code was compiled using gfortran version 7.3.0, and simulations were run on a desktop computer running Ubuntu 18.04 on an Intel i5-7400 CPU with 4 cores. Both methods use the same shared memory parallelization technique (OpenMP). We compare the results for interpolation (and derivative) order 8 and $40\times 40$ grid points.

The Galerkin-type approach needs 118 seconds for the computation, thereby doing 21016 time steps and computing 27751 time derivatives. This means that this implementation computes 178 time steps per second, and 234 time-derivative evaluations per second.

On the other hand, the finite-difference approach needs 14 seconds for the total computation, in which it computes 26901 time steps and 31199 time derivatives. This implementation is able to compute 1922 time steps per second, and 2230 evaluations per second. We clearly see that the finite-difference approach is about ten times faster than the Galerkin-type approach. This also corresponds to the theory (Sections 5.1.2 and 5.2.1): for this 2D problem with 8th order interpolation, the Galerkin-type approach results in 480 nonzeros per matrix row, whereas the finite-difference approach only uses 48 nonzeros per grid point (a factor 10 difference).

6.2 Compressible-wave equations

Unlike the solutions in the previous numerical examples, the solution of the compressible-wave equations has variable propagation speed. This means that the crests of the wave travel faster than the troughs. This effect causes the wave to steepen over time, and finally to become discontinuous. In Figure 6, the initial solution and the numerical approximation at the time $t_{N}$ (when the solution becomes discontinuous) are shown. Obviously, the solution becomes more difficult to approximate for a numerical method as the moment $t_{N}$ approaches. However, we see that away from the shock, the approximations are still accurate. Because the solution becomes increasingly difficult to approximate as time $t=t_{N}$ approaches, results are shown for $t=1.09=t_{N}/2$ and for $t=2.18=t_{N}$ in Table 8. Halfway the simulation, the convergence order approaches the theoretical order; for the final time, only order 1 can be achieved.

Table 9 shows that even the most inaccurate solutions at the final time conserve discrete mass, momentum and energy up to machine precision. However, the physically-relevant solution should have lower total energy after the shock has appeared. Therefore, numerical viscosity can be added, for instance in the form of flux limiters, but such techniques are beyond the scope of this paper.

6.3 Shallow-water equations

Like the compressible-wave equations, the variable propagation speed leads to a steepening wave in the exact solution. The results at time $t=t_{N}$ (the moment that the solution becomes discontinuous) are visualized in Figure 7. The accuracies at time $t=t_{N}/2$ are shown in Table 10: the accuracy orders are approaching the theoretical order. At the final time, only an order 1 can be achieved. Mass, momentum and energy are again conserved up to machine precision, as shown in Table 11.

7 Conclusion

This paper describes the construction of an arbitrary-order symmetry-preserving finite-difference technique on structured curvilinear staggered grids, offering flexibility and accuracy of the numerical approximations. The numerical examples presented in this work show that the method leads to results in which a high accuracy can be obtained, while the discrete mass, momentum and energy are all preserved. Also, the computation times are investigated for the linear-wave model, and the method is shown to be faster than when using a symmetry-preserving Galerkin-type approach. We have shown that the finite-difference method is able to handle difficult operators such as the advection operator in the shallow-water equations.

This paper does not address the application of energy-conserving (symplectic) time integration. This is left outside the scope of this paper, and so are the discussion of boundary conditions and the use of numerical viscosity (for instance using flux limiters), which is necessary when shocks occur in the solution. Finally, future work includes handling local grid refinements.

Acknowledgments The authors gratefully wish to acknowledge the useful comments provided by Jennifer Ryan (University of East Anglia) and Kees Vuik (Delft University of Technology) that helped to shape this work. We thank the reviewers for the suggestions they provided for revising our paper. We also thank Shell and VORtech that allowed us to develop the research.