$$Date$$

Abstract

Calculate ﬁrst derivates of grid functions using ﬁnite diﬀerence stencils that satisfy summation by parts.

Given a discretization ${x}_{0}\dots {x}_{N}$ of a computational domain $x\in \left[a,b\right]$ with gridspacing $h$ a one dimensional ﬁnite diﬀerence operator approximation to a ﬁrst derivative, $D$, is said to satisfy summation by parts (SBP) with respect to a scalar product (deﬁned by its coeﬃcients ${\sigma}_{ij}$)

$${\u27e8u,v\u27e9}_{h}=h\sum _{i=0}^{N}{u}_{i}{v}_{j}{\sigma}_{ij}$$ | (1) |

if the property

$${\u27e8u,Dv\u27e9}_{h}+{\u27e8Du,v\u27e9}_{h}={\left.\left(uv\right)\right|}_{a}^{b}$$ | (2) |

is satisﬁed for all possible gridfunctions $u$ and $v$.

At a given ﬁnite diﬀerence order, there are several diﬀerent ways of doing this depending on the structure of the scalar product. The three commonly considered cases are the diagonal norm, the restricted full norm and the full norm (see ﬁgure 1 for the structure).

$$\sigma =\left(\begin{array}{cccccc}\hfill x\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill x\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill x\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill x\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill 1\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill \end{array}\right),\text{}\sigma =\left(\begin{array}{ccccccc}\hfill x\hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill 1\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill \end{array}\right),\text{}\sigma =\left(\begin{array}{cccccc}\hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill \hfill & \hfill \hfill \\ \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill \hfill & \hfill \hfill \\ \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill \hfill & \hfill \hfill \\ \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill x\hfill & \hfill \hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill 1\hfill & \hfill \hfill \\ \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \hfill & \hfill \ddots \hfill \end{array}\right)$$

In the following we denote the order of accuracy at the boundary by
$\tau $, the order in the interior
by $s$ and the width of
the boundary region^{1}
by $r$.

For the diagonal norm case it turns out that with $r=2\tau $ it is possible to ﬁnd SBP operators with $s=2\tau $ (at least when $\tau \le 4$), i.e. the order of accuracy at the boundary is half the order in the interior. For the restricted full norm case $\tau =s-1$ when $r=\tau +2$ and for the full norm case $\tau =s-1$ when $r=\tau +1$.

The operators are named after their norm and their interior and boundary orders. Thus, for example, we talk about diagonal norm 6-3 operators and restricted full norm 4-3 operators.

In the diagonal case the 2-1 and 4-3 operators are unique whereas the 6-3 and 8-4 operators have 1 and 3 free parameters, respectively. In the restricted full norm case the 4-3 operators have 3 free parameters and the 6-5 operators have 4, while in the full norm case the number of free parameters is 1 less than the restricted full case.

Currently this thorn implements only diagonal and restricted full norm SBP operators. The
diagonal norm 2-1, 4-2 and 6-3 and the restricted full norm 4-3 are the ones listed in [1] where
in the presence of free parameters the set of parameters giving a minimal bandwidth have been
chosen. For the diagonal norm 8-4 case the minimal bandwidth choice are to restrictive with respect
to the Courant factor and in this case parameters are chosen so as to maximize the Courant
factor^{2} .
In addition the restricted full norm 6-5 operator has been implemented. This was calculated using a
Mathematica script kindly provided by José M. Martín-García.

For the restricted full norm 6-5 SBP operators a compatible dissipation operator has been implemented as well based on [3].

The basic functionality of calculating ﬁnite diﬀerence approximations are performed by a set of aliased functions. In this way things are kept ﬂexible in the sense that diﬀerent thorns each can provide their own way of calculating derivatives but if they use the same calling interface the user can switch between diﬀerent derivative schemes by simply activating the appropriate ﬁnite diﬀerencing thorn. Also the same derivative routines can be called from both C and Fortran.

There are currently two diﬀerent calling interfaces corresponding to diﬀerent levels of ﬂexibility. In the ﬁrst case the fact that derivatives are approximated by ﬁnite diﬀerencing is completely hidden for the user. Using this interface the same physics thorn could in principle be used with a ﬁnite diﬀerence driver, a ﬁnite elements driver or a spectral driver. There is however a price to pay for this ﬂexibility due to the fact that the operations are performed on whole grid functions at a time so that storage for the derivatives has to be provided. This can signiﬁcantly increase the amount of memery required for the evolution. The alternative interface returns instead the ﬁnite diﬀerence coeﬃcients and allows the user to calculate derivatives on a pointwise basis. This can signiﬁcantly save memory but limits the physics thorn to use ﬁnite diﬀerences.

The thorn is currently still in development and so is not generally available. Access can be requested by contacting the thorn maintainer.

In order to use the derivative routines the appropriate declarations have to be added to your thorns interface.ccl ﬁle. To use the interface that works on whole grid functions you have to add:

SUBROUTINE Diff_gv ( CCTK_POINTER_TO_CONST IN cctkGH, \

CCTK_INT IN dir, \

CCTK_REAL IN ARRAY var, \

CCTK_REAL OUT ARRAY dvar, \

CCTK_INT IN table_handle )

USES FUNCTION Diff_gv

CCTK_INT IN dir, \

CCTK_REAL IN ARRAY var, \

CCTK_REAL OUT ARRAY dvar, \

CCTK_INT IN table_handle )

USES FUNCTION Diff_gv

Here cctkGH is the pointer to the cactus GH, dir is the direction of the derivative (0 for $x$, 1 for $y$ and 2 for $z$), var is the grid function to calculate derivatives of and dvar is the grid function where the derivative is to be stored and table_handle is an integer that contains a handle for a keyword table with optional parameters. If there are no optional paramters just pass in a negative value.

To use the interface that returns the ﬁnite diﬀerence coeﬃcients you have to add:

SUBROUTINE Diff_coeff ( CCTK_POINTER_TO_CONST IN cctkGH, \

CCTK_INT IN dir, \

CCTK_INT IN nsize, \

CCTK_INT OUT ARRAY imin, \

CCTK_INT OUT ARRAY imax, \

CCTK_REAL OUT ARRAY q, \

CCTK_INT IN table_handle )

CCTK_INT IN dir, \

CCTK_INT IN nsize, \

CCTK_INT OUT ARRAY imin, \

CCTK_INT OUT ARRAY imax, \

CCTK_REAL OUT ARRAY q, \

CCTK_INT IN table_handle )

Here, again, cctkGH is the pointer to the cactus GH, dir is the direction of the derivative, nsize is the grid size in the direction you are interested in, imin and imax are 1D integer arrays of size nsize that on return will contain the minimum and maximum index of the stencil in direction dir, q is a 2D real array of size nsize$\times $nsize that on return will contain the coeﬃcients for the derivatives. The table_handle serves the same purpose as in the previous case.

The following piece of Fortran code, shows how to calculate derivatives of a grid function, f, and store the derivatives in the x-, y- and z-directions in the grid functions dxf, dyf and dzf:

call Diff_gv (cctkGH, 0, f, dxf, -1)

call Diff_gv (cctkGH, 1, f, dyf, -1)

call Diff_gv (cctkGH, 2, f, dzf, -1)

call Diff_gv (cctkGH, 1, f, dyf, -1)

call Diff_gv (cctkGH, 2, f, dzf, -1)

In order to use the interface for doing the pointwise derivatives, the following Fortran90 example shows the necessary declarations and an example of how to calculate the derivatives:

CCTK_REAL, dimension(:,:), allocatable :: qx, qy, qz

CCTK_INT, dimension(:), allocatable :: iminx, imaxx, iminy, &

imaxy, iminz, imaxz

CCTK_REAL :: idelx, idely, idelz

CCTK_INT :: i, j, k

allocate ( qx(cctk_lsh(1),cctk_lsh(1)), &

qy(cctk_lsh(2),cctk_lsh(2)), &

qz(cctk_lsh(3),cctk_lsh(3)), &

iminx(cctk_lsh(1)), imaxx(cctk_lsh(1)), &

iminy(cctk_lsh(2)), imaxy(cctk_lsh(2)), &

iminz(cctk_lsh(3)), imaxz(cctk_lsh(3)) )

call Diff_Coeff ( cctkGH, 0, cctk_lsh(1), iminx, imaxx, qx, -1 )

call Diff_Coeff ( cctkGH, 1, cctk_lsh(2), iminy, imaxy, qy, -1 )

call Diff_Coeff ( cctkGH, 2, cctk_lsh(3), iminz, imaxz, qz, -1 )

idelx = 1.0d0 / CCTK_DELTA_SPACE(1)

idely = 1.0d0 / CCTK_DELTA_SPACE(2)

idelz = 1.0d0 / CCTK_DELTA_SPACE(3)

do k = 1, cctk_lsh(3)

do j = 1, cctk_lsh(2)

do i = 1, cctk_lsh(1)

dxf(i,j,k) = idelx * sum ( qx(iminx(i):imaxx(i),i) * &

f(iminx(i):imaxx(i),j,k) )

dyf(i,j,k) = idely * sum ( qy(iminy(j):imaxy(j),j) * &

f(i,iminy(j):imaxy(j),k) )

dzf(i,j,k) = idelz * sum ( qz(iminz(k):imaxz(k),k) * &

f(i,j,iminz(k):imaxz(k)) )

end do

end do

end do

deallocate ( qx, qy, qz, iminx, imaxx, iminy, imaxy, iminz, imaxz )

CCTK_INT, dimension(:), allocatable :: iminx, imaxx, iminy, &

imaxy, iminz, imaxz

CCTK_REAL :: idelx, idely, idelz

CCTK_INT :: i, j, k

allocate ( qx(cctk_lsh(1),cctk_lsh(1)), &

qy(cctk_lsh(2),cctk_lsh(2)), &

qz(cctk_lsh(3),cctk_lsh(3)), &

iminx(cctk_lsh(1)), imaxx(cctk_lsh(1)), &

iminy(cctk_lsh(2)), imaxy(cctk_lsh(2)), &

iminz(cctk_lsh(3)), imaxz(cctk_lsh(3)) )

call Diff_Coeff ( cctkGH, 0, cctk_lsh(1), iminx, imaxx, qx, -1 )

call Diff_Coeff ( cctkGH, 1, cctk_lsh(2), iminy, imaxy, qy, -1 )

call Diff_Coeff ( cctkGH, 2, cctk_lsh(3), iminz, imaxz, qz, -1 )

idelx = 1.0d0 / CCTK_DELTA_SPACE(1)

idely = 1.0d0 / CCTK_DELTA_SPACE(2)

idelz = 1.0d0 / CCTK_DELTA_SPACE(3)

do k = 1, cctk_lsh(3)

do j = 1, cctk_lsh(2)

do i = 1, cctk_lsh(1)

dxf(i,j,k) = idelx * sum ( qx(iminx(i):imaxx(i),i) * &

f(iminx(i):imaxx(i),j,k) )

dyf(i,j,k) = idely * sum ( qy(iminy(j):imaxy(j),j) * &

f(i,iminy(j):imaxy(j),k) )

dzf(i,j,k) = idelz * sum ( qz(iminz(k):imaxz(k),k) * &

f(i,j,iminz(k):imaxz(k)) )

end do

end do

end do

deallocate ( qx, qy, qz, iminx, imaxx, iminy, imaxy, iminz, imaxz )

This thorn is maintained by Peter Diener. Any questions and comments should be directed by e-mail to diener@cct.lsu.edu.

This thorn grew out of the needs of a multipatch relativity code and as such was initially designed to those needs. The addition of the possibility of passing in a handle for a keyword table was done at the request of Jonathan Thornburg and should make it easy to extend the thorn with additional features. The addition of the coeﬃcient interface was done at the request of Bela Szilagyi who felt that the storage overhead was to large for his code. In the near future the thorn will be extended with second derivatives SBP operators.

We thank José M. Martín-García for very kindly providing us with a very well designed mathematica script to calculate the ﬁnite diﬀerence and scalar product coeﬃcients.