## 1 Introduction

Pseudospectral methods are a powerful numerical approach, offering exponential convergence to the true solution across the entire domain for smooth functions (see Boyd2013 for a review). Thanks to their exquisite accuracy and efficiency they have been successfully applied in, for example, large scale relativistic astrophysics simulations (e.g. SPeC; Boyle2008; Scheel2009). A wealth of astrophysical phenomena, however, involve discontinuities such as fluid shocks in supernova explosions or galactic spiral arms. For this class of problem, pseudospectral schemes perform poorly due to sublinear convergence properties and spurious oscillations induced in the solution known as the Gibbs phenomenon Wilbraham1848; Gibbs1898; Michelson1898.

A substantial body of research has been conducted in the attempt to alleviate the Gibbs phenomenon and recover the exponential convergence of pseudospectral schemes. Some proposed solutions are reprojection onto combinations of different basis sets Gottlieb1992; Gottlieb2011, filtering Vandeven1991, artificial viscosity Tadmor1990, modelling and removal of oscillations Krylov2006; Lipman2010 and mollification Gottlieb1985. In Piotrowska2019 we explored and numerically tested the last approach, using an optimal mollifier as prescribed by Tanner2006. We showed successful recovery of convergence properties away from discontinuities and offered improvements to preserve the discontinuous character of the solution. We also solved a toy problem advecting a discontinuous function, demonstrating greater robustness of mollification as compared to the more popular Gegenbauer reconstruction methods.

In this work we take the next step towards applying pseudospectral methods to solve physically relevant shock problems. This time we allow the discontinuities to develop dynamically, instead of introducing them by construction and apply previously explored methods to detect the created shocks. Robust identification of evolving discontinuities is crucial in pseudospectral shock capturing, regardless of the choice of Gibbs phenomenon removal method. We discuss challenges associated with the detection of discontinuity onset and capturing of steep gradients at constant spatial resolution. We also present our heuristic solution to those issues, offering practical improvements to shock detection methods in case of shocks dynamically developing during the simulation. Finally, we solve a toy problem of the inviscid Burgers’ equation to demonstrate the performance of our methods when applied to shocks developing smoothly over time.

## 2 Background

In this section we review existing methods which provide basis for our improvements presented further in Section 4.

### 2.1 Overview of pseudospectral methods used in this work

The pseudospectral decomposition is a means of discretizing a function of interest with a finite sum of orthogonal basis functions :

(1) |

where is the partial sum of and the basis functions are denoted . The expansion coefficients are chosen such that

(2) |

and are given by:

(3) |

where

(4) |

is the inner product between and on a domain . The integral can be approximated via Gauss quadrature in which case the are the quadrature collocation points and their associated weights, such that:

(5) |

Thus, one can conveniently transform the
collocation-point values into expansion
coefficients through^{1}^{1}1where in the Einstein
notation:

(6) |

where is the inverse of the Vandermonde matrix. One then can also approximate the partial sum of the derivative of : with the following differentiation matrices:

(7) |

such that differentiation conveniently reduces to the following matrix operations:

(8) |

In our work we use a basis of Chebyshev polynomials of the
first kind^{2}^{2}2other choices are possible, depending on the
problem of interest e.g. a Legendre polynomial or Fourier basis
defined on the domain and
the Chebyshev-Gauss-Lobatto quadrature.
The associated grid of collocation points allows us to probe the
solution values at the boundaries with and .
For a detailed
description of Gaussian quadrature and transformations between
the spectral and collocation-point bases we refer the interested
reader to Grandclement2006 and Boyd2013.

### 2.2 Stabilization techniques for non-linear PDEs

In this research we focus on one-dimensional (1D) problems which smoothly develop discontinuities over time, such as the non-linear 1D inviscid Burgers’ equation:

(9) |

In our investigation we assume smooth initial conditions, choosing a Gaussian centered on :

(10) |

where and .

Once we discretize the initial conditions with a pseudospectral projection onto the Chebyshev polynomial basis, our hyperbolic partial differential equation (PDE) reduces to a set of

ordinary differential equations (ODEs) at the collocation points. Thus, we solve Eq. (9) using the method of lines, keeping time as a continuous variable and assuming periodic boundary conditions such that at each time step is set to . Note that our initial condition is non-smooth at the boundaries, since is non-periodic. Our periodic boundary conditions serve as a simplistic proxy for a multi-domain spectral method, with the left and right boundaries serving as inflow and outflow into other spectral domains respectively. This allows us to run our simulation for arbitrary long time, despite the finite-sized domain.In general, any linearly stable numerical scheme does not guarantee stability for non-linear problems. In a pseudospectral method nonlinear terms can produce artificially large short-wavelength coefficients which need to be subdued to avoid a numerical catastrophe. In this work we choose a modification to the spectral viscosity technique Tadmor1990, namely the right hand side filtering introduced by Miller2016 in Eq. (73-75). In this framework filtering is applied to the derivative operator modes via the following diagonal matrix:

(11) |

where is the so-called order of dissipation and we chose in the evolution of our system which proved optimal in our numerical tests. The temporal derivative at each grid point then takes the form:

(12) |

where the constant depends on the particular problem and performed best when tested with the inviscid Burgers’ equation.

### 2.3 Removing the Gibbs phenomenon

The main advantage of pseudospectral methods is their global exponential convergence

for smooth functions. More specifically, the exponential decrease of the interpolation error

integrated across the computational domain as a function of expansion order :(13) |

where , is the highest order of derivative well-defined for , is a constant factor and denotes the norm. In the smooth case, arbitrarily high-order derivatives are well-defined and the inequality holds for all Grandclement2009.

When applied to discontinuous problems, pseudospectral interpolants lose their convergence properties with the norm exhibiting at best convergence Thomson2008. Furthermore, the pseudospectral representations of the solution suffer from spurious oscillations. The oscillation amplitude decreases away from the discontinuities and their frequency varies depending on the position within the domain. The left panel of Fig. 1 illustrates this phenomenon by overlaying a Gibbs-affected reconstruction on the input discontinuous piecewise linear ‘tophat’ function:

(14) |

where and .

In order to mitigate these issues and recover a Gibbs-free solution one needs to post-process the affected spectral representation. This procedure consist of two main parts: edge detection and Gibbs oscillations removal. Edge detection is a step necessary for many Gibbs phenomenon removal techniques, allowing construction of adaptive smoothing kernels in the mollification approach or the correct domain decomposition in the case of the Gegenbauer reconstruction Gottlieb2011a.

Following our work in Piotrowska2019, we choose to locate the edges using the technique introduced by Gelb2008. It relies on applying a range of different concentration factors Gelb1999 to the partial sum of the function derivative , which accelerate the rate at which converges onto the jump function,

(15) |

In this prescription, each choice of an admissible generates a jump function approximation given by the following partial sum:

(16) |

The jump function approximation then converges onto the true jump function at the following rate Gelb2008:

(17) |

and each exhibits a different oscillation pattern around the discontinuity associated with a given Gelb2000. The prescription then takes advantage of the differences in the oscillatory patters, combining a range of jump functions such that the oscillations are removed Gelb2008:

(18) |

The middle panel of Fig. 1 demonstrates three jump function approximations from trigonometric, polynomial and exponential forms as well as a resulting for . The location of discontinuities is indicated by the positions of the peaks in the function represented by the black solid curve. In our work we construct the using defined in Eq. (8-10) in Gelb2008 and their filtered versions constructed using Lanczos Lanczos1966 filters of orders 0, 1, 2 and 3.

Once we identify the discontinuities and their positions we then remove the oscillations through mollification with adaptive one-sided mollifiers Tanner2006; Piotrowska2019. The right panel of Fig. 1 demonstrates the result of post processing by overplotting the mollified result (black solid) on the Gibbs-affected reconstruction (red dashed). The post-processed solution is successfully freed from spurious oscillations, simultaneously preserving its discontinuous character.

Finally, it is important to mention, that even in the ideal cases of well-defined discontinuities the function is polluted by residual oscillations of small amplitude which can be misidentified as jumps in an automated search for extrema. Hence one always has to make a choice of a peak height qualifying as a true extremum. As a result, one then sets the range of discontinuity heights captured by the edge detection framework, which depends on the resolution and the problem at hand.

## 3 New shock detection challenges

The edge detection techniques illustrated in Fig. 1 are an implementation of prescriptions developed by A. Gelb and collaborators Gelb2000; Gelb2008 for stationary discontinuities. Realistic numerical applications, however, involve dynamical development of discontinuities in time and thus require the pseudospectral framework to capture transitions between the smooth and discontinuous (shocked) regimes in an automated fashion. Thus, one requires a robust means of automatically recognizing discontinuous functions to deploy post-processing techniques at appropriate time steps.

Solutions in which shocks develop smoothly over time experience additional problems caused by insufficient resolution of the truncated pseudospectral expansion basis. An example of such is the inviscid Burgers’ equation (eq. 9) where the spatial gradient in the function steepens continuously until the point of wave breaking. In such cases the solution is populated by additional oscillations around steep gradients which result from the lack of sufficient resolution to capture significant variation at very small scales. The left panel of Fig. 2 demonstrates this issue by showing the affected solution of the 1D inviscid Burgers’ equation 9 with smooth Gaussian initial conditions (Eq. 10 with and

), subject to periodic boundary conditions. The panel shows a skewed Gaussian at time step

which has a steep gradient at . The solution in the vicinity of the steep slope in is populated by spurious oscillations of small amplitude which introduce undesired, incorrect, information. These oscillations are an artefact of insufficient resolution and their amplitude converges away with increasing expansion order .What is more, the resolution limit also demonstrates itself in shock detection, altering the function such that it exhibits different convergence properties and features additional peaks. The right panel of Fig. 2 highlights the latter by showing the irregular shape of the to the left of the clear peak corresponding to a steep gradient in . The black dotted lines indicate the location of two spurious peaks identified in an automated search, which would be wrongly interpreted as positions of two discontinuities within the domain. This flavor of deformations present in the due to insufficient resolution causes the framework to detect discontinuities before they have fully formed and wrongly identify locations of shocks which are not present in the solution. As in the case of spurious oscillations in the solution, the artefacts depend on the expansion order and disappear at sufficiently high values of .

In summary, dynamically developing shocks interfere with the automated detection of discontinuities through the stencil’s inability to capture steep gradients. This results in premature shock identification and gives rise to localized Gibbs-like oscillations in the vicinity of unresolved gradients.

## 4 Shock detection improvements

In order to address the described issues we first check whether at a given time step the solution is continuous within our resolution limit. If that is not the case, we then determine whether the solution requires treatment appropriate for a genuine discontinuity or whether the resolution is not sufficient to capture steep gradients in the function.

### 4.1 Identification of resolved continuous cases

As described in Gelb2008 and Section 2.3 the function is designed to approximate the jump function of a discontinuous function . Thus, for smooth solutions vanishes for number of modes . The left panel of Fig. 3 demonstrates this property by showing how the function changes as a function of the number of modes for a smooth Gaussian (Eq. 10, , ) and a discontinuous ‘tophat’ function (Eq. 14, , ). As the number of modes increases, the Gaussian peaks decrease in amplitude with the curve approaching a flat line. In the ‘tophat’ case, however, the peaks do not vanish, retaining their amplitude and becoming narrower in shape as the number of modes increases.

The right panel of Fig 3 further quantifies these behaviors by plotting the variation in the peak height as a function of number of modes for the previously introduced tophat, Gaussian and a sine in the log-linear space. The plot clearly illustrates the exponential decay of the peak height for the smooth functions contrasted with the constant amplitude present in the discontinuous case. The slopes resulting from a linear fit in this space are negative for the smooth functions and vanishing for the discontinuous one.

In order to create Fig. 3 we initially decompose functions in Eq. (10) and (14) with N pseudospectral basis of Chebyshev polynomials. We then ‘downsample’ the initial representation to lower , obtaining spectral expansion coefficients in the new basis through the following transformation:

(19) |

where

(20) | ||||

(21) |

and , for . The matrix applied to expansion coefficients probes the spectral reprojection of order on the collocation points of the new basis of order . These collocation-point values are then translated to spectral coefficients via the inverse of the Vandermonde matrix, .

The simple test cases presented in Fig. 3 indicate that different smooth functions can potentially exhibit different decay slopes for the peak heights with increasing . Thus, in order to find a robust demarcation between purely smooth and marginally smooth or discontinuous solutions we investigate a set of different analytic functions and compute convergence slopes of the peak heights for different initial expansion orders . We then repeat this exercise for a test set of discontinuous problems, compare the results and choose a slope value which would allow us to robustly select smooth cases.

Fig. 4 presents the results of our investigations with the smooth functions tested in the left panel and the discontinuous ones in the right panel. The set of 25 smooth tests contains a variety of exponential, trigonometric, Gaussian-like and polynomial functions color-coded by their respective type labelled in the legend. As presented in the figure all slopes for the analytic test cases are while majority of the discontinuous ones are . The variation in slope values of the continuous functions comes from the variation in the amplitude of the input functions with higher amplitudes resulting in increasingly negative slopes. The polynomial case is characterized by small magnitudes of discontinuities, which escape the edge detection prescription and result in different slope values. After examining our results we choose to use the slope value of to robustly label a solution at a given time step as smooth.

### 4.2 Identification and treatment of unresolved steep gradients

When the solution fails the smoothness test at a given time step we then search for the potential shock locations by identifying local extrema in the function. As described earlier in Section 3 the unresolved gradients create spurious extrema revealed by the automated detection algorithm. Fig. 5 illustrates this issue by plotting a comparison between the functions resulting from a resolution limit and a genuine discontinuity. The former is shown in the left panel, where three spurious peaks in the function are visible, marked in blue dotted lines. Location of the steep gradient is marked in red dash-dotted vertical line. This case is in visible contrast with the right panel, where the function clearly exhibits a single, distinct peak and no spurious instances are identified.

In order to remove unwanted extrema we developed a three-step heuristic procedure performed as an additional part of the edge detection routine. It begins with an initial search for extrema as the first step, which identifies all potential peaks, including those marked with dotted lines in the left panel of Fig. 5.

In the second step, each of the initially found extrema are treated separately to verify their identification. For a given suggested jump location we construct a Gaussian smoothing kernel of width equal to the distance between the two collocation points and , neighboring the extremum:

(23) |

where .

We then convolve the with this kernel:

(24) |

to wash away any small peaks resulting from the inability of the fixed grid of collocation points to capture steep gradients in smooth functions.

As a final step we repeat the search for extrema restricting their width to no more than twice the distance between the adjacent collocation points. If a given extremum is identified in the repeated search it is accepted as a genuine feature, otherwise labelled as a spurious one. As illustrated in Fig. 5, the narrow features removed by smoothing do not show up in the search while the wide extremum is discarded by the maximum width criterion. Thus, only the extremum corresponding to the location of the steepest gradient in the solution is retained in this example.

In the edge detection routine, once the artefacts are recognized, the solution at a given time step is labelled as resolution-limited case. Such cases also require mollification treatment in order to remove Gibbs oscillations and we choose to use a continuous optimal mollifier as introduced in Eq. (4.2) of Tanner2006. This operation allows us to recover the solution without Gibbs oscillations, simultaneously preserving the steep gradient without excessive smoothing.

If, on the other hand no spurious extrema are confirmed in the repeated search, the case is labelled as genuinely discontinuous and the edge location found in the first step is correct (e.g. right panel of Fig. 5). Such a solution is then mollified using a prescription for one-sided mollifiers from Eq. (28) in Piotrowska2019 to recover the Gibbs-free solution along with its sharp discontinuity (shock) feature.

## 5 Application: 1D inviscid Burgers’ equation

We apply our improved shock detection techniques to a toy problem involving dynamical development of discontinuities from smooth solutions. We solve the 1D inviscid Burgers’ equation from Eq. (9) with a pseudospectral Chebyshev basis under periodic boundary conditions. As our initial conditions we choose a unit height Gaussian function in Eq. (10) with and . To evolve the nonlinear system in time we use the filter defined in Eq. (75) of Miller2016 applied to the time derivative to ensure stability of the numerical scheme. We allow the system to evolve until and then treat the Gibbs-affected solution with the edge detection and mollification techniques as described earlier in Sec. 4.

Fig. 6 presents snapshots of our solution at different times with the spectral reconstruction plotted in red dashed line and the mollified reconstruction overlaid as a black solid curve. The time progression is from left to right and from top to bottom with the initial conditions shown in the upper left corner and a snapshot in the right bottom one. The black curve is missing in the first two snapshots because the solution was recognized as smooth within the resolution limit and didn’t require post-processing to remove Gibbs-like oscillations. The third snapshot demonstrates how the resolution affects the vicinity of steep gradients and how well continuous mollifies cope with removing spurious information introduced locally into the solution. Our improved edge detection techniques are also capable of correctly identifying steep gradients located across the domain boundaries (middle bottom panel) allowing for robust treatment of marginally continuous cases located at the least favorable location within the computational domain. The last panel shows a solution with a clearly developed discontinuity where the Gibbs oscillations pollute the entire domain. The edge detection identified a single extremum in the function and hence one-sided mollification was successfully applied to recover the Gibbs-free solution while preserving its discontinuous character.

## 6 Summary and conclusions

Pseudospectral methods are a powerful means of solving PDEs which enjoy global exponential convergence for smooth problems. When applied to treat discontinuous functions, however, they lose their main advantages and suffer spurious oscillations in the solution which deems them unsuitable for astrophysical applications involving for e.g. fluid shocks. Over the years, multiple authors developed post-processing techniques allowing for a successful removal of the Gibbs phenomenon which we have previously tested in practice and applied to a toy advection problem.

In this work we focused on discontinuities which develop dynamically from smooth solutions, highlighting the challenges associated with their correct detection and treatment. We showed how steepening gradients lead to Gibbs-like fluctuations concentrated close to the discontinuities and how they influence the edge detection procedure causing problems in automated identification of the shock creation.

We then propose a heuristic solution to those problems, which allow for a robust automated identification of solutions requiring post-processing. This step is an inseparable part of any (pseudo)spectral shock capturing scheme, regardless of the choice of oscillation removal technique. Thus, in this work we offer solutions which can potentially benefit a wider community, encouraging further application of pseudospectral schemes to solving discontinuous problems.

Our improvements to jump identification utilize the rate of convergence of the function to differentiate perfectly smooth cases from the marginally smooth ones. We also suggest the use of Gaussian smoothing kernels on the function to identify solutions affected by the resolution limit. For these cases we then suggest treatment with symmetric adaptive mollifiers, while for the genuine shocks we recommend the one-sided kind. We tested this prescription on a range of smooth functions involving steep gradients located at random positions within the domain. The procedure was very successful at tackling steep gradients and genuine discontinuities across the whole domain, and suffered less problems close to the domain boundaries.

The combination of Gaussian smoothing and repeated peak search allows us to avoid problems akin to those highlighted in Piotrowska2019 such as spurious, narrow peaks originating close to the domain boundaries. Nevertheless, the method still suffered from occasional misidentification of genuine discontinuities as steep gradients, especially in cases when the discontinuity fell exactly at the domain boundary. Aware of this numerical weakness, we recommend exercising caution when dealing with discontinuities located close to the domain boundaries.

Finally, we solved a 1D inviscid Burgers’ equation to test the performance of our implemented methods in more realistic settings. Our framework was able to automatically determine when the gradients in the solution become steep enough to require post processing and treat the genuine discontinuity once it develops. Although we only tackled the inviscid Burger’s equation with a single set of initial conditions, we argue that our approach extends to arbitrary systems. In particular, one must introduce length and height scales, which need to be incorporated into the convergence rates shown in figure 4.

Encouraged by our progress with dynamically developing discontinuities for test cases we will apply our pseudospectral Chebyshev framework to solving real shock problems in the future. Benefiting from the apparent robustness of the scheme we would also like to extend it to two dimensions in order to tackle astrophysical problems with accuracy higher than the popular finite volume methods.

## Acknowledgements

We would like to thank Erik Schnetter, Josh Dolence and Chris Fryer for their support and constructive comments. JMP would also like to thank Roberto Maiolino and Asa Bluck for their excellent mentorship and valuable suggestions which increased the quality of this paper.

We gratefully acknowledge support from the U.S. Department of Energy (DOE) Office of Science and the Office of Advanced Scientific Computing Research via the SciDAC4 program and Grant DE-SC0018297, from the U.S. NSF grant AST-174267, and from the U.S. DOE through Los Alamos National Laboratory (LANL). This work used resources provided by the LANL Institutional Computing Program. Additional funding was provided by the LDRD Program and the Center for Nonlinear Studies at LANL under project number 20170508DR. LANL is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. DOE (Contract No. 89233218CNA000001).

This article is cleared for unlimited release, LA-UR-19-28857.

We are grateful to the countless developers contributing to open source projects on which we relied in this work, including Python (rossumPythonWhitePaper), numpy and scipy (numpy; scipyLib), and Matplotlib (hunterMatplotlib).