A Fast Numerical Scheme for the Godunov-Peshkov-Romenski Model of Continuum Mechanics

TL;DR

This paper introduces a computationally efficient second-order numerical scheme for the Godunov-Peshkov-Romenski model, balancing accuracy and speed for continuum mechanics simulations.

Contribution

It presents a novel operator splitting scheme combining finite volume WENO reconstruction with analytic solutions for ODEs, offering a practical alternative to higher-order methods.

Findings

The scheme achieves sufficient accuracy for practical purposes.

It is computationally cheaper than existing ADER-WENO schemes.

Convergence studies confirm the expected order of accuracy.

Abstract

A new second-order numerical scheme based on an operator splitting is proposed for the Godunov-Peshkov-Romenski model of continuum mechanics. The homogeneous part of the system is solved with a finite volume method based on a WENO reconstruction, and the temporal ODEs are solved using some analytic results presented here. Whilst it is not possible to attain arbitrary-order accuracy with this scheme (as with ADER-WENO schemes used previously), the attainable order of accuracy is often sufficient, and solutions are computationally cheap when compared with other available schemes. The new scheme is compared with an ADER-WENO scheme for various test cases, and a convergence study is undertaken to demonstrate its order of accuracy.

Tables

Table 1. Table 1: Initial conditions for the slow opposing shear flow test

	$\rho$	$p$	$𝒗$	$A$	$𝑱$
	$1$	$1/\gamma$	$(0,-0.1,0)$	$I_3$	$\mathrm{𝟎}$
$x\ge 0$	$1$	$1/\gamma$	$(0,0.1,0)$	$I_3$	$\mathrm{𝟎}$

Table 2. Table 2: Initial conditions for the heat conduction test

	$\rho$	$p$	$𝒗$	$A$	$𝑱$
	$2$	$1$	$\mathrm{𝟎}$	$\sqrt[3]{2}\cdot I_3$	$\mathrm{𝟎}$
$x\ge 0$	$0.5$	$1$	$\mathrm{𝟎}$	$\frac{1}{\sqrt[3]{2}}\cdot I_3$	$\mathrm{𝟎}$

Table 3. Table 3: Wall time for various tests (all with 200 cells) under the ADER-WENO method and the Split-WENO method

	ADER-WENO	Split-WENO	Speed-up
Stokes’ First Problem ($\mu ={10}^{-2}$)	265s	38s	7.0
Stokes’ First Problem ($\mu ={10}^{-3}$)	294s	38s	7.7
Stokes’ First Problem ($\mu ={10}^{-4}$)	536s	38s	14.1
Viscous Shock	297s	56s	5.3
Heat Conduction in a Gas	544s	94s	5.8

Table 4. Table 4: Time steps taken for various tests (all with 200 cells) under the ADER-WENO method and the Split-WENO method

	Timesteps (ADER-WENO)	Timesteps (Split-WENO)
Stokes’ First Problem ($\mu ={10}^{-2}$)	385	442
Stokes’ First Problem ($\mu ={10}^{-3}$)	386	443
Stokes’ First Problem ($\mu ={10}^{-4}$)	385	442
Viscous Shock	562	645
Heat Conduction in a Gas	942	1077

Table 5. Table 5: Convergence rates for the Split-WENO method ( N = 2 𝑁 2 N=2 )

Grid Size	$ϵ\left(L_1\right)$	$ϵ\left(L_2\right)$	$ϵ\left(L_{\infty }\right)$	$𝒪\left(L_1\right)$	$𝒪\left(L_2\right)$	$𝒪\left(L_{\infty }\right)$
20	$2.87\times {10}^{-3}$	$7.15\times {10}^{-3}$	$6.21\times {10}^{-2}$
40	$5.81\times {10}^{-4}$	$1.62\times {10}^{-3}$	$1.73\times {10}^{-2}$	2.30	2.14	1.85
60	$1.98\times {10}^{-4}$	$5.39\times {10}^{-4}$	$5.94\times {10}^{-3}$	2.65	2.70	2.63
80	$1.23\times {10}^{-4}$	$3.47\times {10}^{-4}$	$3.41\times {10}^{-3}$	1.67	1.52	1.92

Table 6. Table 6: Convergence rates for the Split-WENO method ( N = 3 𝑁 3 N=3 )

Grid Size	$ϵ\left(L_1\right)$	$ϵ\left(L_2\right)$	$ϵ\left(L_{\infty }\right)$	$𝒪\left(L_1\right)$	$𝒪\left(L_2\right)$	$𝒪\left(L_{\infty }\right)$
10	$1.01\times {10}^{-2}$	$2.58\times {10}^{-2}$	$1.27\times {10}^{-1}$
20	$1.68\times {10}^{-3}$	$4.02\times {10}^{-3}$	$2.93\times {10}^{-2}$	2.59	2.68	2.11
30	$5.34\times {10}^{-4}$	$1.57\times {10}^{-3}$	$1.70\times {10}^{-2}$	2.83	2.32	1.34
40	$3.32\times {10}^{-4}$	$8.94\times {10}^{-4}$	$7.55\times {10}^{-3}$	1.65	1.95	2.82

Table 7. Table 7: Convergence rates for the ADER-DG PNPM method ( N , M = 2 𝑁 𝑀 2 N,M=2 )

Grid Size	$ϵ\left(L_1\right)$	$ϵ\left(L_2\right)$	$ϵ\left(L_{\infty }\right)$	$𝒪\left(L_1\right)$	$𝒪\left(L_2\right)$	$𝒪\left(L_{\infty }\right)$
20	$9.44\times {10}^{-3}$	$2.20\times {10}^{-3}$	$2.16\times {10}^{-3}$
40	$1.95\times {10}^{-3}$	$4.50\times {10}^{-4}$	$4.27\times {10}^{-4}$	2.27	2.29	2.34
60	$7.52\times {10}^{-4}$	$1.74\times {10}^{-4}$	$1.48\times {10}^{-4}$	2.35	2.35	2.61
80	$3.72\times {10}^{-4}$	$8.66\times {10}^{-5}$	$7.40\times {10}^{-5}$	2.45	2.42	2.41

Equations

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho v_k\right)}{\partial x_k}$$

$$\frac{\partial \left(\rho v_i\right)}{\partial t}+\frac{\partial (\rho v_i{v}_k+p{\delta }_{ik}-{\sigma }_{ik})}{\partial x_k}$$

$$\frac{\partial A_{ij}}{\partial t}+\frac{\partial \left(A_{ik}{v}_k\right)}{\partial x_j}+v_k\left(\frac{\partial A_{ij}}{\partial x_k}-\frac{\partial A_{ik}}{\partial x_j}\right)$$

$$\frac{\partial \left(\rho J_i\right)}{\partial t}+\frac{\partial \left(\rho J_i{v}_k+T{\delta }_{ik}\right)}{\partial x_k}$$

$$\frac{\partial \left(\rho E\right)}{\partial t}+\frac{\partial \left(\rho E{v}_k+\left(p{\delta }_{ik}-{\sigma }_{ik}\right)v_i+q_k\right)}{\partial x_k}$$

$$p$$

$$\sigma$$

$$T$$

$$\mathbf{q}$$

$$E=E_1\left(\rho ,p\right)+E_2\left(A,\mathbf{J}\right)+E_3\left(\mathbf{v}\right)$$

$$E_1=\frac{p+\gamma p_{\infty }}{\left(\gamma -1\right)\rho }$$

$$E_2=\frac{c_s^2}{4}{\Vert \mathrm{dev}\left(G\right)\Vert }_F^2+\frac{{\alpha }^2}{2}{\Vert \mathbf{J}\Vert }^2$$

$$c_h=\frac{\alpha }{\rho }\sqrt{\frac{T}{c_v}}$$

$$\mathrm{dev}\left(G\right)=G-\frac{1}{3}\mathrm{tr}\left(G\right)I$$

$$E_3=\frac{1}{2}{\Vert \mathbf{v}\Vert }^2$$

$${\theta }_1\left({\tau }_1\right)$$

$${\theta }_2\left({\tau }_2\right)$$

$${\tau }_1$$

$${\tau }_2$$

$$\sigma$$

$$\mathbf{q}$$

$$-\frac{\psi }{{\theta }_1({\tau }_1)}$$

$$-\frac{\rho \mathbf{H}}{{\theta }_2\left({\tau }_2\right)}$$

$$\det \left(A\right)=\frac{\rho }{{\rho }_0}$$

$${\chi }^i=\frac{1}{\Delta x}\left(x-x_i\right)$$

$${\psi }_i\left({\chi }_j\right)={\delta }_{ij}$$

$$\{\begin{array}{ll}{S}_1& =\left\{C_{i-\frac{N}{2}},\dots ,C_{i+\frac{N}{2}}\right\}\\ S_2& =\left\{C_{i-N},\dots ,C_i\right\}\\ S_3& =\left\{C_i,\dots ,C_{i+N}\right\}\end{array}$$

$$\{\begin{array}{ll}{S}_1& =\left\{C_{i-\lfloor \frac{N}{2}\rfloor },\dots ,C_{i+\lceil \frac{N}{2}\rceil }\right\}\\ S_2& =\left\{C_{i-\lceil \frac{N}{2}\rceil },\dots ,C_{i+\lfloor \frac{N}{2}\rfloor }\right\}\\ S_3& =\left\{C_{i-N},\dots ,C_i\right\}\\ S_4& =\left\{C_i,\dots ,C_{i+N}\right\}\end{array}$$

$$\sum _p{\psi }_p\left({\chi }^i\left(x\right)\right){\hat{w}}_p^{ij}$$

$$\frac{1}{\Delta x}{\int }_{x_k}^{x_{k+1}}\sum _p{\psi }_p\left({\chi }^k\left(x\right)\right){\hat{w}}_p^{ij}dx=u_k\forall C_k\in S_j$$

$${\Sigma }_{mn}=\sum _{\alpha =1}^N{\int }_0^1{\psi }_m^{\left(\alpha \right)}{\psi }_n^{\left(\alpha \right)}d\chi$$

$$o_j={\Sigma }_{mn}{\hat{w}}_m^{ij}{\hat{w}}_n^{ij}$$

$$w_i\left(x\right)=\sum _p{\psi }_p\left({\chi }^i\left(x\right)\right){\bar{w}}_p^i$$

$${\omega }_j=\frac{{\tilde{\omega }}_j}{{\sum }_k{\tilde{\omega }}_k}{\tilde{\omega }}_j=\frac{{\zeta }_j}{{\left(o_j+\epsilon \right)}^r}$$

$${\upsilon }^i=\frac{1}{\Delta y}\left(y-y_i\right)$$

$$\sum _{p,q}{\psi }_p\left({\chi }^i\left(x\right)\right){\psi }_q\left({\upsilon }^i\left(x\right)\right){\tilde{w}}_{pq}^{ij}$$

Full text

††journal: Journal of Computational Physics

1. 1 Background

2. 1.1 Motivation

3. 1.2 The GPR Model

4. 1.3 The ADER-WENO Method

5. 2 An Alternative Numerical Scheme

6. 2.1 The Homogeneous System

7. 2.2 The Temporal ODEs

8. 3 GPR-Specific Performance Improvements

9. 3.1 The Thermal Impulse ODEs

10. 3.2 The Distortion ODEs

11. 3.2.1 Reduced Distortion ODEs

12. 3.2.2 Bounds on Reduced Distortion ODEs

13. 3.2.3 Analytical Approximation

14. 4 Numerical Results

15. 4.1 Strain Relaxation

16. 4.2 Stokes’ First Problem

17. 4.3 Viscous Shock

18. 4.4 Heat Conduction in a Gas

19. 4.5 Speed

20. 4.6 Convergence

21. 5 Conclusions

22. 6 References

23. 7 Acknowledgments

24. 8 Appendix

25. 8.1 Jacobian of Distortion ODEs

26. 8.2 Jacobian of Thermal Impulse ODEs

27. 8.3 WENO Matrices for $N=2$

1 Background

1.1 Motivation

The Godunov-Peshkov-Romenski model of continuum mechanics (as described in 1.2) presents an exciting possibility of being able to describe both fluids and solids within the same mathematical framework. This has the potential to streamline development of simulation software by reducing the number of different systems of equations that require solvers, and cutting down on the amount of theoretical work required, for example in the treatment of interfaces in multimaterial problems. In addition to this, the hyperbolic nature of the GPR model ensures that the nonphysical instantaneous transmission of information appearing in certain non-hyperbolic models (such as the Navier-Stokes equations) cannot occur. Parallelization also tends to be easier with hyperbolic models, allowing us to leverage the great advances that have been made in parallel computing architectures in recent years.

At the time of writing, the GPR model has been solved for a variety of fluid and solid problems using the ADER-WENO method (Dumbser et al. [8], Boscheri et al. [4]). ADER-WENO methods (described in 1.3) are extremely effective in producing arbitrarily-high order solutions to hyperbolic systems of PDEs, but in some situations their accompanying computational cost may prove burdensome. A new method is presented in this study that is simple to implement and computationally cheaper than a corresponding ADER-WENO method if only second order accuracy is required. This may prove useful in the design of simulation software addressing problems in which not just accuracy but also speed and usability are of paramount importance.

1.2 The GPR Model

The GPR model, first introduced in Peshkov and Romenski [23], has its roots in Godunov and Romenski’s 1970s model of elastoplastic deformation (see Godunov and Romenski [14]). It was expanded upon in Dumbser et al. [8] to include thermal conduction. This expanded model takes the following form:

[TABLE]

[TABLE]

$\rho$,$\mathbf{v}$,$p$,$\delta$,$\sigma$,$T$,$E$,$\boldsymbol{q}$ retain their usual meanings. $\theta_{1}$ and $\theta_{2}$ are positive scalar functions, chosen according to the properties of the material being modeled. $A$ is the distortion tensor (containing information about the deformation and rotation of material elements), $\mathbf{J}$ is the thermal impulse vector (a thermal analogue of momentum), $\tau_{1}$ is the strain dissipation time, and $\tau_{2}$ is the thermal impulse relaxation time. $\psi=\frac{\partial E}{\partial A}$ and $\boldsymbol{H}=\frac{\partial E}{\partial\boldsymbol{J}}$.

The following definitions are given:

[TABLE]

[TABLE]

To close the system, the equation of state (EOS) must be specified, from which the above quantities and the sources can be derived. $E$ is the sum of the contributions of the energies at the molecular scale (microscale), the material element111The concept of a material element corresponds to that of a fluid parcel from fluid dynamics, applied to both fluids and solids. scale (mesoscale), and the flow scale (macroscale):

[TABLE]

The EOS used in this study (and described in the following passages) is taken from Dumbser et al. [8]. It should be noted, however, that this is just one particular choice, and there are many others that may be used.

For an ideal or stiffened gas, $E_{1}$ is given by:

[TABLE]

where $p_{\infty}=0$ for an ideal gas.

$E_{2}$ is chosen to have the following quadratic form:

[TABLE]

$c_{s}$ is the characteristic velocity of propagation of transverse perturbations. $\alpha$ is a constant related to the characteristic velocity of propagation of heat waves:

[TABLE]

$G=A^{T}A$ is the Gramian matrix of the distortion tensor, and $\operatorname{dev}\left(G\right)$ is the deviator (trace-free part) of $G$:

[TABLE]

$E_{3}$ is the usual specific kinetic energy per unit mass:

[TABLE]

The following forms are chosen:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The justification of these choices is that classical Navier–Stokes–Fourier theory is recovered in the stiff limit $\tau_{1},\tau_{2}\rightarrow 0$ (see Dumbser et al. [8]). This results in the following relations:

[TABLE]

[TABLE]

The following constraint also holds (see Peshkov and Romenski [23]):

[TABLE]

The GPR model and Godunov and Romenski’s 1970s model of elastoplastic deformation in fact relie upon the same equations. The realization of Peshkov and Romenski was that these are the equations of motion for an arbitrary continuum - not just a solid - and so the model can be applied to fluids too. Unlike in previous continuum models, material elements have not only finite size, but also internal structure, encoded in the distortion tensor.

The strain dissipation time $\tau_{1}$ of the HPR model is a continuous analogue of Frenkel’s “particle settled life time” Frenkel [12]; the characteristic time taken for a particle to move by a distance of the same order of magnitude as the particle’s size. Thus, $\tau_{1}$ characterizes the time taken for a material element to rearrange with its neighbors. $\tau_{1}=\infty$ for solids and $\tau_{1}=0$ for inviscid fluids. It is in this way that the HPR model seeks to describe all three major phases of matter, as long as a continuum description is appropriate for the material at hand.

The evolution equation for $\boldsymbol{J}$ and its contribution to the energy of the system are derived from Romenski’s model of hyperbolic heat transfer, originally proposed in Malyshev and Romenskii [19], Romenski [26], and implemented in Romenski et al. [25, 24]. In this model, $\boldsymbol{J}$ is effectively defined as the variable conjugate to the entropy flux, in the sense that the latter is the derivative of the specific internal energy with respect to $\boldsymbol{J}$. Romenski remarks that it is more convenient to evolve $\boldsymbol{J}$ and $E$ than the heat flux or the entropy flux, and thus the equations take the form given here. $\tau_{2}$ characterizes the speed of relaxation of the thermal impulse due to heat exchange between material elements.

1.3 The ADER-WENO Method

The ADER-WENO method was used in Dumbser et al. [8], Boscheri et al. [4] to solve the GPR system. It produces arbitrarily high-order solutions to hyperbolic systems of PDEs and has been shown to be particularly effective for a wide range of systems (e.g. the classical Euler equations of gas dynamics, the special relativistic hydrodynamics and ideal magnetohydrodynamics equations, and the Baer-Nunziato model for compressible two-phase flow - see Balsara et al. [1], Zanotti and Dumbser [28]). The first step in the process - the WENO method - will be used later in this study and is therefore discussed in detail here. The remaining steps are described qualitatively, with references for further information given.

WENO (Weighted Essentially Non-Oscillatory) methods are used to produce high order polynomial approximations to piece-wise constant data. Many variations exist. In this study, the method of Dumbser et al. [11] is used.

Consider the domain $\left[0,L\right]$. Take $K,N\in\mathbb{N}$. The order of accuracy of the resulting method will be $N+1$. Take the set of grid points $x_{i}=\frac{i\cdot L}{K}$ for $i=0,\ldots,K$ and let $\Delta x=\frac{L}{K}$. Denote cell $\left[x_{i},x_{i+1}\right]$ by $C_{i}$. Given cell-wise constant data $u$ on $\left[0,L\right]$, an order $N$ polynomial reconstruction of $u$ in $C_{i}$ will be performed. Define the scaled space variable:

[TABLE]

Denoting the Gauss-Legendre abscissae on $\left[0,1\right]$ by $\left\{\chi_{0},\ldots,\chi_{N}\right\}$, define the nodal basis of order $N$: the Lagrange interpolating polynomials $\left\{\psi_{0},\ldots,\psi_{N}\right\}$ with the following property:

[TABLE]

If $N$ is even, take the stencils:

[TABLE]

If $N$ is odd, take the stencils:

[TABLE]

The data is reconstructed on $S_{j}$ as:

[TABLE]

where the $\hat{w}_{p}^{ij}$ are solutions to the following linear system:

[TABLE]

where $u_{k}$ is the value of $u$ in $C_{k}$. This can be written as $M_{j}\boldsymbol{\hat{w}^{ij}}=\boldsymbol{u_{\left[j_{0}:j_{N}\right]}}$ where $\left\{j_{0},\ldots,j_{N}\right\}$ indexes the cells in $S_{j}$. In this study reconstructions with $N=2$ are used. The matrices of these linear systems are given in 8.3, along with their inverses, which are precomputed to accelerate the solution of these systems.

Define the oscillation indicator matrix:

[TABLE]

and the oscillation indicator for each stencil:

[TABLE]

The full reconstruction in $C_{i}$ is:

[TABLE]

where $\bar{w}_{p}^{i}=\omega_{j}\hat{w}_{p}^{ij}$ is the weighted coefficient of the $p$th basis function, with weights:

[TABLE]

In this study, $r=8$, $\varepsilon=10^{-14}$, $\zeta_{j}=10^{5}$ if $S_{j}$ is a central stencil, and $\zeta_{j}=1$ if $S_{j}$ is a side stencil, as in Dumbser et al. [7].

The reconstruction can be extended to two dimensions by taking:

[TABLE]

and defining stencils in the y-axis in an analogous manner. The data in $C_{i}$ is then reconstructed using stencil $S_{j}$ as:

[TABLE]

where the coefficients of the weighted 1D reconstruction are used as cell averages:

[TABLE]

The oscillation indicator is calculated for each $p$ in the same manner as the 1D case. The reconstruction method is easily further extensible to three dimensions, now using the coefficients $\bar{w}_{pq}$ of the weighted 2D reconstruction as cell averages.

The next process in the ADER-WENO method is to perform a Continuous Galerkin or Discontinuous Galerkin spatio-temporal polynomial reconstruction of the data in each cell, using the WENO reconstruction as initial data at the start of the time step (see Balsara et al. [1] and Dumbser et al. [5] respectively for implementations of these two variations). The order of this reconstruction in time is usually taken to be the same as the spatial order, and the same basis polynomials are used. The process involves finding the root of a non-linear system, and this process is guaranteed to converge in exact arithmetic for certain classes of PDEs (see Jackson [16]). This root finding can be computationally expensive relative to the WENO reconstruction, especially if the source terms of the PDE system are stiff.

The final step in the ADER-WENO method is to perform a finite volume update of the data in each cell, using the boundary-extrapolated values of the cell-local Galerkin reconstructions to calculate the flux terms, and the interior values of the Galerkin reconstructions to calculate the interior volume integrals. See Dumbser et al. [7] for more details.

2 An Alternative Numerical Scheme

Note that (1a), (1b), (1c), (1d), (1e) can be written in the following form:

[TABLE]

As described in Toro [27], a viable way to solve inhomogeneous systems of PDEs is to employ an operator splitting. That is, the following subsystems are solved:

[TABLE]

[TABLE]

The advantage of this approach is that specialized solvers can be employed to compute the results of the different subsystems. Let $H^{\delta t},S^{\delta t}$ be the operators that take data $\boldsymbol{Q}\left(x,t\right)$ to $\boldsymbol{Q}\left(x,t+\delta t\right)$ under systems (27a) and (27b) respectively. A second-order scheme (in time) for solving the full set of PDEs over time step $\left[0,\Delta t\right]$ is obtained by calculating $\boldsymbol{Q_{\Delta t}}$ using a Strang splitting:

[TABLE]

In the scheme proposed here, the homogeneous subsystem will be solved using a WENO reconstruction of the data, followed by a finite volume update, and the temporal ODEs will be solved with appropriate ODE solvers. This new scheme will be referred to here as the Split-WENO method.

2.1 The Homogeneous System

A WENO reconstruction of the cell-averaged data is performed at the start of the time step (as described in 1.3). Focusing on a single cell $C_{i}$ at time $t_{n}$, we have $\boldsymbol{w^{n}}\left(\boldsymbol{x}\right)=\boldsymbol{w^{n}}_{p}\Psi_{p}\left(\boldsymbol{\chi}\left(\boldsymbol{x}\right)\right)$ in $C_{i}$ where $\Psi_{p}$ is a tensor product of basis functions in each of the spatial dimensions. The flux in $C$ is approximated by $\boldsymbol{F}\left(\boldsymbol{x}\right)\approx\boldsymbol{F}\left(\boldsymbol{w}_{p}\right)\Psi_{p}\left(\boldsymbol{\chi}\left(\boldsymbol{x}\right)\right)$. $\boldsymbol{w}_{p}$ are stepped forwards half a time step using the update formula:

[TABLE]

i.e.

[TABLE]

where $\boldsymbol{\chi_{p}}$ is the node corresponding to $\Psi_{p}$. This evolution to the middle of the time step is similar to that used in the second-order MUSCL and SLIC schemes (see Toro [27]) and, as with those schemes, it is integral to giving the method presented here its second-order accuracy.

Integrating (27a) over $C$ gives:

[TABLE]

where

[TABLE]

[TABLE]

where $V$ is the volume of $C$ and $\boldsymbol{Q^{-},Q^{+}}$ are the interior and exterior extrapolated states at the boundary of $C$, respectively.

Note that (27a) can be rewritten as:

[TABLE]

where $\boldsymbol{M}=\frac{\partial\boldsymbol{F}}{\partial\boldsymbol{Q}}+\boldsymbol{B}$. Let $\boldsymbol{n}$ be the normal to the boundary at point $\boldsymbol{s}\in\partial C$. For the GPR model, $\hat{M}=\boldsymbol{M}\left(\boldsymbol{Q}\left(\boldsymbol{s}\right)\right)\cdot\boldsymbol{n}$ is a diagonalizable matrix with decomposition $\hat{M}=\hat{R}\hat{\Lambda}\hat{R}^{-1}$ where the columns of $\hat{R}$ are the right eigenvectors and $\hat{\Lambda}$ is the diagonal matrix of eigenvalues. Define also $\boldsymbol{\hat{F}}=\boldsymbol{F}\cdot\boldsymbol{n}$ and $\hat{B}=\boldsymbol{B}\cdot\boldsymbol{n}$. Using these definitions, the interface terms arising in the FV formula have the following form:

[TABLE]

$\tilde{M}$ is chosen to either correspond to a Rusanov/Lax-Friedrichs flux (see Toro [27]):

[TABLE]

or a simplified Osher–Solomon flux (see Dumbser and Toro [9, 10]):

[TABLE]

where

[TABLE]

$\tilde{B}$ takes the following form:

[TABLE]

It was found that the Osher-Solomon flux would often produce slightly less diffusive results, but that it was more computationally expensive, and also had a greater tendency to introduce numerical artefacts.

$\boldsymbol{P_{i}^{n+\frac{1}{2}}},\boldsymbol{D_{i}^{n+\frac{1}{2}}}$ are calculated using an $N+1$-point Gauss-Legendre quadrature, replacing $\boldsymbol{Q}\left(\boldsymbol{x},t_{n+\frac{1}{2}}\right)$ with $\boldsymbol{w^{n+\frac{1}{2}}}\left(\boldsymbol{x}\right)$.

2.2 The Temporal ODEs

Noting that $\frac{d\rho}{dt}=0$ over the ODE time step, the operator $S$ entails solving the following systems:

[TABLE]

[TABLE]

These systems can be solved concurrently with a stiff ODE solver. The Jacobians of these two systems to be used in an ODE solver are given in 8.1 and 8.2. However, these systems can also be solved separately, using the analytical results presented in Section 3, under specific assumptions. The second-order Strang splitting is then:

[TABLE]

where $D^{\delta t},T^{\delta t}$ are the operators solving the distortion and thermal impulse ODEs respectively, over timestep $\delta t$. This allows us to bypass the relatively computationally costly process of solving these systems numerically.

3 GPR-Specific Performance Improvements

3.1 The Thermal Impulse ODEs

Taking the EOS for the GPR model (3) and denoting by $E_{2}^{\left(A\right)},E_{2}^{\left(J\right)}$ the components of $E_{2}$ depending on $A$ and $\boldsymbol{J}$ respectively, we have:

[TABLE]

where:

[TABLE]

[TABLE]

Over the time period of the ODE (39b), $c_{1},c_{2}>0$ are constant. We have:

[TABLE]

Therefore:

[TABLE]

where

[TABLE]

[TABLE]

Note that this is a generalized Lotka-Volterra system in $\left\{J_{1}^{2},J_{2}^{2},J_{3}^{2}\right\}$. It has the following analytical solution:

[TABLE]

3.2 The Distortion ODEs

3.2.1 Reduced Distortion ODEs

Let $k_{0}=\frac{3}{\tau_{1}}\left(\frac{\rho}{\rho_{0}}\right)^{\frac{5}{3}}>0$ and let $A$ have singular value decomposition $U\Sigma V^{T}$. Then:

[TABLE]

[TABLE]

Therefore:

[TABLE]

It is a common result (see Giles [13]) that:

[TABLE]

and thus:

[TABLE]

Using a fast $3\times 3$ SVD algorithm (such as in McAdams et al. [20]), $U,V,\Sigma$ can be obtained, after which the following procedure is applied to $\Sigma$, giving $A\left(t\right)=U\Sigma\left(t\right)V^{T}$.

Denote the singular values of $A$ by $a_{1},a_{2},a_{3}$. Then:

[TABLE]

Letting $x_{i}=\frac{a_{i}^{2}}{\det\left(A\right)^{\frac{2}{3}}}=\frac{a_{i}^{2}}{\left(\frac{\rho}{\rho_{0}}\right)^{\frac{2}{3}}}$ we have:

[TABLE]

where $\tau=\frac{2}{\tau_{1}}\left(\frac{\rho}{\rho_{0}}\right)^{\frac{7}{3}}t$ and $\bar{x}$ is the arithmetic mean of $x_{1},x_{2},x_{3}$. This ODE system travels along the surface $\Psi=\left\{x_{1},x_{2},x_{3}>0,x_{1}x_{2}x_{3}=1\right\}$ to the point $x_{1},x_{2},x_{3}=1$. This surface is symmetrical in the planes $x_{1}=x_{2}$, $x_{1}=x_{3}$, $x_{2}=x_{3}$. As such, given that the system is autonomous, the paths of evolution of the $x_{i}$ cannot cross the intersections of these planes with $\Psi$. Thus, any non-strict inequality of the form $x_{i}\geq x_{j}\geq x_{k}$ is maintained for the whole history of the system. By considering (53) it is clear that in this case $x_{i}$ is monotone decreasing, $x_{k}$ is monotone increasing, and the time derivative of $x_{j}$ may switch sign.

Note that we have:

[TABLE]

Thus, an ODE solver can be used on these two equations to effectively solve the ODEs for all 9 components of $A$. Note that:

[TABLE]

This has solution:

[TABLE]

where

[TABLE]

In the case that $x_{i,0}=x_{j,0}$, we have $x_{i}=x_{j}$ for all time. Thus, the ODE system for $A$ has been reduced to a single ODE, as $x_{j}\left(x_{i}\right)$ can be inserted into the RHS of the equation for $\frac{dx_{i}}{d\tau}$. However, it is less computationally expensive to evolve the system presented in (54).

3.2.2 Bounds on Reduced Distortion ODEs

If any of the relations in $x_{i}\geq x_{j}\geq x_{k}$ are in fact equalities, equality is maintained throughout the history of the system. This can be seen by noting that the time derivatives of the equal variables are in this case equal. If $x_{j}=x_{k}$ then $x_{i}=\frac{1}{x_{j}^{2}}$. Combining these results, the path of the system in $\left(x_{i},x_{j}\right)$ coordinates is in fact confined to the curved triangular region:

[TABLE]

This is demonstrated in Figure 1 on page 1. By (54), the rate of change of $x_{i}$ at a particular value $x_{i}=x_{i}^{*}$ is given by:

[TABLE]

Note that:

[TABLE]

[TABLE]

Thus, $x_{i}$ decreases fastest on the line $x_{i}=\frac{1}{x_{j}^{2}}$ (the bottom boundary of the region given in Figure 1 on page 1), and slowest on the line $x_{i}=x_{j}$. The rates of change of $x_{i}$ along these two lines are given respectively by:

[TABLE]

[TABLE]

These have implicit solutions:

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

As (53) is an autonomous system of ODEs, it has the property that its limit $x_{1}=x_{2}=x_{3}=1$ is never obtained in finite time, in precise arithmetic. In floating point arithmetic we may say that the system has converged when $x_{i}-1<\epsilon$ (machine epsilon) for each $i$. This happens when:

[TABLE]

This provides a quick method to check whether it is necessary to run the ODE solver in a particular cell. If the following condition is satisfied then we know the system in that cell converges to the ground state over the time interval in which the ODE system is calculated:

[TABLE]

If the fluid is very inviscid, resulting in a stiff ODE, the critical time is lower, and there is more chance that the ODE system in the cell reaches its limit in $\Delta t$. This check potentially saves a lot of computationally expensive stiff ODE solves. The same goes for if the flow is slow-moving, as the system will be closer to its ground state at the start of the time step and is more likely to converge over $\Delta t$. Similarly, if the following condition is satisfied then we know for sure that an ODE solver is necessary, as the system certainly will not have converged over the timestep:

[TABLE]

3.2.3 Analytical Approximation

We now explore cases when even the reduced ODE system (54) need not be solved numerically. Define the following variables:

[TABLE]

[TABLE]

It is a standard result that $m\geq\sqrt[3]{x_{1}x_{2}x_{3}}$. Thus, $m\geq 1$. Note that $u$ is proportional to the internal energy contribution from the distortion. From (53) we have:

[TABLE]

[TABLE]

Combining these equations, we have:

[TABLE]

Therefore:

[TABLE]

We make the following assumption, noting that it is true in all physical situations tested in this study:

[TABLE]

Thus, we have the linearized ODE:

[TABLE]

This is a Sturm-Liouville equation with solution:

[TABLE]

Thus, we also have:

[TABLE]

Once $m_{\Delta t}=1+\eta\left(\frac{2}{\tau_{1}}\left(\frac{\rho}{\rho_{0}}\right)^{\frac{7}{3}}\Delta t\right)$ and $u_{\Delta t}=u\left(\frac{2}{\tau_{1}}\left(\frac{\rho}{\rho_{0}}\right)^{\frac{7}{3}}\Delta t\right)$ have been found, we have:

[TABLE]

[TABLE]

This gives:

[TABLE]

[TABLE]

where

[TABLE]

Note that taking the real parts of the above expression for $x_{i}$ gives:

[TABLE]

[TABLE]

At this point it is not clear which values of $\left\{x_{i},x_{j},x_{k}\right\}$ are taken by $x_{1},x_{2},x_{3}$. However, this can be inferred from the fact that any relation $x_{i}\geq x_{j}\geq x_{k}$ is maintained over the lifetime of the system. Thus, the stiff ODE solver has been obviated by a few arithmetic operations.

4 Numerical Results

4.1 Strain Relaxation

In this section, the approximate analytic solver for the distortion ODEs, presented in 3.2.3, is compared with a numerical ODE solver. Initial data was taken from Barton and Drikakis [2]:

[TABLE]

Additionally, the following parameter values were used: $\rho_{0}=1,c_{s}=1,\mu=10^{-2}$, giving $\tau_{1}=0.06$. As can be seen in Figure 2 on page 2, Figure 3 on page 3, and Figure 4 on page 4, the approximate analytic solver compares well with the numerical solver in its results for the distortion tensor $A$, and thus also the internal energy and stress tensor. The numerical ODE solver was the odeint solver from SciPy 0.18.1, based on the LSODA solver from the FORTRAN library ODEPACK (see Oliphant [22]).

4.2 Stokes’ First Problem

This problem is one of the few test cases with an analytic solution for the Navier-Stokes equations. It consists of two ideal gases in an infinite domain, meeting at the plane $x=0$, initially flowing with equal and opposite velocity $\pm 0.1$ in the $y$-axis. The initial conditions are given in Table 1 on page 1.

The flow has a low Mach number of 0.1, and this test case is designed to demonstrate the efficacy of the numerical methods in this flow regime. The exact solution to the Navier-Stokes equations is given by222In this problem, the Navier-Stokes equations reduce to $v_{t}=\mu v_{xx}$. Defining $\eta=\frac{x}{2\sqrt{\mu t}}$, and assuming $v=f\left(\eta\right)$, this becomes $f^{{}^{\prime\prime}}+2\eta f^{{}^{\prime}}=0$. The result follows by solving this equation with the boundary conditions $v\left(\pm\infty\right)=\pm v_{0}$.:

[TABLE]

Heat conduction is neglected, and $\gamma=1.4$, $c_{v}=1$, $\rho_{0}=1$, $c_{s}=1$. The viscosity is variously taken to be $\mu=10^{-2}$, $\mu=10^{-3}$, $\mu=10^{-4}$ (resulting in $\tau_{1}=0.06$, $\tau_{1}=0.006$, $\tau_{1}=0.0006$, respectively). Due to the stiffness of the source terms in the equations governing $A$ in the case that $\mu=10^{-4}$, the step (30) in the WENO reconstruction under the Split-WENO method was not performed, and $\boldsymbol{w_{p}^{n+\frac{1}{2}}}\equiv\boldsymbol{w_{p}^{n}}$ was taken instead. This avoided the numerical diffusion that otherwise would have emerged at the interface at $x=0$.

The results of simulations with 200 cells at time $t=1$, using reconstruction polynomials of order $N=2$, are presented in Figure 5 on page 5. The GPR model solved with both the ADER-WENO and Split-WENO methods closely matches the exact Navier-Stokes solution. Note that at $\mu=10^{-2}$ and $\mu=10^{-3}$, the ADER-WENO and Split-WENO methods are almost indistinguishable. At $\mu=10^{-4}$ the Split-WENO method matches the curve of the velocity profile more closely, but overshoots slightly at the boundaries of the center region. This overshoot phenomenon is not visible in the ADER-WENO results.

4.3 Viscous Shock

This test is designed to demonstrate that the numerical methods used are also able to cope with fast flows. First demonstrated by Becker Becker [3], the Navier-Stokes equations have an analytic solution for $P_{r}=0.75$ (see Johnson Johnson [17] for a full analysis). As noted by Dumbser Dumbser et al. [8], if the wave has nondimensionalised upstream velocity $\bar{v}=1$ and Mach number $M_{c}$, then its nondimensionalised downstream velocity is:

[TABLE]

The wave’s velocity profile $\bar{v}\left(x\right)$ is given by the roots of the following equation:

[TABLE]

[TABLE]

$c_{1},c_{2}$ are constants that affect the position of the center of the wave, and its stretch factor, respectively. Following the analysis of Morduchow and Libby Morduchow and Libby [21], the nondimensional pressure and density profiles are given by:

[TABLE]

[TABLE]

To obtain an unsteady shock traveling into a region at rest, a constant velocity field $v=M_{c}c_{0}$ is imposed on the traveling wave solution presented here (where $c_{0}$ is the adiabatic sound speed). Thus, if $p_{0},\rho_{0}$ are the downstream (reference) values for pressure and density:

[TABLE]

[TABLE]

These functions are used as initial conditions, along with $A=\sqrt[3]{\bar{\rho}}I$ and $\boldsymbol{J}=\boldsymbol{0}$. The downstream density and pressure are taken to be $\rho_{0}=1$ and $p_{0}=\frac{1}{\gamma}$ (so that $c_{0}=1$). $M_{c}=2$ and $R_{e}=100$. The material parameters are taken to be: $\gamma=1.4$, $p_{\infty}=0$, $c_{v}=2.5$, $c_{s}=5$, $\alpha=5$, $\mu=2\times 10^{-2}$, $\kappa=\frac{28}{3}\times 10^{-2}$ (resulting in $\tau_{1}=0.0048$, $\tau_{2}=0.00522\dot{6}$).

The results of a simulation with 200 cells at time $t=0.2$, using reconstruction polynomials of order $N=2$, are presented in Figure 6 on page 6 and Figure 7 on page 7. The shock was initially centered at $x=0.25$, reaching $x=0.65$ at the final time. Note that the density, velocity, and pressure results for both methods match the exact solution well, with the ADER-WENO method appearing to produce a slightly more accurate solution. The results for the two methods for the stress tensor and heat flux are close.

4.4 Heat Conduction in a Gas

This is a simple test case to ensure that the heat transfer terms in the implementation are working correctly. Two ideal gases at different temperatures are initially in contact at position $x=0$. The initial conditions for this problem are given in Table 2 on page 2.

The material parameters are taken to be: $\gamma=1.4$, $c_{v}=2.5$, $\rho_{0}=1$, $p_{0}=1$, $c_{s}=1$, $\alpha=2$, $\mu=10^{-2}$, $\kappa=10^{-2}$ (resulting in $\tau_{1}=0.06$, $\tau_{2}=0.0025$). The results of a simulation with 200 cells at time $t=1$, using reconstruction polynomials of order $N=2$, are presented in Figure 8 on page 8. The ADER-WENO and Split-WENO methods are in perfect agreement for both the temperature and heat flux profiles. As demonstrated in Dumbser et al. [8], this means that they in turn agree very well with a reference Navier-Stokes-Fourier solution.

4.5 Speed

Both the ADER-WENO scheme and the Split-WENO scheme used in this study were implemented in Python3. All array functions were precompiled with Numba’s JIT capabilities and the root-finding procedure in the Galerkin predictor was performed using SciPy’s Newton-Krylov solver, compiled against the Intel MKL. Clear differences in computational cost between the ADER-WENO and Split-WENO methods were apparent, as is to be expected, owing to the lack of Galerkin method in the Split-WENO scheme. The wall times for the various tests undertaken in this study are given in Table 3 on page 3, comparing the combined WENO and Galerkin methods of the ADER-WENO scheme to the combined WENO and ODE methods of the Split-WENO scheme. All computations were performed using an Intel Core i7-4910MQ, on a single core. The number of time steps taken are given in Table 4 on page 4. The differences between the methods in terms of the number of time steps taken in each test result from the fact that, for numerical stability, CFL numbers of 0.8 and 0.7 were required by the ADER-WENO method and the Split-WENO method, respectively.

Note that, unlike with the ADER-WENO scheme, the wall time for the Split-WENO scheme is unaffected by a decrease in the viscosity in Stokes’ First Problem (and the corresponding increase in the stiffness of the source terms). This is because the analytic approximation to the distortion ODEs obviates the need for a stiff solver. The large difference in ADER-WENO solver times between the $\mu=10^{-3}$ and $\mu=10^{-4}$ cases is due to the fact that, in the latter case, a stiff solver must be employed for the initial guess to the root of the nonlinear system produced by the Discontinuous Galerkin method (as described in Hidalgo and Dumbser [15]).

4.6 Convergence

To assess the rate of convergence of the Split-WENO method, the convected isentropic vortex convergence study from Dumbser et al. [8] was performed. The initial conditions are given as $\rho=1+\delta\rho$, $p=1+\delta p$, $\boldsymbol{v}=\left(1,1,0\right)+\delta\boldsymbol{v}$, $A=\sqrt[3]{\rho}I$, $\boldsymbol{J}=\boldsymbol{0}$, where:

[TABLE]

[TABLE]

The 2D domain is taken to be $\left[0,10\right]^{2}$. $\epsilon$ is taken to be 5. The material parameters are taken to be: $\gamma=1.4$, $c_{v}=2.5$, $\rho_{0}=1$, $p_{0}=1$, $c_{s}=0.5$, $\alpha=1$, $\mu=10^{-6}$, $\kappa=10^{-6}$ (resulting in $\tau_{1}=2.4\times 10^{-5}$, $\tau_{2}=10^{-6}$). Thus, this can be considered to be a stiff test case.

The convergence rates in the $L_{1}$, $L_{2}$, $L_{\infty}$ norms for the density variable are given in Table 5 on page 5 and Table 6 on page 6 for WENO reconstruction polynomial orders of $N=2$ and $N=3$, respectively. As expected, both sets of tests attain roughly second order convergence. For comparison, the corresponding results for this test from Dumbser et al. [8]

• solved using a third-order P2P2 scheme - are given in Table 7 on page 7 for comparison.

5 Conclusions

In summary, a new numerical method based on an operator splitting, and including some analytical results, has been proposed for the GPR model of continuum mechanics. It has been demonstrated that this method is able to match current ADER-WENO methods in terms of accuracy on a range of test cases. It is significantly faster than the other currently available methods, and it is easier to implement. The author would recommend that if very high order-of-accuracy is required, and computational cost is not important, then ADER-WENO methods may present a better option, as by design the new method cannot achieve better than second-order accuracy. This new method clearly has applications in which it will prove useful, however.

In a similar manner to the operator splitting method presented in Leveque and Yee [18], the Split-WENO method is second-order accurate and stable even for very stiff problems (in particular, the reader is referred to the results of the $\mu=10^{-4}$ variation of Stokes’ First Problem in 4.2 and the convergence study in 4.6). However, it will inevitably suffer from the incorrect speed of propagation of discontinuities on regular, structured grids. This is due to a lack of spatial resolution in evaluating the source terms, as detailed in Leveque and Yee [18]. This issue can be rectified by the use of some form of shock tracking or mesh refinement, as noted in the cited paper. It is noted in Dumbser et al. [6] that operator splitting-based methods can result in schemes that are neither well-balanced nor asymptotically consistent. The extent to which these two conditions are violated by the Split-WENO method - and the severity in practise of any potential violation - is a topic of further research.

It should be noted that the assumption (72) used to derive the approximate analytical solver may break down for situations where the flow is compressed heavily in one direction but not the others. The reason for this is that one of the singular values of the distortion tensor will be much larger than the others, and the mean of the squares of the singular values will not be close to its geometric mean, meaning that the subsequent linearization of the ODE governing the mean of the singular values fails. It should be noted that none of the situations covered in this study presented problems for the approximate analytical solver, and situations which may be problematic are in some sense unusual. In any case, a stiff ODE solver can be used to solve the system (54) if necessary, utilizing the Jacobians derived in the appendix, and so the Split-WENO method is still very much usable in these situations, albeit slightly slower.

It should be noted that both the ADER-WENO and Split-WENO methods, as described in this study, are trivially parallelizable on a cell-wise basis. Thus, given a large number of computational cores, deficiencies in the Split-WENO method in terms of its order of accuracy may be overcome by utilizing a larger number of computational cells and cores. The computational cost of each time step is significantly smaller than with the ADER-WENO method, and the number of grid cells that can be used scales roughly linearly with number of cores, at constant time per iteration.

6 References

7 Acknowledgments

I acknowledge financial support from the EPSRC Centre for Doctoral Training in Computational Methods for Materials Science under grant EP/L015552/1.

8 Appendix

8.1 Jacobian of Distortion ODEs

The Jacobian of the source function is used to speed up numerical integration of the ODE. It is derived thus:

[TABLE]

Thus:

[TABLE]

Thus:

[TABLE]

where $G^{{}^{\prime}}=AA^{T}$ and $X^{a,b}$ refers to tensor $X$ with indices $a,b$ transposed.

8.2 Jacobian of Thermal Impulse ODEs

As demonstrated in 3.1, we have:

[TABLE]

where

[TABLE]

[TABLE]

Thus, the Jacobian of the thermal impulse ODEs is:

[TABLE]

8.3 WENO Matrices for $N=2$

[TABLE]

[TABLE]

[TABLE]

[TABLE]