2.1. 3D Euler fluid equations

Let us consider the 3D incompressible Euler equations for the velocity field $\boldsymbol{u}(x,y,z,t)\in \mathbb{R}^{3}$ defined for $(x,y,z)\in \mathbb{R}^{3}$ and in a time interval $t\in [0,T)$ :

We define the vorticity field as the three-dimensional vector field ${\bf\omega}\equiv \boldsymbol{{\rm\nabla}}\wedge \boldsymbol{u}.$

2.2. Symmetry plane

We consider a special configuration of the fields and define a symmetry plane at $z=0$ by the following conditions on the velocity and pressure fields:

for arbitrary $(x,y,z)\in \mathbb{R}^{3}$ and $t\in [0,T).$ It is easy to show that these conditions are consistent with the evolution (2.1a,b ). As for the vorticity, these conditions imply

for all $(x,y,z)\in \mathbb{R}^{3}$ and $t\in [0,T).$

At the symmetry plane $z=0$ , the 3D Euler fluid equations will simplify because: (i) $u_{z}(x,y,0,t)\equiv 0$ so the velocity field is parallel to the plane; (ii) ${\it\omega}_{x}(x,y,0,t)={\it\omega}_{y}(x,y,0,t)=0$ so the vorticity field is perpendicular to the plane. This leads to a new system of equations which is ‘almost’ 2D (except for a pressure term depending on the full 3D velocity field). Let us denote the ‘horizontal’ component of the velocity field and the pressure at the symmetry plane by

The horizontal components of (2.1a,b ) become, at $z=0$ ,

where $\boldsymbol{{\rm\nabla}}_{h}$ denotes the ‘horizontal’ gradient operator, $\boldsymbol{{\rm\nabla}}_{h}=(\partial _{x},\partial _{y})$ . The incompressibility condition in (2.1a,b ) allows us to define the stretching-rate scalar on the symmetry plane:

Therefore, even though $u_{z}=0$ at the symmetry plane $z=0$ , we have $u_{z,z}\neq 0$ at $z=0.$ Let us compute the $z$ -derivative of the $z$ -component of (2.1a,b ) and then evaluate at $z=0$ . We obtain $(\partial {\it\gamma}/\partial t)+\partial _{z}(u_{x}u_{z,x}+u_{y}u_{z,y}+u_{z}u_{z,z})|_{z=0}=-{p_{,zz}|}_{z=0}.$ This can be simplified by noticing that the symmetry-plane conditions (2.2) imply $u_{z}=u_{z,zz}=0$ and $u_{x,z}=u_{y,z}=0$ at $z=0$ so that we are left with

Although this equation does not close (the pressure still depends on the full 3D velocity profile), it is remarkable that the equation for the vorticity at the symmetry plane does close. As mentioned before, the vorticity at the symmetry plane has no horizontal component so we can define the vorticity scalar

An evolution equation for vorticity is obtained by taking the curl of the 2D equations (2.5). We readily obtain

which explains the meaning of ${\it\gamma}$ as the stretching rate of vorticity.

Taken together, equations (2.6)–(2.9) would be a closed system in two dimensions, except for the pressure term which, as usual, depends on the full 3D velocity profile. An important consistency condition on the pressure term is derived after integrating spatially (2.7) over the whole horizontal domain, and discarding boundary terms by assuming either periodic or vanishing boundary conditions on the horizontal velocity field. The condition reads

In this paper we propose a consistent family of closure models for the pressure term, based on an exact solution by Gibbon et al. (Reference Gibbon, Fokas and Doering1999). These models will be discussed in the following subsections.

From here on, we will assume periodic boundary conditions in the two spatial directions $(x,y)$ , for the three-dimensional velocity field components $\boldsymbol{u}=(u_{x},u_{y},u_{z})$ and the pressure scalar $p,$ with periodicity box $[0,2{\rm\pi}]\times [0,2{\rm\pi}]$ . As for the $z$ -direction, we do not assume yet any boundary condition.

2.3. Gibbon et al. exact solution of 3D Euler (and Navier–Stokes)

A general class of exact solutions of 3D Euler (and Navier–Stokes) was presented by Gibbon et al. (Reference Gibbon, Fokas and Doering1999). In the case of the 3D Euler equations in the presence of a symmetry plane, this exact solution becomes

where $\boldsymbol{u}_{h}$ and ${\it\gamma}$ satisfy (2.6)–(2.9), along with the following condition on the pressure:

(a function of time only). Due to the periodicity of the horizontal domain, condition (2.10) implies the closure $f(t)=-2\langle {\it\gamma}^{2}\rangle ,$ where

Correspondingly, (2.6)–(2.9) determine the fate of the full 3D flow via the knowledge of the scalars ${\it\gamma}$ and ${\it\omega}$ . Remarkably, along the characteristics of the horizontal velocity field $\boldsymbol{u}_{h}$ the equation for stretching rate ${\it\gamma}$ is ‘decoupled’ from the system and reads

This allowed Constantin (Reference Constantin2000) to solve for ${\it\gamma}$ along characteristics (and for vorticity ${\it\omega}$ , which can be found a posteriori), proving that the stretching rate ${\it\gamma}$ would blow up in a finite time, with explicit formulae for the singularity time which confirmed the accuracy of the numerical blowup predictions in Ohkitani & Gibbon (Reference Ohkitani and Gibbon2000). We also note here for completeness the work of Gibbon, Moore & Stuart (Reference Gibbon, Moore and Stuart2003) who also proved blowup in the axisymmetric case.

2.4. A physically motivated model on the symmetry plane

The main drawback of the exact-solution ansatz $p_{,zz}=f(t)$ is that it has infinite energy (pressure goes like $z^{2}$ ). In particular, the exact-solution ansatz is not suitable for 3D numerical simulations on fully periodic domains, where the $z$ coordinate is also periodic. One can cure this drawback by re-interpreting the ansatz as a closure on the symmetry plane: $p_{,zz}|_{z=0}=f(t)$ and dealing with a model rather than an exact solution. Importantly, condition (2.10) must be met at all times. However, this closure ansatz has, by definition, a spatially uniform pressure curvature $p_{,zz}|_{z=0}$ on the symmetry plane, a feature that is not observed in numerical simulations. Therefore, an improved type of ansatz is required, satisfying condition (2.10) while at the same time showing a spatial dependence of $p_{,zz}|_{z=0}$ on the symmetry plane. We propose a one-parameter family of models which achieves all of this, while still keeping the structure of characteristics, so blowup can be assessed analytically. Looking at (2.6)–(2.9), the model is defined by the closure

where ${\it\lambda}$ is a real (free) parameter. With this closure, the family of models corresponds to the following equations on the symmetry plane:

where $(x,y)\in \mathbb{T}^{2}\equiv [0,2{\rm\pi}]\times [0,2{\rm\pi}]$ and ${\it\lambda}\in \mathbb{R}$ . The case ${\it\lambda}=0$ recovers the equations of Gibbon et al. (Reference Gibbon, Fokas and Doering1999). The horizontal velocity field $\boldsymbol{u}_{h}=(u_{x},u_{y})$ is defined, as usual, by (2.6) and (2.8). The evolution equations (2.16) and (2.17) are supplemented by initial conditions ${\it\gamma}(x,y,0)={\it\gamma}_{0}(x,y)$ and ${\it\omega}(x,y,0)={\it\omega}_{0}(x,y),$ which must have zero mean: $\langle {\it\gamma}_{0}\rangle =\langle {\it\omega}_{0}\rangle =0$ .

We must stress that we do not know whether this family of models could be extended to three dimensions and become a valid solution of the full 3D Euler equations; this is a matter of further investigation. What we do know is that any smooth solution of 3D Euler with discrete mirror symmetry has a 2D plane (the symmetry plane), where vorticity is a scalar satisfying (2.17), even in the case of finite energy. The effective unknown is the governing equation for vorticity stretching rate at the symmetry plane. We motivate our family of models by the physical interpretation of the pressure term at the symmetry plane and the introduction of a tuneable parameter, which provides a range of phenomenology and asymptotic behaviours. This approach provides a valuable framework for developing methods for assessing blowup.

3.1. Blowup solutions

The solutions along characteristics for stretching rate (3.2) and vorticity (3.3) will develop a singularity if the factor $1+({\it\lambda}+1){\it\gamma}_{0}(X_{0},Y_{0})S(t)$ becomes zero for some time $t$ and position $(X_{0},Y_{0})$ . Since $S(0)=0$ and ${\dot{S}}(t)\geqslant 0,$ it follows that $S(t)$ can only grow in time and thus the singularity will occur first at the characteristic starting at the position of the infimum (if ${\it\lambda}>-1$ ) or the supremum (if ${\it\lambda}<-1$ ) of ${\it\gamma}_{0}$ over $\mathbb{T}^{2}.$ Consequently, the singularity time $T^{\ast }$ is defined by the condition $S(T^{\ast })=S^{\ast },$ where $S^{\ast }({>}0)$ is defined by

From the ODE (3.4) an explicit formula for the singularity time $T^{\ast }$ is derived:

We consider briefly the blowup structure for the stretching and vorticity solutions for different values of the model’s free parameter ${\it\lambda}$ (a detailed explanation will be presented in a forthcoming paper). The parameter space is divided in regions of finite-time blowup alternating with regions of infinite-time blowup as illustrated in figure 1. The region of ${\it\lambda}$ where ${\dot{S}}(t)=0$ depends on the initial conditions; if the local profile of initial stretching near the infimum is parabolloidal (a generic situation), then the regions are as depicted as in figure 1.

3.2. Motivation for the choice of parameter ${\it\lambda}=-3/2$ and generic initial conditions

In this paper we will focus on one particular choice of parameter: ${\it\lambda}=-3/2.$ This choice has several advantages.

1. (i) In this case the evolution equation for $S(t)$ is of the form

   (3.10) $$\begin{eqnarray}\displaystyle {\dot{S}}=1+a^{2}S^{2},\quad S(0)=0, & & \displaystyle\end{eqnarray}$$
   which can be solved analytically for any initial condition ( $a$ depends on the initial condition), leading to explicit analytical solutions in terms of trigonometric functions for all of the blowup quantities. The utility of this is that we can perform direct comparisons between theory and numerics.

2. (ii) This case gives finite-time singularity with ${\dot{S}}(T^{\ast })<\infty ,$ leading to simple asymptotic expressions for the blowup quantities. The singularity is controlled by the supremum of ${\it\gamma}_{0}$ and leads to a blowup of the form $\sup {\it\gamma}(x,y,t)\sim 2(T^{\ast }-t)^{-1}$ when $t$ is close to the singularity time $T^{\ast }.$ This feature is analogous to what is normally expected in a 3D Euler fluid simulation in a potential singularity scenario.

3. (iii) In this case there exists a special conserved quantity: $\langle {\it\gamma}^{2}\rangle ,$ which is reminiscent of the ‘energy’ in 2D and 3D ideal models, and provides an opportunity for a simplified analysis of the Fourier spectrum.

We will exploit analytical solutions (3.3) for vorticity and (3.2) for stretching rate in order to validate direct numerical simulations of the system. For example, using the fact that the back-to-labels transformation is bijective for $t<T^{\ast },$ one can use (3.2)–(3.4) to calculate the infimum and supremum of stretching rate and vorticity from the knowledge of the initial conditions. This gives either explicit expressions in terms of simple functions or numerically computable expressions to any desired accuracy. In table 1, we summarise the relevant analytical solutions for the initial condition $u_{x}(x,y,0)=\cos (x)\sin (y)$ , $u_{y}(x,y,0)=\cos (x)+\sin (y)$ or, in terms of stretching rate and vorticity,

Figure 2 shows contour plots of the initial condition for both the vorticity and stretching rate. The discrete symmetry $(x\rightarrow {\rm\pi}+x,y\rightarrow -y,u_{x}\rightarrow u_{x},u_{y}\rightarrow -u_{y})$ of this initial condition is preserved under the time evolution. Thus, we can restrict the analysis to the quadrant ${\rm\pi}\leqslant x\leqslant 2{\rm\pi}$ , ${\rm\pi}\leqslant y\leqslant 2{\rm\pi}.$

Figure 1. Singular/non-singular behaviour of system (3.2)–(3.4) depending on the value of the model’s parameter ${\it\lambda}.$ The limiting case ${\it\lambda}=-1$ (not analysed in this paper) gives an infinite-time singularity.

Figure 2. Stretching rate ( ${\it\gamma}$ ) and vorticity ( ${\it\omega}$ ) plots for ${\it\lambda}=-3/2$ at $t=0$ (a,b), 0.5 (c,d) and 1.0 (e, f) from the original 2D Euler model.

4.1. Using the mapped system to assess blowup of original system

The mapping (4.1), (4.3) and (4.4) is bijective as long as ${\it\tau}<\infty .$ Correspondingly, integration of the mapped system (4.8) and (4.9) should give enough information to assess blowup quantities in the original variables.

The norm $\Vert {\it\gamma}(\cdot ,t)\Vert _{\infty }$ satisfies the ODE (4.6). In terms of ${\it\tau}$ and the mapped stretching rate ${\it\gamma}_{map}(x,y,{\it\tau}),$ this equation reads

where we have used $(\text{d}{\it\tau}/\text{d}t)=\Vert {\it\gamma}(\cdot ,t({\it\tau}))\Vert _{\infty }.$ Correspondingly we obtain, after a simple ${\it\tau}$ integration,

Note that the right-hand side is written entirely in terms of mapped fields. The integrands ${\it\sigma}_{\infty }$ and ${\it\sigma}_{\infty }\langle {\it\gamma}_{map}^{2}\rangle$ are bounded by $1$ so, remarkably, the blowup assessment of the original variables is done in terms of bounded quantities. In particular, this leads to the following general formula for the singularity time $T^{\ast }$ :

Note that this integral converges if and only if the original problem has a finite-time singularity.

From (2.17) it follows that the norm $\Vert {\it\omega}(\cdot ,t)\Vert _{\infty }$ satisfies

where $\boldsymbol{X}_{{\it\omega}}(t)$ is the position at which the maximum of $|{\it\omega}(\boldsymbol{x},t)|$ is attained. In terms of the mapped variables, this reads $(\text{d}\ln \Vert {\it\omega}(\cdot ,t({\it\tau}))\Vert _{\infty }/\text{d}{\it\tau})={\it\gamma}_{map}(\boldsymbol{X}_{{\it\omega}}(t({\it\tau})),{\it\tau}),$ which gives

Remarkably again, the blowup assessment of the original vorticity depends on a bounded quantity ( $-1\leqslant {\it\gamma}_{map}\leqslant 1$ ).

4.2. Analytical solution for the mapped variables in the case ${\it\lambda}=-3/2$

The analytical solution along characteristics (3.2)–(3.4) leads to a connection with the mapped variables. Let us consider the case ${\it\lambda}=-3/2$ and the initial conditions (3.11) and (3.12). Since ${\it\lambda}<-1,$ the supremum of ${\it\gamma}$ blows up and this happens along the characteristic starting at the position of the supremum of ${\it\gamma}_{0}$ , $\boldsymbol{X}_{+}=((3{\rm\pi}/2),(5{\rm\pi}/4))$ (see figure 2 for reference). On the other hand, the infimum of ${\it\gamma}$ (a negative quantity) remains small in size. Thus, in this case the norm $\Vert {\it\gamma}(\cdot ,t)\Vert _{\infty }$ can be identified with $\sup {\it\gamma}(\cdot ,t)$ (i.e. one can set ${\it\sigma}_{\infty }\equiv 1$ ) for all times $0\leqslant t<T^{\ast }.$ Equation (3.2) evaluated at $\boldsymbol{X}_{0}=\boldsymbol{X}_{+}$ is then compared with definition (4.1) of mapped time ${\it\tau}$ to give

The above relation can be inverted to solve for $t$ as a function of ${\it\tau}$ , provided $t<T^{\ast }$ . The supremum of original stretching rate, in terms of ${\it\tau},$ can be obtained after using this inversion along with the formula in table 1, giving

There is a simple analytical expression (in terms of mapped time ${\it\tau}$ ) for vorticity at the position where the maximum of $|{\it\gamma}(\boldsymbol{x},t)|$ is attained, $\boldsymbol{X}_{{\it\gamma}}(t)$ . Equation (3.3) evaluated at $\boldsymbol{X}_{0}=\boldsymbol{X}_{+}$ gives

Note that this is a lower bound for the $L^{\infty }$ norm of the vorticity. The latter norm can be obtained at all times in terms of the initial conditions, by maximising the right-hand side of (3.3) over the quadrant ${\rm\pi}\leqslant x,y\leqslant 2{\rm\pi}.$ Although this rarely leads to an explicit analytical expression for $\Vert {\it\omega}(\cdot ,t({\it\tau}))\Vert _{\infty },$ it can always be computed numerically to any desired accuracy.

Figure 3. The errors $Q_{{\it\gamma}}$ (a) and $Q_{{\it\omega}}$ (b) for both the original 2D Euler model and the mapped 2D Euler model at various resolutions. Labels are (a) $N=256$ mapped system, (b) $N=256$ original system, (c) $N=512$ mapped system, (d) $N=512$ original, (e) $N=1024$ mapped, (f) $N=1024$ original, (g) $N=2048$ mapped, (h) $N=2048$ original. There is a clear accuracy gain by performing the mapping; the errors remain small for later times. Grey dotted vertical line is at ${\it\tau}=4.34$ , representing the reliability time at which the $N=2048$ case becomes unresolved, see § 5.3 and table 2.

5.1. Numerical solutions of the original and mapped systems

We turn our attention to the original evolution (2.16) and (2.17) and the mapped evolution (4.8) and (4.9). Both systems were solved numerically using a standard pseudospectral method written in CUDA, using CUFFT library and implemented on NVIDIA GPUs. To remove the usual aliasing errors, Hou’s high-order exponential filter (Hou & Li Reference Hou and Li2007) is used, where we multiply the spectrum at each time step by the factor $\exp (36(2k/N)^{36}),$ where $k$ is the modulus of the wavevector and $N$ is the spatial resolution. We checked that the $2/3$ dealiasing rule (in which the last $1/3$ of the high-frequency modes are set to zero) gives similar results, but Hou’s filter provides sensible spectra for a slightly broader range of wavevectors.

In both systems time marching was carried out using the fourth-order Runge–Kutta method. In the mapped system, a uniform time step $\text{d}{\it\tau}$ is used since ${\it\tau}$ stretches the temporal domain such that the singularity is at ${\it\tau}=\infty .$ This uniform $\text{d}{\it\tau}$ would imply that in the original system the time step $\text{d}t$ should be adaptive, via $\text{d}t=\text{d}{\it\tau}/\Vert {\it\gamma}(\cdot ,t)\Vert _{\infty }.$ This adaptive method is normally used in 3D Euler blowup assessment studies for reasons of accuracy close to singularity. In the results below adaptive time stepping in original variables is used for this reason, with the added advantage that the data from the original and the mapped systems are more comparably spaced. A striking result of this paper, to be shown and discussed thoroughly below in this section, is that even within this fair-comparison scenario, the mapped system gains a significant amount of accuracy in the estimation of blowup quantities (cf. figures 3, 10 and 11).

Finally, unlike the integration of the original system, the numerical integration of the mapped system requires a special method: a normalisation so that the mapped system satisfies the constraint $\Vert {\it\gamma}_{map}(\cdot ,{\it\tau})\Vert _{\infty }=1$ , for all ${\it\tau}$ . The accuracy of this normalisation is essential for the performance of the mapped system’s numerical integration. To apply the normalisation, $\Vert {\it\gamma}_{map}(\cdot ,{\it\tau})\Vert _{\infty }$ is computed using a $P_{6,k}$ quarter-section interpolation. This $P_{6,k}$ interpolation is an efficient procedure to compute the field’s maximum and its position using an iterative application of cubic splines at progressively finer resolution. This was tested against several other interpolation methods. See the discussion in appendix A.

The initial conditions were chosen as in (3.11) and (3.12). Figure 2 shows contour plots for both the vorticity and stretching rate at various times.

5.2. Errors in local quantities near blowup: $\Vert {\it\gamma}(\cdot ,t)\Vert _{\infty }$ and ${\it\omega}(\boldsymbol{X}_{{\it\gamma}}(t),t)$

A sensible definition of the error of the numerical simulation of the original 2D Euler model is the relative difference between the supremum norm of ${\it\gamma}$ obtained from the numerical simulation and the exact analytical formula, given in table 1. We can also define the error in ${\it\omega}$ by evaluating it at $\boldsymbol{X}_{{\it\gamma}}(t),$ the location of the supremum of $|{\it\gamma}|$ , also given analytically in table 1. We define the relative errors as

where the subscripts ‘ $num$ ’ and ‘ $ana$ ’ stand for ‘numerical’ and ‘analytic’. While these are useful for instantaneous monitoring purposes, given a time series of a numerical solution of ${\it\gamma}$ or ${\it\omega}$ we would like to measure the error using a single number for each variable. The naive measure in terms of the $L^{2}$ norm of the relative error, $\Vert \mathscr{E}\Vert _{2}\equiv \sqrt{\int _{0}^{\text{ T }}[\mathscr{E}(t)]^{2}\,\text{d}t},$ is not the best choice because it is not bounded a priori. Given two signals $f(t)$ and $g(t)$ for comparison, we define an $L^{2}$ norm of the error (not relative error) normalised with the sum of norms of the individual signals (Perlin & Bustamante Reference Perlin and Bustamante2014):

which has three advantageous properties:

1. (i) $Q$ and $\Vert \mathscr{E}\Vert _{2}$ are proportional if they are small: $\Vert \mathscr{E}\Vert _{2}\propto \sqrt{T}Q$ if $\mathscr{E}\ll 1$ ;

2. (ii) $Q$ , being dimensionless, does not explicitly require a time scale ( $T$ ) so it can be applied in a variety of contexts; in practice, though, and for assessment purposes, we will normally plot $Q$ as a function of the total integration time $T$ ;

3. (iii) $Q$ is bounded, $0\leqslant Q\leqslant 1,$ with value $0$ representing perfect match and value $1$ representing perfect mismatch.

Thus, we work with $Q_{{\it\gamma}}=Q(\Vert {\it\gamma}_{num}(\cdot ,t)\Vert _{\infty },\Vert {\it\gamma}_{ana}(\cdot ,t)\Vert _{\infty })$ and $Q_{{\it\omega}}=Q({\it\omega}_{num}(\boldsymbol{X}_{{\it\gamma}},t),{\it\omega}_{ana}(\boldsymbol{X}_{{\it\gamma}},t))$ .

Figure 4. (a) Ratio between the terms $(2+{\it\lambda})\!\int _{0}^{{\it\tau}}\!{\it\sigma}_{\infty }\langle {\it\gamma}_{map}^{2}\rangle \,\text{d}{\it\tau}^{\prime }$ and $-(1+{\it\lambda})\int _{0}^{{\it\tau}}{\it\sigma}_{\infty }\,\text{d}{\it\tau}^{\prime },$ case ${\it\lambda}=-3/2,$ appearing in the exponent in (4.11). (b) Resolution study. Plots of $\langle {\it\gamma}_{map}(\cdot ,{\it\tau})^{2}\rangle$ as a function of ${\it\tau}$ using data from the numerical integration of the mapped system, at different resolutions: from top to bottom, $N=256$ , 512, 1024, 2048. The curves converge to the analytically-obtained asymptotic regime $\langle {\it\gamma}_{map}(\cdot ,{\it\tau})^{2}\rangle \sim (3/11)\exp (-{\it\tau}).$

To gain an appropriately accurate estimate of these errors we must consider the best method for approximating $\Vert {\it\gamma}_{num}(\cdot ,t)\Vert _{\infty },\boldsymbol{X}_{{\it\gamma}}$ and ${\it\omega}_{num}(\boldsymbol{X}_{{\it\gamma}},t)$ . The simplest approach is to apply the maximum value across the collocation points of the discretised field, however this leads to significant spurious oscillations. The more accurate procedure is to perform some post-processing interpolation. We use the same highly-accurate interpolation as in the normalisation procedure in the mapped system. The various interpolation options are described in appendix A.

The numerical solution of the mapped system does not provide direct access to the original variable $\Vert {\it\gamma}(\cdot ,t)\Vert _{\infty }$ so in order to compute it we employ (4.11), namely $\Vert {\it\gamma}(\cdot ,t({\it\tau}))\Vert _{\infty }=\Vert {\it\gamma}_{0}\Vert _{\infty }\exp [-(1+{\it\lambda})\int _{0}^{{\it\tau}}{\it\sigma}_{\infty }\,\text{d}{\it\tau}^{\prime }+(2+{\it\lambda})\int _{0}^{{\it\tau}}{\it\sigma}_{\infty }\langle {\it\gamma}_{map}^{2}\rangle \,\text{d}{\it\tau}^{\prime }],$ where $\int _{0}^{{\it\tau}}{\it\sigma}_{\infty }\langle {\it\gamma}_{map}^{2}\rangle \,\text{d}{\it\tau}$ is computed using Simpson’s rule. We compare this against (4.16). Note that we are using the global quantity $\langle {\it\gamma}_{map}^{2}\rangle$ for the assessment of the local quantity $\Vert {\it\gamma}(\cdot ,t)\Vert _{\infty }.$ Therefore, errors of the global quantity $\langle {\it\gamma}_{map}^{2}\rangle$ might affect the errors of this assessment. To see at which times these errors might be important, figure 4(a) shows the ratio between the terms $(2+{\it\lambda})\int _{0}^{{\it\tau}}{\it\sigma}_{\infty }\langle {\it\gamma}_{map}^{2}\rangle \,\text{d}{\it\tau}^{\prime }$ and $-(1+{\it\lambda})\int _{0}^{{\it\tau}}{\it\sigma}_{\infty }\,\text{d}{\it\tau}^{\prime },$ appearing in the exponent in (4.11). It is clear that at early times the global quantity has more influence on the size of the error $Q_{{\it\gamma}}$ . This is addressed in § 5.2.1. In contrast, at late times the global quantity is not relevant and this explains why the size of the error $Q_{{\it\gamma}}$ remains small and stable.

To evaluate the original variable ${\it\omega}(\boldsymbol{X}_{{\it\gamma}}(t({\it\tau})),t({\it\tau}))$ from the mapped variables we use

We compare this against (4.17).

Figure 3 shows a comparison of the errors $Q_{{\it\gamma}}$ and $Q_{{\it\omega}}$ at various resolutions and demonstrates how the analysis of the mapped system along with its numerical solution serve to improve the accuracy of the blowup quantities $\Vert {\it\gamma}(\cdot ,t)\Vert _{\infty }$ and ${\it\omega}(\boldsymbol{X}_{{\it\gamma}}(t),t)$ near the singularity time.

We produce standard resolution convergence studies of the quantities relevant to this section. Figure 5(a) shows the classical convergence study of the multiplicative inverse of the supremum norm of stretching rate, using the numerical solution of the original system, plotted as a function of original time $t$ . Good convergence to the analytical result is obtained. This is the basis for the method of computing running estimates of singularity time $T^{\ast }$ (method A in § 5.4). Figure 5(b), shows the less classical (but similar in spirit) lin–log convergence study of the supremum norm of stretching rate, again using the numerical solution of the original system, plotted as a function of mapped time ${\it\tau}.$ Again, good convergence to the analytical asymptote is obtained. Finally, figure 4(a) shows the lin–log convergence study of the spatial average of the square of the mapped stretching rate, $\langle {\it\gamma}_{map}^{2}\rangle ,$ using the numerical solution of the mapped system. Good convergence to the analytical asymptote is verified.

Figure 5. Two different resolution studies of the data from original system at different resolutions $N=256$ , 512, 1024, 2048. (a) Classical plot of $1/\Vert {\it\gamma}(\cdot ,t)\Vert _{\infty },$ giving a convergence to the analytically obtained asymptotic regime $1/\Vert {\it\gamma}(\cdot ,t)\Vert _{\infty }\sim (T^{\ast }-t)/2$ (solid line). (b) New plot, in lin–log scaling, of $\Vert {\it\gamma}(\cdot ,t({\it\tau}))\Vert _{\infty }$ as a function of ${\it\tau},$ giving a convergence to the analytically-obtained asymptotic regime $\Vert {\it\gamma}(\cdot ,t({\it\tau}))\Vert _{\infty }\sim \sqrt{11/2}\exp ({\it\tau}/2)$ (solid line).

5.3. Spectra, analyticity strip and BKM theorem

An effective diagnostic for the spatial collapse associated with a typical finite-time blowup scenario presents itself in the form of the analyticity strip (Sulem et al. Reference Sulem, Sulem and Frisch1983; Bustamante & Brachet Reference Bustamante and Brachet2012). Given the spectrum of spatial Fourier coefficients provided in our numerical simulation, the $L^{2}$ spectrum of ${\it\gamma}$ is defined as the sum of the squares of modulus of the Fourier coefficients over circular shells, or in short, the $L^{2}$ stretching-rate spectrum:

While we know analytical solutions for the blowup quantities, so far we have not found a method to obtain analytical expressions regarding the stretching-rate spectrum, other than the formula valid for ${\it\lambda}=-3/2,$ stating $\sum _{k=1}^{\infty }E(k,t)=3/4\text{(const)}.$

The lack of analytical results for spectra is seen as an advantage since it allows us to test the analyticity-strip method and its bridge with the BKM theorem (Bustamante & Brachet Reference Bustamante and Brachet2012) from a purely numerical point of view. The results presented in this section can therefore be contrasted against future analytical developments.

The function ${\it\gamma}(\boldsymbol{x},t)$ remains analytic in the space variables if $E(k,t)$ can be bounded by

where ${\it\delta}_{E}(t)$ is the analyticity strip width, also known as the logarithmic decrement, and $C_{E}(t),n_{E}(t)$ are positive numbers. We assume this approximation holds for our functions. Figure 8 shows snapshots of the $L^{2}$ spectrum $E(k,t)$ , in log–log as well as lin–log scaling, to provide evidence of the feasibility of this approximation. The common procedure is to find the coefficients $C_{E}(t),n_{E}(t),{\it\delta}_{E}(t)$ by performing a least-squares fit, at each time $t,$ on $\ln E(k,t)$ over some interval $k_{i}\leqslant k\leqslant k_{f}$ . The problem becomes linear in the parameters $\ln C_{E}(t),n_{E}(t)$ and ${\it\delta}_{E}(t)$ . More details can be found in, e.g., Bustamante & Brachet (Reference Bustamante and Brachet2012, equation (5)).

It is customary to define a ‘reliability time barrier’ $t_{rel}$ by the condition

where $\text{d}x$ is the grid spacing of the numerical simulation. This barrier represents the obvious requirement that the smallest scales available in the numerical simulation are well resolved. Figure 7 shows the results at resolution $N=2048$ for the fit parameters at several times (in mapped time ${\it\tau}$ ). Figure 7(d) indicates a reliability time ${\it\tau}_{rel}\approx 4.2$ , corresponding to $t_{rel}\approx 1.121.$ In figure 7(a), the $+$ symbol shows $n_{E}(t)$ with the remarkable convergence to $n_{E}=5/3\pm 0.08$ (dotted horizontal line) near reliability time. The error bar $0.08$ stands for the $5\,\%$ error of the least-squares fit procedure.

Figure 7. Results at resolution $N=2048$ (except (a)) for the fit parameters of the $L^{2}$ spectra $E(k,t)$ and $L^{1}$ spectra $F(k,t)$ at several times (as function of mapped time ${\it\tau}$ ). (a) Resolution study ( $N=256,512,1024,2048$ ) for $n_{E}(t({\it\tau})).$ Progressive convergence is observed towards $n_{E}\approx 5/3$ (within a $5\,\%$ error) near reliability time. (b) Results for $n_{F}(t({\it\tau}))$ (circles). Solid and dashed curves represent the upper and lower bounds $n_{E}/2$ and $(n_{E}-1)/2,$ respectively, cf. inequality (5.15). The curve $n_{F}$ is consistently between these bounds. (c) Results for $C_{E}(t({\it\tau}))$ ( $+$ symbols) and $C_{F}(t({\it\tau}))$ (circles). At late times these coincide, thus confirming the isotropy of the 2D spectrum. (d) Results for ${\it\delta}_{E}(t({\it\tau}))$ ( $+$ symbols) and ${\it\delta}_{F}(t({\it\tau}))$ (circles). These coincide, in agreement with inequalities (5.15). The horizontal line is the smallest resolved scale ${\rm\Delta}x=2{\rm\pi}/2048$ and the dotted line corresponds to the numerical fit ${\it\delta}=\exp (0.53{-}1.33{\it\tau})$ obtained using data in the range $3\leqslant {\it\tau}\leqslant 4.34.$

A classical method used in Bustamante & Brachet (Reference Bustamante and Brachet2012) gives rise to the following rigorous inequality:

This inequality is saturated if and only if there is alignment of the phases of the Fourier components $\hat{{\it\gamma}}(\boldsymbol{k},t)$ that carry a significant amplitude. This alignment is expected to happen near the singularity time, as learned from the 1D inviscid Burgers equation. The fact that a $L^{1}$ norm (rather than the more familiar $L^{2}$ norm $E(k,t)$ ) appears in this rigorous inequality, motivates the introduction of a new type of spectrum, the $L^{1}$ stretching-rate spectrum

In terms of this $L^{1}$ spectrum the above inequality becomes

We now make a connection between this new $L^{1}$ spectrum and the more familiar $L^{2}$ spectrum, via a rigorous equivalence:

where

We have the approximate result valid as $k\rightarrow \infty$ :

Inequalities (5.15) allow us to work with the $L^{1}$ spectrum $F(k,t)$ in the same way as we would work with the more familiar $L^{2}$ spectrum $E(k,t).$ In particular, a fit of the form

implies a fit of the form

One could interpret these two fits as working hypotheses, as in Bustamante & Brachet (Reference Bustamante and Brachet2012). In the limit $k\rightarrow \infty$ the inequalities (5.15) imply a sandwich so ${\it\delta}_{E}(t)={\it\delta}_{F}(t)$ is necessary. This is verified numerically in figure 7(d).

In the intermediate- $k$ range, the inequalities imply the following bounds:

These bounds are verified in figure 7(b).

Finally, the proportionality factors $C_{F}(t)$ and $C_{E}(t)$ cannot be easily related unless the above bounds (5.20) for $n_{F}(t)$ become saturated. For example, if the upper bound was saturated, this would correspond rigorously to a very anisotropic $2D$ spectrum $|\hat{{\it\gamma}}(\boldsymbol{k},t)|$ in terms of orientation of $\boldsymbol{k}.$ In this case we would have $C_{F}\approx \sqrt{C_{E}}.$ On the other hand, if the lower bound was saturated, as it happens at late times (see figure 7(b)) this would correspond rigorously to a very isotropic $2D$ spectrum $|\hat{{\it\gamma}}(\boldsymbol{k},t)|.$ In this case we would have $C_{F}\approx \sqrt{2{\rm\pi}C_{E}}.$ Figure 7(c) shows the curves $C_{F}$ and $\sqrt{2{\rm\pi}C_{E}}$ as functions of time, showing equality at late times. So we can conclude that the $2D$ spectrum becomes nearly isotropic at late times. To support this analysis we provide in figure 9 two snapshots of the stretching-rate 2D spectrum, at mapped times ${\it\tau}=2$ (a) and ${\it\tau}=4$ (b). It is evident that the spectrum is strongly anisotropic at early times and evolves towards isotropy at late times, in agreement with the saturation of inequalities found.

Figure 8. (a) and (b) Snapshots of stretching-rate 1D shell spectra $E(k,t)$ at mapped times ${\it\tau}=1,2,3,4,5$ (curves progressing from bottom to top) in log–log scale (a) and lin–log scale (b). The slopes of the straight lines on the left panel are proportional to the exponent $n_{E}(t)$ . The slopes of the straight lines on the right panel are proportional to the logarithmic decrement ${\it\delta}_{E}(t).$ Resolution $N=2048.$

Figure 9. Snapshots of stretching-rate 2D spectra at mapped times ${\it\tau}=2$ (a) and ${\it\tau}=4$ (b). Resolution $N=2048.$

Following an analogous discussion to that in Bustamante & Brachet (Reference Bustamante and Brachet2012), we see that the left-hand-side of inequality (5.14) has a singular behaviour. In fact, a BKM type of theorem can be demonstrated for the left-hand side, namely we can assume

This is obvious from the analytical solution presented in table 1. This will imply a singular behaviour of the right-hand side of inequality (5.14):

Using the above fit for $F(k,t)$ we get the result

We recognise the Jonquiere’s function $\text{Li}(n_{F}(t),\text{e}^{-{\it\delta}_{F}(t)}).$ We get

Now, in the limit as $t\rightarrow T^{\ast }$ we have ${\it\delta}_{F}\rightarrow 0$ so we can approximate the Jonquiere’s function to get

where $n_{0}=\liminf _{t\rightarrow T^{\ast }}n_{F}(t).$ So the asymptotic behaviour ${\it\delta}_{F}(t)\sim (T^{\ast }-t)^{{\rm\Gamma}}$ would be consistent with singularity behaviour if and only if

From figure 7(b), we see that $n_{0}\approx 0.36$ if we consider the data near reliability time ( ${\it\tau}_{rel}\approx 4.34$ ). This gives ${\rm\Gamma}\geqslant 1.56,$ which is consistent with the fits obtained from figure 7(d), that produce ${\rm\Gamma}\approx 2.66$ by virtue of the analytical result $(T^{\ast }-t)\approx \exp (-{\it\tau}/2).$ At the same time this shows that the inequality (5.14) is not saturated by the field ${\it\gamma}(\boldsymbol{x},t).$ The interpretation of this lack of saturation is that the Fourier phases, $\{\arg \hat{{\it\gamma}}(\boldsymbol{k},t)\}_{k=1}^{\infty },$ do not all align near the singularity time. In fact this is evident from the physical-space snapshots in figure 2, where the singular structure is between a filament and a point.

5.4. Estimating the singularity time $T^{\ast }$ efficiently

It would be meaningless to provide results comparing the accuracy of particular numerical methods, without including some quantification of their relative computational expense.

It is possible to obtain two independent estimates of the singularity time $T^{\ast }$ by using the numerical solutions from either the original system or the mapped system.

Method A. For the original system, we fit the stretching rate norm $\Vert {\it\gamma}(\cdot ,t)\Vert _{\infty }$ as a power law locally in time, using a method introduced in Bustamante & Kerr (Reference Bustamante and Kerr2008) to determine an estimate of singularity time $T^{\ast }$ . The power-law fits are of the form

where $f(t)$ stands for $\Vert {\it\gamma}(\cdot ,t)\Vert _{\infty }$ in this case. This ansatz is justified in this case by the analytically obtainable asymptotic formulae in table 1. The local fits are achieved using the function

Instantaneous running estimates ${\it\alpha}$ and $T_{A}^{\ast }$ are then computed by linear fitting the function $g(t)$ instantaneously using adjacent data points, or more generally over a small time window (of size ${\sim}0.2$ ) containing a good number of data points, in order to eliminate spurious oscillations in the running estimates. Note that this method provides an extra quantity: the exponent ${\it\alpha}$ , which serves as an extra measure of validation. In the case ${\it\lambda}=-3/2$ one should get ${\it\alpha}=-1$ (see table 1). This validation is consistently held throughout the computation (figure not shown).

Method B. For the mapped system, the situation relies on the explicit formula (4.12). There, the only relevant numerical quantity is $\int _{0}^{{\it\tau}}\langle {\it\gamma}_{map}(\cdot ,{\it\tau}^{\prime })^{2}\rangle \,\text{d}{\it\tau}^{\prime }.$ In analogy to the previous case, we will fit the integrand. But, unlike the previous case, we cannot rely on local fits because doing this would lead to accumulation of errors in the estimation of the integral.

Using (4.11) and (4.12) we obtain the estimate

where in both terms we compute the integrand using (4.11). The subscript ‘ $num$ ’ means that we use the numerical solution of the mapped system to compute the time integrals, using Simpson’s rule. As for the subscript ‘ $extrap$ ’, $\widehat{{\it\tau}}$ is a very big number chosen so that the numerical integral converges ( ${\sim}1000$ in practice) and we compute the integrand, using (4.11), as follows:

where we set ${\it\sigma}_{\infty }({\it\tau}^{\prime \prime })={\it\sigma}_{\infty }({\it\tau})$ in the above exponent, and we model the functions appearing in the above exponent using the following fit ansatz that is motivated by the generic asymptotic behaviour of $\langle {\it\gamma}_{map}^{2}\rangle$ as ${\it\tau}^{\prime \prime }\rightarrow \infty$ :

where ${\it\beta}$ and $c$ are two positive fit parameters. In order to obtain these fit parameters we use the numerical data for $\langle {\it\gamma}_{map}^{2}(\cdot ,{\it\tau}^{\prime \prime })\rangle$ in the range $0\leqslant {\it\tau}^{\prime \prime }\leqslant {\it\tau}$ and use a least-squares fit. The integral $\int _{{\it\tau}}^{{\it\tau}^{\prime }}\!{\it\sigma}_{\infty }\langle {\it\gamma}_{map}^{2}\rangle \,\text{d}{\it\tau}^{\prime \prime }$ is done analytically in terms of the fit parameters. Finally, the resulting integral $\int _{{\it\tau}}^{\widehat{{\it\tau}}}\!(1/\Vert {\it\gamma}(\cdot ,t({\it\tau}^{\prime }))\Vert _{\infty })\,\text{d}{\it\tau}^{\prime }$ is done using Simpson’s rule.

It is important to stress that in both original and mapped systems the estimation of the singularity time $T^{\ast }$ depends on two fit parameters. In general, we tried to lever as much accuracy as possible from both methods. For example, we used adaptive time stepping in the original system to get a distribution of data points that is comparable with that of the mapped system, so that a more accurate estimate for $T^{\ast }$ could be obtained in the original system.

We present results from original and mapped systems regarding the assessment of estimates of singularity times. First, in figure 10 we plot the relative error of the running estimates

where $T_{ana}^{\ast }=(4/\sqrt{3})\arctan \,(\sqrt{6}/4)\approx 1.26894$ is the analytically computed singularity time and $T_{A,B}^{\ast }$ stands for the running estimate obtained either from method A (for the original system) or from method B (for the mapped system). It is observed from the figure that: (i) for each method, there is good resolution convergence in the assessment of $T^{\ast }$ ; (ii) method B (for the mapped system) produces much better results as compared with method A (for the original system), with an improvement of about three orders of magnitude at any given resolution.

Figure 10. Errors in the running estimates of singularity time $T^{\ast }$ using data from original system’s numerical integration (lines) and mapped system’s numerical integration (symbols) at different resolutions: from top to bottom in each case, $N=256$ , 512, 1024, 2048.

The second set of results is the following. We produce, from each method, a single estimate (not a running estimate) $T_{A,B}^{0}$ of the singularity time, computed using the running estimates already obtained. It is important to stress the perhaps obvious fact that the procedure to find this single estimate is completely independent of any previous knowledge of the singularity time $T^{\ast }$ . The procedure is simple: since we have a reliability time $t_{rel}$ (or, in mapped time, ${\it\tau}_{rel}$ ) well defined for each resolution, in terms of the analysis of spectra done in § 5.3, we evaluate our running estimates at the reliability time. So we define, for the mapped system, $T_{B}^{0}=T_{B}^{\ast }|_{{\it\tau}_{rel}}.$ As for the original system, setting the single estimate to $T_{A}^{\ast }|_{t_{rel}}$ would be possible, however we found that a better estimate is obtained by averaging the running estimates between $t_{rel}$ and $t_{min},$ where $t_{min}\,({>}t_{rel})$ depends on the resolution and is defined by the time at which the running estimate has a global minimum.

Figure 11 shows the CPU time versus relative error of the estimated singularity time,

for the numerical solutions of both the original (method A) and mapped systems (method B) at various resolutions. It is clear that while the mapping incurs some additional expense in evaluating the extra terms, computing the interpolated supremum and applying the normalisation, it is far outweighed by the positive effect on the errors (three orders of magnitude in this study). In these measures one can make a significant improvement, saving not only CPU time, but also the memory cost of high-resolution runs.

We present in table 2 a useful summary of the reliability times and relative errors $\mathscr{E}_{A,B}$ of the estimated singularity time, for a range of resolutions, showing that the errors are dramatically reduced in the case of the mapped equations.

Figure 11. The singularity-time error $\mathscr{E}_{A,B}$ (error between estimated singularity time and analytically computed singularity time $T^{\ast }$ , (5.33)), as a function of CPU time for both the original system, method A (filled symbols) and the mapped system, method B (open symbols) at various resolutions: $N=256$ (squares, red online), $N=512$ (circles, green online), $N=1024$ (triangles, blue online), $N=2048$ (diamonds, magenta online). The CPU overhead of applying the mapping is shown to be more than covered by the accuracy gain.