Spectrally-tuned compact finite-difference schemes with domain decomposition and applications to numerical relativity
Boris Daszuta
TL;DR
This paper introduces spectrally-tuned compact finite-difference schemes that improve accuracy and efficiency in wave propagation and numerical relativity applications, especially under domain decomposition.
Contribution
It presents a novel method for numerically tuning spectral error in compact FD schemes, enhancing accuracy and robustness in wave and relativity simulations.
Findings
Error reduced by at least an order of magnitude over standard FD
Hybrid strategy maintains accuracy under domain decomposition
Improved phase accuracy in numerical relativity simulations
Abstract
Compact finite-difference (FD) schemes specify derivative approximations implicitly, thus to achieve parallelism with domain-decomposition suitable partitioning of linear systems is required. Consistent order of accuracy, dispersion, and dissipation is crucial to maintain in wave propagation problems such that deformation of the associated spectra of the discretized problems is not too severe. In this work we consider numerically tuning spectral error, at fixed formal order of accuracy to automatically devise new compact FD schemes. Grid convergence tests indicate error reduction of at least an order of magnitude over standard FD. A proposed hybrid matching-communication strategy maintains the aforementioned properties under domain-decomposition. Under evolution of linear wave-propagation problems utilizing exponential integration or explicit Runge-Kutta methods improvement is found to…
| Scheme | ||||
|---|---|---|---|---|
| Scheme | Order | ||||
|---|---|---|---|---|---|
| 4 | |||||
| 6 | |||||
| 8 | |||||
| 10 |
| Scheme | Order | |||
|---|---|---|---|---|
| 3 | ||||
| 5 | ||||
| 7 | ||||
| 9 |
| Scheme | ||
|---|---|---|
| Scheme | ||
|---|---|---|
| Scheme | ||
|---|---|---|
| Scheme | ||
|---|---|---|
| Scheme | |||
|---|---|---|---|
| Scheme | |||
|---|---|---|---|
| 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | |
| 1 | 1 | 1 | 1 | 1 | 1 | |
| 1 | 1 | 1 | 1 | 1 | 1 | |
| 1 | 1 | 1 | 1 | 1 | 1 | |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsNumerical methods for differential equations · Advanced Numerical Methods in Computational Mathematics · Electromagnetic Simulation and Numerical Methods
Spectrally-tuned compact finite-difference schemes with domain decomposition and applications to numerical relativity
Boris Daszuta1
1Theoretisch-Physikalisches Institut, Friedrich-Schiller-Universität Jena, 07743, Jena, Germany
Abstract
Compact finite-difference (FD) schemes specify derivative approximations implicitly, thus to achieve parallelism with domain-decomposition suitable partitioning of linear systems is required. Consistent order of accuracy, dispersion, and dissipation is crucial to maintain in wave propagation problems such that deformation of the associated spectra of the discretized problems is not too severe. In this work we consider numerically tuning spectral error, at fixed formal order of accuracy to automatically devise new compact FD schemes. Grid convergence tests indicate error reduction of at least an order of magnitude over standard FD. A proposed hybrid matching-communication strategy maintains the aforementioned properties under domain-decomposition. Under evolution of linear wave-propagation problems utilizing exponential integration or explicit Runge-Kutta methods improvement is found to remain robust. A first demonstration that compact FD methods may be applied to the c formulation of numerical relativity is provided where we couple our header-only, templated C++ implementation to the highly performant GR-Athena++ code. Evolving c on test-bed problems shows at least an order in magnitude reduction in phase error compared to FD for propagated metric components. Stable binary-black-hole evolution utilizing compact FD together with improved convergence is also demonstrated.
pacs:
04.25.D-, 04.30.Db, 95.30.Sf,
I Introduction
Finite difference (FD) methods provide a well-known and flexible approach on structured grids furnishing derivative approximants at target grid-points. Allowing for an implicit, linear relation between known function samples and sought derivative values leads to compact finite difference (CFD) schemes Lele (1992). At fixed formal order of accuracy and total number of coupled points, CFD allows for construction of narrower stencils Hirsh (1975); Carpenter et al. (1994) when compared with explicit FD together with improved resolution characteristics over a wider range of (spatial) scales Lele (1992); Mehra and Patel (2017). Numerical solution of the associated linear systems (typically banded tri- or penta-diagonal Hirsh (1975); Lele (1992); Fu and Ma (1997); Wang et al. (2013); Qin et al. (2014); Chen et al. (2021)) is more expensive than the direct evaluation of explicit FD. Solution of such banded systems is of linear algorithmic complexity Higham (2002); Golub and Van Loan (2013); Askar and Karawia (2015) and this computational overheard is mitigated through improved resolving efficiency. Indeed this efficiency at resolving widely disparate length-scales has led CFD-based techniques to be extensively widely in areas such as e.g. computational aeroacoustics Kim (2007), direct numerical simulations of Navier-stokes without turbulence modeling Jagannathan and Donzis (2016), and large eddy simulations Bodony and Lele (2005); see also Song et al. (2022).
Given problem-specific requirements it may be preferable to quantify and minimize error at specific spatial scales or alternatively over a desired wavenumber range when viewed in the frequency domain. Indeed increasing resolution capability of CFD further rather than trunction order can enchance overall accuracy more effectively Kim and Lee (1996); Kim (2007). One approach in derivation of CFD schemes is to impose a specific (implicit) stencil arrangement and by Taylor matching derive coefficients that yield a scheme of specific formal order of accuracy Hirsh (1975); Lele (1992). Unfortunately this does not give direct access to control over desired spectral behaviour at particular wavenumbers. Instead rather than matching all possible coefficients a degree of freedom may be left underdetermined. This can be exploited through minimization against a suitably chosen cost function characterizing an error profile which allows for spectral-tuning Kim and Lee (1996); Liu et al. (2008); Zhang and Yao (2013). A unified approach to this for construction of CFD schemes with arbitrary derivative degree, order, and stencil size (and bias) was recently introduced in Deshpande et al. (2019). The advantage of this latter approach in addition to the aforementioned generalizations is that various properties such as coefficient (skew-)symmetries do not need to be imposed a priori but are a consequence of the optimization procedure. While quite general, the framework of Deshpande et al. (2019) does not immediately embed derivation of multi-derivative such as those of Fu and Ma (1997); Qin et al. (2014) or Hermite-FD methods Fornberg (2020).
In realistic problems that significantly deviate from non-ideal conditions parallelizing simulations can be extremely demanding both from the point of view of problem reformulation aspects together with implementation such that high performance computing (HPC) infrastructure may be efficiently utilized Huerta et al. (2019). Implicit CFD schemes introduce data-dependency which must be treated when one seeks to introduce parallelism. Two broad categories of approaches can be envisaged Kim (2013); Chen et al. (2021): the algorithmic approach and the boundary approximation approach. In the algorithmic approach the implicit system is solved in parallel utilizing techniques such as the pipeline Thomas algorithm Povitsky and Morris (2000) or parallel diagonal dominance algorithm Sun (1995); Terekhov (2016). For the boundary approximation approach (BAA) a computational domain is partitioned and closures are imposed on the CFD scheme such that inversion of linear systems may be performed in a decoupled fashion on each sub-domain. Unfortunately this can lead to changes in the resolving efficiency in the vicinity of boundaries and potentially introduce artifacts. Examples of treating partitioned sub-domains include partial overlapping with implicit closures supplemented by one-sided filtering as in Sengupta et al. (2007) or sub-domain ghost layer extension with closures prescribed based on a combination of optimization Kim (2013). It is known that regardless of the style of approach when fluid flow or wave-propagation problems are treated numerical dispersion relation preservation (DRP) is crucial Tam and Webb (1993). It was recently demonstrated in Chen et al. (2021) that for first degree derivative, upwind CFD schemes DRP can be achieved exactly. This was illustrated in a order scheme that showed excellent consistency during numerically evolved flow problems under domain partitioning preserving accuracy particularly as sub-domain sampling was increased.
As in the hydrodynamical context an important concern for numerical relativity (NR) investigations of the binary black hole (BBH) merger problem is careful treatment of the widely-varying range of length and time-scales involved Ashtekar et al. (2015). The numerical solution of the BBH merger problem crucially complements observational gravitational wave (GW) based detection efforts Abbott et al. (2016a, b, 2017a). A variety of successful NR based approaches for the BBH merger problem are known. Time discretization is generally achieved through a method of lines prescription whereas the most common spatial treatment has utilized Cartesian grid coordinatization and been based on a combination of FD and adaptive mesh refinement (AMR) in Pollney et al. (2011); Reisswig et al. (2013); Brown et al. (2009); Sperhake (2007); Zlochower et al. (2005); Herrmann et al. (2007) built upon Goodale et al. (2003) – see also the independent codes of Brügmann et al. (2008); Cao et al. (2008); Clough et al. (2015). Pseudo-spectral methods involving multi-patch decomposition of the computational domain through a combination of topological spheres and cylinders have also found excellent success Szilagyi et al. (2009) and the related discontinuous Galerkin based approaches also show future promise Hilditch et al. (2016); Bugner et al. (2016); Kidder et al. (2017). In light of recent HPC trends towards massively parallel systems a particular emphasis has been placed within new codes on domain-decomposition strategies compatible with excellent scaling properties. In particular domain-decomposition described through oct-tree based grids has recently been explored in the FD-based codes of Fernando et al. (2018); Daszuta et al. (2021). Curiously for NR while AMR based FD approaches, together with (pseudo)-spectral and discontinuous Galerkin methods have been pursued it does not appear that CFD methods have been previously utilized for the evolution problem. On the GW observational front the continuing effort towards improving operating sensitivity of current detectors Abbott et al. (2020) and developing new detectors Akutsu et al. (2020); Amaro-Seoane et al. (2017); Punturo et al. (2010); Abbott et al. (2017b) allows probing of ever more extreme regions of the underlying parameter space and consequently motivates a search for alternative numerical methods and hence we seek to take a first step towards bridging this gap here.
The first goal of this work is to obviate the need for repeated, hand construction of CFD schemes by automatizing the procedure of Deshpande et al. (2019) numerically and extending it to multi-derivative schemes thereby allowing for rapid experimentation on model problems prior to general application. A second goal is to combine this framework with the BAA technique proposed in Chen et al. (2021), extend it to centered schemes and consider how it may be further refined through an iterative procedure. A third goal is to demonstrate that the CFD methods ensuing from this lead to improved accuracy in solution of wave-propagation problems while remaining robust in the parallel context. Finally we seek to demonstrate that CFD shows future promise in application to the evolution problem of numerical relativity in simulation of binary mergers. The rest of this paper is organized as follows: In §II we describe the overall method for construction of general FD schemes, providing a characterization of spectral error and description of how it may be numerically tuned subject to problem requirements. We describe our hybrid-communication strategy that builds upon DRP and can be leveraged during problems involving domain-decomposition. In §III we consider application to a variety of wave-propagation problems whereupon the c system as implemented in Daszuta et al. (2021) is investigated in §III.3. Section IV concludes. For convenience we have also collected a selection of stencils investigated in the text in §A.
II Method
The aim of this section is to provide a method to optimize finite-difference based numerical evaluation of an unknown derivative of a function based on known (derivative) data. Our approach builds upon the framework of Deshpande et al. (2019) and extends it to allow incorporation of data from multiple derivatives in addition to function data in specification of stencils.
Consider a discretization of the interval of uniform spacing where sampled points are denoted . Suppose then function (derivative) samples are denoted . A generalized finite-difference stencil may be viewed as a relation between linear combinations of unknown and specified samples. We therefore introduce the stencil coefficient weights111Typically we suppress indicial ranges when evident from context.222In this work we assume that coefficients take values in .:
[TABLE]
together with function samples:
[TABLE]
The product is interpreted through summation . As ansatz for the -derivative of at the base-point denoted we form:
[TABLE]
where is an indexing set with number of elements denoted and is ordered such that . We assume that . Notice that we have not restricted the number of elements in individual coefficient weights which allows for consideration of biased schemes. Our setup will allow us to form explicit or implicit stencils. Prescriptions that involve function data and also derivatives as known samples such as the Hermite methods of Fornberg (2020) together with their implicit extensions may also be constructed. As a schematic of multi-derivative, implicit stencil node coupling see Fig.1.
In order to determine the -coefficients of a particular scheme we need to impose constraints. Suppose that and are known then Eq. (3) may be viewed as a linear system of equations for unknowns specifying the stencil. To achieve this notice that in Eq.(3) we also have a freedom to fix an overall scaling and so we impose . Next note that for a stencil of formal order we require the derivative approximant to have truncation error . Therefore we consider Taylor series expansion of about to order :
[TABLE]
Substitution of Eq.(4) into Eq.(3) and matching coefficients order by order leads to additional linear constraints on the -coefficients. If it is the case that is selected such that then we have a fully determined linear system that may be solved uniquely resulting in a scheme for numerical calculation of a derivative approximant of fixed formal order. For a general, implicit, multi-derivative scheme this can be a potentially tedious procedure and consequently we make use of a computer algebra system (CAS) to automate construction of these constraints not . As a example consider Fig.1 for , , and . In this case we have , , and can find linear conditions:
[TABLE]
These together with the normalization condition immediately yields:
[TABLE]
which can be verified to be of formal order for computing implicitly. Such order matching strategies are often employed in the construction of the classical Padé schemes Lele (1992).
Of more interest is the case of which yields an underdetermined system. This property can be exploited. We can consider further introduction and constrained minimization against a well-defined objective function (e.g. characterizing some error metric of interest) so as to yield at some desired fixed formal order tuned against Deshpande et al. (2019). One can also impose constraints that fix the values of constants that appear in the truncation error directly which is of use during domain-decomposition (§II.4 and §II.5).
II.1 Characterizing and tuning derivative spectral error
In addition to imposing control on formal order of accuracy in this work we are also concerned with reducing spectral error which can be characterized based on Fourier analysis Lele (1992); Deshpande et al. (2019); LeVeque (2007). This will furnish us with a functional that can be directly optimized numerically for the stencil coefficients subject to the linear constraints previously discussed.
Consider where the (periodic) interval serves as a model for . Recall that family of plane-waves constitute an orthonormal basis for Hesthaven et al. (2007). Suppose that is smooth and consider the truncated expansion:
[TABLE]
Derivative approximants may be constructed utilizing Eq.(7) i.e. by evaluating
[TABLE]
where for later convenience we have defined . An important property is that for and its truncated expansion we have Hesthaven et al. (2007). Suppose now that is discretized through introduction of a sequence of samples of uniform spacing . The discrete analogue of Eq.(7):
[TABLE]
shares similar approximation properties to the continuous expansion for smooth functions provided that is selected large enough such that aliasing error is suppressed Hesthaven et al. (2007); Boyd (2001). With discrete orthogonality focus may be restricted to a single mode where :
[TABLE]
where we have defined the normalized wave number . Set:
[TABLE]
The selected form of together with Eq.(10), Eq.(11) and Eq.(2) allows for Eq.(3) to be rewritten as:
[TABLE]
Thus Eq.(12) provides for a characterization of the spectral error associated with a finite-difference scheme as we may compare the modified, normalized wave number to the analytically expected . To this end we introduce:
[TABLE]
Which allows for consideration of an optimization problem based on the functional:
[TABLE]
where is a non-negative weight function that allows to preferentially tune over prescribed ranges of . The choices of Eq.(13) and Eq.(14) are motivated by the requirement of a convex optimization problem (see discussion in Deshpande et al. (2019) for and prescribed values ). In order to construct when we solve (numerically) subject to the normalization condition and linear constraints arising from Taylor series matching using IPOPT Wächter and Biegler (2006). To achieve this we refactor the functional such that quadrature may be performed numerically so as to yield an objective function involving only the sought after . In brief, some choice of , , and is made. Then we set:
[TABLE]
and take . The enlarged stencil coefficients we denote . Next we construct with respect to in factored form:
[TABLE]
where is a row-vector assembled from the unknown of size and is a square matrix formed from the enlarged to be integrated element-wise. We use CAS to automate this factorization process. To seek a solution based on the initial target stencil sizes further auxiliary linear constraints are imposed:
[TABLE]
Working as above allows for rapid derivation and experimentation with a variety of schemes – we have prepared a public notebook not that treats construction of the Taylor constraints, associated error functional construction and factorization, numerical quadrature evaluation, and solution for the resulting scheme coefficients at arbitrary precision.
II.2 Dispersion and dissipation in wave propagation
Our intention is to apply tuned, derivative approximant schemes to wave propagation problems as described by discretized hyperbolic partial differential equations and consequently we briefly recall dispersion and dissipation properties. This will allow for providing an assessment on potential phase accuracy achieved. For in-depth treatments of this together with stability and convergence properties see LeVeque (2007); Gustafsson et al. (2013). For simplicity we focus on the linear advection equation as it is sufficient to fix conventions and illustrate concepts required later. To this end define with as in §II.1 and consider smooth satisfying:
[TABLE]
with or respectively representing a constant speed right-ward or left-ward propagation of the initial condition (IC) . Due to linearity it is sufficient to focus on a single Fourier mode such as which in turn leads to the time-dependent, translated profile . In order to pass to the semi-discrete setting consider the equidistant discretization where . From Eq.(18) the usual method of lines description now follows:
[TABLE]
where the components of represent values of a discrete derivative stencil in matricial form based on Eq.(3) that embeds the periodic boundary conditions. Formal solution of Eq.(19) can be immediately provided through direct matrix exponentiation . Whether the amplitude of a mode remains bounded under this prescription as is expected analytically depends sensitively on the spectrum of . To see this directly one can work in the modal representation which for an initial condition comprised of a fixed Fourier mode becomes:
[TABLE]
with the discrete derivative appearing in Eq.(19) now described through the associated modified wavenumber of the selected scheme Deshpande et al. (2019). It follows that where . Comparing to the analytical we see that modifies the speed of propagation i.e. the dispersion whereas non-zero introduces amplitude attentuation (dissipation) or amplification depending on sign.
Following Hesthaven et al. (2007) a characterization of the number of points required to achieve an error tolerance for the phase associated with a given Fourier mode in solution of Eq.(18) may now be provided. If it is the case that for a selected spatial derivative approximant then there is no difference in the amplitude of the two solutions and at the semi-discrete level. Indeed we take the phase error as the leading contribution in the relative error:
[TABLE]
If we introduce the relative error of the modified, normalized wavenumber as:
[TABLE]
and we define the number of periods (in time) the solution has propagated as then Eq.(21) becomes:
[TABLE]
Given an IC we have waves in the domain and consequently the number of points per wavelength is . In order to resolve without aliasing we require that . Of interest is the number of samples required per wavelength to attain a specified phase error . For a scheme with formal order of accuracy we may expand and consequently:
[TABLE]
provides an estimate on the number of points required per wavelength required satisfy a phase error tolerance of .
In practice numerical solutions of systems such as Eq.(18) are often constructed based on time-marching schemes involving discretization in time. Explicit Runge-Kutta (ERK) methods are the most widely used example of this. Linear stability and accuracy for an s-stage ERK method is characterized completely by the so-called stability polynomial defined by the method Ketcheson and Ahmadia (2012); Butcher (2008); Hairer and Wanner (2010). Given a constant-coefficient linear system subject to initial conditions, full discretization can be reduced to an iteration involving . We have that where is a time-step, is an approximation to , and is a discretized approximation to Ketcheson and Ahmadia (2012). Based on the iteration it is clear that the propagation of errors is controlled by . For ERK methods defines a closed curve in the complex plane the interior of which forms the so-called absolute stability region . The iteration in time is considered absolutely stable if for all eigenvalues . Tuning to be sufficiently small thus implies remains bounded under iteration and forms a necessary333This is also a sufficient condition for stability if is a normal matrix Hesthaven et al. (2007). For the more general case see e.g. Ketcheson and Ahmadia (2012) and references therein. condition for stable propagation of errors. Typically in numerical solution of a fully-discretized system the Courant-Friedrich-Lewy (CFL) number is taken as one necessary criterion for stability, where depends on the details of and the ERK scheme applied.
II.3 Comparison of select schemes
The purpose of this section is to illustrate spectral properties of derivative approximant schemes that arise upon full Taylor matching in Eq.(3) and furthermore tuning based on exploiting underdeterminedness as described in §II.1. For the latter we make use of a cutoff in Eq.(14):
[TABLE]
which is motivated by seeking to reduce spectral error at low to moderate wavenumber. Scheme stencils we have generated that are later utilized for wave-propagation problems are summarized in §A. A general scheme is here denoted where indices up correspond to the left-hand-side of Eq.(3). Without ambiguity we can condense the notation such that for centered explicit schemes we write . In the case of explicit, one point lop-siding we write . The classical (fully order-matched) Padé Lele (1992) we denote with . When tuning based on Eq.(25) is performed is explicitly indicated. Thus takes , , and in Eq.(25). Hermite schemes involve multiple derivatives which are indicated explicitly.
Using our CAS notebook not we generate a variety of example stencils satisfying a formal order of accuracy of . In order to assess spectral properties we plot the real and imaginary parts of the modified wavenumber as defined by Eq.(12) over in Fig.2. These plots may be extended to through the parity conditions whereas .
As can be seen from Fig.2a the standard explicit scheme has a more accurate dispersion than at low , however the situation reverses as is increased. This is consistent with Chirvasa and Husa (2010) which analyzes the use of lop-sided stencils when treating advective terms. We also observe such an advantage in numerical experiments based on setup of the aforementioned work performed in §III.2. On the other hand it is clear that implicit schemes have better dispersion properties than their explicit counterparts at equivalent formal order of accuracy as can be seen comparing and in the middle panels of Fig.2. Enlarging the number of function samples coupled by a stencil and leveraging underdeterminedness as described in §II.1 allows us to construct which further improves over . In Deshpande et al. (2019) it was observed that increasing the implicit bandwidth further is also advantageous which we also observe based on e.g. . If there is freedom to incorporate function derivative data we find that the implicit Hermite schemes outperform all other choices. Unfortunately a scheme such as requires knowledge of the second derivative of a function to compute the first which may not be available without reformulation of a given problem of interest Fornberg (2020). It has been observed in Lele (1992) that when first degree derivatives are required on a staggered grid a stencil may be devised that does not suffer as . We verify this directly by generating where CC and VC have the meaning of Eq.(26).
To provide another comparison of first degree schemes we consider the semi-discretized form of the one-dimensional advection equation (see §II.2) and the number of points required to resolve a mode at a fixed tolerance on phase error in Tab.1.
II.4 Domain decomposition: dispersion relation preservation
For large-scale problems a high degree of parallelism is required and the implicit nature of compact stencils may seem as a potential obstruction to achieve this. In this work we consider a strategy based on domain-decomposition inspired by Chen et al. (2021). Suppose we have a decomposed and discretized domain . Instead of solving a global, implicit system sequentially where data is communicated between neighboring and in succession a specification of a tuned complementary scheme on is made such that decoupled calculations may be performed in parallel that are local to each .
To illustrate this consider a one-dimensional that has been discretized with uniformly spaced samples. Suppose that is periodically identified. A parameter is selected which controls the number of samples globally. Under domain-decomposition we fix on each such that exactly divides . For we can consider cell-centered (CC) and vertex-centered (VC) sampling that are respectively defined through:
[TABLE]
and the grid spacing we denote by . In the case of the sub-domain shares points with its nearest neighbours . To the exterior of a thin layer of ghost nodes with the same uniform spacing are appended (where ). This facilitates communication of data between sub-domains of the decomposition and evaluation of e.g. centered stencils for all in . We illustrate schematically in Fig.3 the grid structure together with boundary communication strategy.
Consider now evaluation of a centered derivative approximant for a scheme that implicitly couples nearest-neighbor nodes such as . Such a scheme provides the approximation at a point by coupling and relating to a linear combination of function samples . In the case that we treat the domain globally then the derivative approximant under discussion may rewritten as a tridiagonal system which can be directly inverted utilizing the well-known Thomas Algorithm (TDMA) at an algorithmic complexity Higham (2002); Golub and Van Loan (2013). On the other hand when working with individual sub-domains from Fig.3 we see that a closure relation is required at the points with local indices if we are to apply TDMA with data locally available to .
Motivated by solution of wave-like propagation problems it is crucial to match both the formal order of accuracy of derivative approximants together with (approximate) preservation of the modified wavenumber §II.2. Given the relative error of the modified, normalized wavenumber of Eq.(22) recall that and are related to the dispersion and dissipation properties of a scheme respectively. For centered stencils such as the coefficient tuples are (skew)-symmetric about the central element for (odd)-even respectively Deshpande et al. (2019). It therefore follows from Eq.(12) that and consequently . Suppose that is then:
[TABLE]
where only even powers of appear, and the constants are real. On as the nodes with indices are approached from the interior the approximant can eventually no longer be applied if we wish to treat the evaluation locally in a decoupled fashion using a banded linear solver. Instead we construct a (sequence of) tuned, centered scheme(s) with fixed formal order of accuracy where the implicit bandwidth (i.e. ) is reduced as is approached. These closing schemes are instead imposed on nodes where it is not possible to impose . If is reduced but simultaneously may increase then there is sufficient freedom to match the lowest order error constants of the closing schemes to those appearing in Eq.(27). This procedure can be tuned to approximately preserve dispersion. In particular for we can apply on with and at impose explicit, centered closures. For an explicit scheme can be constructed to match the formal order of accuracy. In general the constants appearing in the expansion of for each scheme will not match. Therefore we increase the explicit stencil size and fix the number of equations arising from inserting the Taylor approximant of Eq.(4) into Eq.(3) which fixes the order. This however results in an underdetermined system for the scheme coefficients. Due to underdeterminedness we may impose additional constraints that tune . We have:
[TABLE]
where the error constants are real and all odd powers of must vanish. Therefore to construct we specify Taylor conditions, a unit normalization condition on of , and also impose:
[TABLE]
This gives us linear conditions on the unknown coefficients which may be solved for directly. Hence we find a composite scheme over at consistent formal order of approximation together with consistent lowest order expansion constants. This allows us to treat each independently once the ghost nodes have been populated with known function data. Once this is the case we can solve the banded system:
[TABLE]
where are fixed by the explicit closure whereas the remaining are to be solved for. For the present closure we require that on where in Eq.(30) data at node indices and is supplied through communication from neighbouring sub-domains (see Fig.3). Summarized in Tab.8 and Tab.9 of the appendix §A are coefficients for with matched, explicit closures for and . The numerically optimized, centered schemes of §II.3 may be treated in a similar fashion and are summarized in Tab.12 and Tab.13.
We now turn our attention to degree one, biased schemes as it is possible to close stencils on sub-domain boundaries in a fashion that preserves the dispersion relation exactly Chen et al. (2021). We demonstrate directly that this can be done for higher order derivative approximants by extending the order CCU scheme of Qin et al. (2014). To this end we particularize Eq.(3) as:
[TABLE]
where we seek even-order schemes with which will be denoted . We will also impose Eq.(31) with to construct odd-order schemes with . We refer to schemes with stencils as in Eq.(31) involving the term as left biased and as right biased. Given a left biased scheme of the form of Eq.(31) then the analogous right biased scheme is given by , , and Chen et al. (2021); Deshpande et al. (2019).
In Fu and Ma (1997); Qin et al. (2014) the so-called UCD-style schemes and have been presented however here we use the flexibility of our method in allowing for automatic generation and solution of Taylor constraints to easily extend to higher order approximants. Stencil coefficient tuples together with expansions of the modified wavenumber for a variety of orders are provided in Tab.2 and Tab.3.
An interesting observation made in Chen et al. (2021) that is compatible with the expansions we have provided and verified is that switching between left and right bias does not alter the dispersion for and schemes derived based on Eq.(31) and Taylor matching. In particular we have:
[TABLE]
where the error constants and are real. Due to the bias of these schemes closure is only required on a single point and therefore may be closed at utilizing and similarly for without modifying dispersion. However stability properties for a wave propagation problem depend on the dissipation. To modify this latter without affecting dispersion the closure for is instead based on a modification to where at one imposes a stencil of the form Chen et al. (2021):
[TABLE]
where the values of can be approximated using of Tab.7 and once are determined we have closures for the schemes to be applied at . Incorporating the high-degree derivative as in Eq.(34) means that Eq.(32) is modified to:
[TABLE]
where for the constants are left unmodified whereas now depend on . In order to determine we match the next to leading order error term:
[TABLE]
which gives rise to a linear relation for . The result of this is summarized for the schemes of Tab.2 in Tab.4.
In the case of the closure is based on an analogous modification to where we consider:
[TABLE]
and may again be approximated with of Tab.7. Once is determined we have closures for to be applied at . For the closure specification of Eq.(37) it is the case that Eq.(33) is modified to:
[TABLE]
where for the constants are left unmodified whereas now depend on . In order to determine we match the leading order error term:
[TABLE]
The result of this is summarized for the schemes of Tab.3 in Tab.5.
The schemes and that arise are termed so-called slow and fast schemes due to their respective under- and over-estimation of and its role in wave propagation (see §II.2). We illustrate in Fig.4. Consequently it has been suggested Qin et al. (2014); Chen et al. (2021) to linearly combine schemes with a weighting parameter. To this end we consider:
[TABLE]
where , and analogously for the closure relations. In prescribing Eq.(40) we have formal order of accuracy of at least . The resultant schemes that combine and are generalizations of CCU Qin et al. (2014); Chen et al. (2021).
In many problems of interest may not be available therefore this term must be approximated. One may consider utilizing e.g. unfortunately within the closure approach described this only approximately preserves dispersion under domain decomposition. Such CCU schemes we will denote . Construction of is more involved as the now embedded, implicit specification of requires not only a closure at at for the first derivative approximant but also at for the second. Alternatively one can make use of explicit, derivative approximants such as for (see Tab.6). This results in a deformed for and which exactly preserves dispersion between them. Such CCU schemes we will denote . In both approaches we select with a simple strategy by starting at and reducing until over . In the case of when an embedded, explicit, second degree derivative approximation is made for it turns out that overestimates and so instead we set .
For convenience scheme coefficients and closures are provided in Tab.10 whereas for see Tab.11. The normalized, modified wavenumber for is shown in Fig.5. On a sub-domain in order to compute derivative data for the left biased schemes such as as closed by we solved the banded linear system (cf. Eq.(30)):
[TABLE]
similarly for the right biased cases with coefficients now those of as closed by :
[TABLE]
and in both cases we take the size of the ghost node layer to be .
II.5 Domain decomposition: iteration of closures
Having investigated construction of dispersion relation preservation (DRP) in §II.4 we now turn our attention to a combined, hybrid strategy which features a DRP prescription and further iteratively updates values employed in the implicit system closure through additional communication of points.
Under domain-decomposition application of derivative scheme closures on a sub-domain leads to numerical approximants with error that accumulates towards . While this error is convergent at the constructed formal order, when compared to a unpartitioned, global scheme, an artifact of the decomposition can often be observed. This motivates us to consider construction of a sequence of approximants that iteratively refine as found using the closure prescription of §II.4 and consequently better capture the global scheme.
To this end, consider again with the matched closure . We solve Eq.(41) which provides us with . For we then solve:
[TABLE]
where during discretization on the number of ghosts is selected such that can be imposed and we can select such that . In particular allows us to take . In practice the relation of Eq.(43) may be iterated until a desired tolerance is attained. While this requires additional communication only a single value needs updating prior to solution of Eq.(43) as the sampled function remains unchanged. Furthermore the linear system when considered as a matrix equation appears with diagonal and single sub-diagonal and consequently can be cheaply solved for at algorithmic complexity.
In a similar manner we may treat centered implicit schemes such as with matched closures which fashion us with based on Eq.(30). For one then solves:
[TABLE]
where for each additional iteration derivative values at must now be communicated. For we utilize and impose Eq.(44) with and .
We now turn to tests where prescribed functions are numerically differentiated globally and compared against the decomposed, local method. Subsequently we verify through grid convergence tests that compare against numerical sampling of the analytically expected result that error converges away at the design order of the underlying schemes. This will allow us to assess the performance of the explicit closures with (approximate) dispersion relation preservation under domain decompopsition introduced in §II.4 and the iterated closures described by Eq.(43) and Eq.(44). To this end, define the smooth functions:
[TABLE]
As a first direct numerical probe of the effect of decomposition we fix the number of samples on as and partition to sub-domains with samples each. To inspect the effect on the spectrum directly upon numerical differentiation with our proposed compact finite-difference schemes we consider Fourier decomposition of the result and quantify error in the amplitude and phase respectively according to:
[TABLE]
where is a vector of coefficients (appearing in Eq.(9)) as found from the global scheme, are the coefficients of the analogous domain-decomposed approach, and we select coefficients according to a relative amplitude threshold444In the modal representation smooth functions as discussed here have coefficients with amplitudes that decay exponentially Boyd (2001); Hesthaven et al. (2007) which for a sufficiently high band-limit yields negligible contribution to description of the underlying function. Thus we impose a cut when computing to avoid polluting the error analysis with spurious values arising from coefficients with exponentially small amplitude. taken with respect to the global scheme. We select of Eq.(45) with parameters , , , and such that the amplitude of the normalized, modal coefficients of the Fourier decomposition (Eq.(9)) are at numerical round-off i.e. for . In Fig.6 we investigate first degree numerical differentiation of with the (iterated) schemes: and with iteration on , whereas our CCU(6,7) generalization and are iterated with and respectively. It is clear that the implicit, interior compact finite-difference and the complimentary, explicit, tuned closures lead to error that has comparable order of magnitude (see Fig.6a). The decomposition induced error (near the boundaries of the sub-domains) is mitigated through the iterative prescription of this section. This is particularly effective for the biased schemes as can be seen in Fig.6b.
In a similar fashion we take of Eq.(46) and select parameters , , , , , , and such that for . We investigate second degree numerical differentiation of with the (iterated) schemes: and with iteration on , whereas is iterated with . Results are depicted in Fig.7.
The number of samples (i.e. size) of a sub-domain has an influence on the accuracy when compared with the equivalent global scheme. At fixed sampling for where taking more samples per that is cf. is favoured which is particularly apparent in the degree one case. This is compatible with the expectation that introducing fewer sub-domain boundaries leads to fewer artifacts due to them; however, it may not always be practical to select larger due to e.g. concerns involving partitioning a problem to utilize parallelism efficiently. During tests we have found that CC and VC sampling (Eq.(26)) performs comparably. Due to the higher rate at with which biased schemes converge under iteration – for the schemes we have investigated a single iteration is required to reach saturation – a symmetrization procedure may be considered. To this end we compute derivatives utilizing and and then take as the final result after iteration as the average of the two. Such symmetrized schemes we denote as .
In Fig.8 we verify that the error of our proposed schemes converges as resolution is increased at the formal order anticipated from design choices. We find that symmetrization of the biased scheme leads to a reduction in error at intermediate and as is increased. To demonstrate the generality of our approach we also generate a variety of schemes tailored to the full domain ; it can be seen that for higher implicit bandwidth error is more efficiently reduced at equal order, which is in agreement with general conclusions of Deshpande et al. (2019). For reference with more well-known schemes we also compare explicit finite differencing together with a Fourier based spectral method.
II.6 Embedded operations
As discussed in §II.2 given a discretized wave propagation problem the spectrum of the associated discretized operator governing the time-evolution controls stability properties. With the introduction of domain decomposition and tuned closures for derivative schemes the aforementioned spectrum will be modified cf. the analogous global scheme. It is thus useful to have a practical method for inspecting the spectrum. Consequently we summarize a simple strategy that embeds into matrices a representation of sub-domain communication and imposition of boundary conditions together with the iterated closures discussed in §II.5.
Consider a partitioned, periodic, one-dimensional domain where we select as the sampling parameter over and for where divides and the number of sub-domains is . Thus . Suppose that is sampled on each using of Eq.(26) and we extend each by ghost nodes in each direction. Collectively this data may be concatenated into the partitioned vector:
[TABLE]
where each sub-domain (fixed ) has associated to it the elements . To capture communication operations between sub-domains together with imposition of boundary conditions we consider the introduction of a block-partitioned matrix with blocks where each block is of shape . The diagonal blocks are comprised of unit entries555In the case of VC the approach is similar but for consistency shared nodes are averaged Daszuta et al. (2021). which in indices local to the block satisfy666The standard definition of the Kronecker delta if and [math] otherwise is extended to logical conditionals where is taken if satisfied. . To populate the ghost layer of we utilize data data from where . Thus has entries with local indices whereas has entries . As an example suppose the overall structure of becomes:
[TABLE]
where , , and are , , and zero matrices, and is an identity matrix. In practice, during calculation of an implicit, centered numerical derivative on a sub-domain we form a relation as in Eq.(30) which would yield on nodes through the use of a banded solver. Alternatively we may form the inverse matrix to solve Eq.(30) directly. In order to work with square matrices that have shapes compatible with the blocks appearing in and we define:
[TABLE]
where and have shapes , and zero entries we now suppress. In Eq.(57) and Eq.(65) elements are partitioned according to whether they act on data of or the ghost layer. On each we may thus evaluate derivatives using . For collective application to each sub-domain described by we need copies of embedded as: . This describes an initial communication populating ghost layer data, followed by differentiation, and finalized by another communication. To provide with we need to incorporate the iterated scheme of Eq.(44). To this end define by replacing with in of Eq.(57) whereas of Eq.(65) becomes:
[TABLE]
By introducing and matrix multiplication will allow for extraction and fixing of values at nodes of . We introduce the block-diagonal matrices , , , and . Finally this allows for:
[TABLE]
Given and domain-decomposed function data assembled as in Eq.(48) we may perform numerical experiments involving construction of derivative approximants based on implicit finite-difference schemes under domain-decomposition with closures iterated as described in §II.5. Equation (74) also allows for direct inspection of the eigenvalues which may be used to rapidly gain insight on stability properties for wave propagation problems. If biased schemes are instead considered then the iteration matrix is constructed analogously and therefore we omit details. While we detailed the construction of for a one-dimensional grid this may be extended to higher dimensions through use of Kronecker products.
III Applications
We have considered construction of a variety of compact finite-difference schemes together with their properties under domain-decomposition in the numerical approximation of function derivatives. The goal of this section is to demonstrate their application to wave propagation problems. In particular, numerical solution of the standard, two dimensional homogeneous advection equation Evans (2010) where sensitivity to characteristics LeVeque (2007) will allow for careful examination of biased schemes. Thereafter the shifted wave equation Chirvasa and Husa (2010) will be considered as it provides a simple model Calabrese (2005); Babiuc et al. (2006) for first order in time, second order in space partial differential equations governing numerical relativity formulations such as BSSNOK Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1999) and Z4c Bernuzzi and Hilditch (2010); Ruiz et al. (2011); Weyhausen et al. (2012); Hilditch et al. (2013). These toy problems we take as an initial stepping-stone allowing for verification of implementation and feasibility of the approach thereafter attention is turned to the Z4c system in §III.3.
III.1 Two dimensional advection
We begin with the two-dimensional, spatially periodic, initial value problem for the advection equation. Define with as in §II.1. Consider smooth satisfying:
[TABLE]
with subject to the unit-speed constraint . The direction of propagation of an initial advected profile is parametrized through selection of the angle that forms with the axis. Under single domain, uniformly spaced, spatial discretization with and samples in and directions the system of Eq.(75) becomes:
[TABLE]
where enforcement of the periodic boundary conditions is considered embedded within the discretized derivatives and initial conditions are compatible with periodicity. For simplicity in this section we will restrict attention to cell-centered sampling. We may also view Eq.(76) as where acts on the state vector comprised of elements . Formal solution is provided through matrix exponentiation . Consequently to describe propagating solutions to the semi-discrete problem with bounded amplitude we require that satisfies (see discussion of §II.2). Due to this when utilizing biased schemes for care needs to be taken to appropriately upwind. This can be achieved by ensuring that for we select a biased scheme777These choices can be understood through considering the propagating, single mode solution to Eq.(20) and the modified wavenumber. with whereas for we require . Suppose and where are the previously introduced biased schemes. During calculations based on Eq.(76) involving approximants of we replace:
[TABLE]
In the case of decomposition of into sub-domains the form of Eq.(76) does not change if indices are instead viewed as local to a sub-domain. Furthermore in light of Eq.(74) a global description of the domain decomposed problem with iterated closures for implicitly specified derivative approximants is provided through . This is of particular utility during consideration of the fully discretized problem where is uniformly partitioned into time-steps as stability properties may be assessed based on . For numerical calculation of the eigenvalues we do not need to assemble the full matrix explicitly. It is the case that if then for and we have that Schäcke (2004). We fix and compare the spectrum for constructed with respect to a single domain where with that of where and corresponding to the domain-decomposed problem for a variety of schemes in Fig.9.
We find that the spectrum is deformed during domain-decomposition on account of decoupling the implicit specification of derivative approximants with explicit closures however as is evident a single iteration is sufficient to almost entirely mitigate this effect. Additionally we observe that for the biased stencils of imposing Eq.(77) leads to all schemes investigated satisfying to numerical round-off. From Fig.9 we also see that selecting based on a CFL of satisfies a neccessary (and in this case sufficient) condition for stability of the fully discrete problem evolved with the classical ERK4 as all eigenvalues are contained with the stability polynomial of the method. Furthermore we have verified that these properties are robust under changes to the direction of together with changes in the number of samples and sub-domains.
Having confirmed stability properties through numerical spectra we now turn to numerical solution of the evolution problem. In the semi-discrete case we make use of Eq.(76) supplemented by Eq.(77) for biased schemes and consider formal exponential integration. Recall that and consequently the commutator vanishes. Thus through the use of the Zassenhaus formula Suzuki (1977) we have and therefore:
[TABLE]
where matrix exponentials we evaluate numerically based on Al-Mohy and Higham (2010).
In the fully-discrete case we make use of classical ERK4 which provides a time-integrator of formal order of accuracy Butcher (2008); Hairer and Wanner (2010). In the context of hyperbolic evolution for a spatial discretization of it is common to add dissipation involving derivatives of degree through the standard Kreiss-Oliger prescription Gustafsson et al. (2013) on each field component and in each spatial direction:
[TABLE]
where regulates the strength of the added dissipation. The derivative product may be evaluated through of Tab.7. The choice of order and degree is made such that in the case of non-linear hyperbolic PDE (as tested in §III.3) stability in appropriate norm may be demonstrated and attained for various classes of problems Gustafsson et al. (2013). While our numerical experiments show that addition of dissipation does not appear strictly necessary for full discretization of Eq.(75) for the investigated our purpose here is to ensure that conventions are consistently selected in the context of this simple problem.
For initial conditions we form with is defined in Eq.(45) and parameters selected as whereas is defined in Eq.(46) and we take . As the numerical solutions approximate Eq.(75) we compare the sampled, advected initial condition pointwise to at for the methods described. Given fixed spatial resolution taken to be uniform in and directions Fig.10 depicts the error associated with the result of exponential integration (Eq.(78)) and similarly that of ERK4 based solution for the single domain and domain-decomposed approaches at a variety of angles for .
From Fig.10a we see that schemes tend to have a more pronounced error as the direction of propagation tends towards alignment along the axis. This can be understood from the comparing the spectral content of and . We find that for whereas for at and consequently higher resolution would be required in the direction to achieve a more uniform error as is varied. Crucially we observe that biased schemes remain stable (and indeed error is symmetric under reflection about axes) when upwinding in accordance with Eq.(77). As can be seen in Fig.10b by judiciously selecting CFL when evolving with ERK4 the temporal error may be made comparable to that of the spatial discretization and consequently converges to the error associated with exponential integration (EI) of the semi-discrete problem indicating consistency between approaches. Interestingly we find that for some schemes an intermediate regime of exists where ERK4 outperforms EI. One possible explanation for this is that while in this work we exclusively focus on tuning modified wavenumber for derivative approximants based on arguments involving semi-discretization; ERK schemes propagating the fully discrete system may further modify dispersion and dissipation Hu et al. (1996). For the present setup we find that utilizing implicit schemes for specification of spatial derivatives reduces maximum error by a factor of for the scheme or for when compared with the standard explicit finite-difference approach . This maintains the trend observed in the grid convergence study of §II.5.
III.2 Shifted wave equation
As another example we consider the shifted wave equation Chirvasa and Husa (2010). This will also allow us to test the previously introduced second degree implicit derivative schemes. Suppose is a Lorentzian manifold endowed with metric . For simplicity, suppose and introduce the projector888Geometric quantities may feature space-time or spatial indices respectively. Juxtaposition of an index that appears raised and lowered implies summation on that index. where corresponds to a shift vector. The homogeneous, scalar, wave equation may then be written as:
[TABLE]
We reduce Eq.(80) to a first order in time system by defining the auxiliary field:
[TABLE]
The shifted wave equation can now be written as:
[TABLE]
subject to supplementation with suitable initial and boundary conditions. As our goal here is to provide a dynamical test of the second degree derivative schemes we simplify the problem by working in dimensions, impose spatial periodicity, and freeze to be constant with selected as a flat background. Under these assumptions we thus seek smooth and satisfying:
[TABLE]
where . In particular reduces Eq.(83) to the standard un-shifted case. For the single domain, uniformly spaced, spatial discretization with samples on we define:
[TABLE]
such that the semi-discrete formulation of Eq.(83) is given by:
[TABLE]
where enforcement of the periodic boundary conditions is considered embedded within the discretized derivatives appearing in and initial conditions are compatible with periodicity. If biased derivative approximation schemes are selected in Eq.(84) we consider suitable upwinding based on Eq.(77) with replaced by . For simplicity in this section we assume cell-centered sampling. As in the case of the advection problem we may view the formal solution to the semi-discrete problem as provided through matrix exponentiation. In passing to the domain-decomposed problem over with sub-domains the general form of Eq.(84) and Eq.(85) remains unchanged if discretized derivative operators are understood in the sense of Eq.(74).
Full-discretization with sampled with uniform time-steps entails that the eigenvalues are to be assessed. However of Eq.(84) is non-normal and satisfying for all together with containment within the ERK4 stability region only provide a necessary condition for stability. On the other hand it is known Calabrese et al. (2006); Chirvasa and Husa (2010) that for this system fully-discrete stability can be established for explicit finite differencing based on a modified norm involving additional derivative terms. We do not seek to extend this analytical result here for implicit derivative schemes but rather take it as a guide. In order to gain insight on the properties of the implicit derivative scheme closures and hybrid iteration for this system featuring multiple derivative degrees we select , , and and investigate in Fig.11.
As in the case of the advection problem of §III.1 we find that the spectrum of is deformed when compared with the single domain approach. However as can be seen in Fig.11b this can be mitigated through the use of a single hybrid iteration.
We now consider propagating the initial condition:
[TABLE]
with , , and , which describes a left-ward propagating Gaussian. With regard to the continuum problem Eq.(83) after a crossing time of the profile described by Eq.(86) is reconstructed. This feature must appear in the discretized solution if it is accurate and consequently allows us to characterize error stroboscopically at integer multiples of the crossing time through direct comparison with the initial profile. We perform convergence testing based on solution to the semi-discrete problem posed in single and domain-decomposed form together with verification of convergence of the latter in the fully-discretized context with ERK4 in Fig.12.
As is clear from Fig.12a utilizing implicit schemes for specification of spatial derivatives at can reduce maximum error by a factor of for the scheme or for when compared with the standard explicit finite-difference approach . From Fig.12b we find agreement with observations made in Chirvasa and Husa (2010) in that for non-zero there exists a CFL regime where utilizing a combination of and reduces error when compared to full centering (i.e. using instead for both derivative degrees).
III.3 Z4c: system description
A primary target application of this work is numerical solution of the Cauchy problem for the Einstein field equations (EFE). For problems in the absence of symmetries, this requires considerable computational infrastructure and highly performant code. We therefore utilize the octree-based, adaptive mesh refinement (AMR) infrastructure offered by GR-Athena++ Daszuta et al. (2021) where hybrid MPI-OMP provides parallelism at scale. The generalized finite-difference schemes investigated in prior sections we have coupled via a header-only, templated C++ library. Before describing our numerical tests we briefly recall some formulation details. In the context of the Cauchy problem for the EFE the conformal formulations of BSSNOK Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1999) and c Bernuzzi and Hilditch (2010); Ruiz et al. (2011); Weyhausen et al. (2012); Hilditch et al. (2013), and the generalized harmonic gauge (GHG) approach Friedrich (1985); Pretorius (2005); Lindblom et al. (2006) have found success in simulation of a wide variety of physical problems. The former two formulations leverage consideration of a globally hyperbolic space-time as foliated by a family of non-intersecting spatial slices; initial data is provided on a selected slice of the foliation and well-posed evolution equations compatible with the EFE must be prescribed and solved999For an introductory account of geometric, analytical, and numerical considerations see the textbooks Gourgoulhon (2012); Alcubierre (2008); Baumgarte and Shapiro (2010); Shibata (2016). over . In particular the formulation Bona et al. (2003) augments the EFE through introduction of an auxiliary, dynamical vector field and first-order covariant derivatives thereof. This results in evolution equations involving the variables where and are the induced metric and extrinsic curvature associated with respectively whereas and are normal and spatial projections of . Furthermore Hamiltonian, momentum, and auxiliary vector constraints must also be satisfied such that a numerical space-time is faithful to a solution of the standard EFE. Importantly for a space-time without boundary if for some element of the foliation (e.g. ) then analytically this property extends for all (Bona et al., 2003). This strategy leads to a framework wherein certain strengths of BSSNOK and GHG can be blended Bona et al. (2010). Isolating a spatial conformal degree of freedom leads to c which is a conformal, free-evolution scheme featuring prescribable constraint damping. One defines101010Existence of a global chart with Cartesian coordinatization is assumed throughout.:
[TABLE]
with and where is the determinant of . Define:
[TABLE]
The c system has dynamical variables \big{(}\chi,\,\tilde{\gamma}{}_{ij},\,\hat{K},\tilde{A}{}_{ij},\,\Theta,\,\tilde{\Gamma}{}^{i}\big{)} which in vacuum are governed by the evolution equations:
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
where is the covariant derivative compatible with , and are constraint damping parameters, and in Eq.(93) the trace-free (tf) operation is computed with respect to . The intrinsic curvature is split as and utilizing the conformal connection compatible with allows us to write:
[TABLE]
and:
[TABLE]
Evolved variables must satisfy the dynamical constraints which in terms of transformed variables :
[TABLE]
[TABLE]
[TABLE]
The transformation of Eq.(87) also implies the algebraic constraints \mathcal{C}_{A}:=\big{(}\ln(\tilde{\gamma}),\,\tilde{\gamma}{}^{ij}\tilde{A}{}_{ij}\big{)}=0 which are enforced during a numerical evolution for consistency111111In particular, coupling to the puncture gauge with enforcement of results in a strongly hyperbolic and well-posed system Cao and Hilditch (2012); Bernuzzi and Nagar (2010). Consequently the fully-discrete evolution enforces this condition at each time sub-step..
The c system must be further supplemented by gauge conditions where the lapse and shift describe how the elements of the foliation piece together. In this work we make use of the moving puncture gauge which consists of the Bona-Másso lapse Bona et al. (1995) and the gamma-driver shift Alcubierre et al. (2003):
[TABLE]
where the lapse variant is selected through together with , and is a specifiable damping parameter. During a subset of numerical tests we also make use of the harmonic gauge condition which sets in the dynamical relation for of Eq.(101) and the shift evolution becomes:
[TABLE]
Semi-discretization proceeds as in prior sections however a few remarks are in order. In the evolution equations fields to be sampled over a domain are sampled at points assembled from a tensor product grid of . Points are of the form where are to be understood as grid indices for a given axis. Sampled fields thus carry suppressed grid indices e.g. and similarly . Derivatives are approximated according to:
[TABLE]
Some care is required with mixed partial derivatives such as as they commute when applied to functions. Consequently we explicitly symmetrize the discrete approximants for and explicitly replace for . In the case of shift-advective terms derivative approximants are treated as in Eq.(77) where the value of is replaced by the pointwise value of the relevant shift vector component sampled on the underlying grid. In passing to the fully-discrete setting the implementation within GR-Athena++ performs time-evolution using the order RK low-storage method of Ketcheson (2010). To ensure numerical stability Kreiss-Oliger dissipation is incorporated according to Eq.(79) and is applied to each field component, in each spatial direction.
III.4 Z4c: numerical tests - gauge wave evolution
Armed with the c system (§III.3) our first goal is to ensure that coupling our header-only, templated, generalized finite-difference code to GR-Athena++ leads to successful evolution on simple test problems. For this we make use of suitably modified Apples with Apples (AwA) test-beds Alcubierre et al. (2004); Babiuc et al. (2008); Daverio et al. (2018). The intention here is to quantify solution quality on small scale idealized problems and to probe for any sources of potential instability that may have been introduced through modifying derivative approximants away from the well-known properties of standard finite-difference (FD) explored in Daszuta et al. (2021). Furthermore it allows for evolution of c while gradually bridging the complexity gap from linear propagation problems of §III.1 and §III.2 towards binary black hole merger discussed later. Due to the task-based infrastructure and sophisticated treatment of sub-domain communication as a first step, compact FD approximants computed by the code we have coupled to GR-Athena++ relies solely on the closures described in §II.4 and does not feature the hybrid strategy procedure described in §II.5.
The AwA test-beds are specified for of topology hence is considered as periodic in each spatial direction. The effective dynamics occur over one (or two) spatial dimensions depending on the details of the test. In these directions the grid is taken as such that . For verification of the full system the remaining direction(s) fix this spacing and take the sampling parameter as . During domain-decomposition partitioning is performed as discussed previously over directions with effective dynamics where sampled sub-domains have and are extended by ghost points to facilitate communication and derivative stencil evaluation. Overall we set with serving to adjust resolution during convergence tests as required. This choice is motivated by the formal spatial order of the schemes we employ. For the tests presented here numerical evolution is performed over with CFL of and we select constraint damping parameters and together with Kreiss-Oliger dissipation unless otherwise stated.
The initial AwA test we performed was that of robust stability. An initial spatial slice of Minkowski space-time and to each sampled grid point an independent uniform random value drawn from with is added. This choice of effectively linearizes the system. We utilized the moving puncture gauge of Eq.(101) with initial conditions and . The shift-damping parameter is set as . The dynamics are considered to be effectively one-dimensional. The quantity together with constraints such as were monitored over the course of a calculation. As in the case of finite-difference tests made in Daszuta et al. (2021) we found that when adopting a wide variety of combinations of compact stencils as selected in the toy problems of §III.1 and §III.2 leads to contraint quantities decaying towards a plateau in norm with values comparable to the FD case. This indicates that error associated with numerical evolution of the principal part of the c system does not appear to induce spurious growth of unstable exponential modes when utilizing compact finite-difference (CFD) approximants. The second AwA test is the linearized wave test. Effective one-dimensional dynamics are induced through where together with with remaining components zero. Gauge is chosen as in the robust stability test. An amplitude forces non-linear terms to numerical round-off thus linearizing the c system when numerical calculations are performed in double-precision arithmetic. While these choices lead to a numerical solution that is well described as a simple travelling (i.e. advected) wave, the puncture gauge is not necessarily compatible with pure advection, and furthermore the initial data are constraint violating Cao and Hilditch (2012). We thus focus instead on the gauge wave tests as they describe propagation of simple constraint satisfying data to the full non-linear c system.
Consider the AwA aligned, unshifted, gauge wave test in the form presented in Daverio et al. (2018) with components permuted for propagation along :
[TABLE]
[TABLE]
[TABLE]
[TABLE]
where:
[TABLE]
Evaluating Eqs. 104, 105, 106 and 107 at and setting together with yields a one-parameter family of initial data parametrized by amplitude . We select as it is known that large values (e.g. ) can lead to issues with stability in a variety of formulations and regardless of puncture or harmonic gauge choice (Daverio et al., 2018; Cao and Hilditch, 2012; Boyle et al., 2007). For compatiblity with the analytical gauge we make use of the harmonic prescription of Eq.(101) for the lapse evolution with and Eq.(102) for the shift. To assess solution quality we consider numerical evolution repeated at a triplet of resolutions where . Convergence rates of the overall approximation of the corresponding field data may be examined by comparing differences of solutions at distinct resolutions. For an approximation of order one finds based on Taylor expansion that where we have introduced the so-called convergence factor:
[TABLE]
As we know the space-time metric over the full foliation we may directly compare the RMS error of the numerical solution at any sampled . Additionally, we may inspect the associated phase error as suggested in Daverio et al. (2018). To do this evolved, field data on sub-domains is reassembled on a single, discretized domain (e.g. ) with respect to which we define:
[TABLE]
where we denote the phase of each complex coefficient . Phase error can thus be quantified as . We also also compute the offset of the numerical profile relative to the amplitude through . To simultaneously assess convergence and absolute error define the normalized error:
[TABLE]
In the convergent regime the resolution triplet induced by satisfies with absolute scale given by . Results of a calculation involving are shown in Fig.13.
In the case of the aligned, unshifted, gauge wave test we see that shift-advective terms are analytically zero. Thus we restrict discussion to replacement of standard FD approximants with centered CFD. As can be seen in Fig.13a we find clean order convergence for the derivative schemes investigated. Furthermore we see that all tested CFD schemes outperform the FD scheme. In particular, for the phase error at we find that reduces error by a factor of when compared with whereas the spectrally-tuned reduces error by a factor of . The best improvement in phase error is found for the approximant which when compared to FD results in a reduction in error of a factor of . Considering instead leads to qualitatively similar conclusions. Constraints are also well satisfied and better preserved when using CFD schemes – indeed for the RMS of we find a reduction of a factor of when compared with FD.
In order to test shift-advective terms we consider the aligned, shifted gauge wave test where Daverio et al. (2018):
[TABLE]
[TABLE]
[TABLE]
[TABLE]
[TABLE]
with and are defined as in Eq.(108). The method of setup and quantities analyzed are as in the unshifted case and we once again select . During this test for standard FD resolution in induced through the triplet whereas for CFD we set . As is common we may compare again with the normalization prescription of Eq.(111). We depict the result of numerical evolution in Fig.14.
As shift-advective terms are now non-zero we additionally make use of upwinded stencils. In Fig.14a we again find clean order convergence for the derivative schemes investigated (cf. Fig.13a). In a similar vein as the unshifted test errors are reduced when making use of CFD schemes. Comparing phase error at between schemes shows that replacing only the centered derivatives with reduces error by a factor of when compared with the fully explicit . If is utilized then instead find a reduction in phase error of a factor of when compared with FD. The best improvement in phase error for schemes tested here is found for the approximants which when compared to FD results in a reduction in error of a factor of .
III.5 Z4c: numerical tests - binary black hole evolution
We close our numerical tests with a preliminary investigation of binary black hole (BBH) evolution utilizing CFD schemes. This test departs from those presented earlier in this work as the underlying computational domain is no longer periodic and non-trivial boundary conditions (BC) must be applied on evolved field components. In particular, the c dynamical equations supplemented by gauge conditions populate on whereas Sommerfeld BC are applied to the field components . In addition, due to the range of spatial scales, we make use of adaptive mesh refinement (AMR) for computational efficiency. Suppose is an element of a domain-decomposition of an of interest. Within GR-Athena++ one can prescribe a conditional (i.e. a target resolution over a region of within a given distance of some feature described by the evolved fields) which controls the AMR. The sub-domain is then recursively (de)refined to satisfy the conditional under the further restriction the nearest-neighbour sub-domains can differ in resolution by at most a ratio. Extensive details on the treatment of BC and AMR made in GR-Athena++ which we utilize for this problem are described in Daszuta et al. (2021).
To model the BBH evolution itself initial data compatible with the constraints must first be provided. To this end we consider the initial geometry as modelled by Brill-Lindquist wormhole topology describing black holes with disconnected, asymptotically flat ends. Each disconnected end is diffeomorphic to minus a compact ball Dain (2002). An end is compactified and identified with a point on . The coordinate singularity that occurs at a given is a so-called puncture which describes the location of a black hole. This allows the constraints to be solved based on Ansorg et al. (2004) thus providing initial data. Gauge conditions are initialized based on a “precollapsed” lapse and zero-shift Campanelli et al. (2006). The damping parameter in Eq.(101) is now taken as which is fixed in terms of the ADM mass Arnowitt et al. (2008) of the underlying system. The AMR criterion is based on a mock “box-in-box” oct-tree structure which adapts resolution based on tracking puncture centers during the course of a simulation. Given an overall that has been domain-decomposed and refined the resolution at both punctures is controlled by the maximum number of refinement levels as . Unless otherwise stated we make use of for the number of samples along each direction taken on a sub-domain.
As we would like to inspect convergence for a variety of derivative approximants we investigate an equal-mass initial configuration leading to a short evolution. The BBH system has two non-spinning punctures, initially centered on-axis at with initial momenta , and with bare-masses . This configuration results in orbits before merger at evolution time . For the overall grid extent we select such that is causally disconnected from the interior strong-field dynamics during the course of the initial inspiral through merger. As a diagnostic for assessing numerical simulation quality we consider extracted gravitational wave (GW) content associated with the strong-field dynamics. This is done by first assembling the four-dimensional Weyl tensor from the evolved c variables. Subsequent projection over a suitable null tetrad Brügmann et al. (2008); Daszuta et al. (2021) yields the complex, out-going Weyl scalar . A mode-decomposition with respect to spherical harmonics of spin-weight at extraction radius based on numerical quadrature over geodesic spheres Daszuta et al. (2021) furnishes us with radiated GW content in the mode from .
We compute the dominant mode for simulations involving a variety of derivative approximants and choices of and show the result in Fig.15.
We find that stable BBH evolution is possible utilizing the CFD schemes discussed in this work with resulting phase error of extracted compatible with the order of accuracy design of the underlying spatial derivative approximant schemes. In Fig.15b we observe that replacing centered FD with CFD leads to a reduction in the associated phase error (when comparing fixed ) of a factor . A variety of effects influence this factor. As shown in §III.4 it may be important to adequately treat shift-advective terms for maximum improvement. Another delicate matter is transferring field data between between sub-domains at differing levels of refinement. In GR-Athena++ this is achieved through use of prolongation and restriction operations based on centered polynomial interpolation at formal order of accuracy matched to the underlying (C)FD scheme. Without additional care this may potentially degrade properties of the modified wavenumber discussed in §II.3. Additionally as observed during grid convergence tests (see e.g. Fig.6a) when CFD stencils are utilized in the context of domain-decomposition error tends to accumulate at sub-domain boundaries. Transferring data between differing levels of refinement also occurs in this region and it is not entirely evident as to whether polynomial interpolation will amplify or diminish this source of error. Nonetheless we have described how error in the vicinity of can be mitigated through usage of a hybrid-communication strategy in §II.5. We have not yet implemented this strategy as this would require involved modification to core GR-Athena++ functionality. We aim to address this in future.
IV Summary and conclusion
In this work we have shown that the unified compact finite difference (CFD) framework of Deshpande et al. (2019) may be extended in numerical generation of new schemes which may be biased or centered, of arbitrary extent, and involve not only function data but also prescribed function derivative data. Upon fixing formal order of accuracy Taylor matching yields a linear system specifying a stencil. When the aforementioned system is underdetermined we may minimize a functional characterizing spectral error. This allows us to further extend the approach to construct implicit extensions to the Hermite methods described in Fornberg (2020). Facilitating simpler construction of such schemes allows for rapid experimentation on practical problems and consequently we also have open-sourced our notebook not .
Large-scale problems crucially depend on exploiting parallelism for efficiency and consequently in order to treat solution of the implicit problem specifying a given CFD scheme we investigated modification of a dispersion-relation-preserving (DRP) method due to Chen et al. (2021). This featured treating domain-decomposition with tailored closures for decoupling the implicit, linear system specifying a derivative approximant to individual, decoupled sub-domains. We applied the method to generalize the biased, first degree order CCU scheme Chen et al. (2021) to CCU and CCU which have formal order of accuracy of and order respectively. We also show directly that the DRP method may be exactly employed for biased second degree schemes, and furthermore, applied approximately in decoupling centered CFD methods. Grid convergence tests involving smooth functions revealed that edge artifacts can be induced under domain-decomposition with this strategy. We therefore proposed a hybrid-communication strategy that iterates upon results from decoupled sub-domains and allows for mitigation of error. Efficacy was verified through grid convergence testing on smooth functions. As a further test geared towards wave-propagation problems the numerical spectra of the semi-discretized two-dimensional advection equation, and shifted wave equation were directly inspected. This was achieved through an embedded description of sub-domain communication, iterated closure, and differentiation as block-partitioned system. Indeed domain-decomposition and DRP leads to deformed spectra which we demonstrate can be mitigated with our proposed hybrid-communication approach.
We numerically verified anticipated stability properties for the aforementioned toy-problems through solution of semi-discretized and fully-discretized formulations based on exponential integration and explicit Runge-Kutta methods respectively. Spatial discretization was fixed to have formal order of accuracy . In the case of two-dimensional advection we found that when compared to standard finite-difference (FD) the well-known Padé scheme allowed for a reduction in maximum error by a factor of whereas our new spectrally-tuned scheme attained a factor of both of which remained robust under domain-decomposition. Comparable improvements were observed for the shifted wave equation.
As a first novel application we considered evolution of the c formulation of numerical relativity. To this end we wrote a stand-alone, header-only, templated, generalized finite difference library which we coupled to the GR-Athena++ Daszuta et al. (2021) code. This is a dramatic increase in complexity of the underlying system being evolved as not only is the system now quasi-linear, with significantly more involved stability criteria, the number of independent field components is larger. Indeed this precluded a numerical investigation of spectra associated with semi-discretization. As a first step we coupled our code to GR-Athena++ and investigated introducing CFD with the (approximate) DRP prescription for closure under domain-decomposition. In the context of the AwA aligned, unshifted, gauge wave evolution test (without shift-advective terms) this allowed for reducing the phase error of propagated metric components when compared with standard FD by a factor of for and for . For the case of the shifted variant of this test we found that it is important to replace not only centered FD terms with an appropriate centered CFD prescription but to also treat lop-sided FD terms with biased CFD. The best improvement here we found for the combinations test was which reduced phase error by a factor of when compared with standard FD.
We also considered a preliminary application to binary-black-hole evolution. For this class of problem the discrete treatment of the underlying domain in GR-Athena++ features adaptive mesh refinement (AMR) for computational efficiency. During investigation of a orbit, equal mass, non-spinning BBH simulation we found that utilizing CFD leads to stable evolution. By inspecting the phase of the dominant mode of the extracted gravitational waves we verified convergence compatible with a order trend. A factor reduction in phase error at fixed resolution when compared with standard FD was also observed. This reduced efficiency in error reduction can potentially be attributed to a combination of: insufficiently aggressive CFL, approximate DRP and lack of closure iteration based on hybrid-communication, and the overly simple treatment of level-to-level transfer of field data between sub-domains at differing levels of refinement. The influence and precise tuning of dissipation together with potential replacement by compact filters Kim (2013) may also significantly affect error reduction efficiency. We leave a thorough investigation to future work.
Acknowledgements.
The author thanks David Radice and Sebastiano Bernuzzi for discussions and constructive input during several stages of preparing this work. Special thanks to Martin Bernreuther for assistance with high performance computing support. The author is indebted to Beppe Starnazza. The author acknowledges support by the EU H2020 under ERC Starting Grant, no. BinGraSp-714626. Simulations with GR-Athena++ were performed on the ARA cluster at Friedrich Schiller University Jena, SuperMUC-NG at the Leibniz-Rechenzentrum (LRZ, www.lrz.de) Munich, and HPE Apollo Hawk at the High Performance Computing Center Stuttgart (HLRS, www.hlrs.de). The ARA cluster is funded in part by DFG grants INST 275/334-1 FUGG and INST 275/363-1 FUGG, and ERC starting Grant, grant agreement no. BinGraSp-714626. The author acknowledges the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS Supercomputer SuperMUC-NG at LRZ (allocation pn68wi). The author acknowledges HLRS for funding this project by providing access to the supercomputer HPE Apollo Hawk under the grant number INTRHYGUE/44215.
Appendix A Collection of derivative approximants
For convenience we collect a variety of derivative approximants as specified through the coefficients that enter Eq.(3) that have been constructed utilizing the method outlined in §II. We tabulate expansions of the associated relative error in the normalized wavenumber where is defined in Eq.(12).
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Lele (1992) S. K. Lele, Journal of Computational Physics 103 , 16 (1992) . · doi ↗
- 2Hirsh (1975) R. S. Hirsh, Journal of Computational Physics 19 , 90 (1975) . · doi ↗
- 3Carpenter et al. (1994) M. H. Carpenter, D. Gottlieb, and S. Abarbanel, Journal of Computational Physics 111 , 220 (1994) . · doi ↗
- 4Mehra and Patel (2017) M. Mehra and K. S. Patel, ACM Transactions on Mathematical Software 44 , 23:1 (2017) . · doi ↗
- 5Fu and Ma (1997) D. Fu and Y. Ma, Journal of Computational Physics 134 , 1 (1997) . · doi ↗
- 6Wang et al. (2013) Z. Wang, K. Fidkowski, R. Abgrall, F. Bassi, D. Caraeni, A. Cary, H. Deconinck, R. Hartmann, K. Hillewaert, H. Huynh, N. Kroll, G. May, P.-O. Persson, B. van Leer, and M. Visbal, International Journal for Numerical Methods in Fluids 72 , 811 (2013) , _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/fld.3767. · doi ↗
- 7Qin et al. (2014) Q. Qin, Z. A. Xia, and Z. F. Tian, International Journal of Heat and Mass Transfer 71 , 405 (2014) . · doi ↗
- 8Chen et al. (2021) J. Chen, P. Yu, H. Ouyang, and Z. F. Tian, Journal of Scientific Computing 86 , 5 (2021) . · doi ↗