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 x0xN of a computational domain x [a,b] with gridspacing h a one dimensional finite difference operator approximation to a first derivative, D, is said to satisfy summation by parts (SBP) with respect to a scalar product (defined by its coefficients σij)

u,vh = h i=0Nu ivjσij (1)

if the property

u,Dvh + Du,vh = (uv)ab (2)

is satisfied for all possible gridfunctions u and v.

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 figure 1 for the structure).


σ = x x x x 1 , σ = x x x x x xxxx x x x x xxxx 1 , σ = xxxx x x x x xxxx x x x x 1


Figure 1: 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 s and the width of the boundary region1 by r.

For the diagonal norm case it turns out that with r = 2τ it is possible to find SBP operators with s = 2τ (at least when τ 4), i.e. the order of accuracy at the boundary is half the order in the interior. For the restricted full norm case τ = s 1 when r = τ + 2 and for the full norm case τ = s 1 when r = τ + 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.

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

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:

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

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 finite difference coefficients 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 )

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×nsize that on return will contain the coefficients for the derivatives. The table_handle 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:

  call Diff_gv (cctkGH, 0, f, dxf, -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 )

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 José M. Martín-García for very kindly providing us with a very well designed mathematica script to calculate the finite difference and scalar product coefficients.

References

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

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

[3]   Ken Mattsson, Magnus Svärd and Jan Nordström, 2004, Journal of Scientific Computing, 21, 57–79.

6 Parameters




check_grid_sizes
Scope: restricted  BOOLEAN



Description: Should we check grid sizes and ghost zones



  Default: yes






diss_fraction
Scope: restricted  REAL



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



Range   Default: 0.2
0:0.5






dissipation_type
Scope: restricted  KEYWORD



Description: Type of dissipation operator



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



[1]

Mattson-Svard-Nordstrom




epsdis
Scope: restricted  REAL



Description: Dissipation strength



Range   Default: 0.2
*:*
Values typical between 0 and 1






h_scaling
Scope: restricted  REAL



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



Range   Default: 1.0
0:*
Positive please






norm_type
Scope: restricted  KEYWORD



Description: Type of norm



Range   Default: Diagonal
Diagonal
Diagonal norm
Full restricted
Full restricted norm






onesided_interpatch_boundaries
Scope: restricted  BOOLEAN



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



  Default: yes






onesided_outer_boundaries
Scope: restricted  BOOLEAN



Description: Evaluate derivatives within ghost zones of the outer boundary



  Default: yes






operator_type
Scope: restricted  KEYWORD



Description: Type of operator



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






order
Scope: restricted  INT



Description: Order of accuracy



Range   Default: 2
2:8:2






poison_derivatives
Scope: restricted  BOOLEAN



Description: Should we poison Dvar at boundary_shiftout perimeter when taking derivatives



  Default: no






poison_dissipation
Scope: restricted  BOOLEAN



Description: Should we poison rhs at boundary_shiftout perimeter when applying dissipation



  Default: no






poison_value
Scope: restricted  REAL



Description: Degree of intoxication



Range   Default: 666.0
*:*
Anything you want






sbp_1st_deriv
Scope: restricted  BOOLEAN



Description: Should the 1st derivative operator be SBP



  Default: yes






sbp_2nd_deriv
Scope: restricted  BOOLEAN



Description: Should the 2nd derivative operator be SBP



  Default: yes






sbp_upwind_deriv
Scope: restricted  BOOLEAN



Description: Should the upwind derivative operator be SBP



  Default: yes






scale_with_h
Scope: restricted  BOOLEAN



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



  Default: no






use_dissipation
Scope: restricted  BOOLEAN



Description: Should we add dissipation



  Default: no






use_shiftout
Scope: restricted  BOOLEAN



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



  Default: no






use_variable_deltas
Scope: restricted  BOOLEAN



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



  Default: no






vars
Scope: restricted  STRING



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



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






zero_derivs_y
Scope: restricted  BOOLEAN



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



  Default: no






zero_derivs_z
Scope: restricted  BOOLEAN



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



  Default: no



7 Interfaces

General

Implements:

summationbyparts

Grid Variables

7.0.1 PUBLIC GROUPS




  Group Names    Variable Names    Details   




normmask   compact0
nmask   descriptionMask for the norm calculation
  dimensions3
  distributionDEFAULT
  group typeGF
  tagstensortypealias=”scalar” Prolongation=”None” checkpoint=”no”
  timelevels1
 variable typeREAL




deltas   compact0
sbp_dx   descriptionDissipation deltas
sbp_dy   dimensions3
sbp_dz   distributionDEFAULT
  group typeGF
  tagstensortypealias=”U” Prolongation=”None”
  timelevels1
 variable typeREAL




Provides:

Diff_gf to

Diff_gv to

Diff_up_gv to

Diff2_gv to

Diff_coeff to

Diff_up_coeff to

Diff2_coeff 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

 

Always: Conditional:
normmaskdeltas
   

Scheduled Functions

CCTK_BASEGRID (conditional)

  sbp_setnormmask

  setup the mask for the calculation of the norm

 

 Language:fortran
 Type: function

CCTK_BASEGRID (conditional)

  sbp_deltainitial

  initialize dissipation deltas

 

 Language:fortran
 Type: function

CCTK_POSTINITIAL (conditional)

  sbp_checkgridsizes

  check grid sizes and ghost zones

 

 Language:fortran
 Type: function

MoL_PostRHS (conditional)

  sbp_dissipationadd

  add sbp compatible dissipation to the right hand sides

 

 Language:c
 Type: function