Implementation of 6th order dissipation operator compatible with the 6-5

finite difference operators. Not fully tested yet. The MoL dissipation
interface was shamelessly stolen from Eriks Dissipation thorn (in AEIThorns).

Things to do:

Add a test for the presence of MoLQueryEvolvedRHS (since the use of
dissipation is optional, I didn't want to require the presence of
this routine).

Add compatible dissipation operators to the other finite difference operators.

More testing.


git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@21 f69c4107-0314-4c4f-9ad4-17e986b73f4a

6 files changed, 688 insertions, 2 deletions

diff --git a/interface.ccl b/interface.ccl
index b8fd0f9..8b4304d 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -46,7 +46,8 @@ REQUIRES FUNCTION SymmetryTableHandleForGrid
 CCTK_INT FUNCTION SymmetryHandleOfName (CCTK_STRING IN sym_name)
 REQUIRES FUNCTION SymmetryHandleOfName
 
-
+CCTK_INT FUNCTION MoLQueryEvolvedRHS (CCTK_INT IN EvolvedIndex)
+USES FUNCTION MoLQueryEvolvedRHS
 
 public
 
diff --git a/param.ccl b/param.ccl
index bcfa006..fc70283 100644
--- a/param.ccl
+++ b/param.ccl
@@ -11,3 +11,17 @@ INT order "Order of accuracy" STEERABLE=always
 {
   2:8:2 :: ""
 } 2
+
+BOOLEAN use_dissipation "Should we add dissipation"
+{
+} "no"
+
+REAL epsdis "Dissipation strength" STEERABLE=always
+{
+  *:* :: "Values typical between 0 and 1"
+} 0.2
+
+STRING vars "List of evolved grid functions that should have dissipation added" STEERABLE=always
+{
+  .* :: "Must be a valid list of grid functions"
+} ""
diff --git a/schedule.ccl b/schedule.ccl
index 7575d69..77148a8 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -13,3 +13,8 @@ SCHEDULE SBP_CheckGridSizes AT postinitial
 {
   LANG: Fortran
 } "Check grid sizes and ghost zones"
+
+SCHEDULE SBP_DissipationAdd IN MoL_PostRHS
+{
+  LANG: C
+} "Add SBP compatible dissipation to the right hand sides"
diff --git a/src/Dissipation_6_5.F90 b/src/Dissipation_6_5.F90
new file mode 100644
index 0000000..34f39f1
--- /dev/null
+++ b/src/Dissipation_6_5.F90
@@ -0,0 +1,547 @@
+! $Header$
+
+#include "cctk.h"
+
+subroutine dissipation_6_5 (var, ni, nj, nk, bb, gsize, delta, epsdis, rhs)
+
+  implicit none
+
+  integer ::  ni, nj, nk
+  CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var
+  CCTK_REAL, dimension(ni,nj,nk), intent(inout) :: rhs
+  CCTK_INT, dimension(6), intent(in) :: bb
+  CCTK_INT, dimension(3), intent(in) :: gsize
+  CCTK_REAL, dimension(3), intent(in) :: delta
+  CCTK_REAL, intent(in) :: epsdis 
+
+  CCTK_REAL :: zero = 0.0
+  integer, parameter :: wp = kind(zero)
+  integer :: i, j, k
+  CCTK_REAL, dimension(4), save :: ac = (/ -20.0_wp, 15.0_wp, -6.0_wp, 1.0_wp /)
+  CCTK_REAL, dimension(10,7) :: a
+  CCTK_REAL :: idel
+
+  CCTK_INT :: il, ir, jl, jr, kl, kr
+
+!  print*,"dissipation_6_5 called"
+!  print*,ni,nj,nk
+!  print*,bb
+!  print*,gsize
+!  print*,delta
+!  print*,epsdis
+!  stop
+  call set_coeff ( a, delta(1) )
+
+  idel = epsdis / ( 64 * delta(1) )
+
+  if ( bb(1) == 0 ) then
+    il = 1 + gsize(1)
+  else
+    rhs(1,:,:) = rhs(1,:,:) + &
+                 ( a(1,1) * var(1,:,:) + a(2,1) * var(2,:,:) + &
+                   a(3,1) * var(3,:,:) + a(4,1) * var(4,:,:) ) * idel
+    rhs(2,:,:) = rhs(2,:,:) + &
+                 ( a(1,2) * var(1,:,:) + a(2,2) * var(2,:,:) + &
+                   a(3,2) * var(3,:,:) + a(4,2) * var(4,:,:) + &
+                   a(5,2) * var(5,:,:) + a(6,2) * var(6,:,:) + &
+                   a(7,2) * var(7,:,:) + a(8,2) * var(8,:,:) ) * idel
+    rhs(3,:,:) = rhs(3,:,:) + &
+                 ( a(1,3) * var(1,:,:) + a(2,3) * var(2,:,:) + &
+                   a(3,3) * var(3,:,:) + a(4,3) * var(4,:,:) + &
+                   a(5,3) * var(5,:,:) + a(6,3) * var(6,:,:) + &
+                   a(7,3) * var(7,:,:) + a(8,3) * var(8,:,:) + &
+                   a(9,3) * var(9,:,:) ) * idel
+    rhs(4,:,:) = rhs(4,:,:) + &
+                 ( a(1,4) * var(1,:,:) + a(2,4) * var(2,:,:) + &
+                   a(3,4) * var(3,:,:) + a(4,4) * var(4,:,:) + &
+                   a(5,4) * var(5,:,:) + a(6,4) * var(6,:,:) + &
+                   a(7,4) * var(7,:,:) + a(8,4) * var(8,:,:) + &
+                   a(9,4) * var(9,:,:) + a(10,4) * var(10,:,:) ) * idel
+    rhs(5,:,:) = rhs(5,:,:) + &
+                 ( a(1,5) * var(1,:,:) + a(2,5) * var(2,:,:) + &
+                   a(3,5) * var(3,:,:) + a(4,5) * var(4,:,:) + &
+                   a(5,5) * var(5,:,:) + a(6,5) * var(6,:,:) + &
+                   a(7,5) * var(7,:,:) + a(8,5) * var(8,:,:) + &
+                   a(9,5) * var(9,:,:) + a(10,5) * var(10,:,:) ) * idel
+    rhs(6,:,:) = rhs(6,:,:) + &
+                 ( a(1,6) * var(1,:,:) + a(2,6) * var(2,:,:) + &
+                   a(3,6) * var(3,:,:) + a(4,6) * var(4,:,:) + &
+                   a(5,6) * var(5,:,:) + a(6,6) * var(6,:,:) + &
+                   a(7,6) * var(7,:,:) + a(8,6) * var(8,:,:) + &
+                   a(9,6) * var(9,:,:) + a(10,6) * var(10,:,:) ) * idel
+    rhs(7,:,:) = rhs(7,:,:) + &
+                 ( a(1,7) * var(1,:,:) + a(2,7) * var(2,:,:) + &
+                   a(3,7) * var(3,:,:) + a(4,7) * var(4,:,:) + &
+                   a(5,7) * var(5,:,:) + a(6,7) * var(6,:,:) + &
+                   a(7,7) * var(7,:,:) + a(8,7) * var(8,:,:) + &
+                   a(9,7) * var(9,:,:) + a(10,7) * var(10,:,:) ) * idel
+    il = 8
+  end if
+  if ( bb(2) == 0 ) then
+    ir = ni - gsize(1)
+  else
+    rhs(ni-6,:,:) = rhs(ni-6,:,:) + &
+                    ( a(1,7) * var(ni,:,:) + a(2,7) * var(ni-1,:,:) + &
+                      a(3,7) * var(ni-2,:,:) + a(4,7) * var(ni-3,:,:) + &
+                      a(5,7) * var(ni-4,:,:) + a(6,7) * var(ni-5,:,:) + &
+                      a(7,7) * var(ni-6,:,:) + a(8,7) * var(ni-7,:,:) + &
+                      a(9,7) * var(ni-8,:,:) + a(10,7) * var(ni-9,:,:) ) * idel
+    rhs(ni-5,:,:) = rhs(ni-5,:,:) + &
+                    ( a(1,6) * var(ni,:,:) + a(2,6) * var(ni-1,:,:) + &
+                      a(3,6) * var(ni-2,:,:) + a(4,6) * var(ni-3,:,:) + &
+                      a(5,6) * var(ni-4,:,:) + a(6,6) * var(ni-5,:,:) + &
+                      a(7,6) * var(ni-6,:,:) + a(8,6) * var(ni-7,:,:) + &
+                      a(9,6) * var(ni-8,:,:) + a(10,6) * var(ni-9,:,:) ) * idel
+    rhs(ni-4,:,:) = rhs(ni-4,:,:) + &
+                    ( a(1,5) * var(ni,:,:) + a(2,5) * var(ni-1,:,:) + &
+                      a(3,5) * var(ni-2,:,:) + a(4,5) * var(ni-3,:,:) + &
+                      a(5,5) * var(ni-4,:,:) + a(6,5) * var(ni-5,:,:) + &
+                      a(7,5) * var(ni-6,:,:) + a(8,5) * var(ni-7,:,:) + &
+                      a(9,5) * var(ni-8,:,:) + a(10,5) * var(ni-9,:,:) ) * idel
+    rhs(ni-3,:,:) = rhs(ni-3,:,:) + &
+                    ( a(1,4) * var(ni,:,:) + a(2,4) * var(ni-1,:,:) + &
+                      a(3,4) * var(ni-2,:,:) + a(4,4) * var(ni-3,:,:) + &
+                      a(5,4) * var(ni-4,:,:) + a(6,4) * var(ni-5,:,:) + &
+                      a(7,4) * var(ni-6,:,:) + a(8,4) * var(ni-7,:,:) + &
+                      a(9,4) * var(ni-8,:,:) + a(10,4) * var(ni-9,:,:) ) * idel
+    rhs(ni-2,:,:) = rhs(ni-2,:,:) + &
+                    ( a(1,3) * var(ni,:,:) + a(2,3) * var(ni-1,:,:) + &
+                      a(3,3) * var(ni-2,:,:) + a(4,3) * var(ni-3,:,:) + &
+                      a(5,3) * var(ni-4,:,:) + a(6,3) * var(ni-5,:,:) + &
+                      a(7,3) * var(ni-6,:,:) + a(8,3) * var(ni-7,:,:) + &
+                      a(9,3) * var(ni-8,:,:) ) * idel
+    rhs(ni-1,:,:) = rhs(ni-1,:,:) + &
+                    ( a(1,2) * var(ni,:,:) + a(2,2) * var(ni-1,:,:) + &
+                      a(3,2) * var(ni-2,:,:) + a(4,2) * var(ni-3,:,:) + &
+                      a(5,2) * var(ni-4,:,:) + a(6,2) * var(ni-5,:,:) + &
+                      a(7,2) * var(ni-6,:,:) + a(8,2) * var(ni-7,:,:) ) * idel
+    rhs(ni,:,:)   = rhs(ni,:,:) + &
+                    ( a(1,1) * var(ni,:,:) + a(2,1) * var(ni-1,:,:) + &
+                      a(3,1) * var(ni-2,:,:) + a(4,1) * var(ni-3,:,:) ) * idel
+    ir = ni - 7
+  end if
+  rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
+                   ( -20.0_wp * var(il:ir,:,:) + &
+                      15.0_wp * ( var(il-1:ir-1,:,:) + &
+                                  var(il+1:ir+1,:,:) ) - &
+                       6.0_wp * ( var(il-2:ir-2,:,:) + &
+                                  var(il+2:ir+2,:,:) ) + &
+                                ( var(il-3:ir-3,:,:) + &
+                                  var(il+3:ir+3,:,:) ) ) * idel
+
+  call set_coeff ( a, delta(2) )
+
+  idel = epsdis / ( 64 * delta(2) )
+
+  if ( bb(3) == 0 ) then
+    jl = 1 + gsize(2)
+  else
+    rhs(:,1,:) = rhs(:,1,:) + &
+                 ( a(1,1) * var(:,1,:) + a(2,1) * var(:,2,:) + &
+                   a(3,1) * var(:,3,:) + a(4,1) * var(:,4,:) ) * idel
+    rhs(:,2,:) = rhs(:,2,:) + &
+                 ( a(1,2) * var(:,1,:) + a(2,2) * var(:,2,:) + &
+                   a(3,2) * var(:,3,:) + a(4,2) * var(:,4,:) + &
+                   a(5,2) * var(:,5,:) + a(6,2) * var(:,6,:) + &
+                   a(7,2) * var(:,7,:) + a(8,2) * var(:,8,:) ) * idel
+    rhs(:,3,:) = rhs(:,3,:) + &
+                 ( a(1,3) * var(:,1,:) + a(2,3) * var(:,2,:) + &
+                   a(3,3) * var(:,3,:) + a(4,3) * var(:,4,:) + &
+                   a(5,3) * var(:,5,:) + a(6,3) * var(:,6,:) + &
+                   a(7,3) * var(:,7,:) + a(8,3) * var(:,8,:) + &
+                   a(9,3) * var(:,9,:) ) * idel
+    rhs(:,4,:) = rhs(:,4,:) + &
+                 ( a(1,4) * var(:,1,:) + a(2,4) * var(:,2,:) + &
+                   a(3,4) * var(:,3,:) + a(4,4) * var(:,4,:) + &
+                   a(5,4) * var(:,5,:) + a(6,4) * var(:,6,:) + &
+                   a(7,4) * var(:,7,:) + a(8,4) * var(:,8,:) + &
+                   a(9,4) * var(:,9,:) + a(10,4) * var(:,10,:) ) * idel
+    rhs(:,5,:) = rhs(:,5,:) + &
+                 ( a(1,5) * var(:,1,:) + a(2,5) * var(:,2,:) + &
+                   a(3,5) * var(:,3,:) + a(4,5) * var(:,4,:) + &
+                   a(5,5) * var(:,5,:) + a(6,5) * var(:,6,:) + &
+                   a(7,5) * var(:,7,:) + a(8,5) * var(:,8,:) + &
+                   a(9,5) * var(:,9,:) + a(10,5) * var(:,10,:) ) * idel
+    rhs(:,6,:) = rhs(:,6,:) + &
+                 ( a(1,6) * var(:,1,:) + a(2,6) * var(:,2,:) + &
+                   a(3,6) * var(:,3,:) + a(4,6) * var(:,4,:) + &
+                   a(5,6) * var(:,5,:) + a(6,6) * var(:,6,:) + &
+                   a(7,6) * var(:,7,:) + a(8,6) * var(:,8,:) + &
+                   a(9,6) * var(:,9,:) + a(10,6) * var(:,10,:) ) * idel
+    rhs(:,7,:) = rhs(:,7,:) + &
+                 ( a(1,7) * var(:,1,:) + a(2,7) * var(:,2,:) + &
+                   a(3,7) * var(:,3,:) + a(4,7) * var(:,4,:) + &
+                   a(5,7) * var(:,5,:) + a(6,7) * var(:,6,:) + &
+                   a(7,7) * var(:,7,:) + a(8,7) * var(:,8,:) + &
+                   a(9,7) * var(:,9,:) + a(10,7) * var(:,10,:) ) * idel
+    jl = 8
+  end if
+  if ( bb(4) == 0 ) then
+    jr = nj - gsize(2)
+  else
+    rhs(:,nj-6,:) = rhs(:,nj-6,:) + &
+                    ( a(1,7) * var(:,nj,:) + a(2,7) * var(:,nj-1,:) + &
+                      a(3,7) * var(:,nj-2,:) + a(4,7) * var(:,nj-3,:) + &
+                      a(5,7) * var(:,nj-4,:) + a(6,7) * var(:,nj-5,:) + &
+                      a(7,7) * var(:,nj-6,:) + a(8,7) * var(:,nj-7,:) + &
+                      a(9,7) * var(:,nj-8,:) + a(10,7) * var(:,nj-9,:) ) * idel
+    rhs(:,nj-5,:) = rhs(:,nj-5,:) + &
+                    ( a(1,6) * var(:,nj,:) + a(2,6) * var(:,nj-1,:) + &
+                      a(3,6) * var(:,nj-2,:) + a(4,6) * var(:,nj-3,:) + &
+                      a(5,6) * var(:,nj-4,:) + a(6,6) * var(:,nj-5,:) + &
+                      a(7,6) * var(:,nj-6,:) + a(8,6) * var(:,nj-7,:) + &
+                      a(9,6) * var(:,nj-8,:) + a(10,6) * var(:,nj-9,:) ) * idel
+    rhs(:,nj-4,:) = rhs(:,nj-4,:) + &
+                    ( a(1,5) * var(:,nj,:) + a(2,5) * var(:,nj-1,:) + &
+                      a(3,5) * var(:,nj-2,:) + a(4,5) * var(:,nj-3,:) + &
+                      a(5,5) * var(:,nj-4,:) + a(6,5) * var(:,nj-5,:) + &
+                      a(7,5) * var(:,nj-6,:) + a(8,5) * var(:,nj-7,:) + &
+                      a(9,5) * var(:,nj-8,:) + a(10,5) * var(:,nj-9,:) ) * idel
+    rhs(:,nj-3,:) = rhs(:,nj-3,:) + &
+                    ( a(1,4) * var(:,nj,:) + a(2,4) * var(:,nj-1,:) + &
+                      a(3,4) * var(:,nj-2,:) + a(4,4) * var(:,nj-3,:) + &
+                      a(5,4) * var(:,nj-4,:) + a(6,4) * var(:,nj-5,:) + &
+                      a(7,4) * var(:,nj-6,:) + a(8,4) * var(:,nj-7,:) + &
+                      a(9,4) * var(:,nj-8,:) + a(10,4) * var(:,nj-9,:) ) * idel
+    rhs(:,nj-2,:) = rhs(:,nj-2,:) + &
+                    ( a(1,3) * var(:,nj,:) + a(2,3) * var(:,nj-1,:) + &
+                      a(3,3) * var(:,nj-2,:) + a(4,3) * var(:,nj-3,:) + &
+                      a(5,3) * var(:,nj-4,:) + a(6,3) * var(:,nj-5,:) + &
+                      a(7,3) * var(:,nj-6,:) + a(8,3) * var(:,nj-7,:) + &
+                      a(9,3) * var(:,nj-8,:) ) * idel
+    rhs(:,nj-1,:) = rhs(:,nj-1,:) + &
+                    ( a(1,2) * var(:,nj,:) + a(2,2) * var(:,nj-1,:) + &
+                      a(3,2) * var(:,nj-2,:) + a(4,2) * var(:,nj-3,:) + &
+                      a(5,2) * var(:,nj-4,:) + a(6,2) * var(:,nj-5,:) + &
+                      a(7,2) * var(:,nj-6,:) + a(8,2) * var(:,nj-7,:) ) * idel
+    rhs(:,nj,:)   = rhs(:,nj,:) + &
+                    ( a(1,1) * var(:,nj,:) + a(2,1) * var(:,nj-1,:) + &
+                      a(3,1) * var(:,nj-2,:) + a(4,1) * var(:,nj-3,:) ) * idel
+    jr = nj - 7
+  end if
+  rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
+                   ( -20.0_wp * var(:,jl:jr,:) + &
+                      15.0_wp * ( var(:,jl-1:jr-1,:) + &
+                                  var(:,jl+1:jr+1,:) ) - &
+                       6.0_wp * ( var(:,jl-2:jr-2,:) + &
+                                  var(:,jl+2:jr+2,:) ) + &
+                                ( var(:,jl-3:jr-3,:) + &
+                                  var(:,jl+3:jr+3,:) ) ) * idel
+
+  call set_coeff ( a, delta(3) )
+
+  idel = epsdis / ( 64 * delta(3) )
+
+  if ( bb(5) == 0 ) then
+    kl = 1 + gsize(3)
+  else
+    rhs(:,:,1) = rhs(:,:,1) + &
+                 ( a(1,1) * var(:,:,1) + a(2,1) * var(:,:,2) + &
+                   a(3,1) * var(:,:,3) + a(4,1) * var(:,:,4) ) * idel
+    rhs(:,:,2) = rhs(:,:,2) + &
+                 ( a(1,2) * var(:,:,1) + a(2,2) * var(:,:,2) + &
+                   a(3,2) * var(:,:,3) + a(4,2) * var(:,:,4) + &
+                   a(5,2) * var(:,:,5) + a(6,2) * var(:,:,6) + &
+                   a(7,2) * var(:,:,7) + a(8,2) * var(:,:,8) ) * idel
+    rhs(:,:,3) = rhs(:,:,3) + &
+                 ( a(1,3) * var(:,:,1) + a(2,3) * var(:,:,2) + &
+                   a(3,3) * var(:,:,3) + a(4,3) * var(:,:,4) + &
+                   a(5,3) * var(:,:,5) + a(6,3) * var(:,:,6) + &
+                   a(7,3) * var(:,:,7) + a(8,3) * var(:,:,8) + &
+                   a(9,3) * var(:,:,9) ) * idel
+    rhs(:,:,4) = rhs(:,:,4) + &
+                 ( a(1,4) * var(:,:,1) + a(2,4) * var(:,:,2) + &
+                   a(3,4) * var(:,:,3) + a(4,4) * var(:,:,4) + &
+                   a(5,4) * var(:,:,5) + a(6,4) * var(:,:,6) + &
+                   a(7,4) * var(:,:,7) + a(8,4) * var(:,:,8) + &
+                   a(9,4) * var(:,:,9) + a(10,4) * var(:,:,10) ) * idel
+    rhs(:,:,5) = rhs(:,:,5) + &
+                 ( a(1,5) * var(:,:,1) + a(2,5) * var(:,:,2) + &
+                   a(3,5) * var(:,:,3) + a(4,5) * var(:,:,4) + &
+                   a(5,5) * var(:,:,5) + a(6,5) * var(:,:,6) + &
+                   a(7,5) * var(:,:,7) + a(8,5) * var(:,:,8) + &
+                   a(9,5) * var(:,:,9) + a(10,5) * var(:,:,10) ) * idel
+    rhs(:,:,6) = rhs(:,:,6) + &
+                 ( a(1,6) * var(:,:,1) + a(2,6) * var(:,:,2) + &
+                   a(3,6) * var(:,:,3) + a(4,6) * var(:,:,4) + &
+                   a(5,6) * var(:,:,5) + a(6,6) * var(:,:,6) + &
+                   a(7,6) * var(:,:,7) + a(8,6) * var(:,:,8) + &
+                   a(9,6) * var(:,:,9) + a(10,6) * var(:,:,10) ) * idel
+    rhs(:,:,7) = rhs(:,:,7) + &
+                 ( a(1,7) * var(:,:,1) + a(2,7) * var(:,:,2) + &
+                   a(3,7) * var(:,:,3) + a(4,7) * var(:,:,4) + &
+                   a(5,7) * var(:,:,5) + a(6,7) * var(:,:,6) + &
+                   a(7,7) * var(:,:,7) + a(8,7) * var(:,:,8) + &
+                   a(9,7) * var(:,:,9) + a(10,7) * var(:,:,10) ) * idel
+    kl = 8
+  end if
+  if ( bb(6) == 0 ) then
+    kr = nk - gsize(3)
+  else
+    rhs(:,:,nk-6) = rhs(:,:,nk-6) + &
+                    ( a(1,7) * var(:,:,nk) + a(2,7) * var(:,:,nk-1) + &
+                      a(3,7) * var(:,:,nk-2) + a(4,7) * var(:,:,nk-3) + &
+                      a(5,7) * var(:,:,nk-4) + a(6,7) * var(:,:,nk-5) + &
+                      a(7,7) * var(:,:,nk-6) + a(8,7) * var(:,:,nk-7) + &
+                      a(9,7) * var(:,:,nk-8) + a(10,7) * var(:,:,nk-9) ) * idel
+    rhs(:,:,nk-5) = rhs(:,:,nk-5) + &
+                    ( a(1,6) * var(:,:,nk) + a(2,6) * var(:,:,nk-1) + &
+                      a(3,6) * var(:,:,nk-2) + a(4,6) * var(:,:,nk-3) + &
+                      a(5,6) * var(:,:,nk-4) + a(6,6) * var(:,:,nk-5) + &
+                      a(7,6) * var(:,:,nk-6) + a(8,6) * var(:,:,nk-7) + &
+                      a(9,6) * var(:,:,nk-8) + a(10,6) * var(:,:,nk-9) ) * idel
+    rhs(:,:,nk-4) = rhs(:,:,nk-4) + &
+                    ( a(1,5) * var(:,:,nk) + a(2,5) * var(:,:,nk-1) + &
+                      a(3,5) * var(:,:,nk-2) + a(4,5) * var(:,:,nk-3) + &
+                      a(5,5) * var(:,:,nk-4) + a(6,5) * var(:,:,nk-5) + &
+                      a(7,5) * var(:,:,nk-6) + a(8,5) * var(:,:,nk-7) + &
+                      a(9,5) * var(:,:,nk-8) + a(10,5) * var(:,:,nk-9) ) * idel
+    rhs(:,:,nk-3) = rhs(:,:,nk-3) + &
+                    ( a(1,4) * var(:,:,nk) + a(2,4) * var(:,:,nk-1) + &
+                      a(3,4) * var(:,:,nk-2) + a(4,4) * var(:,:,nk-3) + &
+                      a(5,4) * var(:,:,nk-4) + a(6,4) * var(:,:,nk-5) + &
+                      a(7,4) * var(:,:,nk-6) + a(8,4) * var(:,:,nk-7) + &
+                      a(9,4) * var(:,:,nk-8) + a(10,4) * var(:,:,nk-9) ) * idel
+    rhs(:,:,nk-2) = rhs(:,:,nk-2) + &
+                    ( a(1,3) * var(:,:,nk) + a(2,3) * var(:,:,nk-1) + &
+                      a(3,3) * var(:,:,nk-2) + a(4,3) * var(:,:,nk-3) + &
+                      a(5,3) * var(:,:,nk-4) + a(6,3) * var(:,:,nk-5) + &
+                      a(7,3) * var(:,:,nk-6) + a(8,3) * var(:,:,nk-7) + &
+                      a(9,3) * var(:,:,nk-8) ) * idel
+    rhs(:,:,nk-1) = rhs(:,:,nk-1) + &
+                    ( a(1,2) * var(:,:,nk) + a(2,2) * var(:,:,nk-1) + &
+                      a(3,2) * var(:,:,nk-2) + a(4,2) * var(:,:,nk-3) + &
+                      a(5,2) * var(:,:,nk-4) + a(6,2) * var(:,:,nk-5) + &
+                      a(7,2) * var(:,:,nk-6) + a(8,2) * var(:,:,nk-7) ) * idel
+    rhs(:,:,nk)   = rhs(:,:,nk) + &
+                    ( a(1,1) * var(:,:,nk) + a(2,1) * var(:,:,nk-1) + &
+                      a(3,1) * var(:,:,nk-2) + a(4,1) * var(:,:,nk-3) ) * idel
+    kr = nk - 7
+  end if
+  rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
+                   ( -20.0_wp * var(:,:,kl:kr) + &
+                      15.0_wp * ( var(:,:,kl-1:kr-1) + &
+                                  var(:,:,kl+1:kr+1) ) - &
+                       6.0_wp * ( var(:,:,kl-2:kr-2) + &
+                                  var(:,:,kl+2:kr+2) ) + &
+                                ( var(:,:,kl-3:kr-3) + &
+                                  var(:,:,kl+3:kr+3) ) ) * idel
+contains
+
+  subroutine set_coeff ( a, dx )
+
+    implicit none
+
+    CCTK_REAL, dimension(10,7), intent(OUT) :: a
+    CCTK_REAL, intent(in) :: dx
+    CCTK_REAL :: zero = 0.0
+    integer, parameter :: wp = kind(zero)
+    CCTK_REAL :: b1, b2, b3, b4, b5, b6
+
+    b1 = dx**2; b2 = (2+25*b1)/27; b3 = (7+20*b1)/27;
+    b4 = (1+b1)/2; b5 = (20+7*b1)/27; b6 = (25+2*b1)/27;
+
+    a(1,1) = -4.566666666666666666666666666666666666666666666667_wp*b1 &
+            - 4.566666666666666666666666666666666666666666666667_wp*b2 &
+            - 4.566666666666666666666666666666666666666666666667_wp*b3
+    a(2,1) = 13.7_wp*b1 + 13.7_wp*b2 + 13.7_wp*b3
+    a(3,1) = -13.7_wp*b1 - 13.7_wp*b2 - 13.7_wp*b3
+    a(4,1) = 4.566666666666666666666666666666666666666666666667_wp*b1 &
+           + 4.566666666666666666666666666666666666666666666667_wp*b2 &
+           + 4.566666666666666666666666666666666666666666666667_wp*b3
+    a(5,1) = zero
+    a(6,1) = zero
+    a(7,1) = zero
+    a(8,1) = zero
+    a(9,1) = zero
+    a(10,1) = zero
+    
+    a(1,2) = 0.797987233606384804466079266492184563887676357853_wp*b1 &
+           + 0.797987233606384804466079266492184563887676357853_wp*b2 &
+           + 0.797987233606384804466079266492184563887676357853_wp*b3
+    a(2,2) = -2.393961700819154413398237799476553691663029073559_wp*b1 &
+            - 2.393961700819154413398237799476553691663029073559_wp*b2 &
+            - 2.393961700819154413398237799476553691663029073559_wp*b3 &
+            + 0.307359194675104101198317748428244417750742420031_wp*b4
+    a(3,2) = 2.393961700819154413398237799476553691663029073559_wp*b1 &
+           + 2.393961700819154413398237799476553691663029073559_wp*b2 &
+           + 2.393961700819154413398237799476553691663029073559_wp*b3 &
+           - 0.922077584025312303594953245284733253252227260092_wp*b4 &
+           - 0.0336105345682464004522484486235834134456512678_wp*b5
+    a(4,2) = -0.797987233606384804466079266492184563887676357853_wp*b1 &
+            - 0.797987233606384804466079266492184563887676357853_wp*b2 &
+            - 0.797987233606384804466079266492184563887676357853_wp*b3 &
+            + 0.922077584025312303594953245284733253252227260092_wp*b4 &
+            + 0.1008316037047392013567453458707502403369538034_wp*b5 &
+            - 0.2929089778363370916823203340749080588489943555_wp*b6
+    a(5,2) = 0.09404587787607681409837743471424421680278209869_wp &
+           - 0.307359194675104101198317748428244417750742420031_wp*b4 &
+           - 0.1008316037047392013567453458707502403369538034_wp*b5 &
+           + 0.878726933509011275046961002224724176546983066499_wp*b6
+    a(6,2) = -0.282137633628230442295132304142732650409647604264_wp &
+            + 0.0336105345682464004522484486235834134456512678_wp*b5 &
+            - 0.878726933509011275046961002224724176546983066499_wp*b6
+    a(7,2) = 0.282137633628230442295132304142732650412591061208_wp &
+           + 0.2929089778363370916823203340749080588489943555_wp*b6
+    a(8,2) = -0.094045877876076814098377434714244216807708544941_wp
+    a(9,2) = zero
+    a(10,2) = zero
+    
+    a(1,3) = -0.570050756575490634640420716139511539225905452034_wp*b1 &
+            - 0.570050756575490634640420716139511539225905452034_wp*b2 &
+            - 0.570050756575490634640420716139511539225905452034_wp*b3
+    a(2,3) = 1.710152269726471903921262148418534617677716356102_wp*b1 &
+           + 1.710152269726471903921262148418534617677716356102_wp*b2 &
+           + 1.710152269726471903921262148418534617677716356102_wp*b3 &
+           + 0.65012922388697115376186819507780973585422870837_wp*b4
+    a(3,3) = -1.710152269726471903921262148418534617677716356102_wp*b1 &
+            - 1.710152269726471903921262148418534617677716356102_wp*b2 &
+            - 1.710152269726471903921262148418534617677716356102_wp*b3 &
+            - 1.950387671660913461285604585233429207562686125111_wp*b4 &
+            + 0.222968626602313512338688373037925009289627989014_wp*b5
+    a(4,3) = 0.570050756575490634640420716139511539225905452034_wp*b1 &
+           + 0.570050756575490634640420716139511539225905452034_wp*b2 &
+           + 0.570050756575490634640420716139511539225905452034_wp*b3 &
+           + 1.950387671660913461285604585233429207562686125111_wp*b4 &
+           - 0.668905879806940537016065119113775027868883967041_wp*b5 &
+           - 0.10042075802449997038657833546993392801822845339_wp*b6
+    a(5,3) = -0.223539235539512073697159351113375565174639048593_wp &
+            - 0.65012922388697115376186819507780973585422870837_wp*b4 &
+            + 0.668905879806940537016065119113775027868883967041_wp*b5 &
+            + 0.30126227407349991115973500640980178405468536017_wp*b6
+    a(6,3) = 0.74867234655369702626509041115338883998616371461_wp &
+           - 0.222968626602313512338688373037925009289627989014_wp*b5 &
+           - 0.30126227407349991115973500640980178405468536017_wp*b6
+    a(7,3) = -0.904781626424018636612315126779913128909937434049_wp &
+            + 0.10042075802449997038657833546993392801822845339_wp*b6
+    a(8,3) = 0.45770315534499448921799642455316199855922050042_wp
+    a(9,3) = -0.078054639935160805173612357813262144460088314167_wp
+    a(10,3) = zero
+    
+    a(1,4) = -0.432850298736358291696660164313057377022426239952_wp*b1 &
+            - 0.432850298736358291696660164313057377022426239952_wp*b2 &
+            - 0.432850298736358291696660164313057377022426239952_wp*b3
+    a(2,4) = 1.298550896209074875089980492939172131067278719856_wp*b1 &
+           + 1.298550896209074875089980492939172131067278719856_wp*b2 &
+           + 1.298550896209074875089980492939172131067278719856_wp*b3 &
+           - 0.13602180681849030343231048896787005715681064672_wp*b4
+    a(3,4) = -1.298550896209074875089980492939172131067278719856_wp*b1 &
+            - 1.298550896209074875089980492939172131067278719856_wp*b2 &
+            - 1.298550896209074875089980492939172131067278719856_wp*b3 &
+            + 0.40806542045547091029693146690361017147043194016_wp*b4 &
+            + 0.139708407134619223127669197053218276788706051412_wp*b5
+    a(4,4) = 0.432850298736358291696660164313057377022426239952_wp*b1 &
+           + 0.432850298736358291696660164313057377022426239952_wp*b2 &
+           + 0.432850298736358291696660164313057377022426239952_wp*b3 &
+           - 0.40806542045547091029693146690361017147043194016_wp*b4 &
+           - 0.419125221403857669383007591159654830366118154235_wp*b5 &
+           + 0.422237850757834010213265187828443211383102779457_wp*b6
+    a(5,4) = -0.231917524240705800245249424344181664391294057734_wp &
+            + 0.13602180681849030343231048896787005715681064672_wp*b4 &
+            + 0.419125221403857669383007591159654830366118154235_wp*b5 &
+            - 1.26671355227350203063979556348532963414930833837_wp*b6
+    a(6,4) = 0.632304366383038264396056896345266141026261929548_wp &
+           - 0.139708407134619223127669197053218276788706051412_wp*b5 &
+           + 1.26671355227350203063979556348532963414930833837_wp*b6
+    a(7,4) = -0.484284312078557755005703255606118914056056082961_wp &
+            - 0.422237850757834010213265187828443211383102779457_wp*b6
+    a(8,4) = -0.021798019655498318906737367811423460076462751066_wp
+    a(9,4) = 0.126819131218045846472604038781047420172516321492_wp
+    a(10,4) = -0.021123641626322236710970887364589522674965359279_wp
+    
+    a(1,5) = 0.04525671266969037217439909305498103903084448236_wp*b1 &
+           + 0.04525671266969037217439909305498103903084448236_wp*b2 &
+           + 0.04525671266969037217439909305498103903084448236_wp*b3
+    a(2,5) = -0.135770138009071116523197279164943117092533447081_wp*b1 &
+            - 0.135770138009071116523197279164943117092533447081_wp*b2 &
+            - 0.135770138009071116523197279164943117092533447081_wp*b3 &
+            - 0.356074970503407919432108002084434090165232256258_wp*b4
+    a(3,5) = 0.135770138009071116523197279164943117092533447081_wp*b1 &
+           + 0.135770138009071116523197279164943117092533447081_wp*b2 &
+           + 0.135770138009071116523197279164943117092533447081_wp*b3 &
+           + 1.068224911510223758296324006253302270495696768774_wp*b4 &
+           - 0.161252732964068670989303604260577833997844116999_wp*b5
+    a(4,5) = -0.04525671266969037217439909305498103903084448236_wp*b1 &
+            - 0.04525671266969037217439909305498103903084448236_wp*b2 &
+            - 0.04525671266969037217439909305498103903084448236_wp*b3 &
+            - 1.068224911510223758296324006253302270495696768774_wp*b4 &
+            + 0.483758198892206012967910812781733501993532350997_wp*b5 &
+            + 0.000913145598627662444167435877738793583486988286_wp*b6
+    a(5,5) = 0.778188619899941335162434896805956583938706965447_wp &
+           + 0.356074970503407919432108002084434090165232256258_wp*b4 &
+           - 0.483758198892206012967910812781733501993532350997_wp*b5 &
+           - 0.002739436795882987332502307633216380750460964857_wp*b6
+    a(6,5) = -2.888060650467623805586479517528702488708606593183_wp &
+            + 0.161252732964068670989303604260577833997844116999_wp*b5 &
+            + 0.002739436795882987332502307633216380750460964857_wp*b6
+    a(7,5) = 4.064047176584057140151337032365458641610672353212_wp &
+           - 0.000913145598627662444167435877738793583486988286_wp*b6
+    a(8,5) = -2.645663825945841938559482959983726831967447155003_wp
+    a(9,5) = 0.760485624510301003198698408956104774243768795869_wp
+    a(10,5) = -0.068996944580833734366507860615090679117094366342_wp
+    
+    a(1,6) = 0.234241201265594841727615788033296717516774995874_wp*b1 &
+           + 0.234241201265594841727615788033296717516774995874_wp*b2 &
+           + 0.234241201265594841727615788033296717516774995874_wp*b3
+    a(2,6) = -0.702723603796784525182847364099890152550324987621_wp*b1 &
+            - 0.702723603796784525182847364099890152550324987621_wp*b2 &
+            - 0.702723603796784525182847364099890152550324987621_wp*b3 &
+            + 0.112108192839478902858478028045243413779536843362_wp*b4
+    a(3,6) = 0.702723603796784525182847364099890152550324987621_wp*b1 &
+           + 0.702723603796784525182847364099890152550324987621_wp*b2 &
+           + 0.702723603796784525182847364099890152550324987621_wp*b3 &
+           - 0.336324578518436708575434084135730241338610530087_wp*b4 &
+           - 0.311645621737838510439662177349969575726926106823_wp*b5
+    a(4,6) = -0.234241201265594841727615788033296717516774995874_wp*b1 &
+            - 0.234241201265594841727615788033296717516774995874_wp*b2 &
+            - 0.234241201265594841727615788033296717516774995874_wp*b3 &
+            + 0.336324578518436708575434084135730241338610530087_wp*b4 &
+            + 0.934936865213515531318986532049908727180778320468_wp*b5 &
+            - 0.814776445474468482860158311580599018777969910384_wp*b6
+    a(5,6) = 1.318650558110101765365290411968234643910731822102_wp &
+           - 0.112108192839478902858478028045243413779536843362_wp*b4 &
+           - 0.934936865213515531318986532049908727180778320468_wp*b5 &
+           + 2.444329336423405448580474934741797056333909731151_wp*b6
+    a(6,6) = -4.321946659420587264930126779432660223055921154843_wp &
+            + 0.311645621737838510439662177349969575726926106823_wp*b5 &
+            - 2.444329336423405448580474934741797056333909731151_wp*b6
+    a(7,6) = 4.959408369623725459954942028713585060789912972866_wp &
+           + 0.814776445474468482860158311580599018777969910384_wp*b6
+    a(8,6) = -2.133050733448670443936969529227140283141530210559_wp
+    a(9,6) = 0.082410205158004740903168030202993056583347011382_wp
+    a(10,6) = 0.094528259977425742643695837774987744913459559052_wp
+    
+    a(1,7) = -0.021123641626322236710970887364589522673829626284_wp*b1 &
+            - 0.021123641626322236710970887364589522673829626284_wp*b2 &
+            - 0.021123641626322236710970887364589522673829626284_wp*b3
+    
+    a(2,7) = 0.063370924878966710132912662093768568021488878853_wp*b1 &
+           + 0.063370924878966710132912662093768568021488878853_wp*b2 &
+           + 0.063370924878966710132912662093768568021488878853_wp*b3 &
+           - 0.005626019701867024233595198521322111094015702611_wp*b4
+    a(3,7) = -0.063370924878966710132912662093768568021488878853_wp*b1 &
+            - 0.063370924878966710132912662093768568021488878853_wp*b2 &
+            - 0.063370924878966710132912662093768568021488878853_wp*b3 &
+            + 0.016878059105601072700785595563966333282047107834_wp*b4 &
+            + 0.23814816884096023561030675752649121424056599846_wp*b5
+    a(4,7) = 0.021123641626322236710970887364589522673829626284_wp*b1 &
+           + 0.021123641626322236710970887364589522673829626284_wp*b2 &
+           + 0.021123641626322236710970887364589522673829626284_wp*b3 &
+           - 0.016878059105601072700785595563966333282047107834_wp*b4 &
+           - 0.71444450652288070683092027257947364272169799538_wp*b5 &
+           + 0.504229942386996567949909229607651725484488940042_wp*b6
+    a(5,7) = -2.568464018793247324511052938299838510846083027331_wp &
+            + 0.005626019701867024233595198521322111094015702611_wp*b4 &
+            + 0.71444450652288070683092027257947364272169799538_wp*b5 &
+            - 1.512689827160989703849727688822955176453466820125_wp*b6
+    a(6,7) = 10.531909539708674517698111289364420212328345593768_wp &
+            - 0.23814816884096023561030675752649121424056599846_wp*b5 &
+            + 1.512689827160989703849727688822955176453466820125_wp*b6
+    a(7,7) = -17.15862642080199236829756567570752704681_wp &
+             - 0.504229942386996567949909229607651725484488940042_wp*b6
+    a(8,7) = 13.96906221208640324381455867393444497492_wp
+    a(9,7) = -5.747563226635290830973600786704797104493652582603_wp
+    a(10,7) = 0.973681914435452762269549437413297474901185356942_wp
+  end subroutine set_coeff
+
+end subroutine dissipation_6_5
diff --git a/src/dissipation.c b/src/dissipation.c
new file mode 100644
index 0000000..50d99a8
--- /dev/null
+++ b/src/dissipation.c
@@ -0,0 +1,119 @@
+/* $Header$ */
+
+#include <assert.h>
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+void CCTK_FCALL CCTK_FNAME(dissipation_6_5) (const CCTK_REAL *var,
+                                             const CCTK_INT *ni,
+                                             const CCTK_INT *nj,
+                                             const CCTK_INT *nk,
+                                             const CCTK_INT *bbox,
+                                             const CCTK_INT *gsize,
+                                             const CCTK_REAL *dx,
+                                             const CCTK_REAL *epsdis,
+                                             CCTK_REAL *rhs);
+
+static void
+apply (int const varindex, char const * const optstring, void * const arg);
+
+void
+SBP_DissipationAdd (CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  CCTK_TraverseString (vars, apply, cctkGH, CCTK_GROUP_OR_VAR);
+}
+
+void
+apply (int const varindex, char const * const optstring, void * const arg)
+{
+  cGH const * const cctkGH = (cGH const *) arg;
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  int rhsindex;
+  int vargroup, rhsgroup;
+  cGroup vardata, rhsdata;
+  CCTK_REAL const * varptr;
+  CCTK_REAL       * rhsptr;
+  CCTK_INT ni, nj, nk, bbox[6], gsize[3];
+  CCTK_REAL dx[3];
+  int n;
+  int d;
+  int ierr;
+  
+  assert (varindex >= 0);
+  
+  for (d=0; d<3; ++d) {
+    dx[d] = CCTK_DELTA_SPACE(d);
+    gsize[d] = cctk_nghostzones[d];
+  }
+  
+  for (d=0; d<6; ++d) {
+    bbox[d] = cctk_bbox[d];
+  }
+
+  rhsindex = MoLQueryEvolvedRHS (varindex);
+  if (rhsindex < 0) {
+    char * const fullvarname = CCTK_FullName (varindex);
+    assert (fullvarname);
+    CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                "There is no RHS variable registered with MoL for the evolved variable \"%s\"",
+                fullvarname);
+    free (fullvarname);
+  }
+  assert (rhsindex >= 0);
+  
+/*  if (verbose) {
+    char * const fullvarname = CCTK_FullName (varindex);
+    char * const fullrhsname = CCTK_FullName (rhsindex);
+    assert (fullvarname);
+    assert (fullrhsname);
+    CCTK_VInfo (CCTK_THORNSTRING,
+                "Applying dissipation to \"%s\" (RHS \"%s\")",
+                fullvarname, fullrhsname);
+    free (fullvarname);
+    free (fullrhsname);
+  } */
+  
+  vargroup = CCTK_GroupIndexFromVarI (varindex);
+  assert (vargroup >= 0);
+  rhsgroup = CCTK_GroupIndexFromVarI (rhsindex);
+  assert (rhsgroup >= 0);
+  
+  ierr = CCTK_GroupData (vargroup, &vardata);
+  assert (!ierr);
+  ierr = CCTK_GroupData (rhsgroup, &rhsdata);
+  assert (!ierr);
+  
+  assert (vardata.grouptype == CCTK_GF);
+  assert (vardata.vartype == CCTK_VARIABLE_REAL);
+  assert (vardata.dim == cctk_dim);
+  assert (rhsdata.grouptype == CCTK_GF);
+  assert (rhsdata.vartype == CCTK_VARIABLE_REAL);
+  assert (rhsdata.dim == cctk_dim);
+  
+  varptr = CCTK_VarDataPtrI (cctkGH, 0, varindex);
+  assert (varptr);
+  rhsptr = CCTK_VarDataPtrI (cctkGH, 0, rhsindex);
+  assert (rhsptr);
+  
+  if ( CCTK_Equals(norm_type,"Diagonal") ) {
+    assert(0);
+  } else {
+    switch(order) {
+    case 6: {
+      CCTK_FNAME(dissipation_6_5)
+        (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+                            bbox, gsize, dx, &epsdis, rhsptr);
+      break;
+    }
+    default:
+      assert(0);
+    }
+  }
+}
diff --git a/src/make.code.defn b/src/make.code.defn
index cc9407d..e835d0f 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
 # $Header$
 
 # Source files in this directory
-SRCS =  Derivatives_2_1.F90 Derivatives_4_2.F90 Derivatives_6_3.F90 Derivatives_8_4.F90 call_derivs.c Derivatives_4_3.F90 Derivatives_6_5.F90 call_derivs_name.c Get_Coeff.F90 set_norm_mask.F90 CheckGridSizes.F90
+SRCS =  Derivatives_2_1.F90 Derivatives_4_2.F90 Derivatives_6_3.F90 Derivatives_8_4.F90 call_derivs.c Derivatives_4_3.F90 Derivatives_6_5.F90 call_derivs_name.c Get_Coeff.F90 set_norm_mask.F90 CheckGridSizes.F90 dissipation.c Dissipation_6_5.F90
 
 # Subdirectories containing source files
 SUBDIRS =