Spectral shock detection for dynamically developing discontinuities
SSpectral shock detection for dynamically developing discontinuities
Joanna Piotrowska a,b,c, ∗ , Jonah M. Miller a,d,e a Center for Nonlinear Studies, Los Alamos National Laboratory, Los Alamos, NM 87545, USA b Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK c Cavendish Laboratory, University of Cambridge, 9 JJ Thompson Avenue, Cambridge, CB3 0HE, UK d CCS-2, Los Alamos National Laboratory, Los Alamos, NM 87545, US e Center for Theoretical Astrophysics, Los Alamos National Laboratory, Los Alamos, NM 87545, US
Abstract
Pseudospectral schemes are a class of numerical methods capable of solving smooth problemswith high accuracy thanks to their exponential convergence to the true solution. Whenapplied to discontinuous problems, such as fluid shocks and material interfaces, due to theGibbs phenomenon, pseudospectral solutions lose their superb convergence and suffer fromspurious oscillations across the entire computational domain. Luckily, there exist theoreticalremedies for these issues which have been successfully tested in practice for cases of welldefined discontinuities. We focus on one piece of this procedure—detecting a discontinuityin spectral data. We show that realistic applications require treatment of discontinuitiesdynamically developing in time and that it poses challenges associated with shock detection.More precisely, smoothly steepening gradients in the solution spawn spurious oscillations dueto insufficient resolution, causing premature shock identification and information loss. Weimprove existing spectral shock detection techniques to allow us to automatically detect truediscontinuities and identify cases for which post-processing is required to suppress spuriousoscillations resulting from the loss of resolution. We then apply these techniques to solve aninviscid Burgers’ equation in 1D, demonstrating that our method correctly treats genuineshocks caused by wave breaking and removes oscillations caused by numerical constraints.
Keywords: numerical methods, shock detection, pseudospectral schemes
1. Introduction
Pseudospectral methods are a powerful numerical approach, offering exponential conver-gence to the true solution across the entire domain for smooth functions (see [1] for a review).Thanks to their exquisite accuracy and efficiency they have been successfully applied in, forexample, large scale relativistic astrophysics simulations (e.g. [2, 3, 4]). A wealth of as-trophysical phenomena, however, involve discontinuities such as fluid shocks in supernova ∗ I am the corresponding author
Email addresses: [email protected] (Joanna Piotrowska), [email protected] (Jonah M. Miller)
Preprint submitted to Journal of Computational Physics October 3, 2019 a r X i v : . [ m a t h . NA ] O c t xplosions or galactic spiral arms. For this class of problem, pseudospectral schemes per-form poorly due to sublinear convergence properties and spurious oscillations induced in thesolution known as the Gibbs phenomenon [5, 6, 7].A substantial body of research has been conducted in the attempt to alleviate the Gibbsphenomenon and recover the exponential convergence of pseudospectral schemes. Someproposed solutions are reprojection onto combinations of different basis sets [8, 9], filtering[10], artificial viscosity [11], modelling and removal of oscillations [12, 13] and mollification[14]. In [15] we explored and numerically tested the last approach, using an optimal mol-lifier as prescribed by [16]. We showed successful recovery of convergence properties awayfrom discontinuities and offered improvements to preserve the discontinuous character ofthe solution. We also solved a toy problem advecting a discontinuous function, demon-strating greater robustness of mollification as compared to the more popular Gegenbauerreconstruction methods.In this work we take the next step towards applying pseudospectral methods to solvephysically relevant shock problems. This time we allow the discontinuities to develop dynam-ically, instead of introducing them by construction and apply previously explored methodsto detect the created shocks. Robust identification of evolving discontinuities is crucialin pseudospectral shock capturing, regardless of the choice of Gibbs phenomenon removalmethod. We discuss challenges associated with the detection of discontinuity onset andcapturing of steep gradients at constant spatial resolution. We also present our heuristicsolution to those issues, offering practical improvements to shock detection methods in caseof shocks dynamically developing during the simulation. Finally, we solve a toy problem ofthe inviscid Burgers’ equation to demonstrate the performance of our methods when appliedto shocks developing smoothly over time.
2. Background
In this section we review existing methods which provide basis for our improvementspresented further in Section 4.
The pseudospectral decomposition is a means of discretizing a function of interest witha finite sum of N orthogonal basis functions T ( x ) n : u ( x ) ≈ S N [ u ]( x ) = N (cid:88) n =0 ˆ u n T n ( x ) , (1)where S N [ u ] is the partial sum of u and the basis functions are denoted T n ( x ). The expansioncoefficients ˆ u n are chosen such that lim N →∞ S N [ u ]( x ) = u ( x ) (2)2nd are given by: ˆ u n = ( u, T n )( T n , T n ) (3)where ( u, T n ) = (cid:90) Ω u ( x ) T n ( x ) w ( x ) dx (4)is the inner product between u and T n on a domain Ω. The integral can be approximatedvia Gauss quadrature in which case the x i are the quadrature collocation points and w i theirassociated weights, such that:lim N →∞ N (cid:88) i =0 u ( x i ) T n ( x i ) w i = (cid:90) Ω u ( x ) T n ( x ) w ( x ) dx . (5)Thus, one can conveniently transform the collocation-point values u ( x i ) ≡ u i into expan-sion coefficients ˆ u n through : ˆ u n = V ni u i , V ni ≡ T n ( x i ) w i ( T n , T n ) (6)where V ni is the inverse of the Vandermonde matrix . One then can also approximate thepartial sum of the derivative of u : S N [ ∂ x u ] with the following differentiation matrices: M mn = ( ∂ x T m , T n )( T n , T n ) , N ji = V − jn M nm V mi , (7)such that differentiation conveniently reduces to the following matrix operations:lim N →∞ dd x ˆ u m = M mn ˆ u n , lim N →∞ dd x u j = N ji u i . (8)In our work we use a basis of Chebyshev polynomials of the first kind defined on thedomain Ω = [ − ,
1] and the
Chebyshev-Gauss-Lobatto quadrature. The associated grid ofcollocation points allows us to probe the solution values at the boundaries with x = − x N = 1. For a detailed description of Gaussian quadrature and transformations between thespectral and collocation-point bases we refer the interested reader to [17] and [1]. where (cid:80) i a i b i = a i b i in the Einstein notation other choices are possible, depending on the problem of interest e.g. a Legendre polynomial or Fourierbasis .2. Stabilization techniques for non-linear PDEs In this research we focus on one-dimensional (1D) problems which smoothly developdiscontinuities over time, such as the non-linear
1D inviscid Burgers’ equation : ∂ t u + u ∂ x u = 0 . (9)In our investigation we assume smooth initial conditions, choosing a Gaussian g ( x ) cen-tered on x = 0: g ( x ) = exp (cid:18) − ( x − x ) σ (cid:19) , (10)where x = 0 and σ = 0 . N + 1 ordinary differential equations (ODEs) at the N + 1 collocation points. Thus, wesolve Eq. (9) using the method of lines, keeping time t as a continuous variable and assumingperiodic boundary conditions such that at each time step u = u ( −
1) is set to u N = u (1).Note that our initial condition is non-smooth at the boundaries, since g ( x ) is non-periodic.Our periodic boundary conditions serve as a simplistic proxy for a multi-domain spectralmethod, with the left and right boundaries serving as inflow and outflow into other spectraldomains respectively. This allows us to run our simulation for arbitrary long time, despitethe 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 largeshort-wavelength coefficients which need to be subdued to avoid a numerical catastrophe. Inthis work we choose a modification to the spectral viscosity technique [11], namely the righthand side filtering introduced by [18] in Eq. (73-75). In this framework filtering is appliedto the derivative operator modes via the following diagonal matrix: F mn := (cid:16) nN (cid:17) s δ mn , (11)where s is the so-called order of dissipation and we chose s = 2 in the evolution of oursystem which proved optimal in our numerical tests. The temporal derivative at each gridpoint then takes the form: ∂ t u i = u i N ij u j − c V − im F mn V nj u j , (12)where the constant c depends on the particular problem and c ∼ .
01 performed best whentested with the inviscid Burgers’ equation. 4 .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 u − S N [ u ] integrated across the computational domain as a function of expansion order N : (cid:107) u − S N [ u ] (cid:107) ≤ αN p p (cid:88) k =0 (cid:13)(cid:13) u ( k ) (cid:13)(cid:13) , (13)where 0 ≤ p ≤ p max , p max is the highest order of derivative well-defined for u , α is a constantfactor and (cid:107) f ( x ) (cid:107) ≡ (cid:82) Ω | f ( x ) | dx denotes the L norm. In the smooth case, arbitrarily high-order derivatives are well-defined and the inequality holds for all p > L norm exhibiting at best O (cid:0) N − (cid:1) convergence [20]. Further-more, the pseudospectral representations of the solution suffer from spurious oscillations.The oscillation amplitude decreases away from the discontinuities and their frequency variesdepending on the position within the domain. The left panel of Fig. 1 illustrates this phe-nomenon by overlaying a Gibbs-affected N = 60 reconstruction on the input discontinuouspiecewise linear ‘tophat’ function: h ( x ) = (cid:40) x < x < x , (14)where x = − . x = − . edgedetection and Gibbs oscillations removal . Edge detection is a step necessary for many Gibbsphenomenon removal techniques, allowing construction of adaptive smoothing kernels in the mollification approach or the correct domain decomposition in the case of the Gegenbauerreconstruction [22].Following our work in [15], we choose to locate the edges using the minmod techniqueintroduced by [21]. It relies on applying a range of different concentration factors µ ( kN ) [23]to the partial sum of the function derivative ˜ S N u (cid:48) ( x ), which accelerate the rate at which˜ S N u (cid:48) ( x ) converges onto the jump function ,[ u ]( x ) = u ( x + ) − u ( x − ) . (15)In this prescription, each choice of an admissible µ ( kN ) generates a jump function approxi-mation j µ ( x ) given by the following partial sum: j µ ( x ) = π √ − x N N (cid:88) k =1 µ (cid:18) kN (cid:19) ˆ u k T (cid:48) k ( x ) . (16)The jump function approximation then converges onto the true jump function at thefollowing rate [21]: | j µ ( x ) − [ u ]( x ) | ≤ Const · log NN (17)5 . − . . . . x . . . f ( x ) Pseudospectralreconstruction originalreconstructed − . − . . . . x − . − . . . Jump functions j − n and theresulting minmod ( j · · · j n ) minmod ( j , j , j ) j j j − . − . . . . x . . . f ( x ) Mollified reconstruction reconstructedmollified
Figure 1: Left panel: demonstration of the Gibbs oscillations affecting the pseudospectral reconstruction (reddashed line) of a discontinuous ‘tophat’ function (Eq. (14), blue solid line). Middle panel: exemplary jumpfunctions and the resulting minmod function from [21]. The extrema indicate locations of the discontinuitiesin the ‘tophat’. Right panel: mollified reconstruction (black solid line) post-processed using one-sidedmollifiers from [15]. The reconstruction is free from the oscillations while preserving the sharp character ofits discontinuities. This illustration was made for expansion order N = 60. and each j µ ( x ) exhibits a different oscillation pattern around the discontinuity associatedwith a given µ ( kN ) [24]. The minmod prescription then takes advantage of the differences inthe oscillatory patters, combining a range of jump functions such that the oscillations areremoved [21]: minmod ( j , . . . , j n )( x ) := min( j ( x ) , . . . , j n ( x )) , if j ( x ) , . . . , j n ( x ) > j ( x ) , . . . , j n ( x )) , if j ( x ) , . . . , j n ( x ) < , otherwise. (18)The middle panel of Fig. 1 demonstrates three jump function approximations j µ ( x ) fromtrigonometric, polynomial and exponential µ ( kN ) forms as well as a resulting minmod ( x ) for N = 60. The location of discontinuities is indicated by the positions of the peaks in thefunction represented by the black solid curve. In our work we construct the minmod using µ ( kN ) defined in Eq. (8-10) in [21] and their filtered versions constructed using Lanczos [25]filters of orders 0, 1, 2 and 3.Once we identify the discontinuities and their positions we then remove the oscillationsthrough mollification with adaptive one-sided mollifiers [16, 15]. The right panel of Fig. 1demonstrates the result of post processing by overplotting the mollified result (black solid) onthe Gibbs-affected reconstruction (red dashed). The post-processed solution is successfullyfreed from spurious oscillations, simultaneously preserving its discontinuous character.Finally, it is important to mention, that even in the ideal cases of well-defined discon-tinuities the minmod function is polluted by residual oscillations of small amplitude whichcan be misidentified as jumps in an automated search for extrema. Hence one always has tomake a choice of a peak height qualifying as a true extremum. As a result, one then sets therange of discontinuity heights captured by the edge detection framework, which depends onthe resolution and the problem at hand. 6 . − . . . . x . . . . . u ( x ) Solution of Burgers’ equation at t = 0 . − . − . . . . x − . − . − . . . . m i n m o d ( x ) Minmod function at t = 0 . spurious peaksteepest gradient Figure 2: Issues arising in discontinuities developing dynamically from smooth solutions. Left panel: mildoscillations in the solution present in the vicinity of steep gradients. The pseudospectral basis fails to captureoverly steep gradients due to its resolution limit. Right panel: corresponding minmod function polluted byspurious peaks. These are misidentified as shock locations in an automated search.
3. New shock detection challenges
The edge detection techniques illustrated in Fig. 1 are an implementation of prescriptionsdeveloped by A. Gelb and collaborators [24, 21] for stationary discontinuities. Realistic nu-merical applications, however, involve dynamical development of discontinuities in time andthus require the pseudospectral framework to capture transitions between the smooth anddiscontinuous (shocked) regimes in an automated fashion. Thus, one requires a robust meansof automatically recognizing discontinuous functions to deploy post-processing techniquesat appropriate time steps.Solutions in which shocks develop smoothly over time experience additional problemscaused by insufficient resolution of the truncated pseudospectral expansion basis. An ex-ample of such is the inviscid Burgers’ equation (eq. 9) where the spatial gradient in thefunction steepens continuously until the point of wave breaking. In such cases the solutionis populated by additional oscillations around steep gradients which result from the lack ofsufficient resolution to capture significant variation at very small scales. The left panel ofFig. 2 demonstrates this issue by showing the affected solution u ( x ) of the 1D inviscid Burg-ers’ equation 9 with smooth Gaussian initial conditions (Eq. 10 with x = 0 and σ = 0 . t = 0 .
12 which has a steep gradient at x ≈ .
45. The solution in the vicinity of the steepslope in u ( x ) is populated by spurious oscillations of small amplitude which introduce un-desired, incorrect, information. These oscillations are an artefact of insufficient resolutionand their amplitude converges away with increasing expansion order N .What is more, the resolution limit also demonstrates itself in shock detection, altering7he minmod function such that it exhibits different convergence properties and featuresadditional peaks. The right panel of Fig. 2 highlights the latter by showing the irregularshape of the minmod to the left of the clear peak corresponding to a steep gradient in u ( x ).The black dotted lines indicate the location of two spurious peaks identified in an automatedsearch, which would be wrongly interpreted as positions of two discontinuities within thedomain. This flavor of deformations present in the minmod due to insufficient resolutioncauses the framework to detect discontinuities before they have fully formed and wronglyidentify locations of shocks which are not present in the solution. As in the case of spuriousoscillations in the solution, the minmod artefacts depend on the expansion order N anddisappear at sufficiently high values of N .In summary, dynamically developing shocks interfere with the automated detection ofdiscontinuities through the stencil’s inability to capture steep gradients. This results inpremature shock identification and gives rise to localized Gibbs-like oscillations in the vicinityof unresolved gradients.
4. Shock detection improvements
In order to address the described issues we first check whether at a given time step thesolution is continuous within our resolution limit. If that is not the case, we then determinewhether the solution requires treatment appropriate for a genuine discontinuity or whetherthe resolution is not sufficient to capture steep gradients in the function.
As described in [21] and Section 2.3 the minmod function is designed to approximatethe jump function [ f ]( x ) = f ( x + ) − f ( x − ) of a discontinuous function f ( x ). Thus, forsmooth solutions minmod vanishes for number of modes N → ∞ . The left panel of Fig. 3demonstrates this property by showing how the minmod function changes as a functionof the number of modes N for a smooth Gaussian (Eq. 10, x = − . σ = 0 .
15) and adiscontinuous ‘tophat’ function (Eq. 14, x = − . x = − . mindmod peaks decrease in amplitude with the curve approachinga flat line. In the ‘tophat’ case, however, the minmod peaks do not vanish, retaining theiramplitude 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 inthe minmod peak height as a function of number of modes N for the previously introducedtophat, Gaussian and a sine f ( x ) = sin(6 x ) in the log-linear space. The plot clearly illus-trates the exponential decay of the peak height for the smooth functions contrasted with theconstant amplitude present in the discontinuous case. The slopes resulting from a linear fitin 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) withN= 60 pseudospectral basis of Chebyshev polynomials. We then ‘downsample’ the initial8 . − . . x − . − . − . . . . . m i n m o d ( x ) Gaussian − . − . . x tophat 30405060 N
30 40 50 60 N . . . l og ( h p e a k / h m p e a k , i n i t i a l ) slope -0.04slope -0.03 slope 0.00 minmod peak height fordownsampled functions gaussiansinetophat Change in the minmod functionfor downsampled orders N Figure 3: Left panel: behavior of the minmod function under a change of the expansion order N . Thecolor of each line corresponds to an order labelled in the colorbar. For a smooth function the minmod peaks converge to 0, while they retain constant amplitude for the discontinuous one. Right panel: change ofminmod peak height with the number of modes sampled downward from N=60. Continuous functions showexponential convergence of the minmod peak height with linear slopes in the log-linear space of slope ≈ − . − .
03. The discontinuous minmod height stays approximately constant with a vanishing slope in thelog-linear space. representation to lower N = K , obtaining spectral expansion coefficients in the new basis a k through the following transformation: a k = D km a m , (19)where D km = V kj P jm , (20) P jm = M (cid:88) m =0 T m ( x j ) (21)(22)and [ m ] = { , · · · M } , [ j ] = { , · · · K } for K < M . The P jm matrix applied to expansioncoefficients a m probes the spectral reprojection of order M on the K + 1 x j collocationpoints of the new basis of order K < M . These collocation-point values are then translatedto spectral coefficients a k via the inverse of the Vandermonde matrix, V .The simple test cases presented in Fig. 3 indicate that different smooth functions canpotentially exhibit different decay slopes for the minmod peak heights with increasing N .Thus, in order to find a robust demarcation between purely smooth and marginally smoothor discontinuous solutions we investigate a set of different analytic functions and computeconvergence slopes of the minmod peak heights for different initial expansion orders N . Wethen repeat this exercise for a test set of discontinuous problems, compare the results andchoose a slope value which would allow us to robustly select smooth cases.9 N − . − . − . − . . s l o p e Smooth downsamplingslopes cutoffexponentialtrigonometricgaussianpolynomial 30 40 50 60 N − . − . − . . . . s l o p e Discontinuousdownsampling slopes cutofftophatsawtoothdiscontinuous cosinediscontinuous polynomial
Figure 4: Convergence slopes in the log-linear space for smooth (left panel) and discontinuous (right panel)functions explored for different initial expansion orders N . The investigated range of continuous cases isalways characterized by a negative slope, while majority of discontinuous cases by a positive one. Anexperimentally suggested demarcation to distinguish purely smooth functions is plotted as a black dashedline. Fig. 4 presents the results of our investigations with the smooth functions tested in theleft panel and the discontinuous ones in the right panel. The set of 25 smooth tests containsa variety of exponential, trigonometric, Gaussian-like and polynomial functions color-codedby their respective type labelled in the legend. As presented in the figure all minmod slopesfor the analytic test cases are < − .
015 while majority of the discontinuous ones are (cid:38) − .
125 to robustly label a solution at a giventime step as smooth.
When the solution fails the smoothness test at a given time step we then search forthe potential shock locations by identifying local extrema in the minmod function. Asdescribed earlier in Section 3 the unresolved gradients create spurious extrema revealed bythe automated detection algorithm. Fig. 5 illustrates this issue by plotting a comparisonbetween the minmod functions resulting from a resolution limit and a genuine discontinuity.10 . . x − . − . . . m i n m o d ( x ) Resolution-limited case minmod smoothed minmod correct peakfalse peaks − . − . . x − . − . − . . . m i n m o d ( x ) Real discontinuity minmod smoothed minmod correct peak minmod functions with applied Gaussian smoothing
Figure 5: Smoothing procedure which aids in distinguishing clear discontinuous cases from resolution-drivenoscillations. Left panel: resolution limited case, where the minmod function is polluted by spurious peaks.Initially, multiple peaks are recognised, which are then rejected in secondary search after Gaussian smoothingis applied. Right panel: true discontinuity case where no spurious peaks are identified in which case thenumber of detected peaks does not change after applying gaussian smoothing.
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-dottedvertical line. This case is in visible contrast with the right panel, where the minmod functionclearly exhibits a single, distinct peak and no spurious instances are identified.In order to remove unwanted extrema we developed a three-step heuristic procedureperformed as an additional part of the edge detection routine. It begins with an initialsearch for extrema as the first step, which identifies all potential peaks, including thosemarked with dotted lines in the left panel of Fig. 5.In the second step, each of the initially found extrema are treated separately to verifytheir identification. For a given suggested jump location we construct a Gaussian smooth-ing kernel of width equal to the distance between the two collocation points x i and x i +1 ,neighboring the extremum: k ( x − x (cid:48) ) = exp (cid:18) − ( x − x (cid:48) ) ω (cid:19) , (23)where ω = ( x i +1 − x i ) / minmod with this kernel: minmod ( x ) smoothed = (cid:90) − minmod ( x (cid:48) ) k ( x − x (cid:48) ) dx (cid:48) (24)to wash away any small peaks resulting from the inability of the fixed grid of collocationpoints to capture steep gradients in smooth functions.11s a final step we repeat the search for extrema restricting their width to no more thantwice the distance between the adjacent collocation points. If a given extremum is identifiedin the repeated search it is accepted as a genuine feature, otherwise labelled as a spuriousone. As illustrated in Fig. 5, the narrow features removed by smoothing do not show up inthe 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 isretained in this example.In the edge detection routine, once the artefacts are recognized, the solution at a giventime step is labelled as resolution-limited case. Such cases also require mollification treat-ment in order to remove Gibbs oscillations and we choose to use a continuous optimalmollifier as introduced in Eq. (4.2) of [16]. This operation allows us to recover the solutionwithout Gibbs oscillations, simultaneously preserving the steep gradient without excessivesmoothing.If, on the other hand no spurious extrema are confirmed in the repeated search, thecase is labelled as genuinely discontinuous and the edge location found in the first step iscorrect (e.g. right panel of Fig. 5). Such a solution is then mollified using a prescription forone-sided mollifiers from Eq. (28) in [15] to recover the Gibbs-free solution along with itssharp discontinuity (shock) feature.
5. Application: 1D inviscid Burgers’ equation
We apply our improved shock detection techniques to a toy problem involving dynam-ical development of discontinuities from smooth solutions. We solve the 1D inviscid Burg-ers’ equation from Eq. (9) with a N = 60 pseudospectral Chebyshev basis under periodicboundary conditions. As our initial conditions we choose a unit height Gaussian function inEq. (10) with x = 0 and σ = 0 .
15. To evolve the nonlinear system in time we use the filterdefined in Eq. (75) of [18] applied to the time derivative to ensure stability of the numericalscheme. We allow the system to evolve until t = 3 . t = 0 .
48 snapshot in the right bottom one.The black curve is missing in the first two snapshots because the solution was recognized assmooth within the resolution limit and didn’t require post-processing to remove Gibbs-likeoscillations. The third snapshot demonstrates how the resolution affects the vicinity of steepgradients and how well continuous mollifies cope with removing spurious information intro-duced locally into the solution. Our improved edge detection techniques are also capable ofcorrectly identifying steep gradients located across the domain boundaries (middle bottompanel) allowing for robust treatment of marginally continuous cases located at the least fa-vorable location within the computational domain. The last panel shows a solution with12 .02.53.03.5 φ ( x ) t = 0 . t = 0 . t = 0 . x φ ( x ) t = 0 . reconstructedmollified -0.5 0.0 0.5 1.0 x t = 0 .
34 -0.5 0.0 0.5 1.0 x t = 0 . Real space projection of the solution
Figure 6: Snapshots from pseudospectral solution of the 1D inviscid Burgers’ equation post processed withour improved shock detection and appropriate choice of mollification technique. The red dashed lines showthe raw solution while the black solid ones our post-processed result. The time progression is left to rightand top to bottom. The first two snapshots do not feature mollification because the solution was recognisedas smooth. The solid curves show a solution free from spurious oscillations, smoothly evolving to form adiscontinuity captured by our method. a clearly developed discontinuity where the Gibbs oscillations pollute the entire domain. Theedge detection identified a single extremum in the minmod function and hence one-sidedmollification was successfully applied to recover the Gibbs-free solution while preserving itsdiscontinuous character.
6. Summary and conclusions
Pseudospectral methods are a powerful means of solving PDEs which enjoy global expo-nential convergence for smooth problems. When applied to treat discontinuous functions,however, they lose their main advantages and suffer spurious oscillations in the solutionwhich deems them unsuitable for astrophysical applications involving for e.g. fluid shocks.Over the years, multiple authors developed post-processing techniques allowing for a suc-cessful removal of the Gibbs phenomenon which we have previously tested in practice andapplied to a toy advection problem.In this work we focused on discontinuities which develop dynamically from smooth so-lutions, highlighting the challenges associated with their correct detection and treatment.We showed how steepening gradients lead to Gibbs-like fluctuations concentrated close to13he discontinuities and how they influence the edge detection procedure causing problemsin automated identification of the shock creation.We then propose a heuristic solution to those problems, which allow for a robust au-tomated identification of solutions requiring post-processing. This step is an inseparablepart of any (pseudo)spectral shock capturing scheme, regardless of the choice of oscilla-tion removal technique. Thus, in this work we offer solutions which can potentially benefita wider community, encouraging further application of pseudospectral schemes to solvingdiscontinuous problems.Our improvements to jump identification utilize the rate of convergence of the mimod function to differentiate perfectly smooth cases from the marginally smooth ones. We alsosuggest the use of Gaussian smoothing kernels on the minmod function to identify solutionsaffected by the resolution limit. For these cases we then suggest treatment with symmetricadaptive mollifiers, while for the genuine shocks we recommend the one-sided kind. Wetested this prescription on a range of smooth functions involving steep gradients located atrandom positions within the domain. The procedure was very successful at tackling steepgradients and genuine discontinuities across the whole domain, and suffered less problemsclose to the domain boundaries.The combination of Gaussian smoothing and repeated peak search allows us to avoidproblems akin to those highlighted in [15] such as spurious, narrow peaks originating closeto the domain boundaries. Nevertheless, the method still suffered from occasional misidenti-fication of genuine discontinuities as steep gradients, especially in cases when the discontinu-ity fell exactly at the domain boundary. Aware of this numerical weakness, we recommendexercising 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 im-plemented methods in more realistic settings. Our framework was able to automaticallydetermine when the gradients in the solution become steep enough to require post process-ing and treat the genuine discontinuity once it develops. Although we only tackled theinviscid Burger’s equation with a single set of initial conditions, we argue that our approachextends 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 caseswe will apply our pseudospectral Chebyshev framework to solving real shock problems inthe future. Benefiting from the apparent robustness of the scheme we would also like toextend it to two dimensions in order to tackle astrophysical problems with accuracy higherthan the popular finite volume methods.
Acknowledgements
We would like to thank Erik Schnetter, Josh Dolence and Chris Fryer for their supportand constructive comments. JMP would also like to thank Roberto Maiolino and Asa Bluckfor their excellent mentorship and valuable suggestions which increased the quality of thispaper. 14e gratefully acknowledge support from the U.S. Department of Energy (DOE) Officeof Science and the Office of Advanced Scientific Computing Research via the SciDAC4program and Grant DE-SC0018297, from the U.S. NSF grant AST-174267, and from theU.S. DOE through Los Alamos National Laboratory (LANL). This work used resourcesprovided by the LANL Institutional Computing Program. Additional funding was providedby the LDRD Program and the Center for Nonlinear Studies at LANL under project number20170508DR. LANL is operated by Triad National Security, LLC, for the National NuclearSecurity 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 whichwe relied in this work, including Python [26], numpy and scipy [27, 28], and Matplotlib [29].
References [1] J. Boyd, Chebyshev and Fourier Spectral Methods: Second Revised Edition, Dover Books on Mathe-matics, Dover Publications, 2013.URL https://books.google.com/books?id=b4TCAgAAQBAJ [2] L. Kidder, H. Pfeiffer, M. Scheel, et al., Spectral Einstein Code (2000).URL [3] L. Boyle, M. Kesden, S. Nissanke, Binary black-hole merger: Symmetry and the spin expansion, Phys.Rev. Lett. 100 (2008) 151101. doi:10.1103/PhysRevLett.100.151101 .URL http://link.aps.org/doi/10.1103/PhysRevLett.100.151101 [4] M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. Matthews, H. P. Pfeiffer, High-accuracy waveformsfor binary black hole inspiral, merger, and ringdown, Phys. Rev. D 79 (2009) 024003. doi:10.1103/PhysRevD.79.024003 .URL http://link.aps.org/doi/10.1103/PhysRevD.79.024003 [5] H. Wilbraham, On a certain periodic function, Cambridge and Dublin Math. J 3 (198) (1848) 1848.[6] J. W. Gibbs, Fouriers series, Nature 59 (1522) (1898) 200.[7] A. A. Michelson, Fourier’s Series, Nature 59 (1898) 200. doi:10.1038/059200a0 .[8] D. Gottlieb, C.-W. Shu, A. Solomonoff, H. Vandeven, On the Gibbs phenomenon 1: Recovering ex-ponential accuracy from the Fourier partial sum of a non-periodic analytic function 43 (1992) 81–98. doi:10.1016/0377-0427(92)90260-5 .URL http://hdl.handle.net/2060/19920015730 [9] S. Gottlieb, J. H. Jung, S. Kim, A review of David Gottlieb’s work on the resolution of the Gibbsphenomenon, Communications in Computational Physics 9 (3) (2011) 497–519. doi:10.4208/cicp.301109.170510s .[10] H. Vandeven, Family of spectral filters for discontinuous problems, Journal of Scientific Computing6 (2) (1991) 159–192. doi:10.1007/BF01062118 .[11] E. Tadmor, Shock capturing by the spectral viscosity method, Computer Methods in Applied Mechanicsand Engineering 80 (1-3) (1990) 197–208. doi:10.1016/0045-7825(90)90023-F .[12] V. Krylov, A. Stroud, Approximate Calculation of Integrals, Dover Books on Mathematics, DoverPublications, 2006.URL https://books.google.com/books?id=lswsAwAAQBAJ [13] Y. Lipman, D. Levin, Approximating piecewise-smooth functions, IMA Journal of Numerical Analysis30 (4) (2010) 1159–1183. doi:10.1093/imanum/drn087 .URL http://dx.doi.org/10.1093/imanum/drn087 [14] D. Gottlieb, E. Tadmor, Recovering pointwise values of discontinuous data within spectral accuracy,in: E. M. Murman, S. S. Abarbanel (Eds.), Progress and Supercomputing in Computational FluidDynamics, Progress in Scientific Computing, Springer, Birkh¨auser Boston, 1985.
15] J. Piotrowska, J. M. Miller, E. Schnetter, Spectral methods in the presence of discontinuities, Journal ofComputational Physics 390 (2019) 527 – 547. doi:https://doi.org/10.1016/j.jcp.2019.03.048 .URL [16] J. Tanner, Optimal filter and mollifier for piecewise smooth spectral data, Mathematics of Computation75 (254) (2006) 767–790.URL [17] P. Grandclement, Introduction to spectral methods, EAS Publications Series 21. doi:10.1051/eas:2006112 .[18] J. M. Miller, E. Schnetter, An operator-based local discontinuous galerkin method compatible withthe BSSN formulation of the einstein equations, Classical and Quantum Gravity 34 (1) (2016) 015003. doi:10.1088/1361-6382/34/1/015003 .URL https://doi.org/10.1088%2F1361-6382%2F34%2F1%2F015003 [19] P. Grandcl´ement, J. Novak, Spectral methods for numerical relativity, Living Reviews in Relativity12 (1) (2009) 1. doi:10.12942/lrr-2009-1 .URL https://doi.org/10.12942/lrr-2009-1 https://books.google.com/books?id=vA9d57GxCKgC [21] A. Gelb, D. Cates, Detection of edges in spectral data iii—refinement of the concentration method,Journal of Scientific Computing 36 (1) (2008) 1–43. doi:10.1007/s10915-007-9170-8 .[22] S. Gottlieb, J. H. Jung, S. Kim, A review of David Gottlieb’s work on the resolution of the Gibbsphenomenon, Communications in Computational Physics 9 (3) (2011) 497–519. doi:10.4208/cicp.301109.170510s .[23] A. Gelb, E. Tadmor, Detection of Edges in Spectral Data, Applied and Computational HarmonicAnalysis 7 (1999) 101–135. arXiv:0112016 , doi:10.1137/S0036142999359153 .URL http://arxiv.org/abs/math/0112016 [24] A. Gelb, E. Tadmor, Detection of edges in spectral data ii. nonlinear enhancement, SIAM Journal onNumerical Analysis 38 (4) (2000) 1389–1408. arXiv:https://doi.org/10.1137/S0036142999359153 , doi:10.1137/S0036142999359153 .URL https://doi.org/10.1137/S0036142999359153 [25] C. Lanczos, Discourse on Fourier series, University mathematical monographs, Oliver & Boyd, 1966.URL https://books.google.co.uk/books?id=yVDvAAAAMAAJ [26] G. Rossum, Python reference manual, Tech. rep., Amsterdam, The Netherlands, The Netherlands(1995).[27] S. van der Walt, S. C. Colbert, G. Varoquaux, The numpy array: A structure for efficient numericalcomputation, Computing in Science Engineering 13 (2) (2011) 22–30. doi:10.1109/MCSE.2011.37 .[28] E. Jones, et al., SciPy: Open source scientific tools for Python, (2001–).[29] J. D. Hunter, Matplotlib: A 2D graphics environment, Computing In Science & Engineering 9 (3)(2007) 90–95.(2001–).[29] J. D. Hunter, Matplotlib: A 2D graphics environment, Computing In Science & Engineering 9 (3)(2007) 90–95.