pith. sign in

arxiv: 2606.05155 · v1 · pith:I4N3FM2Knew · submitted 2026-06-03 · 🌀 gr-qc · cs.NA· math.NA

High-Order Summation-By-Parts Schemes for First-Order Hyperbolic Systems in Curvilinear Coordinates with Singularities

Stamatis Vretinaris , Erik Schnetter This is my paper

Pith reviewed 2026-06-28 05:03 UTC · model grok-4.3

classification 🌀 gr-qc cs.NAmath.NA
keywords summation-by-partsfinite differencesspherical coordinateshyperbolic systemsnumerical stabilitycurvilinear coordinateshigh-order schemesorigin singularity
0
0 comments X

The pith

Discrete gradient and divergence operators can obey summation-by-parts on spherical grids that place a point at the origin.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper constructs high-order finite-difference operators for first-order hyperbolic systems in spherical coordinates that remain energy stable despite the coordinate singularity at r=0. It achieves this by defining discrete gradient and divergence operators whose summation-by-parts identity mirrors the continuous integration-by-parts relation even when a 1/r^p singularity is present. Explicit constructions are given through sixth order, and the operators place a grid point directly at the origin rather than straddling it. The work extends the 2013 framework of Gundlach et al. and demonstrates the operators on the scalar wave equation. A reader cares because such operators enable stable, high-accuracy long-time evolutions of waves and fluids near compact objects without artificial coordinate patches.

Core claim

We present a method for constructing high-order accurate, energy-stable finite difference operators satisfying the Summation-by-Parts (SBP) property on spherical domains. We define discrete gradient and divergence operators that mirror the continuous integration-by-parts principle, even though there is a 1/r^p coordinate singularity present at the origin. We explicitly construct such operators up to order six. Our operators place a grid point directly on the origin.

What carries the argument

Summation-by-parts (SBP) finite-difference operators for the gradient and divergence that satisfy a discrete integration-by-parts identity on grids containing r=0.

If this is right

  • Energy stability follows directly for any first-order hyperbolic system once the SBP property holds.
  • Accuracy reaches sixth order at every point, including the origin, without special limiting procedures.
  • Time-step restrictions are set by the spectral radius of the resulting spatial operator.
  • The same construction applies to systems with variable coefficients or nonlinear terms provided the continuous problem satisfies an energy estimate.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • The operators could reduce the need for multiple coordinate patches in spherical simulations.
  • Similar constructions might extend to other singular coordinate systems such as cylindrical or toroidal grids.
  • Long-term accuracy in wave propagation near the origin becomes limited mainly by dispersion rather than stability loss.

Load-bearing premise

Finite-difference operators satisfying the summation-by-parts property can be constructed on a spherical domain containing the origin by placing a grid point directly at r=0 while still mirroring the continuous integration-by-parts identity.

What would settle it

Evolve the scalar wave equation with these sixth-order operators on a spherical grid including the origin and check whether the discrete energy remains bounded over many light-crossing times.

Figures

Figures reproduced from arXiv: 2606.05155 by Erik Schnetter, Stamatis Vretinaris.

Figure 1
Figure 1. Figure 1: FIG. 1: Discretization error when applying 4th-order discrete divergence operator to [PITH_FULL_IMAGE:figures/full_fig_p011_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: FIG. 2: Discretization error when applying 6th-order discrete divergence operator to [PITH_FULL_IMAGE:figures/full_fig_p011_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: FIG. 3: Time evolution for the 6th-order non-diagonal SBP operator in spherical symmetry ( [PITH_FULL_IMAGE:figures/full_fig_p014_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: FIG. 4: Error in energy conservation for the non-diagonal-norm SBP wave evolutions. The left and right subfigures [PITH_FULL_IMAGE:figures/full_fig_p014_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: FIG. 5: Error in energy conservation for the staggered SBP wave evolutions. The left and right subfigures show the [PITH_FULL_IMAGE:figures/full_fig_p015_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: FIG. 6: Convergence diagnostics for the 4th-order and 6th-order non-diagonal-norm SBP operators at [PITH_FULL_IMAGE:figures/full_fig_p016_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: FIG. 7: Convergence diagnostics for the 4th-order and 6th-order staggered SBP operators at [PITH_FULL_IMAGE:figures/full_fig_p017_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: FIG. 8: Convergence diagnostics for the 4th-order and 6th-order non-diagonal-norm SBP operator at [PITH_FULL_IMAGE:figures/full_fig_p018_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: FIG. 9: Convergence diagnostics for the 4th-order and 6th-order staggered SBP operator at [PITH_FULL_IMAGE:figures/full_fig_p019_9.png] view at source ↗
read the original abstract

Formulating stable numerical methods for hyperbolic systems in curvilinear coordinate with singularities, e.g. spherical coordinates, is complicated by the presence of these singularities. We present a method for constructing high-order accurate, energy-stable finite difference operators satisfying the Summation-by-Parts (SBP) property on spherical domains, extending ideas presented by [C. Gundlach, J. M. Mart\'in-Garc\'ia, and D. Garfinkle, CQG 30, 145003 (2013)]. We define discrete gradient and divergence operators that mirror the continuous integration-by-parts principle, even though there is a $1/r^p$ coordinate singularity present at the origin. We explicitly construct such operators up to order six. Our operators place a grid point directly on the origin. We also review how to construct stable SBP operators that straddle the origin. We analyze the accuracy and spectral radii of these operators, and we show example evolutions of the scalar wave equation to demonstrate the advantages of such operators.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

0 major / 3 minor

Summary. The manuscript presents a construction of high-order (up to sixth-order) summation-by-parts (SBP) finite-difference operators for discrete gradient and divergence on spherical domains that include the coordinate origin. Extending the 2013 Gundlach et al. framework, the authors define operators that satisfy the discrete integration-by-parts identity exactly despite the 1/r^p singularity by placing a grid point at r=0, supply explicit stencil coefficients, review operators that straddle the origin, analyze truncation errors and spectral radii, and demonstrate stable evolutions of the scalar wave equation.

Significance. If the constructions and verifications hold, the work supplies a practical, explicit route to energy-stable high-order discretizations in curvilinear coordinates without artificial regularization at the origin. The explicit operator coefficients, machine-precision confirmation of the SBP identity, and parameter-free character of the construction (no free parameters or ad-hoc entities introduced) are concrete strengths that facilitate adoption in numerical relativity and related hyperbolic PDE applications.

minor comments (3)
  1. [Abstract] Abstract: the statement that operators are constructed 'up to order six' would be clearer if it indicated which specific orders receive fully explicit coefficient tables.
  2. [§4] §4 (accuracy analysis): the reported convergence rates in the numerical tests are slightly lower than the formal order for some resolutions; a brief remark on whether this is due to the origin treatment or boundary closure would aid interpretation.
  3. [§5] §5 (spectral radii): adding a direct comparison of the spectral radii of the new spherical operators against equivalent Cartesian SBP operators would quantify any additional stiffness introduced by the curvilinear formulation.

Simulated Author's Rebuttal

0 responses · 0 unresolved

We thank the referee for their positive summary of the manuscript, recognition of its strengths (explicit coefficients, machine-precision SBP verification, and parameter-free construction), and recommendation of minor revision. No specific major comments appear in the report, so we have no individual points requiring rebuttal or clarification.

Circularity Check

0 steps flagged

No significant circularity; explicit constructions are self-contained

full rationale

The paper's central contribution is the explicit construction of discrete gradient and divergence operators (up to order 6) on spherical grids with a point at r=0 that satisfy the SBP identity exactly, verified to machine precision, extending the independent 2013 Gundlach et al. framework. No steps reduce by definition, fitted parameters renamed as predictions, or load-bearing self-citations; the operators, modified weights, and wave-equation tests are presented as direct, falsifiable outputs independent of the claims.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The work extends existing summation-by-parts theory without introducing new free parameters, axioms beyond standard finite-difference consistency, or invented entities.

axioms (1)
  • domain assumption Discrete operators can be defined to satisfy a summation-by-parts identity that mirrors continuous integration-by-parts on a domain containing a coordinate singularity.
    Invoked when the authors state they define gradient and divergence operators that mirror the continuous principle despite the 1/r^p singularity.

pith-pipeline@v0.9.1-grok · 5719 in / 1193 out tokens · 27203 ms · 2026-06-28T05:03:14.709284+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

42 extracted references · 3 canonical work pages

  1. [1]

    4th order accurate non-diagonal norm operators For a fourth order accurate vector norm V we use the ansatz V = 1 0 0 0 0 0 · · · 0 0 v2 v2,3 0 0 0 · · · 0 0 v2,3 v3 v3,4 0 0 · · · 0 0 0 v3,4 v4 v4,5 0 · · · 0 0 0 0 v4,5 v5 0 · · · 0 0 0 0 0 0 s6 · · · 0 ... ... ... ... ... ... ... ... 0 0 0 0 0 0 · · · sn     (34) We ad...

  2. [2]

    We cannot demand this near the outer boundary because the underlying Cartesian SBP gradient loses accuracy near the outer boundary

    Everywhere, except close to the outer boundary, we demand that the divergence is exact when applied to the linear functionf(r) =r, i.e.Dr=p+1. We cannot demand this near the outer boundary because the underlying Cartesian SBP gradient loses accuracy near the outer boundary

  3. [3]

    We do this to increase the number of conditions to match the number of unknowns (see below)

    Near the origin we also demand that the divergence is exact forf(r) =r 3. We do this to increase the number of conditions to match the number of unknowns (see below)

  4. [4]

    that R f(r)r p dris exact forf(r) = 1

    Finally we demand that the volume of the domain is exact, i.e. that R f(r)r p dris exact forf(r) = 1. We do so because the first and third conditions are special cases, and it seems advantageous to handle them exactly. We find a sufficient number of off-diagonal non-zero elements by trial and error. The ansatz here leads to a fourth- order operator with g...

  5. [5]

    multiply by (p/r)

    6th order accurate non-diagonal norm operators To construct a sixth-order operator, we start with a similar ansatz for the vector norm V, with more non-zero off-diagonal elements: V =   1 0 0 0 0 0 0 0 0 0 0· · ·0 0v 2,2 v2,3 v2,4 0 0 0 0 0 0 0· · ·0 0v 2,3 v3,3 v3,4 v3,5 0 0 0 0 0 0· · ·0 0v 2,4 v3,4 v4,4 v4,5 0 0 0 0 0 0· · ·0 0 0v...

    2000

  6. [6]

    Gundlach, J

    C. Gundlach, J. M. Mart´ ın-Garc´ ıa, and D. Garfinkle, Summation by parts methods for spherical harmonic decompositions of the wave equation in any dimensions, Classical and Quantum Gravity30, 145003 (2013)

    2013

  7. [7]

    H. O. Kreiss and G. Scherer,Mathematical Aspects of Finite Elements in Partial Differential Equations, edited by C. D. Boor(Academica Press, 1974)

    1974

  8. [8]

    H. O. Kreiss and G. Scherer,On the Existence of Energy Estimates for Difference Approximations for Hyperbolic Systems, Tech. Rep. (Department of Scientific Computing, Uppsala University, 1977)

    1977

  9. [9]

    M. W. Choptuik, Universality and scaling in gravitational collapse of a massless scalar field, Phys. Rev. Lett.70, 9 (1993)

    1993

  10. [10]

    Gundlach and J

    C. Gundlach and J. M. Mart´ ın-Garc´ ıa, Critical phenomena in gravitational collapse, Living Reviews in Relativity10, 5 (2007)

    2007

  11. [11]

    Gundlach, D

    C. Gundlach, D. Hilditch, and J. M. Mart´ ın-Garc´ ıa, Critical phenomena in gravitational collapse, arXiv preprint arXiv:2507.07636 (2025), arXiv:2507.07636 [gr-qc]

    arXiv 2025

  12. [12]

    M¨ uller, Hydrodynamics of core-collapse supernovae and their progenitors, Living Reviews in Computational Astrophysics 6, 3 (2020)

    B. M¨ uller, Hydrodynamics of core-collapse supernovae and their progenitors, Living Reviews in Computational Astrophysics 6, 3 (2020)

    2020

  13. [13]

    S. L. Liebling and C. Palenzuela, Dynamical boson stars, Living Reviews in Relativity20, 5 (2017)

    2017

  14. [14]

    Lehner and F

    L. Lehner and F. Pretorius, Black strings, low viscosity fluids, and violation of cosmic censorship, Phys. Rev. Lett.105, 101102 (2010)

    2010

  15. [15]

    Prochnow, O

    B. Prochnow, O. O’Reilly, E. M. Dunham, and N. A. Petersson, Treatment of the polar coordinate singularity in axisym- metric wave propagation using high-order summation-by-parts operators on a staggered grid, Computers & Fluids149, 138 (2017)

    2017

  16. [16]

    O’Reilly and N

    O. O’Reilly and N. A. Petersson, Energy conservative sbp discretizations of the acoustic wave equation in covariant form on staggered curvilinear grids, Journal of Computational Physics411, 109386 (2020)

    2020

  17. [17]

    Srinivasan, M

    A. Srinivasan, M. Dumett, C. Paolini, G. F. Miranda, and J. E. Castillo, Mimetic finite difference operators and higher order quadratures, GEM - International Journal on Geomathematics14, 19 (2023)

    2023

  18. [18]

    Mohseni and T

    K. Mohseni and T. Colonius, Numerical treatment of polar coordinate singularities, Journal of Computational Physics 157, 787 (2000)

    2000

  19. [19]

    M. D. Griffin, E. Jones, and J. D. Anderson, A computational fluid dynamic technique valid at the centerline for non- axisymmetric problems in cylindrical coordinates, Journal of Computational Physics30, 352 (1979)

    1979

  20. [20]

    Constantinescu and S

    G. Constantinescu and S. Lele, A highly accurate technique for the treatment of flow equations at the polar axis in cylindrical coordinates using series expansions, Journal of Computational Physics183, 165 (2002)

    2002

  21. [21]

    Bhole, B

    A. Bhole, B. Nkonga, S. Pamela, G. Huijsmans, and M. Hoelzl, Treatment of polar grid singularities in the bi-cubic hermite-b´ ezier approximations: Isoparametric finite element framework, Journal of Computational Physics471, 111611 (2022)

    2022

  22. [22]

    K. O. Friedrichs, Symmetric hyperbolic linear differential equations, Communications on Pure and Applied Mathematics 7, 345 (1954), https://onlinelibrary.wiley.com/doi/pdf/10.1002/cpa.3160070206

    work page doi:10.1002/cpa.3160070206 1954

  23. [23]

    Friedrichs and H

    K. Friedrichs and H. Lewy, ¨Uber die eindeutigkeit und das abh¨ angigkeitsgebiet der l¨ osungen beim anfangswertproblem linearer hyperbolischer differentialgleichungen, Mathematische Annalen98, 192 (1928)

    1928

  24. [24]

    K. O. Friedrichs, Symmetric positive linear differential equations, Communications on Pure and Applied Mathematics11, 333 (1958), https://onlinelibrary.wiley.com/doi/pdf/10.1002/cpa.3160110306

    work page doi:10.1002/cpa.3160110306 1958

  25. [25]

    Gustafsson, H

    B. Gustafsson, H. Kreiss, and J. Oliger,Time-Dependent Problems and Difference Methods(Wiley, 2013)

    2013

  26. [26]

    Strand, Summation by parts for finite difference approximations for d/dx, Journal of Computational Physics110, 47 (1994)

    B. Strand, Summation by parts for finite difference approximations for d/dx, Journal of Computational Physics110, 47 (1994)

    1994

  27. [27]

    D. C. Del Rey Fern´ andez, P. D. Boom, and D. W. Zingg, A generalized framework for nodal first derivative summation- by-parts operators, Journal of Computational Physics266, 214 (2014). 21

    2014

  28. [28]

    D. C. Del Rey Fern´ andez and D. W. Zingg, Generalized summation-by-parts operators for the second derivative, SIAM Journal on Scientific Computing37, A2840 (2015), https://doi.org/10.1137/140992205

    work page doi:10.1137/140992205 2015

  29. [29]

    Mattsson and J

    K. Mattsson and J. Nordstr¨ om, Summation by parts operators for finite difference approximations of second derivatives, Journal of Computational Physics199, 503 (2004)

    2004

  30. [30]

    Diener, E

    P. Diener, E. N. Dorband, E. Schnetter, and M. Tiglio, Optimized high-order derivative and dissipation operators satisfying summation by parts, and applications in three-dimensional multi-block evolutions, J. Sci. Comput.32, 109 (2007)

    2007

  31. [31]

    Mattsson, M

    K. Mattsson, M. Sv¨ ard, and J. Nordstr¨ om, Stable and accurate artificial dissipation, Journal of Scientific Computing21, 57 (2004)

    2004

  32. [32]

    Mattsson, Summation by parts operators for finite difference approximations of second-derivatives with variable coeffi- cients, J

    K. Mattsson, Summation by parts operators for finite difference approximations of second-derivatives with variable coeffi- cients, J. Sci. Comput.51, 650 (2012)

    2012

  33. [33]

    Mattsson, M

    K. Mattsson, M. Sv¨ ard, and M. Shoeybi, Stable and accurate schemes for the compressible navier–stokes equations, Journal of Computational Physics227, 2293 (2008)

    2008

  34. [34]

    Hairer and G

    E. Hairer and G. Wanner,Solving Ordinary Differential Equations II(Springer Berlin Heidelberg, 1996)

    1996

  35. [35]

    Shampine, Solving odes and ddes with residual control, Applied Numerical Mathematics52, 113 (2005)

    L. Shampine, Solving odes and ddes with residual control, Applied Numerical Mathematics52, 113 (2005)

    2005

  36. [36]

    Hairer, S

    E. Hairer, S. P. Nørsett, and G. Wanner,Solving Ordinary Differential Equations I: Nonstiff Problems, 2nd ed., Springer Series in Computational Mathematics, Vol. 8 (Springer Berlin, Heidelberg, 1993)

    1993

  37. [37]

    Rackauckas and Q

    C. Rackauckas and Q. Nie, DifferentialEquations.jl–a performant and feature-rich ecosystem for solving differential equa- tions in Julia, Journal of Open Research Software5(2017)

    2017

  38. [38]

    Vretinaris, SphericalSBPOperators.jl (2026)

    S. Vretinaris, SphericalSBPOperators.jl (2026)

    2026

  39. [39]

    Danisch and J

    S. Danisch and J. Krumbiegel, Makie.jl: Flexible high-performance data visualization for Julia, Journal of Open Source Software6, 3349 (2021)

    2021

  40. [40]

    D. K. Zhang and A. Aiken, High-performance branch-free algorithms for extended-precision floating-point arithmetic, in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’25 (Association for Computing Machinery, New York, NY, USA, 2025) p. 695–710

    2025

  41. [41]

    R. A. Smith, GenericSchur.jl: Julia package for Schur decomposition of matrices with generic element types (2025)

    2025

  42. [42]

    Ranocha, SummationByPartsOperators.jl: A Julia library of provably stable semidiscretization techniques with mimetic properties, Journal of Open Source Software6, 3454 (2021)

    H. Ranocha, SummationByPartsOperators.jl: A Julia library of provably stable semidiscretization techniques with mimetic properties, Journal of Open Source Software6, 3454 (2021)

    2021