xx Summation By Parts

Summation By Parts

Peter Diener <diener@cct.lsu.edu>

Date

Abstract

Calculate first derivates of grid functions using finite difference stencils that satisfy summation by parts.

1 Introduction

Given a discretization of a computational domain with gridspacing a one dimensional finite difference operator approximation to a first derivative, , is said to satisfy summation by parts (SBP) with respect to a scalar product (defined by its coefficients ) if the property is satisfied for all possible gridfunctions and .

At a given finite difference order, there are several different 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 figurex1 for the structure).

Figurex1: The structure of the scalar product matrix in the diagonal case (left), the restricted full case (middle) and the full case (right) for the 4th order interior operators. Only non-zero elements are shown.

In the following we denote the order of accuracy at the boundary by , the order in the interior by and the width of the

boundary region1 ^{1The width of the region where standard centered finite differences are not used.} by .

For the diagonal norm case it turns out that with it is possible to find SBP operators with (at least when ), i.e.xthe order of accuracy at the boundary is half the order in the interior. For the restricted full norm case when and for the full norm case when .

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.

2 Numerical Implementation

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

factor2 2This is done by choosing parameters that minimizes the largest eigenvalue of the amplification matrix for a simple 1D advection equation with a penalty method implementation of periodic boundary conditions as suggested in [2]. . In addition the restricted full norm 6-5 operator has been implemented. This was calculated using a Mathematica script kindly provided by Josx M. Martxn-Garcxa.

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

3 Using This Thorn The basic functionality of calculating finite difference approximations are performed by a set of aliased functions. In this way things are kept flexible in the sense that different thorns each can provide their own way of calculating derivatives but if they use the same calling interface the user can switch between different derivative schemes by simply activating the appropriate finite differencing thorn. Also the same derivative routines can be called from both C and Fortran.

There are currently two different calling interfaces corresponding to different levels of flexibility. In the first case the fact that derivatives are approximated by finite differencing is completely hidden for the user. Using this interface the same physics thorn could in principle be used with a finite difference driver, a finite elements driver or a spectral driver. There is however a price to pay for this flexibility 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 significantly increase the amount of memery required for the evolution. The alternative interface returns instead the finite difference coefficients and allows the user to calculate derivatives on a pointwise basis. This can significantly save memory but limits the physics thorn to use finite differences.

3.1 Obtaining This Thorn

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

3.2 Basic Usage

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

SUBROUTINExDiff_gvx(xCCTK_POINTER_TO_CONSTxINxcctkGH,x\
xxxxxxxxxxxxxxxxxxxxxCCTK_INTxINxdir,x\
xxxxxxxxxxxxxxxxxxxxxCCTK_REALxINxARRAYxvar,x\
xxxxxxxxxxxxxxxxxxxxxCCTK_REALxOUTxARRAYxdvar,x\
xxxxxxxxxxxxxxxxxxxxxCCTK_INTxINxtable_handlex)
USESxFUNCTIONxDiff_gv

Here cctkGH is the pointer to the cactus GH, dir is the direction of the derivative (0 for , 1 for and 2 for ), var is the grid function to calculate derivatives of and dvar is the grid function where the derivative is to be stored and tablexhandle 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 finite difference coefficients you have to add:

SUBROUTINExDiff_coeffx(xCCTK_POINTER_TO_CONSTxINxcctkGH,x\
xxxxxxxxxxxxxxxxxxxxxxxxCCTK_INTxINxdir,x\
xxxxxxxxxxxxxxxxxxxxxxxxCCTK_INTxINxnsize,x\
xxxxxxxxxxxxxxxxxxxxxxxxCCTK_INTxOUTxARRAYximin,x\
xxxxxxxxxxxxxxxxxxxxxxxxCCTK_INTxOUTxARRAYximax,x\
xxxxxxxxxxxxxxxxxxxxxxxxCCTK_REALxOUTxARRAYxq,x\
xxxxxxxxxxxxxxxxxxxxxxxxCCTK_INTxINxtable_handlex)

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 nsizensize that on return will contain the coefficients for the derivatives. The tablexhandle serves the same purpose as in the previous case.

3.3 Examples

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:

xxcallxDiff_gvx(cctkGH,x0,xf,xdxf,x-1)
xxcallxDiff_gvx(cctkGH,x1,xf,xdyf,x-1)
xxcallxDiff_gvx(cctkGH,x2,xf,xdzf,x-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:

xxCCTK_REAL,xdimension(:,:),xallocatablex::xqx,xqy,xqz
xxCCTK_INT,xdimension(:),xallocatablex::ximinx,ximaxx,ximiny,x&
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxximaxy,ximinz,ximaxz
xxCCTK_REALx::xidelx,xidely,xidelz
xxCCTK_INTx::xi,xj,xk
xxallocatex(xqx(cctk_lsh(1),cctk_lsh(1)),x&
xxxxxxxxxxxxxqy(cctk_lsh(2),cctk_lsh(2)),x&
xxxxxxxxxxxxxqz(cctk_lsh(3),cctk_lsh(3)),x&
xxxxxxxxxxxxximinx(cctk_lsh(1)),ximaxx(cctk_lsh(1)),x&
xxxxxxxxxxxxximiny(cctk_lsh(2)),ximaxy(cctk_lsh(2)),x&
xxxxxxxxxxxxximinz(cctk_lsh(3)),ximaxz(cctk_lsh(3))x)
xxcallxDiff_Coeffx(xcctkGH,x0,xcctk_lsh(1),ximinx,ximaxx,xqx,x-1x)
xxcallxDiff_Coeffx(xcctkGH,x1,xcctk_lsh(2),ximiny,ximaxy,xqy,x-1x)
xxcallxDiff_Coeffx(xcctkGH,x2,xcctk_lsh(3),ximinz,ximaxz,xqz,x-1x)
xxidelxx=x1.0d0x/xCCTK_DELTA_SPACE(1)
xxidelyx=x1.0d0x/xCCTK_DELTA_SPACE(2)
xxidelzx=x1.0d0x/xCCTK_DELTA_SPACE(3)
xxdoxkx=x1,xcctk_lsh(3)
xxxxdoxjx=x1,xcctk_lsh(2)
xxxxxxdoxix=x1,xcctk_lsh(1)
xxxxxxxxdxf(i,j,k)x=xidelxx*xsumx(xqx(iminx(i):imaxx(i),i)x*x&
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxf(iminx(i):imaxx(i),j,k)x)
xxxxxxxxdyf(i,j,k)x=xidelyx*xsumx(xqy(iminy(j):imaxy(j),j)x*x&
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxf(i,iminy(j):imaxy(j),k)x)
xxxxxxxxdzf(i,j,k)x=xidelzx*xsumx(xqz(iminz(k):imaxz(k),k)x*x&
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxf(i,j,iminz(k):imaxz(k))x)
xxxxxxendxdo
xxxxendxdo
xxendxdo
xxdeallocatex(xqx,xqy,xqz,ximinx,ximaxx,ximiny,ximaxy,ximinz,ximaxzx)

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

4 History

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 coefficient 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.

5 Acknowledgements

We thank Josx M.xMartxn-Garcxa for very kindly providing us with a very well designed mathematica script to calculate the finite difference and scalar product coefficients.

References

[1]xxxBo Strand, 1994, Journal of Computational Physics, 110, 47–67.

[2]xxxLuis Lehner, Oscar Reula and Manuel Tiglio, in preparation.

[3]xxxKen Mattsson, Magnus Svxrd and Jan Nordstrxm, 2004, Journal of Scientific Computing, 21, 57–79.

6 Parameters

11checkxgridxsizes 10Scope: restrictedx01 BOOLEAN

11Description: Should we check grid sizes and ghost zones

11                         10  x01  Default: yes

11

11dissxfraction  10Scope: restrictedx01  REAL

11Description: Fractional size of the transition region for the full restricted dissipation operator

11 Range 10 x01 Default: 0.2
110:0.5  10

11

11dissipationxtype  10Scope: restrictedx01  KEYWORD

11Description: Type of dissipation operator

11 Range 10 x01 Default: Mattson-Svard-Nordstrom
11see [1] below  10Mattson, Svaerd and Nordstroem type
11Kreiss-Oliger  10Kreiss-Oliger modified near the boundaries

11

Mattson-Svard-Nordstrom

11epsdis  10Scope: restrictedx01  REAL

11Description: Dissipation strength

11 Range 10 x01 Default: 0.2
11*:*  10Values typical between 0 and 1

11

11hxscaling  10Scope: restrictedx01  REAL

11Description: Scaling factor for the local grid spacing in the dissipation operators

11 Range 10 x01 Default: 1.0
110:*  10Positive please

11

11normxtype  10Scope: restrictedx01  KEYWORD

11Description: Type of norm

11 Range 10 x01 Default: Diagonal
11Diagonal  10Diagonal norm
11Full restricted  10Full restricted norm

11

11onesidedxinterpatchxboundaries 10Scope: restrictedx01  BOOLEAN

11Description: Evaluate derivatives near the local grid boundary if it is an inter-patch boundary

11                         10  x01  Default: yes

11

11onesidedxouterxboundaries 10Scope: restrictedx01  BOOLEAN

11Description: Evaluate derivatives within ghost zones of the outer boundary

11                         10  x01  Default: yes

11

11operatorxtype  10Scope: restrictedx01  KEYWORD

11Description: Type of operator

11 Range 10 x01 Default: Optimized
11Minimal Bandwidth  10Minimal bandwidth (except for 8-4 which is minimal spectral radius)
11Optimized  10Optimized for performance

11

11order  10Scope: restrictedx01  INT

11Description: Order of accuracy

11 Range 10 x01 Default: 2
112:8:2  10

11

11poisonxderivatives  10Scope: restrictedx01  BOOLEAN

11Description: Should we poison Dvar at boundaryxshiftout perimeter when taking derivatives

11                         10  x01  Default: no

11

11poisonxdissipation  10Scope: restrictedx01  BOOLEAN

11Description: Should we poison rhs at boundaryxshiftout perimeter when applying dissipation

11                         10  x01  Default: no

11

11poisonxvalue  10Scope: restrictedx01  REAL

11Description: Degree of intoxication

11 Range 10 x01 Default: 666.0
11*:*  10Anything you want

11sbpx1stxderiv  10Scope: restrictedx01  BOOLEAN

11Description: Should the 1st derivative operator be SBP

11                         10  x01  Default: yes

11

11sbpx2ndxderiv  10Scope: restrictedx01  BOOLEAN

11Description: Should the 2nd derivative operator be SBP

11                         10  x01  Default: yes

11

11sbpxupwindxderiv  10Scope: restrictedx01  BOOLEAN

11Description: Should the upwind derivative operator be SBP

11                         10  x01  Default: yes

11

11scalexwithxh  10Scope: restrictedx01  BOOLEAN

11Description: Should we scale the dissipation with the grid spacing h

11                         10  x01  Default: no

11

11usexdissipation  10Scope: restrictedx01  BOOLEAN

11Description: Should we add dissipation

11                         10  x01  Default: no

11

11usexshiftout  10Scope: restrictedx01  BOOLEAN

11Description: Should we use the boundaryxshiftxout parameters from CoordBase to shift the stencils of derivatives and dissipation

11                         10  x01  Default: no

11

11usexvariablexdeltas  10Scope: restrictedx01  BOOLEAN

11Description: Use extra grid functions to allow for variable delta’s in the dissipation operators

11                         10  x01  Default: no

11

11vars  10Scope: restrictedx01  STRING

11Description: List of evolved grid functions that should have dissipation added

11 Range 10 x01 Default: (none)
11.*  10Must be a valid list of grid functions

11

11zeroxderivsxy  10Scope: restrictedx01  BOOLEAN

11Description: set all derivatives to 0 in the y-direction

11                         10  x01  Default: no

11

11zeroxderivsxz  10Scope: restrictedx01  BOOLEAN

11Description: set all derivatives to 0 in the z-direction

11                         10  x01  Default: no

11

7 Interfaces General Implements: summationbyparts

Grid Variables

7.0.1 PUBLIC GROUPS

11x Group Names x10x Variable Names xx01 Details x11x
11  normmask  10  x01  compact110
11                 10  nmask  x01 description11Mask for the norm calculation
11                 10  x01 dimensions113
11                 10  x01 distribution11DEFAULT
11                 10  x01 group type11GF
11 10 x01 tags11tensortypealias=”scalar” Prolongation=”None” checkpoint=”no”
11                 10  x01 timelevels111
11                 10  x01variable type11REAL
11  deltas  10  x01  compact110
11                 10  sbpxdx  x01 description11Dissipation deltas
11                 10  sbpxdy  x01 dimensions113
11                 10  sbpxdz  x01 distribution11DEFAULT
11                 10  x01 group type11GF
11 10 x01 tags11tensortypealias=”U” Prolongation=”None”
11                 10  x01 timelevels111
11                 10  x01variable type11REAL

11

Provides: Diffxgf to Diffxgv to Diffxupxgv to Diff2xgv to Diffxcoeff to Diffxupxcoeff to Diff2xcoeff to GetScalProdDiag to GetScalProdCoeff to GetLshIndexRanges to GetBoundWidth to

8 Schedule This section lists all the variables which are assigned storage by thorn CactusNumerical/SummationByParts. Storage can either last for the duration of the run (Always means that if this thorn is activated storage will be assigned, Conditional means that if this thorn is activated storage will be assigned for the duration of the run if some condition is met), or can be turned on for the duration of a schedule function.

Storage x

11Always: 11Conditional:
11normmask11deltas

11x  11x

11

Scheduled Functions

CCTKxBASEGRID (conditional) x sbpxsetnormmask x setup the mask for the calculation of the norm x

11x11Language:11fortran
11 x 11Type: 11function
11 x 11Writes: 11summationbyparts::nmask(everywhere)

11

CCTKxBASEGRID (conditional) x sbpxdeltainitial x initialize dissipation deltas x

11x11Language:11fortran
11 x 11Type: 11function
11 x 11Writes: 11summationbyparts::deltas(everywhere)

11

CCTKxPOSTINITIAL (conditional) x sbpxcheckgridsizes x check grid sizes and ghost zones x

11x11Language:11fortran
11 x 11Type: 11function

11

MoLxPostRHS (conditional) x sbpxdissipationadd x add sbp compatible dissipation to the right hand sides

11x11Language:11c
11 x 11Reads: 11summationbyparts::deltas(everywhere)
11 x 11Type: 11function